# Dynamic Sizing of Frequency Control Ancillary Service Requirements for a
Philippine Grid

Elgar John S. del Rosario1, Jordan Rel C. Orillaza2 Electrical and Electronics
Engineering Institute
University of the Philippines Diliman
Quezon City, Philippines
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Sizing frequency control ancillary service (FCAS) requirements is crucial for
the reliable operation of power systems amid a continuous influx of variable
renewable energy (VRE) generation. Reserve sizing is especially pertinent for
the Philippine grids due to an expected transition to new FCAS classifications
established by its Grid Code. In lieu of the existing deterministic
formulation, this work proposes a dynamic approach for sizing secondary and
tertiary reserves that accounts for the stochasticity and variability of load
demand and VRE. We propose a method where historical power imbalances were
calculated and clustered according to the time and day of week they occurred.
The conditional probabilities of forecast and noise errors were characterized
using kernel density estimation. Recursive convolution was performed to obtain
the total reserve requirement probability distribution. The method was tested
on Visayas grid’s historical system operation data and used target reliability
levels on the error distributions to size upward and downward reserve needs.
Finally, the methodology was extended to demonstrate through a numerical
experiment that sizing FCAS at temporal resolutions higher than one-hour,
e.g., five-minute, provides the benefit of shrinking the required capacities
by as much as 86.2% compared to current deterministic FCAS sizing.

###### Index Terms:

frequency control ancillary services, dynamic sizing, secondary reserves,
tertiary reserves, kernel density estimation, convolution, power system
reliability

##  I Introduction

Continuous and precise frequency control is critical in maintaining power
system reliability. Frequency deviations are counteracted by deploying active
power control reserves, also known as frequency control ancillary services
(FCAS) [1]. The adequate sizing of reserves is crucial as, on the one hand,
undersizing can have serious negative impact, including load shedding,
renewable energy curtailment, equipment damage, and in the worst scenario,
blackouts. Reserve procurement costs, on the other hand, will rise if reserves
are oversized. The Philippine electric power industry is also about to
activate a reserve market, hence FCAS sizing is more critical [2].

Owing to the increasing grid integration of variable renewable energy (VRE)
sources, such as wind and solar PV generation, quantifying reserve
requirements considering VREs has received considerable attention in the last
decade [3]. While VRE growth promotes the development of a clean and
sustainable grid, VREs are also nonsynchronous, highly fluctuating and
difficult to predict, and could induce higher levels of FCAS requirements [4].
Dynamic methods of sizing FCAS requirements due to the increasing integration
of VREs using kernel density estimation [5], k-nearest neighbors and k-means
clustering [6, 7], and dynamic Bayesian belief networks [8] have been proposed
in the literature. In contrast with deterministic approaches that are
typically based on rules of thumb, these dynamic probabilistic approaches vary
reserve requirements depending on expected system conditions, considering the
severity and probabilities of a range of potential power imbalances.

In the Philippines, reserves are still sized based on deterministic rules.
Presently the Philippine system operator sizes three types of FCAS on a day-
ahead basis: regulating, contingency, and dispatchable reserves. Regulating
reserves are set at 4% of the hourly forecast demand in each hourly dispatch
interval, which is based on intra-hour load variations in 2010 [9]. Meanwhile,
contingency and dispatchable reserve levels are determined by the highest and
second highest generating power outputs expected to be online, in accordance
with the N-1 and N-1-1 criteria, respectively.

The current reserve classifications lack requirements for primary frequency
response, often resulting in the activation of automatic load dropping (ALD)
schemes as a first line of defense against large disturbances [10]. To address
this gap, a revision of the Philippine Grid Code of 2016 prescribes new FCAS
classifications according to the hierarchy of frequency control. In the new
FCAS classification, primary reserves must provide primary frequency response
via governor control, secondary reserves must operate to restore the system
frequency to 60 Hz either under automatic generation control (AGC) or manually
upon the command of the system operator, and tertiary reserves aim to
replenish depleted secondary reserves [10]. This development in FCAS framework
necessitates a new approach in sizing reserves.

Furthermore, current practices have yet to determine the incremental reserve
needs introduced by the growing VRE penetration. As a country indigenously
rich in VRE resources, the Philippines targets to achieve 35% renewable energy
penetration by 2030 and 50% by 2040 [11]. In NREL’s Greening the Grid report,
high renewable energy scenarios for the Luzon-Visayas system of the
Philippines by 2030 were investigated and the issue of adequate reserve
provision was raised [12]. However, the said study retained an assumption of
deterministic reserve rules that govern system operations until now due to the
absence of detailed reserve requirements.

The increasingly stochastic nature of system operations demand dynamic, data-
driven and probabilistic approaches that consider the additional variability
and uncertainty introduced by VREs. Considering this fact, this work proposes
a dynamic probabilistic approach as an appropriate methodology for sizing
reserves for frequency restoration, particularly, secondary and tertiary
reserves, in the context of the Philippines. The method is based on clustering
imbalances according to the hour and day they happened and kernel density
estimation.

Additionally, a recent enhancement of the Philippine Wholesale Electricity
Spot Market (WESM) shifts the energy markets in the Luzon and Visayas grids
from one-hour to five-minute intervals [13]. While reserve schedules continue
to be based on hourly requirements, the Department of Energy tasks the system
operator to conduct studies on the determination of required reserves vis-à-
vis the 5-minute dispatch interval implementation [2]. This work further
extends the aforementioned dynamic reserve sizing methodology to perform
numerical experiments for subhourly sizing.

In Section 2, the proposed methodology is discussed. Section 3 provides the
results and discussion on scenario simulations done on real Philippine grid
operations data. Section IV concludes with recommendations for further
research.

##  II Proposed FCAS Sizing Methodology

An overview of the proposed reserve sizing process is shown in Fig. 1.
Probability density estimation of various drivers for power imbalance is done
based on historical values of power imbalances, which are mainly due to the
discrepancy between the forecast and actual values of demand, wind and solar
generation, as well as forced (unplanned) outages of power plants or
interconnectors. The estimation process must be continuously done to
incorporate new data and make the FCAS sizing process more accurate. On top of
the estimated probability distribution functions (PDFs), load and VRE
forecasts are employed as basis for sizing. The status quo reserve sizing is
done on the day-ahead for sizing reserves for each hour of the next day, but
sizing and scheduling reserves closer to real time (n-hours ahead) may prove
to benefit forecast accuracy, especially when reserves are sized for subhourly
intervals. The result of FCAS sizing shall then be the basis of FCAS
procurement and the separate upward and downward capacities reserved for each
intra-day interval.

An overview of the methodology for reserve sizing is shown in Fig. 2. The
necessary data are collected for the calculation of historical imbalances. The
historical imbalances are used to estimate the power imbalance PDFs. The
convolution of these power imbalance PDFs produce PDFs for total and secondary
reserve needs. The reserve requirements are determined by applying a target
reliability level, referring to the acceptable percentage of time that there
could be a shortage of reserves. The tertiary reserve requirement is taken as
the difference between the derived total and secondary reserves. The following
sections expound on these steps.

Figure 1: Proposed framework for dynamic FCAS sizing in the Philippines

Figure 2: Overview of Proposed Methodology

###  II-A Estimation of Error Probability Distributions

####  II-A1 Forecast and Noise Errors

Historical time series of load, wind and solar PV forecast and actual values
are used to calculate forecast and noise errors. The load forecast error is
the deviation of the mean value of load in a time interval from the forecast
value, whereas the load noise error (or “load oscillation”) is the deviation
of the actual value from the interval mean [17]. Calculating forecast and
noise errors are depicted in Figs. 3(a) and 3(b).

(a) Forecast Errors

(b) Noise Errors

Figure 3: Calculation of Forecast and Noise Errors

To avoid making strong assumptions on the shape of the underlying probability
distributions, this work uses kernel density estimation to estimate the PDFs
of wind, solar and load forecast and noise errors as in [6]. Kernel density
estimation has been extensively used to characterize wind speed distributions
and wind forecast errors in the literature [14]. Given the forecast and noise
error values, $\epsilon_{i,t}$, the probability distribution function
$f_{i}(u)$ is given by:

$f_{i}(u)=\frac{1}{T{h_{i}}}\sum_{t=1}^{T}K\left(\frac{u-\epsilon_{i,t}}{h_{i}}\right)$
(1)

where $K(\cdot)$ is the kernel smoothing function; $u$, the evaluation point
for function $f_{i}$; $i\in[1,24\times 7]$ denotes an index to a cluster,
which, in this work, is the hour and day of week of the error; $t$, an index
for time step; and $T$, the maximum time step. The kernel smoothing function
is assumed to be a normal probability density function [5, 6]. The optimal
bandwidth is proportional to the standard deviation of the cluster $i$,
$\sigma_{i}$, and can be calculated from the empirical data as [15]:

$h_{i}=\left(\frac{4}{3T}\right)^{0.2}\sigma_{i}.$ (2)

From kernel density estimation, we get six sets of PDFs, representing forecast
and noise errors for load, wind and solar PV generation for each hour of the
week: $f^{x}_{y}$ where $x\in\\{\text{``forecast'', ``noise''}\\}$ and
$y\in\\{\text{``load'', ``wind'', ``solar''}\\}$.

####  II-A2 Power Plant Outages

The forced outage rate, $FOR$, of a single generating plant or unit is
calculated from the number of hours it was on forced outage, $N_{outage}$, and
the total number of hours in the period considered, $N_{hours}$:

$FOR=\frac{N_{outage}}{N_{hours}}.$ (3)

The forced outage probability, $FOP$, reflects the probability of generation
capacity going out, instead of being out [16]:

$FOP=\frac{FOR}{MTTR}$ (4)

where $MTTR$ refers to the mean time to repair. In determining the total
capacity outage distribution, the FOP of each generating plant or unit is
first calculated. The error distribution of a single plant is a piecewise-
defined function

$f_{j}^{\text{outage}}(X)=\begin{cases}FOP_{j},&\text{if }X=0\\\
1-FOP_{j},&\text{if }X=P^{\text{rated}}_{j}\end{cases}$ (5)

where $j$ is an index for a generating station and $P^{\text{rated}}_{j}$ is
the generating capacity of station $j$. Then the error distributions of all
power stations are recursively convolved to create a single error
distribution[5], the total capacity outage distribution,
$f^{\text{outage}}_{\text{total}}$:

$f^{\text{outage}}_{\text{total}}=f^{\text{outage}}_{1}*f^{\text{outage}}_{2}*...*f^{\text{outage}}_{N_{\text{gen}}-1}*f^{\text{outage}}_{N_{\text{gen}}}.$
(6)

where $N_{\text{gen}}$ is the total number of generating stations and $*$ is a
shorthand operator for convolution. Convolution is defined for two random
variables X and Y and their sum $\textbf{Z}=\textbf{X}+\textbf{Y}$, with their
respective PDFs $f_{X},f_{Y},f_{Z}$:

$f_{Z}=f_{X}*f_{Y}$ (7a) $f_{Z}(z)=\sum_{k=-\infty}^{\infty}f_{X}(k)\cdot
f_{Y}(z-k)$ (7b)

###  II-B Determination of Reserve Requirements

####  II-B1 Convolution of Error PDFs

As each power imbalance driver is treated as a random variable, convolution is
employed to determine the PDF of the sum of these imbalance drivers. The total
reserve PDF is derived from convolution of the PDFs of all the defined power
imbalance drivers:

$\begin{split}f^{\text{total}}_{\text{reserve}}=f^{\text{forecast}}_{\text{load}}*f^{\text{noise}}_{\text{load}}*f^{\text{forecast}}_{\text{wind}}*f^{\text{noise}}_{\text{wind}}\\\
*f^{\text{forecast}}_{\text{solar}}*f^{\text{noise}}_{\text{solar}}*f^{\text{outage}}_{\text{total}}.\end{split}$
(8)

Figure 4: Power imbalance driver PDFs convolved to derive the PDF of total
reserve needs. Data based on 2018-2019 Visayas system operations.

Figure 4 depicts the convolution of seven individual PDFs. Next, the secondary
reserve PDF is the convolution of the PDFs of only the power imbalance drivers
that need to be handled by secondary control, namely noise errors (fast
fluctuations of demand and VRE) and capacity outages [17]:

$\begin{split}f^{\text{secondary}}_{\text{reserve}}=f^{\text{noise}}_{\text{load}}*f^{\text{noise}}_{\text{wind}}*f^{\text{noise}}_{\text{solar}}*f^{\text{outage}}_{\text{total}}.\end{split}$
(9)

####  II-B2 Scaling of Error PDFs

The predicted demand, wind and solar PV capacities for future scenarios are
used to determine the reserve requirements. The PDFs of the power imbalance
drivers are scaled by a growth factor equivalent to the ratio of the forecast
quantity, $P^{\text{future}}_{i}$ [MW] to its historical peak value
$P^{\text{base}}_{i}$ [MW]. Furthermore, to account for future improvements in
forecast performance, a forecast improvement factor,
$k^{\text{forecast}}_{i}$, can be multiplied to the error values to get the
scaled error values $\epsilon_{i}^{\text{future}}$.

$\epsilon_{i}^{\text{future}}=\epsilon_{i}^{\text{base}}\cdot
k_{i}^{\text{forecast}}\cdot\frac{P^{\text{future}}_{i}}{P^{\text{base}}_{i}}$
(10)

####  II-B3 Application of Reliability Margins

The reliability margin $\rho_{margin}$, also referred to as the security
level, is a predefined target that represents the fraction of time, typically
over a one-year period, in which reserves are adequate to serve the demand [4,
6]. The deficit ($\rho_{\text{deficit}}$) and surplus
($\rho_{\text{surplus}}$) probabilities correspond to the the fractions of
time that shortage of downward and upward reserves, respectively, are
acceptable [17], and they are set to be equal in this work:

$\rho_{\text{deficit}}=\rho_{\text{surplus}}=\frac{100\%-\rho_{\text{margin}}}{2}.$
(11)

The cumulative distribution functions $F$ of the total and secondary reserve
are derived as follows:

$F_{\text{reserve}}(z)=\int_{-\infty}^{z}f_{\text{reserve}}(z)dz.$ (12)

The reserve requirements $R$ are determined by seeking the values of $F$ that
satisfy the following inequality conditions with respect to the set values of
$\rho_{\text{deficit}}$ and $\rho_{\text{surplus}}$:

$F_{\text{reserve}}\left(R_{\text{up}}\right)\leq 1-\rho_{\text{surplus}}$
(13a)
$F_{\text{reserve}}\left(R_{\text{down}}\right)\geq\rho_{\text{deficit}}$
(13b)

Equations 13a and 13b apply to both total and secondary reserve requirements.
Finally, tertiary reserves are taken as the difference between the total and
secondary reserve requirements:

$R^{\text{tertiary}}=R^{\text{total}}-R^{\text{secondary}}.$ (14)

Equation 14 applies to both upward ($R_{\text{up}}$) and downward
($R_{\text{down}}$) reserve requirements. The proposed separation scheme
between secondary (fast-response) and tertiary (slow-response) reserves is
considered cost-effective because it reduces secondary reserves, which are
generally more expensive than tertiary reserves [18].

##  III Results and Discussion

The Philippines is comprised of three synchronous grids, Luzon, Visayas, and
Mindanao, each with a regional system operator with the National Grid
Corporation of the Philippines (NGCP). Reserve requirements are determined for
each region separately at present. This section demonstrates the reserve
sizing using two years’ worth of data for the Visayas grid, which has the
highest VRE share among the three grids. In 2021, solar accounted for 12.5%
share of total capacity and 4.01% share of total gross generation, whereas
wind accounted for a 2.4% share of total capacity and 1.07% share of total
gross generation in the Visayas grid [19].

A one-week representative scenario is developed from real Visayas system
operations data. The forecast values for each hour of the week are set to be
equivalent to the maximum historical values of forecast demand and wind and
solar PV power generation that were observed in the same hour of the week in
years 2018-2019. $k_{i}^{\text{forecast}}=1$ is assumed, which implies that
the error distributions did not change from the historical period used. The
following historical data from the NGCP for the Visayas grid were used:

  1. 1.

hourly forecast and per-minute actual system demand from 1 January 2018 to 31
December 2019.

  2. 2.

hourly forecast VRE generation from 2 wind and 11 solar PV plants from 1 Jan.
2018 to 1 Dec. 2019.

  3. 3.

per-minute actual VRE generation from 1 Jan. 2018 to 31 Dec. 2019.

  4. 4.

forced outages from 1 Jan. 2016 to 31 Dec. 2019.

###  III-A Static and Dynamic Secondary Reserve vs. Status Quo RR Requirements

In this section, we compare the secondary reserve requirement with status quo
regulating reserve requirement and the results of a static method. Secondary
reserves and regulating reserves (RR) as currently analogous in Philippine
practice as they are both intended to be operated under automatic generation
control to restore the frequency to its nominal value. The RR requirement is
set at 2% of the forecast demand for upward regulation and 2% of the forecast
demand for downward regulation [20].

On the other hand, static and dynamic methods require the same types of error
distributions, but differ in the forecast data used to estimate the PDF and
size the reserves. The static method is based on the peak forecast in the
whole period for which reserves will be sized, and all historical values of an
imbalance driver are included in the PDF estimation, regardless of its time of
occurrence. The dynamic method is based on hourly forecasts. Every hour of the
week has a distinct PDF for each imbalance driver, wherein only historical
values that occurred in a similar hour of the week are included.

The results of dynamic secondary reserve sizing is shown in Figure 5. The mean
dynamic reserve requirements are lower than the static reserve requirements,
as likewise observed in previous studies [5, 6]. Furthermore, both static and
mean dynamic reserve requirements are larger than the current $\pm$2%
requirement for regulating reserves. Throughout the one-week scenario, it
appears that most hours require levels higher than the current regulating
reserve requirements. Hours 1-7 and 20 provide exceptions as they are well
within the status quo requirement on most days. Based on a 99.0% reliability
margin for a one-week Visayas grid scenario, the study reveals a need for
asymmetric upward and downward reserves, which at some hours may be greater
than the levels determined by existing practice.

Figure 5: Dynamic Upward (+) and Downward (-) Secondary Reserves vs. Static
and Regulating Reserve Requirements

###  III-B Subhourly Reserve Sizing

To contribute to the discussions on reserve requirements in view of the
subhourly dispatch regime, we apply the proposed methodology to quantify the
potential effect of varying the reserve lengths from one-hour to higher
resolutions of 30-minute, 15-minute, and 5-minute.

As subhourly forecasts were not yet available in 2018-2019, we synthesized
subhourly forecasts of demand and VRE generation using historical actual
measurements. At the end of each dispatch interval, we assumed that subhourly
forecast values are equivalent to the actual measurements. Thus, only the
power imbalances within the interval are allocated reserves.

Table I shows the results of subhourly sizing for total reserves. Compared to
hourly reserves, more granular reserves are considerably reduced. In
particular, five-minute reserve sizing can reduce the total reserve
requirements by as much as 86.2% on average. The results suggest economic
benefits from the decreased reserve need while still upholding the required
level of system reliability. More research is needed to determine the
necessary adjustments in operations, as well as regulations, to ensure that
dynamic reserve sizing of five-minute lengths, or subhourly lengths, is indeed
feasible.

TABLE I: Average Dynamic Total Reserve Requirements at Subhourly Intervals and 99.0% Reliability Time | Downward | Upward
---|---|---
Res. | Mean[MW] | %Reduction1 | Mean[MW] | %Reduction1
60-min | 356.0 | 0% | 205.1 | 0%
30-min | 201.2 | -43.5% | 122.3 | -40.4%
15-min | 115.6 | -67.5% | 75.3 | -63.3%
5-min | 49.1 | -86.2% | 35.4 | -82.7%
1 =(mean subhourly total - mean hourly total)/(mean hourly total).

##  IV Conclusions and Recommendations

A dynamic reserve sizing methodology is proposed as the Philippine grids
transition to a new set of reserve classification and face increasing levels
of VRE penetration. In comparison with the current allocation method, the
proposed method indicate a larger capacities for reserves to achieve a 99.0
percent reliability. As the reserve requirements are anchored to a
predetermined reliability target, further study is required to determine the
optimal reliability level to aspire to, ideally considering the prices of
reservation and activation.

Numerical experiments also revealed that dynamic five-minute reserve sizing
can reduce reserve requirements by as much as 86.2% compared to hourly
requirements. New five-minute forecasts from recent operations and higher
resolution actual data would be preferred in order to assess subhourly sizing
more precisely. Validating the effectiveness of the proposed method through
pilot implementation and the evaluation of its effect on frequency quality can
also be put forward.

## References

  * [1] E. Ela, M. Milligan, and B. Kirby, “Operating Reserves and Variable Generation,” NREL/TP-5500-51978, 1023095, Aug. 2011.
  * [2] Philippine Department of Energy, DC2021-03-0009: Adopting A General Framework Governing the Operationalization of the Reserve Market in the Wholesale Electricity Spot Market and Providing Further Policies to Supplement DC2019-12-0018, 2019.
  * [3] H. Holttinen et al., ”Methodologies to Determine Operating Reserves Due to Increased Wind Power,” in IEEE Transactions on Sustainable Energy, vol. 3, no. 4, pp. 713-723, Oct. 2012.
  * [4] L. Hirth and I. Ziegenhagen, ”Balancing power and variable renewables: Three links,” Renewable and Sustainable Energy Reviews, vol. 50, pp. 1035-1051, 2015.
  * [5] D. Jost, M. Speckmann, F. Sandau, R. Schwinn, ”A new method for day-ahead sizing of control reserve in Germany under a 100% renewable energy sources scenario,” Electric Power Systems Research, vol. 119, pp. 485-491, 2015.
  * [6] M. Bucksteeg, L. Niesen, and C. Weber, “Impacts of Dynamic Probabilistic Reserve Sizing Techniques on Reserve Requirements and System Costs,” IEEE Transactions on Sustainable Energy, vol. 7, no. 4, pp. 1408–1420, Oct. 2016.
  * [7] K. De Vos, N. Stevens, O. Devolder, A. Papavasiliou, B. Hebb, and J. Matthys-Donnadieu, “Dynamic dimensioning approach for operating reserves: Proof of concept in Belgium,” Energy Policy, vol. 124, pp. 272–285, Jan. 2019.
  * [8] F. Fahiman, S. Disano, S. M. Erfani, P. Mancarella and C. Leckie, ”Data-Driven Dynamic Probabilistic Reserve Sizing Based on Dynamic Bayesian Belief Networks,” in IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 2281-2291, May 2019.
  * [9] Philippine Energy Regulatory Commission, Ancillary Services Procurement Plan, Draft Nov. 2011.
  * [10] Philippine Energy Regulatory Commission, The Philippine Grid Code: 2016 Edition, Resolution No. 22, Series of 2016.
  * [11] Philippine Department of Energy, National Renewable Energy Plan 2020-2040, 2021.
  * [12] C. P. Barrows et al., “Greening the Grid: Solar and Wind Grid Integration Study for the Luzon-Visayas System of the Philippines,” NREL/TP–6A20-68594, 1418740, Jan. 2018.
  * [13] Philippine Wholesale Electricity Spot Market. Electric Power Industry Reform Act of 2001: Wholesale Electricity Spot Market Rules (WESM Rules), May 2022.
  * [14] Y. Zhang, J. Wang, and X. Wang, ”Review on probabilistic forecasting of wind power generation,” Renewable and Sustainable Energy Reviews, vol. 32, pp. 255-270, 2014.
  * [15] A.W. Bowman and A. Azzalini. Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations. Vol. 18. OUP Oxford, 1997.
  * [16] R. Doherty and M. O’Malley, ”A new approach to quantify reserve demand in systems with significant installed wind capacity,” in IEEE Transactions on Power Systems, vol. 20, no. 2, pp. 587-595, May 2005.
  * [17] C. Maurer, S. Krahl, and H. Weber, “Dimensioning of secondary and tertiary control reserve by probabilistic methods,” European Transactions on Electrical Power, vol. 19, no. 4, pp. 544–552, May 2009.
  * [18] K. De Vos, J. Morbee, J. Driesen and R. Belmans, ”Impact of wind power on sizing and allocation of reserve requirements,” IET Renewable Power Generation, 7: 1-9, 2013.
  * [19] Philippine Department of Energy, 2021 Power Statistics as of 31 December 2021, June 2022.
  * [20] Philippine Wholesale Electricity Spot Market. Protocol for Central Scheduling and Dispatch of Energy and Contracted Reserves Issue 3.0, WESM Manual, June 2021.


# An interpretation of spacetimes with expanding impulsive gravitational waves
generated by snapped cosmic strings

D. Kofroň1, M. Karamazov1, and R. Švarc1

1 Institute of Theoretical Physics,
Charles University, Faculty of Mathematics and Physics,
V Holešovičkách 2, 180 00 Prague 8, Czech Republic.

d.kofron@gmail.commichal.karamazov<EMAIL_ADDRESS>
###### Abstract

The geometric properties of spacetimes representing expanding impulsive
gravitational waves, propagating on a flat background and generated by snapped
cosmic strings, are studied. The construction of the line element is reviewed,
and suitable forms of the string-generating complex mapping are derived for
various configurations such as previously studied examples of a pair of
snapping cosmic strings. Moreover, these mappings are related to the topology
of the flat half-space in front of the wave. Their understanding seems to be
crucial for further analysis of the global geometry, the relation between
half-spaces on both sides of the impulse, and the physical interpretation of,
in principle, observable effects. The spacetime structure is connected with
the motion of free test particles crossing the impulse, where the recent
results allow us to discuss their displacement and induced velocities that are
caused by the interaction with the expanding gravitational impulsive wave.

##  1 Introduction

Interestingly, the prediction of gravitational waves (GW) became an immediate
consequence of Einstein’s geometric description of gravity in terms of curved
spacetime. The analysis of weak field regime [1, 2] as well as pioneering
studies of invariants within full theory [3, 4] pointed out their inherent
properties and gave quantitative estimates about the magnitude of typical
amplitudes. Simultaneously, it became clear that their observation is going to
be very challenging or even impossible. In 1974 Hulse and Taylor discovered
the binary pulsar PSR 1913+16 [5]. The systematic measurements of its orbital
period shifts exactly corresponded to the gravitational energy loss radiated
away in the form of GWs. This was in perfect agreement with the general
relativity (GR) prediction and endorsed the hunt for their direct detection.
The unrelenting effort was rewarded in 2015 when the LIGO facility detected a
signal produced by merging pair of black holes more than one billion light-
years away [6]. Until now, dozens of other events have been detected with the
prominent case of merging neutron stars [7], where the GW signal was
supplemented by observations in the whole range of the electromagnetic
spectrum.

These ultimate experimental achievements are based on sophisticated numerical
simulations. However, to better understand the GW properties, or to verify
various models and numerical schemes, even the exact analytical radiative
spacetimes are necessary. The most important wave-like solutions to Einstein’s
general relativity belong to the expanding Robinson–Trautman and non-expanding
Kundt classes of exact spacetimes [8, 9, 10, 11], see [12, 13] for the
comprehensive review. Interestingly, both these classes allow for the
_impulsive_ profiles of the wave, where the propagating curvature is located
only on the singular null hypersurface. Technically, it corresponds to the
presence of distributional terms in the spacetime geometric quantities.
Therefore, the impulsive solutions are interpreted as extremely short, but
simultaneously very intense, bursts of gravitational radiation. Surprisingly,
these geometries are also of a purely mathematical interest since careful
manipulation with inherently non-linear distributional terms is required that
necessarily goes beyond the classical GR, see e.g. [14, 15].

Within this paper, we are interested in particular aspects of expanding
Robinson–Trautman impulses propagating on a flat Minkowski background. It is
worth mentioning that the elegant geometric construction of these models goes
back to the seventies when Roger Penrose proposed the ‘cut and paste’ approach
[16]. Two decades later Penrose and Nutku found the continuous line element
for these solutions [17]. Subsequently, various extensions of such a
construction were presented, e.g, including the cosmological constant or other
additional parameters [18, 19, 20, 21]. The historical content, summary of
other construction methods (e.g., limits of expanding sandwich waves or
infinite acceleration limit of the C-metric), and detailed list of references
can be found, e.g., in [13, 22, 23]. In the original work [17], the snapped
cosmic string serving as a wave source was described and the possibility of
impulses generated by a pair of colliding and snapping strings was outlined.
This scenario was elaborated by Podolský and Griffiths [24]. However, the
complete understanding of the spacetime topology remained an open problem.
Simultaneously, the topology plays a crucial role in the analysis of a
geodesic motion affected by the induced impulsive wave, and vice versa, the
geodesic motion reflects geometry of the impulse. Therefore, our present
contribution aims to fill up this gap and extend discussion of the wave
sources. In section 2, the description of expanding impulsive waves is
summarized. Subsequently, the construction of string-like geometries is
reviewed in section 3 and the description of geodesics is presented in section
4. In the last two sections 5 and 6, we discuss and interpret specific
properties of one and two-string geometries, respectively, analysing the
interaction of corresponding gravitational impulses with free test particles.

##  2 Expanding gravitational impulses on a flat background

It is natural to begin with a description of the Penrose geometric ‘cut and
paste’ construction in the simplest situation of an expanding impulse
propagating on a flat Minkowski background

${\rm d}s^{2}=-{\rm d}t^{2}+{\rm d}x^{2}+{\rm d}y^{2}+{\rm d}z^{2}\,,$ (1)

which can be simply rewritten using the double null coordinates,

$t=\frac{1}{\sqrt{2}}({\cal V}+{\cal U})\,,\qquad z=\frac{1}{\sqrt{2}}({\cal
V}-{\cal U})\,,\qquad x=\frac{1}{\sqrt{2}}(\eta+\bar{\eta})\,,\qquad
y=\frac{1}{{i}\sqrt{2}}(\eta-\bar{\eta})\,,$ (2)

to the form

${\rm d}s^{2}=-2{\rm d}{\cal U}{\rm d}{\cal V}+2{\rm d}\eta{\rm
d}\bar{\eta}\,.$ (3)

Applying further transformation to the Minkowski metric (3),

${\cal V}=\frac{V}{p}-\epsilon U\,,\qquad{\cal
U}=\frac{|Z|^{2}}{p}\,V-U\,,\qquad\eta\,=\frac{Z}{p}\,V$ (4)

with

$p=1+\epsilon Z\bar{Z}\,,\qquad\hbox{where}\qquad\epsilon=-1,0,+1\,,$ (5)

we get the line element

${\rm d}s^{2}=2\,\frac{V^{2}}{p^{2}}\,{\rm d}Z{\rm d}\bar{Z}+2\,{\rm d}U{\rm
d}V-2\epsilon\,{\rm d}U^{2}\,,$ (6)

which explicitly describes the foliation of the flat spacetime by null cones
that are labelled by constant values of the coordinate $U$. The parameter
$\epsilon$ encodes the Gaussian curvature of spatial two-surfaces
${U=\mbox{const}}$ and ${V=\mbox{const}}$, see [13] for the detailed geometric
picture.

However, another more involved transformation of (3) can be performed
(simultaneously leading to the explicit null-cone foliation), namely

$\displaystyle{\cal V}=AV-DU\,,\qquad{\cal U}=BV-EU\,,\qquad\eta\,=CV-FU\,,$
(7)

where

$\displaystyle A=$ $\displaystyle\frac{1}{p|h^{\prime}|}\,,$ $\displaystyle
D=$
$\displaystyle\frac{1}{|h^{\prime}|}\left\\{\frac{p}{4}\left|\frac{h^{\prime\prime}}{h^{\prime}}\right|^{2}+\epsilon\left[1+\frac{Z}{2}\frac{h^{\prime\prime}}{h^{\prime}}+\frac{\bar{Z}}{2}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}\right]\right\\},$
$\displaystyle B=$ $\displaystyle\frac{|h|^{2}}{p|h^{\prime}|}\,,$
$\displaystyle E=$
$\displaystyle\frac{|h|^{2}}{|h^{\prime}|}\bigg{\\{}\frac{p}{4}\left|\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right|^{2}+\epsilon\left[1+\frac{Z}{2}\left(\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right)+\frac{\bar{Z}}{2}\left(\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}-2\frac{\bar{h}^{\prime}}{\bar{h}}\right)\right]\bigg{\\}},$
(8) $\displaystyle C=$ $\displaystyle\frac{h}{p|h^{\prime}|}\,,$
$\displaystyle F=$
$\displaystyle\frac{h}{|h^{\prime}|}\bigg{\\{}\frac{p}{4}\left(\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right)\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}+\epsilon\left[1+\frac{Z}{2}\left(\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right)+\frac{\bar{Z}}{2}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}\right]\bigg{\\}},$

with ${h=h(Z)}$ representing an arbitrary complex holomorphic function (apart
from its singular points) and a prime denoting its derivative. The resulting
line element becomes

${\rm d}s^{2}=2\left|\frac{V}{p}\,{\rm d}Z+U\,p\bar{H}\,{\rm
d}\bar{Z}\right|^{2}+2\,{\rm d}U{\rm d}V-2\epsilon\,{\rm d}U^{2}\,,$ (9)

where $H$ is the Schwarzian derivative of the function $h$ defined as

$H(h(Z))\equiv\frac{1}{2}\left[\frac{h^{\prime\prime\prime}}{h^{\prime}}-\frac{3}{2}\left(\frac{h^{\prime\prime}}{h^{\prime}}\right)^{2}\right].$
(10)

Although the metric (9) still represents the flat space, the non-triviality of
function $H$ leads to the topological defects primarily induced by the choice
of its generator $h(Z)$.

Finally, to construct the expanding impulsive gravitational wave on a flat
background one has to cut the Minkowski spacetime along the null cone and then
re-attached the two half-spaces with an appropriate warp, see figure 1. Using
the coordinates of (6) and taking the null cone ${U=0}$, corresponding to the
expanding sphere ${t^{2}=x^{2}+y^{2}+z^{2}}$, the half-spaces ${{\cal M}^{-}}$
(with ${U\leq 0}$) and ${{\cal M}^{+}}$ (with ${U\geq 0}$) have to be
identified across the null hypersurface ${\cal N}$ as

$\left[Z,\,\bar{Z},\,V,\,U=0_{-}\right]_{{}_{{\cal
M}^{-}}}\equiv\left[h(Z),\,\bar{h}({\bar{Z}}),\,\frac{1+\epsilon\,h\bar{h}}{1+\epsilon
Z\bar{Z}}\frac{V}{|h^{\prime}|},\,U=0_{+}\right]_{{}_{{\cal M}^{+}}}\,.$ (11)

Figure 1: Two parts ${\cal M}^{-}$ and ${\cal M}^{+}$ of the flat Minkowski
background are re-attached along a future-oriented null cone ${\cal N}$ with a
suitable warp given by the junction Penrose conditions (11). Such construction
results in the presence of expanding spherical gravitational impulse located
on ${\cal N}$ whose specific nature is encoded in the (holomorphic) function
${h(Z)}$ identifying points of ${\cal M}^{-}$ and ${\cal M}^{+}$ across ${\cal
N}$.

The Penrose junction conditions directly correspond to the evaluation of
transformations (4) and (7), respectively, on the impulse ${U=0}$, i.e.,
${{\cal U}{\cal V}-\eta\bar{\eta}=0}$. The global _continuous_ line element
can be then written as a combination of (6) and (9), namely

${\rm d}s^{2}=2\left|\frac{V}{p}\,{\rm d}Z+U\Theta(U)\,p\bar{H}\,{\rm
d}\bar{Z}\right|^{2}+2\,{\rm d}U{\rm d}V-2\epsilon\,{\rm d}U^{2}\,,$ (12)

where the product of $U$ and Heaviside step $\Theta(U)$ represents the
continuous kink function which is typically denoted as ${{U_{+}}\equiv
U\Theta(U)}$. The metric (12) then solves the vacuum Einstein field equations
everywhere except at the singular impulse origin (${U=0=V}$), and possible
poles of ${p^{2}H}$ as can be inferred from the curvature invariants, see e.g.
[23].

##  3 Geometry of expanding impulses

In this section, we will discuss the specific effects of the warp function
${h(Z)}$ representing the source of an expanding impulse. Concerning the
Schwarzian derivative (10), the wave-like nature of $h(Z)$ can be
distinguished. In particular, the general Möbius transformation of the complex
plane to itself,

$h(Z):\quad Z\ \mapsto\ \frac{a\,Z+b}{c\,Z+d}\,,$ (13)

leaves $H(Z)$ unchanged and corresponds to the simple Lorentz transformation.
On the other hand, going beyond the linear fractional transformation brings
the non-trivial Schwarzian derivative. In our case, it will be interpreted as
a topological defect related to the presence of cosmic string. To
geometrically describe these effects in the flat Minkowski space, it is useful
to employ stereographic projection, see e.g. [24]. The mapping ${Z\mapsto
h(Z)}$ in the complex Argand plane corresponds to the geometric
identifications of points $P^{-}$ and $P^{+}$ on a Riemann sphere, see figure
2.

Figure 2: Stereographic projection corresponds to the mutual identification
between points of the Riemann sphere and complex Argand plane. In the case of
expanding impulses, it gives direct geometric interpretation to the
holomorphic mapping ${h(Z)}$ and the Penrose junction conditions (11).

The Riemann sphere can be further identified within the background Cartesian
coordinates. In particular, taking the continuous metric (12) and the
coordinate transformations leading to the half-spaces _in front of_ ($U\geq
0$) and _behind_ (${U\leq 0}$) the impulse, respectively, we find that the
impulsive surface ${U=0}$ is a sphere ${(t_{\rm i}^{\pm})^{2}=(x_{\rm
i}^{\pm})^{2}+(y_{\rm i}^{\pm})^{2}+(z_{\rm i}^{\pm})^{2}}$, where the $\pm$
sign corresponds to a specific half-space from which $\mathcal{N}$ is
approached, and the index $\,{}_{\rm i}$ indicates values obtained on the
impulse $\mathcal{N}$, i.e., on ${U=0}$. Then, we get

$Z_{\rm i}=\frac{\eta_{\rm i}^{-}}{{\cal V}_{\rm i}^{-}}=\frac{x_{\rm
i}^{-}+{i}y_{\rm i}^{-}}{t_{\rm i}^{-}+z_{\rm i}^{-}}\,,\qquad h(Z_{\rm
i})=\frac{\eta_{\rm i}^{+}}{{\cal V}_{\rm i}^{+}}=\frac{x_{\rm
i}^{+}+{i}y_{\rm i}^{+}}{t_{\rm i}^{+}+z_{\rm i}^{+}}\,.$ (14)

Inversely, we can write

$\frac{x_{\rm i}^{-}}{t_{\rm i}^{-}}=\frac{Z_{\rm i}+\bar{Z}_{\rm
i}}{1+|Z_{\rm i}|^{2}}\,,\qquad\frac{y_{\rm i}^{-}}{t_{\rm
i}^{-}}=-{i}\frac{Z_{\rm i}-\bar{Z}_{\rm i}}{1+|Z_{\rm
i}|^{2}}\,,\qquad\frac{z_{\rm i}^{-}}{t_{\rm i}^{-}}=\frac{1-|Z_{\rm
i}|^{2}}{1+|Z_{\rm i}|^{2}}\,,$ (15) $\frac{x_{\rm i}^{+}}{t_{\rm
i}^{+}}=\frac{h+\bar{h}}{1+|h|^{2}}\,,\qquad\frac{y_{\rm i}^{+}}{t_{\rm
i}^{+}}=-{i}\frac{h-\bar{h}}{1+|h|^{2}}\,,\qquad\frac{z_{\rm i}^{+}}{t_{\rm
i}^{+}}=\frac{1-|h|^{2}}{1+|h|^{2}}\,,$ (16)

where the function ${h}$ is evaluated at $Z_{\rm i}$. The expressions (15) and
(16) can thus be understood as the stereographic identification between points
in the complex plane and their images on a unit Riemann sphere representing
the re-scaled impulsive surface. Subsequently, in terms of such a unit sphere
endowed with the Cartesian axes, one can interpret the effects of mapping
${h(Z)}$.

In particular, the construction of explicit form of the function $h(Z)$ can be
decomposed into operations representing either pure Lorentz transformations of
the form (13) or mappings inducing non-trivial Schwartzian derivative (10).
Here, let us define elementary operations which will be sequentially applied
within the following discussion:

  * •

the natural starting point is an _identical_ mapping111The subscript
$\,{}_{j}$ of $h_{j}(Z)$ identifies particular step in a sequence of the final
${h(Z)}$ construction., i.e.,

$h_{0}(Z)=Z\,.$ (17)

  * •

the spatial rotations ${\mathcal{R}_{\\{\varphi,\vartheta,\psi\\}}}$
parameterized by the Euler angles ${\\{\varphi,\,\vartheta,\,\psi\\}}$ lead to

$h_{j+1}(Z)=\mathcal{R}_{\\{\varphi,\vartheta,\psi\\}}h_{j}(Z)=e^{{i}\varphi}\frac{-\sin(\vartheta/2)+e^{{i}\psi}\cos(\vartheta/2)\,h_{j}}{\cos(\vartheta/2)+e^{{i}\psi}\sin(\vartheta/2)\,h_{j}}\,,$
(18)

see the left part of figure 3 representing the rotated Riemann sphere,

  * •

the Lorentz boost in the direction of $z$-axis ${\mathcal{B}_{\\{w\\}}}$,
parameterized by the value $w$, gives

$h_{j+1}(Z)=\mathcal{B}_{\\{w\\}}h_{j}(Z)=w\,h_{j}(Z)\,,$ (19)

see the middle plots in figure 3 for a pure boost and its combination with a
rotation (18),

  * •

finally, the simplest string-like structure can be constructed by “cutting”
out the wedge $2\pi\delta$ around $z$-axis, represented by the action of
${\mathcal{S}_{\\{\delta\\}}}$, namely

$h_{j+1}(Z)=\mathcal{S}_{\\{\delta\\}}h_{j}(Z)=[h_{j}(Z)]^{1-\delta}\,,$ (20)

where the wedge is missing symmetrically in the direction of negative
$x$-axis, see the last example in figure 3. To be more precise, this operation
does not only remove a given angle. However, the spherical surface is cut
along the ${z-x}$ plane (for the negative $x$ values) and then the angle
${2\pi\delta}$ is opened, while the surface is “compressed” (as an accordion
or a paper lantern), which may affect already existing defects as we show
later.

These operations have to be understood as the active transformations of the
sphere, while the coordinate system and axes are kept fixed. In general, the
Euler angles are arbitrary and it thus seems to be possible to cut out an
arbitrary number of strings along different axes, which are moving with
different velocities. We will show explicit examples below. Finally, the
elementary operations (17)–(20) could be supplemented with other operations,
e.g., boosts in the $x$ and $y$ directions. However, these additional
operations can be simply understood as their compositions.

Figure 3: Visualisation of specific elementary operations acting on the whole
Riemann sphere. Their combinations enter the construction of a particular form
of the function $h(Z)$ that, subsequently, encodes the geometric nature of a
given expanding impulsive wave described by (12).

##  4 Interaction of geodesics with the expanding impulses

We aim to analyze properties of particular expanding impulses prescribed by
explicit choices of the generating function ${h(Z)}$, see sections 5 and 6.
Such a discussion is closely related to the geodesic motion of test observers
affected by interaction with the gravitational impulse. However, due to the
presence of a kink function in the continuous metric (12), the distributional
terms appear in the geodesic equation and its analysis becomes more tricky.
The particular case of ${Z=\hbox{const}}$ geodesics was studied in [25].
Assuming the ${C^{1}}$-geodesics, the refraction formulas for their
interaction with a generic impulse were derived in [26]. Subsequently, the
_existence_ and _global uniqueness_ of such $C^{1}$-geodesics crossing
expanding impulse, propagating on all constant curvature backgrounds, were
rigorously proved in [27, 23] using the Filippov solution concept [28, 29].
Such refraction formulas connect initial data ${\\{{\cal U}_{\rm
i}^{+},\,{\cal V}_{\rm i}^{+},\,\eta_{\rm i}^{+},\,\dot{{\cal U}}_{\rm
i}^{+},\,\dot{{\cal V}}_{\rm i}^{+},\,\dot{\eta}_{\rm i}^{+}\\}}$ and
${\\{{\cal U}_{\rm i}^{-},\,{\cal V}_{\rm i}^{-},\,\eta_{\rm
i}^{-},\,\dot{{\cal U}}_{\rm i}^{-},\,\dot{{\cal V}}_{\rm
i}^{-},\,\dot{\eta}_{\rm i}^{-}\\}}$ for the straight lines parameterized by
$\tau$, namely

$\displaystyle{\cal U}^{\pm}=$ $\displaystyle\dot{{\cal U}}_{\rm
i}^{\pm}\,\tau+{\cal U}_{\rm i}^{\pm}\,,$ $\displaystyle{\cal V}^{\pm}=$
$\displaystyle\dot{{\cal V}}_{\rm i}^{\pm}\,\tau+{\cal V}_{\rm i}^{\pm}\,,$
(21) $\displaystyle\eta^{\pm}=$ $\displaystyle\dot{\eta}_{\rm
i}^{\pm}\,\tau+\eta_{\rm i}^{\pm}\,,$

i.e., geodesics in the Minkowski half-spaces ${\cal M}^{+}$ and ${\cal
M}^{-}$, starting/ending on the impulse ${\cal N}$ at ${\tau=0}$, see (22) and
(23) below.

Here, let us summarize the main result of [26, 23] important for our further
discussion. The explicit $C^{1}$-matching of geodesics crossing the impulse
can be expressed in the form of the refraction formulas encoding the shift of
positions and change of the velocities with respect to the fiducial
interpretative background. These are derived starting from the fact that the
geodesics in coordinates (12) are unique $C^{1}$-lines across the impulsive
wavefront $\mathcal{N}$ given by ${U=0}$, i.e., components of position and
velocity evaluated on the impulsive boundary $\mathcal{N}$ (denoted by the
subscript i) are the same irrespectively whether $\mathcal{N}$ is approached
from the region ${\cal M}^{+}$ with ${U\geq 0}$ (denoted by the superscript +)
or from the complementary half-space ${\cal M}^{-}$ with ${U\leq 0}$ (denoted
by the superscript -). However, to observe the influence of the impulse on
test particles, it is natural to employ the fiducial background coordinates
(3) for ${\cal M}^{+}$ and ${\cal M}^{-}$, respectively. With respect to the
background space, the global geodesics do not cross the impulse continuously
and the effects of impulse on their motion become explicit. In particular,
evaluation of the transformations (7) and (4) on ${U=0}$, and elimination of
the continuous coordinates, gives the position shift,

${\cal U}_{\rm i}^{-}=|h^{\prime}|\frac{|Z_{\rm i}|^{2}}{|h|^{2}}\,{\cal
U}_{\rm i}^{+}\,,\qquad{\cal V}_{\rm i}^{-}=|h^{\prime}|\,{\cal V}_{\rm
i}^{+}\,,\qquad\eta_{\rm i}^{-}=|h^{\prime}|\frac{Z_{\rm i}}{h}\,\eta_{\rm
i}^{+}\,,$ (22)

while the same procedure for derivatives of (7) and (4) leads to the
refraction of the velocities,

$\displaystyle\dot{{\cal U}}_{\rm i}^{-}$ $\displaystyle=a_{{}_{\cal
U}}\dot{{\cal U}}_{\rm i}^{+}+b_{{}_{\cal U}}\dot{{\cal V}}_{\rm
i}^{+}+\bar{c}_{{}_{\cal U}}\dot{\eta}_{\rm i}^{+}+c_{{}_{\cal
U}}\dot{\bar{\eta}}_{\rm i}^{+}\,,$ $\displaystyle\dot{{\cal V}}_{\rm i}^{-}$
$\displaystyle=a_{{}_{\cal V}}\dot{{\cal U}}_{\rm i}^{+}+b_{{}_{\cal
V}}\dot{{\cal V}}_{\rm i}^{+}+\bar{c}_{{}_{\cal V}}\dot{\eta}_{\rm
i}^{+}+c_{{}_{\cal V}}\dot{\bar{\eta}}_{\rm i}^{+}\,,$ (23)
$\displaystyle\dot{\eta}_{\rm i}^{-}$ $\displaystyle=a_{\eta}\,\dot{{\cal
U}}_{\rm i}^{+}+b_{\eta}\dot{{\cal V}}_{\rm
i}^{+}+\bar{c}_{\eta}\,\dot{\eta}_{\rm i}^{+}+c_{\eta}\,\dot{\bar{\eta}}_{\rm
i}^{+}\,,$

where the coefficients are constants evaluated on $\mathcal{N}$ with ${Z_{\rm
i}}$ obtained via ${h(Z_{\rm i})=\frac{\eta_{\rm i}^{+}}{{\cal V}_{\rm
i}^{+}}}$, namely

$\displaystyle a_{{}_{\cal U}}$
$\displaystyle=\frac{1}{|h^{\prime}|}\left|1+\frac{Z_{\rm
i}}{2}\frac{h^{\prime\prime}}{h^{\prime}}\right|^{2}\,,$ (24) $\displaystyle
a_{{}_{\cal V}}$
$\displaystyle=\frac{1}{4|h^{\prime}|}\left|\frac{h^{\prime\prime}}{h^{\prime}}\right|^{2}\,,$
(25) $\displaystyle a_{\eta}$
$\displaystyle=\frac{1}{2|h^{\prime}|}\bigg{(}1+\frac{Z_{\rm
i}}{2}\frac{h^{\prime\prime}}{h^{\prime}}\bigg{)}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}\,,$
(26) $\displaystyle b_{{}_{\cal U}}$
$\displaystyle=\frac{|h|^{2}}{|h^{\prime}|}\left|1+\frac{Z_{\rm
i}}{2}\bigg{(}\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\bigg{)}\right|^{2}\,,$
(27) $\displaystyle b_{{}_{\cal V}}$
$\displaystyle=\frac{|h|^{2}}{4|h^{\prime}|}\left|\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right|^{2}\,,$
(28) $\displaystyle b_{\eta}$
$\displaystyle=\frac{|h|^{2}}{2|h^{\prime}|}\bigg{[}1+\frac{Z_{\rm
i}}{2}\bigg{(}\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\bigg{)}\bigg{]}\bigg{(}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}-2\frac{\bar{h}^{\prime}}{\bar{h}}\bigg{)}\,,$
(29) $\displaystyle c_{{}_{\cal U}}$
$\displaystyle=-\frac{h}{|h^{\prime}|}\bigg{[}1+\frac{Z_{\rm
i}}{2}\left(\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right)\bigg{]}\bigg{[}1+\frac{\bar{Z}_{\rm
i}}{2}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}\bigg{]}\,,$ (30)
$\displaystyle c_{{}_{\cal V}}$
$\displaystyle=-\frac{h}{4|h^{\prime}|}\left(\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right)\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}\,,$
(31) $\displaystyle c_{\eta}$
$\displaystyle=-\frac{h}{2|h^{\prime}|}\bigg{[}1+\frac{Z_{\rm
i}}{2}\left(\frac{h^{\prime\prime}}{h^{\prime}}-2\frac{h^{\prime}}{h}\right)\bigg{]}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}\,,$
(32) $\displaystyle\bar{c}_{\eta}$
$\displaystyle=-\frac{\bar{h}}{2|h^{\prime}|}\bigg{(}1+\frac{Z_{\rm
i}}{2}\frac{h^{\prime\prime}}{h^{\prime}}\bigg{)}\bigg{(}\frac{\bar{h}^{\prime\prime}}{\bar{h}^{\prime}}-2\frac{\bar{h}^{\prime}}{\bar{h}}\bigg{)}\,.$
(33)

Notice that ${\bar{c}_{{}_{\cal V}}=\overline{c_{{}_{\cal V}}}}$,
${\bar{c}_{{}_{\cal U}}=\overline{c_{{}_{\cal U}}}}$. As shown in [23], the
normalization of velocity is preserved across the impulse, i.e.,
${\dot{\eta}_{\rm i}^{-}\dot{\bar{\eta}}_{\rm i}^{-}-\dot{{\cal U}}_{\rm
i}^{-}\dot{{\cal V}}_{\rm i}^{-}=\dot{\eta}_{\rm i}^{+}\dot{\bar{\eta}}_{\rm
i}^{+}-\dot{{\cal U}}_{\rm i}^{+}\dot{{\cal V}}_{\rm i}^{+}}$. The above
(local) expressions do not depend on the Gaussian curvature $\epsilon$,
however, to construct a global picture the parameter $\epsilon$ encoding the
spacetime foliation has to be considered, see [13, 23]. The refraction
formulas become identical in the trivial case ${h(Z)=Z}$ that implies ${H=0}$
and lacks the impulse. However, one should be careful in the case of non-
trivial $h(Z)$ and still trivial ${H=0}$, where the above expressions identify
two Minkowski half-spaces via Möbius transformation (13), see section 6. This
ambiguity in mutual background Cartesian coordinate identification on both
sides of the impulse, non-physically affecting in the above refraction
formulas, arises from the absence of global Cartesian-like impulsive metric
for the expanding waves in contrast to the non-expanding case, see e.g. [30].

To identify of the pure wave action on the test observers, a specific choice
of the initial data ${\\{{\cal U}_{\rm i}^{+},\,{\cal V}_{\rm
i}^{+},\,\eta_{\rm i}^{+},\,\dot{{\cal U}}_{\rm i}^{+},\,\dot{{\cal V}}_{\rm
i}^{+},\,\dot{\eta}_{\rm i}^{+}\\}}$ has to be made which enters the above
expressions. Due to the time shift given by a combination of (22) and (2), we
can naturally consider either a swarm of test particles which is hit by the
impulsive wavefront ${\cal{N}}$ at the same _constant coordinate time
$t^{+}_{\rm i}$_ in $\cal{M}^{+}$, or vice versa, which appears in the region
$\cal{M}^{-}$ simultaneously at _constant coordinate time $t^{-}_{\rm i}$_.
The second possibility can be also understood as the case of impulse passing
through the continuous dust-like distribution of particles where we observe
one fixed emergence slice given by constant time $t^{-}_{\rm i}$. To emphasize
the geometric effect of the particular impulse realization encoded in the
mapping ${h(Z)}$, we will assume test particles at rest in $\cal{M}^{+}$ (the
rest is defined with respect to the background Minkowskian coordinates, i.e.,
${\dot{{\cal U}}_{\rm i}^{+}=\dot{{\cal V}}_{\rm i}^{+}=\dot{\eta}_{\rm
i}^{+}=0}$), which are spherically distributed with radius ${t^{\pm}_{\rm
i}=\hbox{const}}$. A schematic visualization of the above cases (in the
simplest one-string situation) is given in figure 4. The explicit initial data
constraints are summarized in the following subsections.

Figure 4: A swarm of particles at rest in $\cal{M}^{+}$ hit by the typical
expanding impulse (null cone $\cal{N}$) as seen from the half-space
$\cal{M}^{+}$ (lower part of the schema) and its emergence at null cone
$\cal{N}$ as seen from the half-space $\cal{M}^{-}$ (upper part of the
schema). On the left, the test particles interact with the impulse at
${t^{+}_{\rm i}=\hbox{const}}$ while the right part is the case with
${t^{-}_{\rm i}=\hbox{const}}$. In this particular case, it can be seen that
the particles attain a speed in the radial direction being ‘attracted’ to the
axis at the time they leave the wave which causes a formation of caustics
later.

###  4.1 Spheres in front of the wave

A Cartesian sphere (formed by test particles) of a constant radius $t_{\rm
i}^{+}$ in front of the wave corresponds to the particle’s displacement on a
particular cut of the null impulsive cone ${U=0}$ encoded in specific values
of the global coordinate $V_{\rm i}$ at the instant of interaction. In
particular, evaluation of the transformation (7) with (2) on the null cone
${U=0}$ gives the constraint

$t_{\rm i}^{+}=\frac{V_{\rm
i}}{\sqrt{2}|h^{\prime}|}\,\frac{1+|h|^{2}}{1+\epsilon|Z_{\rm
i}|^{2}}=\hbox{const}\,,$ (34)

which, combined with the initial assumption ${t_{\rm i}^{+}=\hbox{const}}$,
fixes the values $V_{\rm i}$ for specific choices of complex plane positions
$h(Z_{\rm i})$ (and inversely $Z_{\rm i}$). The $h(Z_{\rm i})$ value is
related to the Cartesian parameterisation via stereographic projection (16),

$x_{\rm i}^{+}=\frac{h+\bar{h}}{1+|h|^{2}}\,t_{\rm i}^{+}\,,\qquad y_{\rm
i}^{+}=-{i}\frac{h-\bar{h}}{1+|h|^{2}}\,t_{\rm i}^{+}\,,\qquad z_{\rm
i}^{+}=\frac{1-|h|^{2}}{1+|h|^{2}}\,t_{\rm i}^{+}\,,$ (35)

In the case of non-trivial $h(Z)$ with respect to its Schwartzian derivative
(10), the Cartesian image will suffer for different cutouts. At this moment,
it is useful to adopt natural parameterisation,

$Z_{\rm i}=\tan\frac{\theta}{2}\,e^{{i}\phi}\,,$ (36)

so that ${(x_{\rm i}^{\pm},\,y_{\rm i}^{\pm},\,z_{\rm i}^{\pm})}$ are
functions of ${\theta\in\langle 0,\,\pi)}$ and ${\phi\in\langle-\pi,\pi)}$.
Subsequently, the deformation of such a test ball is explicitly described by
the shift of positions (22), or directly in terms of the continuous
coordinates using stereographic projection (15) of null cone $\cal{N}$ as
viewed from the region behind the impulse. The scale in terms of the original
constant radius $t_{\rm i}^{+}$ is given by

$t_{\rm i}^{-}=\frac{V_{\rm i}}{\sqrt{2}}\,\frac{1+|Z_{\rm
i}|^{2}}{1+\epsilon|Z_{\rm i}|^{2}}=|h^{\prime}|\,\frac{1+|Z_{\rm
i}|^{2}}{1+|h|^{2}}\,t_{\rm i}^{+}\,,$ (37)

while the deformed surface in front of the wave can be plotted as

$\displaystyle x_{\rm i}^{-}=|h^{\prime}|\frac{Z_{\rm i}+\bar{Z}_{\rm
i}}{1+|h|^{2}}\,t_{\rm i}^{+}\,,\qquad y_{\rm
i}^{-}=-{i}|h^{\prime}|\frac{Z_{\rm i}-\bar{Z}_{\rm i}}{1+|h|^{2}}\,t_{\rm
i}^{+}\,,\qquad z_{\rm i}^{-}=|h^{\prime}|\frac{1-|Z_{\rm
i}|^{2}}{1+|h|^{2}}\,t_{\rm i}^{+}\,$ (38)

which still satisfy ${({x_{\rm i}^{-}})^{2}+({y_{\rm i}^{-}})^{2}+({z_{\rm
i}^{-}})^{2}=({t_{\rm i}^{-}})^{2}}$, however, ${t_{\rm i}^{-}}$ is no more a
constant.

###  4.2 Spheres behind the wave

The second very natural choice of the initial data is such that the test
particles form a sphere at a given constant time $t_{\rm i}^{-}$ behind the
wave, i.e., in the region without any strings. The null cone cut is fixed in
terms of values $V_{\rm i}$ of the global coordinate $V$ by the condition

$t_{\rm i}^{-}=\frac{V_{\rm i}}{\sqrt{2}}\,\frac{1+|Z_{\rm
i}|^{2}}{1+\epsilon|Z_{\rm i}|^{2}}=\hbox{const}\,,$ (39)

and the Cartesian positions on a sphere are related to the values $Z_{\rm i}$
by (15),

$x_{\rm i}^{-}=\frac{Z_{\rm i}+\bar{Z}_{\rm i}}{1+|Z_{\rm i}|^{2}}\,t_{\rm
i}^{-}\,,\qquad y_{\rm i}^{-}=-{i}\frac{Z_{\rm i}-\bar{Z}_{\rm i}}{1+|Z_{\rm
i}|^{2}}\,t_{\rm i}^{-}\,,\qquad z_{\rm i}^{-}=\frac{1-|Z_{\rm
i}|^{2}}{1+|Z_{\rm i}|^{2}}\,t_{\rm i}^{-}\,.$ (40)

Viewed from the region ${U>0}$, this corresponds to the deformed initial
displacement given by (16) satisfying ${({x_{\rm i}^{+}})^{2}+({y_{\rm
i}^{+}})^{2}+({z_{\rm i}^{+}})^{2}=({t_{\rm i}^{+}})^{2}}$ with non-constant
scaling

$t_{\rm i}^{+}=\frac{V_{\rm
i}}{\sqrt{2}|h^{\prime}|}\,\frac{1+|h|^{2}}{1+\epsilon|Z_{\rm
i}|^{2}}=\frac{1}{|h^{\prime}|}\,\frac{1+|h|^{2}}{1+|Z_{\rm i}|^{2}}\,t_{\rm
i}^{-}\,,$ (41)

and Cartesian positions given by

$\displaystyle x_{\rm
i}^{+}=\frac{1}{|h^{\prime}|}\,\frac{h+\bar{h}}{1+|Z_{\rm i}|^{2}}\,t_{\rm
i}^{-}\,,\qquad y_{\rm
i}^{+}=-\frac{{i}}{|h^{\prime}|}\,\frac{h-\bar{h}}{1+|Z_{\rm i}|^{2}}\,t_{\rm
i}^{-}\,,\qquad z_{\rm
i}^{+}=\frac{1}{|h^{\prime}|}\,\frac{1-|h|^{2}}{1+|Z_{\rm i}|^{2}}\,t_{\rm
i}^{-}\,.$ (42)

The resulting deformed displacement in $\cal{M}^{+}$ contains defects given by
a particular choice of ${h(Z)}$.

###  4.3 Visualisations and location of strings

Although the deformations (35) and (42), respectively, do not depend on the
value of the Gaussian curvature $\epsilon$, it is natural to assume the choice
$\epsilon=1$. This way, the cuts of the null cone ${U=0}$, ${V_{\rm
i}=\mbox{const}}$, parameterised by the remaining global coordinate values
$Z_{\rm i}$, are manifestly spheres behind the wave. The radius is $t_{\rm
i}^{-}$ proportional only to ${V_{\rm i}}$, namely ${(x_{\rm
i}^{-})^{2}+(y_{\rm i}^{-})^{2}+(z_{\rm i}^{-})^{2}=\frac{1}{2}V_{\rm
i}^{2}}$, see (39). This assumption allows for simpler visualisations in terms
of the schematic cut and paste figure 1, while other choices of $\epsilon$
lead to more complicated sections of the null cone.

To understand the resulting geometry it is natural to investigate where the
string ends are attached to the null cone. Their location is related to the
scaling factors ${t_{\rm i}^{\pm}}$. In particular, in the case of initial
configuration representing a constant sphere in front of the wave, see section
4.1, these points correspond to the divergence of $t_{\rm i}^{-}$ given by
(37), while for the constant sphere behind the impulse, see section 4.2, we
are inversely looking for zeros of $t_{\rm i}^{+}$ given by (41). They thus
represent extremes of the ‘radial’ distance of the deformed ‘spherical’
surface.

Determination of the string positions allows to distinguish qualitatively
different physical situations. For example, we can arrive at the same string
configurations by applying the elementary operations in a different order.
However, this typically leads to the different forms of the generating
function ${h(Z)}$. Analogously, the rotations will also not change the mutual
string configurations, but the functions ${h(Z)}$ will differ. This can be
solved by finding the string ends. Therefore, we map the constant sphere-like
configuration to obtain its deformed image behind the wave and find points of
divergence on such a surface. Then, we can measure mutual angles between all
possible pairs of such points for particular functions ${h(Z)}$. Two
configurations are identical if angles agree in both cases.

Finally, keep in mind that we infer all the information about the strings from
the behaviour on the null cone or properties of its spherical cuts. To proceed
more explicitly, the coordinates (9) in front of the wave, where the strings
are present (and, possibly, moving), should be employed in the whole space.
However, they are extremely complicated and their analytic inversion to the
Cartesian coordinate system, where the topological defects can be directly
interpreted, seems to be impossible since it requires finding the inverse of
(7) with (8).

##  5 One string geometries

Although the one-string case has been frequently studied, we would like to
provide its description to point out the important technical aspects of the
construction that will subsequently appear in more involved two-string cases.
The simplest situation of one string located along the $z$ axis, and therefore
inducing the deficit angle ${2\pi\delta}$ around it, is given by the mapping

$h(Z)=Z^{1-\delta}\,.$ (43)

As we have already mentioned in section 4, test particles forming an initial
spherical shell in the region $\mathcal{M}^{+}$, and interacting with the
impulse simultaneously at a constant time ${t_{\rm i}^{+}}$, will be displaced
in both space and time directions of $\mathcal{M}^{-}$, and vice versa. The
surfaces representing such initial conditions are visualized in figure 5. In
the case of sphere with a constant radius $t_{\rm i}^{+}$ in
$\mathcal{M}^{+}$, we get

$t_{\rm i}^{-}\rightarrow\infty\,,\qquad z_{\rm i}^{-}\rightarrow\pm\infty\,,$
(44)

along the axis, corresponding to the (37) diverges at the places
${Z=\\{0,\,\infty\\}}$ where the string ends are attached to the null cone,
see the left part of the figure 5. This divergence is in reverse translated
into the shape of the initial time-slice in the right part of the figure 5
leading to particles emerging simultaneously at constant time $t_{\rm i}^{-}$.

Figure 5: In part (a$+$), the test particles are chosen to form a sphere of a
constant radius $t_{\rm i}^{+}$ and deficit angle ${2\pi\delta}$ along $z$
axis in the region $\mathcal{M}^{+}$ (in front of the wave). The ‘sphere’ is
embedded in the Euclidean space with Cartesian coordinates $(x_{\rm
i}^{+},y_{\rm i}^{+},z_{\rm i}^{+})$. The thin dashed lines represent lines of
constant $\theta$ and $\phi$ of the stereographic projection (36). This sphere
is distorted behind the wave in $(x_{\rm i}^{-},y_{\rm i}^{-},z_{\rm i}^{-})$
coordinates as visualised in part (a$-$). The particles appear _at different
times_ and different locations, see (22). The time shift is depicted by a
colour change and the lines of constant $t_{\rm i}^{-}$ are shown as thick
blue lines. There is a singularity along the $z^{-}$ (the spikes reach
infinity) while the deficit angle disappears. The part (b$+$) represents such
initial conditions in $\mathcal{M}^{+}$ that the particles emerge as a real
geometric sphere in $\mathcal{M}^{-}$ with radius $t_{\rm i}^{-}=\hbox{const}$
and coordinates $(x_{\rm i}^{-},y_{\rm i}^{-},z_{\rm i}^{-})$, see part
(b$-$).

Moreover, to visualize the effect of the Penrose junction condition (11) as a
null cone mapping, which was schematically illustrated in figure 1, we plot
its explicit realization for the one string case (43), see figure 6.

Figure 6: The initial ${t_{\rm i}^{+}=\hbox{const}}$ sphere is cut by a plane
passing through the origin and including both string ends (a). Since the plane
does not contain the deficit angle, taking specific values $t_{\rm i}^{+}$
corresponds to the full circles of test particles (b) that are lying on the
null cone $\cal{N}$, where the time goes in the direction of the cone axis.
The function ${h(Z)}$ then maps the initial circles (b) along the null cone
$\mathcal{N}$ as in figure 1. In (c) we can see the resulting displacement of
test particles on $\cal{N}$ viewed from the $\mathcal{M}^{-}$ region. The
divergences correspond to the moving string ends related to the infinities of
(37).

The above one-string situation can be non-trivially extended by its boosting
in the perpendicular direction. The corresponding complex mapping ${h(Z)}$ can
be constructed as a series of string creation along $z$ axis, perpendicular
rotation, boost, and final backward rotation, namely

$h(Z)=\mathcal{R}_{\\{-\frac{\pi}{2},0,0\\}}\mathcal{B}_{\\{w\\}}\mathcal{R}_{\\{\frac{\pi}{2},0,0\\}}\mathcal{S}_{\\{\delta\\}}Z\,.$
(45)

It explicitly becomes

$h(Z)=-\frac{w-1-(w+1)Z^{1-\delta}}{w+1-(w-1)Z^{1-\delta}}\,.$ (46)

The boost-induced asymmetry reflected in the deformation of the natural static
spherical initial conditions, with $t_{\rm i}^{+}$ or $t_{\rm i}^{-}$ being
constant, respectively, is visualized in figure 7.

Figure 7: Description of the plot is the same as in figure 5. The differences
in the shape of the particles swarm interacting with the impulse are caused by
the additional string motion in the $x$ direction related to the form of
function $h(Z)$ given by (46) with the boost parameter ${w=0.6}$. In the
second picture in (a$-$) we can observe a sharp ‘fin’ with the time
equipotentials being non-smooth. That is the region through which the string
travels. In these visualisations, the relation between (a) and (b) subfigures
clearly corresponds to the multiplication of the radius by a scaling factor
${|h^{\prime}|\,\frac{1+|Z_{\rm i}|^{2}}{1+|h|^{2}}}$. The string ends are at
the north and south pole behind the wave regardless of its transversal motion.

The example of an explicit null cone mapping for the moving string (46) is
visualized in figure 8.

Figure 8: The planar cut of the initial swarm of particles is visualized in
analogy with figure 6. However, here in part (a), the plane is chosen to
explicitly contain the topological defect caused by the string (b) while the
string ends are out of the plane. The fin-like deformation in part (c)
corresponds to the string motion in the $x$ direction, however, it does not
contain divergences since the string ends are out of our planar cut.

Finally, as a comprehensive picture characterizing the dynamical structure of
boosted one-string impulsive spacetime with (46) we plot a sequence of
deformations of an initially static swarm of particles. Their initial space
and time displacement in the region $\mathcal{M}^{+}$ is chosen in such a way
that they all emerge simultaneously in the region $\mathcal{M}^{-}$ as a
sphere of constant radius ${t_{\rm i}^{-}}$, i.e, it is described by figure 7
(b${+}$). Subsequently, due to the non-trivial impulsive effect, the particles
start moving along geodesics (straight lines) in the Minkowski background (21)
with initial data given by (22) and (23), see figure 9. The vertical
deformation and particle acceleration along the $z$ axis are given by their
attraction by the moving ends of the snapped string. Near poles, where the
string ends are attached to the impulsive sphere, the magnitude of velocity
approaches the speed of light. The attractive effect of the string ends
results in the caustics formed by a mutual crossing of trajectories starting
on opposite sides of the initial sphere with respect to the $z$ axis.
Moreover, the horizontal asymmetry is induced by the boost of the string in
the $x$ direction. To emphasize the boost effect we plot the cut by ${z=0}$
plane in the left part of figure 9.

Figure 9: The effect of snapping boosted string (46) on free test particles
initially at rest in $\mathcal{M}^{+}$ with positions satisfying ${{x_{\rm
i}^{+}}^{2}+{y_{\rm i}^{+}}^{2}+{z_{\rm i}^{+}}^{2}={t_{\rm i}^{+}}^{2}}$ such
that particular choice of ${t_{\rm i}^{+}}$ leads to ${t_{\rm
i}^{-}=\hbox{const}}$, i.e, the data are prepared as in the bottom part of
figure 7 to emerge as a real sphere at given time ${t_{\rm i}^{-}}$. However,
the test observers gain nontrivial velocities. Their time evolution is plotted
in the left column. To visualize the asymmetry induced by the string motion in
the $x$ direction we show the cut by plane ${z=0}$, perpendicular to the
string, in the right column.

The above discussion shows the effects of a moving snapped cosmic string on
free test observers and the way how these effects can be understood and
visualized. Notice that the subcase without string boost was already studied
in [25, 26] considering also cosmological backgrounds. In the subsequent
section, we follow [17, 24] and analyze various impulses generated by a pair
of snapped strings in the same way as in the one-string case.

Naturally, in the case of boosted string one can be worried by the fact that
the cut out planes are not extrinsically flat, as can be seen in figure 7.
However, it can be analytically shown that following situations are
equivalent:

  * •

initially static particles and their interaction with impulse generated by a
string boosted by $w$ in the positive $x$ direction, i.e., ${h(Z)}$ given by
(46),

  * •

initially moving particles in the negative $x$ direction with the velocity
magnitude $v=\frac{1-w^{2}}{1+w^{2}}$ interacting with a static string
described by (43).

##  6 Two strings geometries

Here, we review and extend the results presented in [24] and discuss their
geometric properties via analysis of induced geodesic motion. Subsequently, an
alternative sequences of basic steps entering the two strings complex mapping
$h{(Z)}$ construction will be employed and particular differences in the
resulting motion identified.

###  6.1 Original results by Podolský and Griffiths

The possibility of expanding impulse generated by a colliding and snapping
pair of cosmic strings was originally anticipated by Nutku and Penrose in
[17]. Simultaneously, there were doubts that explicit realization of the
corresponding function $h{(Z)}$ is hard to find, however, its existence should
be guaranteed by the Riemann theorem. Surprisingly, a few years later Podolský
and Griffiths explicitly performed such a construction in [24]. Their simplest
non-boosted two-string formula reads

$h(Z)=\mathcal{S}_{\\{\varepsilon\\}}\mathcal{R}_{\\{\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2}\\}}\mathcal{S}_{\\{\delta\\}}=\left(\frac{iZ^{1-\delta}-1}{Z^{1-\delta}-i}\right)^{1-\varepsilon}\,.$
(47)

In terms of fixed Cartesian coordinates, it can be described in such a way
that the string parameterized by $\delta$ and placed along the $z$ axis is
rotated to take a place along the $y$ axis and the second string (encoded in
$\varepsilon$) is then created along the $z$ axis. However, the construction
should be more precisely understood in terms of the active Lorentz
transformations (13). This can be seen from the no string limit of (47).
Taking both parameters trivial, namely ${\delta=0=\varepsilon}$, does not
provide the identity, but

$h(Z)=\frac{iZ-1}{Z-i}\,,$ (48)

which is exactly the residual rotation
${\mathcal{R}_{\\{\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2}\\}}\,Z}$ employed
within construction of (47). The Schwarzian derivative (10) is vanishing and
the seemingly non-trivial position shift (22) directly shows the unphysical
rotation or the artificial identification of the background Cartesian frame.
However, to directly gain all relevant information about test particles
interacting with the impulse from the refraction formulas (22) and (23), it is
important and useful to remove such a coordinate discrepancy. Before we do
that, let us show the ultimate result of [24] adding a boost to the simplest
interaction of static strings (47), and therefore, interpreted as the
collision of strings that induces their snap and subsequent creation of the
impulse. The particular complex mapping reads

$h(Z)=w_{2}\frac{h_{c}^{1-\varepsilon}-1}{h_{c}^{1-\varepsilon}+1}$ (49)

with

$h_{c}(Z)=-\frac{(w_{1}-i)Z^{1-\delta}+(w_{1}+i)}{(w_{1}+i)Z^{1-\delta}+(w_{1}-i)}\,.$
(50)

Now, to improve the artificial coordinate effect in (47), and subsequently
also in (49), let us perform the following sequence of mappings, namely

$\displaystyle h_{0}(Z)$ $\displaystyle=Z\,,$ (51) $\displaystyle h_{1}(Z)$
$\displaystyle=\mathcal{S}_{\\{\delta\\}}\,h_{0}(Z)\,,$ (52) $\displaystyle
h_{2}(Z)$
$\displaystyle=\mathcal{B}_{\\{w_{2}\\}}\mathcal{R}_{\\{\frac{\pi}{2},\frac{\pi}{2},\pi\\}}\,h_{1}(Z)$
(53) $\displaystyle h_{3}(Z)$
$\displaystyle=\mathcal{B}_{\\{w_{3}\\}}\mathcal{R}_{\\{-\frac{\pi}{2},0,0\\}}\,h_{2}(Z)\,,$
(54) $\displaystyle h_{4}(Z)$
$\displaystyle=\mathcal{S}_{\\{\varepsilon\\}}\,h_{3}(Z)\,,$ (55)
$\displaystyle h_{5}(Z)$
$\displaystyle=\mathcal{B}_{\\{w_{5}\\}}\mathcal{R}_{\\{\frac{\pi}{2},\frac{\pi}{2},\frac{\pi}{2}\\}}\,h_{4}(Z)\,.$
(56)

These operations can be geometrically understood in terms of the Riemann
sphere, see figures 10 and 11 for the special cases, while the fully general
mapping takes the form

$h(Z)=iw_{5}\frac{-1+i\left(w_{3}\frac{(i+w_{2})Z^{\delta}-(i-w_{2})Z}{(i-w_{2})Z^{\delta}-(i+w_{2})Z}\right)^{1-\varepsilon}}{1+i\left(w_{3}\frac{(i+w_{2})Z^{\delta}-(i-w_{2})Z}{(i-w_{2})Z^{\delta}-(i+w_{2})Z}\right)^{1-\varepsilon}}\,,$
(57)

which naturally becomes identical for trivial boosts and deficit angles, i.e.,
${w_{j}=1}$ and ${\delta=0=\varepsilon}$. The simplest static interaction of
the strings is then described by choosing $w_{2}=w_{3}=w_{5}=1$ in (57) which
then corresponds to (47) with the artificial coordinate rotation removed,
namely

$h(Z)=i\frac{-1+i\left(-i\frac{1-iZ^{1-\delta}}{1+iZ^{1-\delta}}\right)^{1-\varepsilon}}{1+i\left(-i\frac{1-iZ^{1-\delta}}{1+iZ^{1-\delta}}\right)^{1-\varepsilon}}\,,$
(58)

while for a specific choice of boosts ${w_{j}}$ we may obtain an improved
version of (49). On a general level, the location of ends of the strings for
(57) is given by

$Z=\left\\{0,\infty,\left(\frac{i+w_{2}}{i-w_{2}}\right)^{\frac{1}{1-\delta}},\left(\frac{i-w_{2}}{i+w_{2}}\right)^{\frac{1}{1-\delta}}\right\\}\,.$
(59)

Here, let us also emphasize that the elementary operation
$\mathcal{S}_{\\{\varepsilon\\}}$, see (55), which creates the second string
placed along $z$-axis, inherently shifts the position in the $x$-direction of
the already existing string. In figure 10, this corresponds to the step
$h_{3}\,\rightarrow\,h_{4}$, where the sphere is distorted and the existing
string is effectively boosted in the $x$-direction. We can counterbalance this
effect by additional boost in the opposite direction, namely

$w_{2}=\tan\frac{\pi}{4(1-\varepsilon)}\,,$ (60)

and ${w_{3}=1=w_{5}}$, see figure 11 and changes in the sequence
$h_{2}\,\rightarrow\,h_{3}\,\rightarrow\,h_{4}$. The resulting nodal points
are then aligned with the axes and there is no transversal velocity (as could
be seen from their action on test particles). These situations are also
compared in figure 13.

This discussion is connected with the real location of the strings viewed from
the region behind the wave. Geometrically, after cutting out two wedges
(deficit angles) we are glueing the corresponding ‘lips’ back together,
however, it is done in a particular order. Here, the string governed by
parameter $\delta$ remains straight along the $z$ axis, while using (59) the
string parts governed by $\varepsilon$, and lying in the ${x}$-${y}$ plane,
can be identified by the polar angle $\phi_{s}$, namely

$\phi_{s}=\pm\frac{\pi-\arg\frac{w_{2}+i}{w_{2}-i}}{1-\delta}\,.$ (61)

In the case with (60) and ${w_{3}=1=w_{5}}$, i.e., the constants are set so
that there is no transversal motion of the strings, pieces of one string are
attached to the north and south pole, respectively, while the second string
ends form the mutual angle

$\Delta\phi=\left(2-\frac{1}{2(1-\delta)(1-\varepsilon)}\right)\pi\,,$ (62)

which is found as a direct application of the general formula (61).

Figure 10: A sequence of elementary mappings (51)–(56) employed within the
construction of the function (57) encoding interaction of two cosmic strings
and their snap. Here we plot the simplest case $w_{2}=w_{3}=w_{5}=1$ directly
improving (47). The function $h(Z)$ then enters the junction conditions (11)
of two half-spaces ‘in front of’ and ‘behind’ the impulsive wave.
Interestingly, the second string creation, step ${h_{3}\,\rightarrow\,h_{4}}$,
inherently induces an additional boost. This can be counterbalanced by a
specific choice of $w_{j}$, see figure 11. Figure 11: The same sequence of
elementary mappings as in the figure 10. However, here we decided to take
${w_{2}=\tan\frac{\pi}{4(1-\varepsilon)}}$ and ${w_{3}=1=w_{5}}$ so that the
additional inherent boost in ${h_{3}\,\rightarrow\,h_{4}}$ is compensated and
the strings are placed along the axes, compare $h_{4}$ steps.

Finally, based on the geodesic motion let us describe the effect of an impulse
generated by (58). We visualize its interaction with initially static test
particles prepared in $\mathcal{M}^{+}$ to emerge synchronously on a sphere in
$\mathcal{M}^{-}$, see figure 12 identifying the initial data and figure 14
depicting the overall deformation caused by the impulse. The shape evolution
corresponds to the straight motion on a flat background. However, there is a
non-trivial distribution of velocities given by (23). The map of velocity
magnitude,

$v^{-}=\sqrt{(v_{x}^{-})^{2}+(v_{y}^{-})^{2}+(v_{y}^{-})^{2}}$ (63)

with

$\left(v_{x}^{-},\,v_{y}^{-},\,v_{z}^{-}\right)-\equiv\left(\frac{\dot{x}_{\rm
i}^{-}}{\dot{t}_{\rm i}^{-}},\,\frac{\dot{y}_{\rm i}^{-}}{\dot{t}_{\rm
i}^{-}},\,\frac{\dot{z}_{\rm i}^{-}}{\dot{t}_{\rm i}^{-}}\right)\,,$ (64)

is plotted in figure 15, which again indicates the string ends since in their
neighbourhood the test particles approach the speed of light.

Figure 12: As in figure 5 two natural sets of static initial data are prepared
for the improved two-string case (57) with ${w_{3}=1=w_{5}}$ and $w_{2}$ given
by (60). The plot of the initial condition for (47) would differs by rotation
and a trivial choice $w_{2}=1$. Figure 13: In part (a) we plot the Riemann
sphere as in figure 11 $h_{5}$, i.e., corresponding to the mapping (57) with
${w_{2}=\tan\frac{\pi}{4(1-\epsilon)}}$ compensating the additional boost. The
cut by $x$-$y$ plane of the $t_{\rm i}^{+}=\mbox{const}$ initial condition
image in $\mathcal{M}^{-}$ is shown in part (b). The string ends follow the
red lines, whose mutual angle is given by (61). The same situations with a
generic boost ${w_{2}=0.4}$ are visualized in figures (c) and (d). In (a) and
(c) we plot the blue surfaces given by parameterization of the wedge edges. If
there is no transversal motion of the strings, the surfaces are planes with
vanishing extrinsic curvature, while in the case of boosted strings, the
situation becomes more complicated. Figure 14: In analogy with figure 9 we
plot a burst of test particles hit by the impulsive wave generated by the
improved mapping (58), which can be directly used in refraction formulas (22)
and (23) entering geodesics (21). This choice is visualized in figure 10. The
employed initial conditions represent static particles in front of the wave
displaced to emerge synchronously at constant time $t_{\rm i}^{-}$. The
typical values ${\delta=\frac{1}{4}}$ and ${\varepsilon=\frac{1}{8}}$ are
used. We observe the dragging of test particles in the directions of four
moving ends of the snapped string pairs. This is combined with inherent (non-
compensated) boost in the $x$ direction induced by the second string creation.
Figure 15: The distribution of velocity magnitudes (63) in the region behind
the impulse. The initially static particles in front of the impulse emerging
simultaneously behind the wave generated by (58) are considered with the
typical values ${\delta=\frac{1}{4}}$ and ${\varepsilon=\frac{1}{8}}$ of the
deficit parameters. The magnitude reaches the speed of light near the four
string ends. Here, we show the Kavrayskiy VII projection (left) and the
spatial scheme (right).

###  6.2 Peculiar form of the simple two-string mapping

In this part, let us show an alternative construction of the function $h(Z)$
in the case of a pair of cosmic strings, which should demonstrate the subtlety
of the construction process.

Let us start with the sequence of elementary steps

$\displaystyle h_{0}(Z)$ $\displaystyle=Z\,,$ (65) $\displaystyle h_{1}(Z)$
$\displaystyle=\mathcal{B}_{\\{w_{3}\\}}\mathcal{S}_{\\{\delta\\}}\,h_{0}(Z)\,,$
(66) $\displaystyle h_{2}(Z)$
$\displaystyle=\mathcal{R}_{\\{\frac{\pi}{2},0,0\\}}\,h_{1}(Z)\,,$ (67)
$\displaystyle h_{3}(Z)$
$\displaystyle=\mathcal{B}_{\\{w_{7}\\}}\mathcal{S}_{\\{\varepsilon\\}}\,h_{2}(Z)\,,$
(68) $\displaystyle h_{4}(Z)$
$\displaystyle=\mathcal{R}_{\\{-\frac{\pi}{2},0,0\\}}\,h_{3}(Z)\,,$ (69)

which leads to the mapping

$h(Z)=-\frac{1-w_{7}\left(\frac{1+w_{3}Z^{1-\delta}}{1-w_{3}Z^{1-\delta}}\right)^{1-\epsilon}}{1+w_{7}\left(\frac{1+w_{3}Z^{1-\delta}}{1-w_{3}Z^{1-\delta}}\right)^{1-\epsilon}}\,.$
(70)

Surprisingly, we identify only three string ends located at

$Z=\left\\{0,\,\infty,\,w_{3}^{-\frac{1}{1-\delta}}\right\\}\,.$ (71)

The sequence (65)–(65) thus cuts out the wedge $2\pi\delta$ along $z$ axis and
performs a boost, makes rotation about $y$ axis, applies another string-like
cut parameterized by $\varepsilon$ along the original $x$ axis with another
boost, and finally rotates backwards about unchanged $y$ axis. The simplified
construction is visualized in figure 16. The mapping then leads to identity
for ${\delta=0=\varepsilon}$ and ${w_{3}=1=w_{7}}$.

Figure 16: Particular steps of the construction of peculiar mapping (70)
simplified by assuming trivial boost parameters ${w_{3}=1=w_{7}}$.

Similarly to in the previous section and the case of mapping (57), the second
string creation induces an additional boost. Then the edges of Riemann sphere
cutouts, parameterized by $V$, represent a generic curved surface. The
corresponding extrinsic curvature vanishes and the surface becomes a plane if
the inherent boost is compensated by a suitable choice of the artificial
boosts in our construction. In particular, taking

$\displaystyle w_{3}$ $\displaystyle=1\,,$ $\displaystyle w_{7}$
$\displaystyle=\sqrt{\frac{1-\sin\frac{\epsilon\pi}{2}\left(\tan\frac{\delta\pi}{2}\right)^{-(1-\epsilon)}}{1-\sin\frac{\varepsilon\pi}{2}\left(\tan\frac{\delta\pi}{2}\right)^{+(1-\varepsilon)}}}$
(72)

guarantees that the strings do not move in the transverse direction at all.
This is illustrated in figure 18.

Based on the construction with two string-like operations, one would expect
that it describes a pair of perpendicular cosmic strings with four moving ends
after their snap. However, inspecting (71) related to the static initial data
choice in figure 17, we may conclude that there are only _three_ string
pieces. This is exactly the above-mentioned subtlety in the ${h(Z)}$
construction. One part of the first string disappears by the creation of the
second string (68) which removes the associated nodal point of the Riemann
sphere, see ${h_{2}\,\rightarrow\,h_{3}}$ in figure 16.

Figure 17: The static initial data choice for the simplified complex mapping
(70) with ${w_{3}=1=w_{7}}$. The infinities in the part (a$-$) indicate the
presence of only three moving string pieces. Figure 18:  In the case of
generic boosts, one surface swept by the Riemann sphere edge has non-vanishing
extrinsic curvature, while the second one is a trivial plane, see two views in
part (a). In the static case, given by the boost choice (72), all the surfaces
are planar. The pictures are reflectively symmetric with respect to the
${y=0}$ plane corresponding to the blue and green axis .

Finally, we can visualize the kinematic effect of the impulse generated by
(70) with ${w_{3}=1=w_{7}}$ on the motion of initially static test observers,
see figure 19 for the evolution picture and figure 20 for the velocity
magnitude distribution. Due to the absence of a string piece along the
negative $x$ axis and non-compensated inherent boost the resulting picture is
asymmetric. Without an appropriate global view, one may be confused. In
particular, restricting the analysis to just a quarter of the picture
including two perpendicular string pieces in the positive $x$ and $z$
directions, and taking the deficit parameters ${\delta=\varepsilon}$ one would
expect a symmetric picture in the presumed case of a pair of complete strings
with four ends. However, we observe an induced motion of test particles that
prefers the single-end direction.

Figure 19: Deformation of an initially static swarm of particles that emerge
instantly as a sphere of constant radius ${t_{\rm i}^{-}}$. The impulse is
generated by (70) with ${w_{3}=1=w_{7}}$ and string parameters
${\delta=\frac{1}{4}}$, ${\varepsilon=\frac{1}{8}}$. The asymmetry in the
resulting motion of test particles is caused by the absence of one semi-
infinite cosmic string along the negative $x$ axis and a boost induced by the
second string creation. Figure 20: In analogy with figure 15, the distribution
of velocity magnitudes (63) of initially static test particles, simultaneously
escaping the impulsive surface, is plotted. Here, the impulse is generated by
(70) with ${w_{3}=1=w_{7}}$ and ${\delta=\frac{1}{4}}$,
${\varepsilon=\frac{1}{8}}$. There are only three regions with a velocity
magnitude approaching the speed of light corresponding to three moving string
pieces, see the Kavrayskiy VII projection (left) and the spatial picture
(right).

###  6.3 Two parallel strings

As a last example let us briefly show ${h(Z)}$ describing _two parallel
strings_. By the definition, they cannot be made standing simultaneously.
Therefore, the straightforward idea is to create the first (generically
boosted) string, apply the rotation and another boost to induce its motion in
the perpendicular direction, rotate back and create the second parallel string
again with a generic boost. The relevant sequence in terms of elementary
operations is

$h(Z)=\mathcal{R}_{\\{\pi/2,0,0\\}}\mathcal{B}_{\\{w_{10}\\}}\mathcal{R}_{\\{-3\pi/2,0,0\\}}\mathcal{S}_{\\{\varepsilon\\}}\mathcal{R}_{\\{3\pi/2,0,0\\}}\mathcal{B}_{\\{w_{5}\\}}\mathcal{R}_{\\{3\pi/2,0,0\\}}\mathcal{B}_{\\{w_{3}\\}}\mathcal{S}_{\\{\delta\\}}Z\,.$
(73)

with the string ends given by

$Z=\left\\{0,\infty,\left(\frac{1}{w_{3}}\frac{1+w_{5}}{1-w_{5}}\right)^{\frac{1}{1-\delta}},\left(\frac{1}{w_{3}}\frac{1-w_{5}}{1+w_{5}}\right)^{\frac{1}{1-\delta}}\right\\}\,.$
(74)

The mapping (73) could be analysed in the same way as in two previous cases.
Here, we only show the effect of such a mapping in terms of the Riemann sphere
and deformation of the related ${t_{\rm i}^{+}=\mbox{const}}$ spherical
initial data, see figure 21.

Figure 21:  For the two parallel string configuration (73) we plot the
resulting unit Riemann sphere mapping in $\cal{M}^{+}$ (left). The string ends
are visualised in light blue. Divergences of its image (right) in
$\cal{M}^{-}$ indicated motion and mutual angles between these ends. In the
right picture, we also plot the unit sphere in dark blue to represent the
scaling.

The physical interpretation of the above construction (73) can be directly
deduced from the resulting motion of initially static test particles, see
figure 22. The deformation is again induced by the dragging of geodesics due
to the motion of the string ends.

Figure 22: The initially static spherical distribution of test particles
interacts with the expanding impulse generated by (73) corresponding to a
snapped pair of parallel cosmic strings. The only one typical evolution step
of such a distribution along geodesics (21) is depicted.

##  7 Conclusions

We studied geometries representing expanding impulsive gravitational waves.
Our main aim was a geometric description and physical interpretation of the
complex mapping $h(Z)$ entering the Penrose junction conditions (11). In
particular, situations related to the string-like nature of the wave source
were elaborated. The $h(Z)$ properties and its refractive effects were
connected with the motion of free test particles crossing the null wave
surface. This was possible due to employing the recent rigorous results on the
geodesic motion in expanding impulses. To clarify the role of elementary steps
in $h(Z)$ construction, we analyzed specific initially static classes of
geodesic congruences. Such an approach was introduced in the case of boosted
one-string as the simplest possibility. Subsequently, it was followed by three
situations with two strings, where we clarified and extended previous results.
In particular, we studied a snapping pair cosmic strings, its degenerate
subclass with only three string pieces generating the impulse, and finally, a
case of two parallel strings. In general, the effects of expanding impulses
generated by snapped cosmic strings acting onto geodesic motion can be
described as a dragging of test particles by the string ends and induced
motion in their directions. This analysis also showed non-trivial inherent
boost-like effects within the construction of the two string scenarios and the
way of its compensation.

## Acknowledgements

DK acknowledges the support of the Czech Science Foundation, Grant 21-11268S.
RŠ was supported by the Czech Science Foundation Grant No. GAČR 22-14791S.

## References

  * [1] Einstein A 1916 Näherungsweise Integration der Feldgleichungen der Gravitation Sitz. Preuss. Akad. Wiss. Berlin 688–96
  * [2] Einstein A 1918 Über Gravitationwellen Sitz. Preuss. Akad. Wiss. Berlin 1 154–67
  * [3] Pirani F A E 1957 Invariant formulation of gravitational radiation theory Phys. Rev. 105 1089–1099
  * [4] Bondi H 1957 Plane gravitational waves in general relativity Nature 179 1072–1073
  * [5] Hulse R A and Taylor J H 1975 Discovery of a pulsar in a binary system Astroph. J. 195 L51–L53
  * [6] Abbott B P et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2016 Ofor bservation of Gravitational Waves from a Binary Black Hole Merger Ptwo-stringett. 116 061102
  * [7] Abbott B P et al. (LIGO Scientific Collaboration and Virgo Collaboration) 2017 GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral Phys. Rev. Lett. 119 161101
  * [8] Robinson I and Trautman A 1960 Spherical gravitational waves Phys. Rev. Lett. 4 431–2
  * [9] Robinson I and Trautman A 1962 Some spherical gravitational waves in general relativity Proc. Roy. Soc. A 265 463–73
  * [10] Kundt W 1961 The plane-fronted gravitational waves Z. Physik 163 77–86
  * [11] Kundt W 1962 Exact solutions of the field equations: twist-free pure radiation fields Proc. Roy. Soc. A 270 328–34
  * [12] Stephani H, Kramer D, MacCallum M, Hoenselaers C, and Herlt E 2003 Exact Solutions of Einstein’s Field Equations 2nd ed. (Cambridge: Cambridge University Press)
  * [13] Griffiths J B and Podolský J 2009 Exact Space-Times in Einstein’s General Relativity (Cambridge: Cambridge University Press)
  * [14] Grosser M, Kunzinger M, Oberguggenberger M, and Steinbauer R 2001 Geometric Theory of Generalized Functions Mathematics and its Applications 537 (Dordrecht: Kluwer)
  * [15] Steinbauer R and Vickers J A 2006 The use of generalized functions and distributions in general relativity Class. Quantum Grav. 23 R91–R114
  * [16] Penrose R 1972 The geometry of impulsive gravitational waves; in General Relativity, ed. L. O’Raifeartaigh (Oxford: Clarendon Press) 101–15
  * [17] Nutku Y and Penrose R 1992 On impulsive gravitational waves Twistor Newsletter 34 9–12
  * [18] Hogan P A 1992 A spherical gravitational wave in the de Sitter universe Phys. Lett. Math. A 171 21–22
  * [19] Podolský J and Griffiths J B 1999 Expanding impulsive gravitational waves Class. Quant. Grav. 16 2937–46
  * [20] Hogan P A 1995 Imploding-exploding gravitational waves Lett. Math. Phys. 35 277–280
  * [21] Aliev A N and Nutku Y 2001 Impulsive spherical gravitational waves Class. Quant. Grav. 18 891–906
  * [22] Podolský J 2002 Exact impulsive gravitational waves in spacetimes of constant curvature; in Gravitation: Following the Prague Inspiration, ed. O. Semerák, J. Podolský and M. Žofka (Singapore: World Scientific) 205–46
  * [23] Podolský J, Sämann C, Steinbauer R, and Švarc R 2016 The global uniqueness and $C^{1}$-regularity of geodesics in expanding impulsive gravitational waves Class. Quantum Grav. 33 195010
  * [24] Podolský J and Griffiths J B 2000 The collision and snapping of cosmic strings generating spherical impulsive gravitational waves Class. Quant. Grav. 17 1401–13
  * [25] Podolský J and Steinbauer R 2003 Geodesics in spacetimes with expanding impulsive gravitational waves Phys. Rev. D 67 064013
  * [26] Podolský J and Švarc R 2010 Refraction of geodesics by impulsive spherical gravitational waves in constant-curvature spacetimes with a cosmological constant Phys. Rev. D 81 124035
  * [27] Steinbauer R 2014 Every Lipschitz metric has $\mathcal{C}^{1}$-geodesics Class. Quantum Grav. 31 057001
  * [28] Filippov A F 1988 Differential Equations with Discontinuous Righthand Sides (Dordrecht: Kluwer Academic Publishers)
  * [29] Cortés J 2008 Discontinuous dynamical systems IEEE Control Systems Magazine 28 36–73
  * [30] Podolský J, Sämann C, Steinbauer R, and Švarc R 2015 The global existence, uniqueness and $C^{1}$-regularity of geodesics in nonexpanding impulsive gravitational waves Class. Quantum Grav. 32 025003