Short-lived flashes of gamma radiation in a quasi-spherical Szekeres metric
Andrzej Krasiński





N. Copernicus Astronomical Centre, Polish Academy of Sciences, 

Bartycka 18, 00 716 Warszawa, Poland akr@camk.edu.pl




In previous papers it was shown that gamma rays with characteristics similar to
those of the gamma-ray bursts (GRBs) observed by astronomers may arise from
suitably shaped nonuniformities in the Big Bang in a quasi-spherical Szekeres
(QSS) model. The gamma radiation arises by blueshifting the light emitted by
hydrogen atoms at the end of the last scattering epoch along preferred
directions that exist in QSS models. However, the durations of the gamma flashes
and of their afterglows implied by the model were much longer than those of the
observed GRBs. In this paper it is shown that for the gamma-ray flash a duration
of correct order results if the blueshifted radiation, on its way to the present
observer, passes through another QSS region where it is deflected. The angle of
deflection changes with time because of the cosmological drift mechanism, so the
high-frequency ray will miss the observer after a while. It is shown by explicit
numerical calculation that a gamma-ray flash will no longer be visible to the
present observer after 10 minutes.


Short-lived flashes of gamma radiation in a quasi-spherical Szekeres metric
    Andrzej Krasiński
    
===========================================================================





§ MOTIVATION AND BACKGROUND




In previous papers by this author <cit.> it was
shown that flashes of gamma radiation with characteristics similar to those of
the gamma-ray bursts (GRBs), now routinely observed by astronomers
<cit.>, may arise from a non-simultaneous Big Bang (BB) in a
quasi-spherical Szekeres (QSS) model <cit.>. The
gamma radiation arises by blueshifting the light emitted by hydrogen atoms at
the end of the last scattering epoch along preferred directions existing in QSS
models <cit.> (for the first discussions of blueshifting see
Refs. <cit.>). With the BB profile t_B(r) chosen suitably,
the blueshift is sufficiently strong to move the observed frequency of the
radiation to the gamma range. A multitude of sources can be created by
distributing many small QSS regions over a Friedmann background <cit.>.
However, the durations of the gamma flashes and of their afterglows implied by
the model were much longer than those measured for the GRBs.

In this paper it is shown that a duration of the gamma-ray flash of correct
order results if the blueshifted radiation, on its way to the present observer,
passes through another QSS region, where it is deflected. The angle of
deflection changes with time in consequence of the cosmic drift
<cit.>, so the high-frequency ray will miss the
observer after a while.

Sections <ref> – <ref> are partly repeated after Ref. <cit.>;
they present the QSS model used in this paper (Sec. <ref>), the null
geodesic equations and properties of redshift along them (Sec. <ref>)
and the definition of the extremum redshift surface (ERS, Sec. <ref>).

In Sec. <ref>, the configuration of the QSS regions is presented, and
the time-dependent deflection of light rays in them is demonstrated on a
numerical example. In Sec. <ref> the rate of angular drift of the
deflected rays is calculated. In Sec. <ref> it is shown by explicit
numerical calculation that in the configuration of Sec. <ref> the
present observer who registered a gamma-radiation flash at time t_o will no
longer see it at t'_o = t_o + 10 minutes. Instead, the radiation coming from
the same direction at t'_o will have its frequency in the ultraviolet range.
(The 10 minutes were chosen as an exemplary duration of an observed GRB
<cit.>, but the gamma-ray flash is, in this model, instantaneous; see
Sec. <ref>.)

In Sec. <ref> it is shown that the brief duration of the gamma-ray flash
is preserved if the initial point and direction of the later-arriving ray are
slightly perturbed. Thus, the result of Sec. <ref> is not overly
sensitive to numerical inaccuracies. Section <ref> is a summary of the
results and of the problems that remain to be solved. The appendices present
some details of the computations.



§ THE QUASISPHERICAL SZEKERES (QSS) SPACETIME USED IN THIS
PAPER




The signature and labelling of coordinates will be (+, -, -, -) and
(x^0, x^1, x^2, x^3) = (t, r, x, y) or (t, r, ϑ,
φ).

The metric of the QSS spacetimes is <cit.>

    d s^2 =  d t^2 - (Φ,_r - Φ E,_r/
    E)^2/1 + 2 E(r) d r^2 - (Φ/ E)^2
    ( d x^2 +  d y^2),

where

    ES/2[(x - P/S)^2 + (y -
    Q/S)^2 + 1],

P(r), Q(r), S(r) and E(r) being arbitrary functions such that S ≠ 0
and E ≥ -1/2 at all r.

The source in the Einstein equations is dust (p = 0) with the velocity field
u^α = δ_0^α. The surfaces of constant t and r are
nonconcentric spheres, and (x, y) are the stereographic coordinates on each
sphere. At a fixed r, they are related to the spherical coordinates by

    x    =    P + S (ϑ/2) cosφ, 
    
    y    =    Q + S (ϑ/2) sinφ.

The functions (P, Q, S) determine the centres of the spheres in the spaces of
constant t (see illustrations in Ref. <cit.>). Because of the
non-concentricity, the QSS spacetimes in general have no symmetry
<cit.>.

With Λ = 0 assumed, Φ(t,r) obeys

    Φ,_t^2 = 2 E(r) + 2 M(r)/Φ,

where M(r) is an arbitrary function. We will consider only models with E >
0, then the solution of (<ref>) is

    Φ(t,r)    =   M/2E (coshη - 1), 
    sinhη - η   =   (2E)^3/2/M[t - t_B(r)],

where t_B(r) is an arbitrary function; t = t_B(r) is the time of the BB
singularity, at which Φ(t_B, r) = 0. We assume Φ,_t > 0 (the Universe
is expanding).

The mass density implied by (<ref>) is

    κρ = 2 (M,_r - 3 M  E,_r /  E)/Φ^2
    (Φ,_r - Φ E,_r /  E),   κ8
    π G/c^2.

This is a mass-dipole superposed on a spherical monopole <cit.>,
<cit.>. The dipole vanishes where E,_r = 0. The density is
minimum where E,_r/ E is maximum and vice versa <cit.>.

The arbitrary functions must be such that no singularities exist after the BB.
This is ensured by <cit.>:

    M,_r/3M   ≥    P/S,     E,_r/2E >  P/S    ∀ r, 
     where   P      √((S,_r)^2 + (P,_r)^2 + (Q,_r)^2).

These inequalities imply <cit.>

    M,_r/3M≥ E,_r/ E,     E,_r/2E >
     E,_r/ E    ∀ r.


The extrema of E,_r/ E with respect to (x, y) are
<cit.>

    . E,_r/ E|_ ex = ε_2 
    P/S,     ε_2 = ± 1,

with + at a maximum and - at minimum; they occur at

    x = P + SP,_r/ P + ε_2 S,      y = Q + SQ,_r/ P + ε_2 S.


The Lemaître <cit.> – Tolman <cit.> (L–T) models are
contained in (<ref>) – (<ref>) as the limit of constant (P, Q, S).
The Friedmann limit is obtained from QSS when E / M^2/3 and t_B are
constant (then (P, Q, S) can be made constant by a coordinate transformation).
QSS and Friedmann spacetimes can be matched at any constant r.

Because of p = 0, the QSS models can describe the evolution of the Universe no
further back in time than to the last scattering hypersurface (LSH); see Sec.
<ref>.

We will consider such QSS spacetimes whose L–T limit is Model 2 of Ref.
<cit.>. The r-coordinate is chosen so that

    M = M_0 r^3,

and M_0 = 1 (kept in formulae for dimensional clarity) <cit.>. The
function E(r), assumed in the form

    2E/r^2  -k = 0.4,

is the same as in the background Friedmann model.

The units used in numerical calculations were introduced and justified in Ref.
<cit.>. Taking <cit.>

    1  pc = 3.086 × 10^13  km,    1  y = 3.156 ×
    10^7  s,

the numerical length unit (NLU) and the numerical time unit (NTU) are defined as
follows:

    1  NTU = 1  NLU = 9.8 × 10^10  y = 3 × 10^4 
    Mpc.




The BB profile is the same as in Ref. <cit.>, see Fig.
<ref>. It consists of two curved arcs and of a straight line
segment joining them. The upper-left arc is a segment of the curve

    r^6/B_1^6 + (t - t_ Bf - A_0)^6/B_0^6 =
    1,

where

    t_ Bf = -0.13945554689046649  NTU≈ -13.67 × 10^9 
    years.

The t_ Bf is the asymptotic value of t_B(r) in the L–T model that
mimicked accelerating expansion <cit.>. This differs by
∼ 1.6 % from (- T), where T is the age of the Universe determined by
the Planck satellite <cit.>

    T = 13.819 × 10^9  y = 0.141  NTU.

The lower-right arc is a segment of the ellipse[In further figures it
looks like a vertical straight line segment because A_1 is extremely small.]

    (r - B_1 - A_1)^2/A_1^2 + (t - t_ Bf -
    A_0)^2/A_0^2 = 1.

The straight segment passes through the point where the full curves would meet;
its slope is determined by x_0.

The free parameters are A_0, A_1, B_0, B_1 and x_0. In Fig.
<ref> the values of x_0 and A_1 are greatly exaggerated to
improve readability. The actual values are <cit.>

    ([ A_0; B_0; A_1; B_1; x_0;     ]) = ([ 0.000026  NTU; 0.000091  NTU;    1 × 10^-10;         0.015;        10^-11;               ]),

and this BB profile is shown[Figure <ref> here is different
from Fig. 2 in Ref. <cit.> – the latter showed an intermediate profile
before final optimization.] in Fig. <ref>.



The QSS model is axially symmetric, with P(r) = Q(r) = 0 and

    S(r) = √(a^2 + r^2)

(the same as in Ref. <cit.>), where a^2 = 0.001. This S(r) obeys
(<ref>), which, using (<ref>) – (<ref>), reduces to

    1/r > S,_r / S.

The axis of symmetry is at x = y = 0. At r > r_F, where

    r_F = A_1 + B_1 = 0.0150000001,

the BB profile becomes horizontal straight, and the geometry of the model
becomes Friedmannian. See Sec. <ref> for remarks on the choice of
coordinates in that region.

The Friedmann background is defined by

    d s^2 =  d t^2 - ℛ^2(t) [ d r^2/1 - k
    r^2 + r^2 ( dϑ^2 + sin^2 ϑ dφ^2)],

where k and t_ Bf are given by (<ref>) and (<ref>), and

    Λ = 0

is assumed; t = 0 is the present time, t_ Bf is the BB time. The scale
factor ℛ(t) is determined by (<ref>) with Φ = r ℛ(t), (<ref>) – (<ref>) and t_B(r) = t_ Bf = constant.

The mass-density at the last scattering time in the now-standard ΛCDM
model is <cit.>

    κρ_ LS = 56.1294161975316 × 10^9  ( NLU)^-2.

In the L–T and Szekeres models it is assumed that a light ray intersects the
LSH when the density calculated along it becomes equal to the ρ_ LS in
(<ref>).

With (<ref>), (<ref>), (<ref>) and (<ref>), the ρ_
LS of (<ref>) occurs at the redshift relative to the present time

    1 + z^b_ LS = 952.611615159.

This differs by ∼ 12.7 % from the ΛCDM value z_ LS = 1090
<cit.>. To bring our model to agreement with this, laborious
re-calculations would be required. Since the model needs improvements anyway, we
will rather stick to (<ref>), (<ref>) and (<ref>), to be able to
compare the present results with the earlier ones.



§ NULL GEODESICS IN THE AXIALLY SYMMETRIC QSS
SPACETIMES




In (<ref>) – (<ref>) x = ∞ and y = ∞ are at the pole of
the stereographic projection of a sphere. This is a coordinate singularity where
numerical integration of geodesics breaks down. So, we introduce the coordinates
(ϑ, φ) by

    x = S_F (ϑ/2) cosφ,      y = S_F (ϑ/2) sinφ,

where

    S_F  S(r_F) = √(a^2 + r_F^2)

is the value of S at the Szekeres/Friedmann boundary. This changes (<ref>)
and (<ref>) to

    d s^2    =    d t^2 -  N^2  d r^2/1 + 2 E(r) -
    (Φ/ F)^2 ( dϑ^2 + sin^2
    ϑ dφ^2), 
              
     F   =   S_F/2 S (1 + cosϑ) + S/2 S_F (1 -
    cosϑ),

where


    NΦ,_r - Φ F,_r/ F.

In general, (ϑ, φ) are not the spherical polar coordinates
because F depends on ϑ. The dipole equator F,_r = 0
is at (ϑ_ eq/2) = S/S_F (so ϑ_ eq = π/2 in
the Friedmann region including the QSS boundary, see further in this section).
At r = r_F we have F = 1 and (ϑ, φ) become the
spherical coordinates with the origin at r = 0.

Along a geodesic we denote

    (k^t, k^r, k^ϑ, k^φ) (t, r, ϑ,
    φ)λ,

where λ is an affine parameter. In the (ϑ, φ)
coordinates, the geodesic equations for (<ref>) – (<ref>) are

    k^tλ +  N N,_t/1 + 2E(k^r )^2
    + ΦΦ,_t/ F^2[(k^ϑ)^2 + sin^2
    ϑ(k^φ)^2] = 0,


    k^rλ   +    2  N,_t/ N k^t k^r 
       +   ( N,_r/ N - E,_r/1 + 2E)
    (k^r)^2 + 2 S,_r sinϑΦ/S  F^2  Nk^r k^ϑ
       -   Φ (1 + 2E)/ F^2  N[(k^ϑ)^2 + sin^2 ϑ(k^φ)^2] = 0, 
    k^ϑλ   +    2 Φ,_t/Φ k^t k^ϑ -
    S,_r sinϑ N/S Φ (1 + 2E) (k^r)^2 + 2
     N/Φ k^r k^ϑ
       +    F,_ϑ/ F [- (k^ϑ)^2
    + sin^2 ϑ(k^φ)^2] 
       -   cosϑsinϑ(k^φ)^2 = 0, 
    k^φλ   +    2 Φ,_t/Φ k^t k^φ + 2
     N/Φ k^r k^φ
       +    2 [cosϑ/sinϑ - 
    F,_ϑ/ F] k^ϑ k^φ = 0.

The geodesics determined by (<ref>) – (<ref>) are null when

    (k^t)^2 -  N^2 (k^r)^2/1 + 2E(r) -
    (Φ/ F)^2 [(k^ϑ)^2 + sin^2
    ϑ(k^φ)^2] = 0.


On past-directed rays k^t < 0, and λ along each of them can be chosen
such that at the observation point

    k^t_o = -1.

(On future-directed rays k^t > 0 and a convenient choice of λ is
k^t_e = +1.)

For correspondence with Ref. <cit.>, in the Friedmann region we choose
the coordinates so that

    S = √(a^2 + r_F^2) = S_F.

Then, throughout the Friedmann region, F = 1 and (ϑ, φ)
are the spherical coordinates. They coincide with the coordinates of the QSS
region at r = r_F.

To calculate k^r on nonradial rays, (<ref>) will be used, which is
insensitive to the sign of k^r. This sign will be changed at each point where
k^r reaches zero by the numerical program integrating {(<ref>),
(<ref>) – (<ref>)}.

Note that ϑ≡ 0 and ϑ≡π are solutions of
(<ref>). These rays intersect every space of constant t on the symmetry
axis; they are called axial rays. As follows from (<ref>), there exist no
null geodesics on which k^φ≡ 0 and ϑ has any constant
value other than 0 or π (because with 0 ≠ϑ≠π,
k^φ≡ 0 and k^ϑ = 0 at a point, Eq. (<ref>)
implies k^ϑλ≠ 0). Consequently, in the axially
symmetric case the only analogues of radial directions are ϑ = 0 and
ϑ = π; along these geodesics φ is undetermined.

Along a ray emitted at P_e and observed at P_o, with k^α being
affinely parametrised, we have

    1 + z = (u_α k^α)_e/(u_α
    k^α)_o,

where u_α are four-velocities of the emitter and of the observer
<cit.>. In our case, both the emitter and the observer comove with the
cosmic matter, so u_α = δ^0_α, and the affine parameter
is chosen so that (<ref>) holds; then

    1 + z = - k_e^t.

Equation (<ref>) has the first integral:

    k^φsin^2 ϑΦ^2 /  F^2 = J_0,

where J_0 is constant along each geodesic. When (<ref>) is substituted in
(<ref>), the following results:

    (k^t)^2 =  N^2 (k^r)^2/1 + 2E + (Φ/ F)^2 (k^ϑ)^2 + (J_0  F/sinϑΦ)^2.

At the observation/emission point, (<ref>)/(<ref>), respectively,
apply. Equations (<ref>) and (<ref>) show that for rays emitted at the
BB, where Φ = 0, the observed redshift is infinite when J_0 ≠ 0. A
necessary condition for infinite blueshift (1 + z_o = 0) is thus J_0 = 0, so

(a) either k^φ = 0,

(b) or ϑ = 0, π along the ray (note that (<ref>) implies
J_0/sinϑ→ 0 when ϑ→ 0, π).

Condition (b) appears to be also sufficient, but so far this has been
demonstrated only numerically in concrete examples of QSS models
(<cit.>).

Condition (a) is not sufficient, and Ref. <cit.> contains
numerical counterexamples: there exist rays that proceed in a surface of
constant φ, but approach the BB with z →∞; the value of
ϑ along them changes and is different from 0, π. For those rays,
(<ref>) with the last term being zero implies one more thing

    If lim_t → t_B z = ∞  and lim_t → t_B|k^r| < ∞
           then lim_t → t_B k^ϑ = ±∞,

i.e., such rays approach the BB tangentially to the surfaces of constant r.
Examples will appear in Sec. <ref>.

Consider a ray proceeding from event P_1 to P_2 and then from P_2 to
P_3. Denote the redshifts acquired in the intervals [P_1, P_2], [P_2, P_3]
and [P_1, P_3] = [P_1, P_2] ∪ [P_2, P_3] by z_12, z_23 and
z_13, respectively. Then, from (<ref>),

    1 + z_13 = (1 + z_12) (1 + z_23).

Thus, for a ray proceeding to the past from P_1 to P_2, and then back to the
future from P_2 to P_1:

    1 + z_12 = 1/1 + z_21.




§ THE EXTREMUM REDSHIFT SURFACE




Let a null geodesic stay in the surface {ϑ, φ} = {π, 
constant}. Then k^r ≠ 0 at all its points, see (<ref>). Assume it
is past-directed and has its initial point at r = 0. Thus, r has to increase
along it and can be used as a parameter. Using (<ref>), we then obtain from
(<ref>)

    z r =  N N,_t/1 + 2E k^r.

Since N≠ 0 from no-shell-crossing conditions <cit.> and
k^r > 0, the extrema of z on such a geodesic occur where

    N,_t ≡Φ,_tr - Φ,_t  F,_r/ F = 0.

Since ϑ = π was assumed, the set defined by (<ref>) is
2-dimensional; it is the Extremum Redshift Surface (ERS) <cit.>.
Equation (<ref>) is equivalent to <cit.>

    χ^4 + χ^3 = - k^3 [r t_B,r/4 M_0 (1 - r S,_r /
    S)]^2,

where

    χsinh^2 (η/2)

With k < 0, (<ref>) is solvable for χ at any r, since its left-hand
side is independent of r, monotonic in χ and varies from 0 to +∞,
while the right-hand side is non-negative. The right half of the ERS profile
with the parameters of (<ref>) is shown in Fig. <ref>.

In the limit S,_r → 0 (which occurs at a →∞), (<ref>)
reproduces the equation of the Extremum Redshift Hypersurface (ERH) of
Ref. <cit.>.

Equation (<ref>) was derived for null geodesics proceeding along ϑ
= π, where the mass dipole is maximum, F,_r/ F = S,_r/S > 0.
With S given by (<ref>) we have

    F_1  1/(1 - rS,_r/S) = (r/a)^2 + 1 > 1,

so, at a given r, the ERS has a greater η (greater t - t_B) than the
ERH of the corresponding L–T model. Moreover, the χ of (<ref>) is
greater when a is smaller.

Conversely, for a ray proceeding along the dipole minimum (where ϑ =
0), the factor F_1 is replaced by

    F_2  1/(1 + rS,_r/S) = a^2 + r^2/a^2 + 2r^2 < 1,

and so the ERS has a smaller t - t_B than the ERH in L–T. Also here, a
smaller a has a more pronounced effect.

Extrema of redshift exist along other directions than ϑ = 0 and
ϑ = π (see examples in Sec. <ref>), but a general equation
defining their loci remains to be derived.



§ THE TIME-DEPENDENT DEFLECTION OF LIGHT RAYS




From here on we will consider only the rays proceeding within the surface of
constant φ.

For the source of the gamma rays we take the BB hump of Fig. <ref>. We
assume that the maximally blueshifted axial ray emitted there (at the big dot in
Fig. <ref> and near to the tip of the arrow in the QSS1 circle in Fig.
<ref>) passes, on its way to the observer, through another axially
symmetric QSS region (QSS2 in Fig. <ref>), where it is deflected. The
parameters of QSS2 are assumed the same as in QSS1. In general, the angle of
deflection changes with time in consequence of the cosmic drift
<cit.>. Consequently, if the maximally blueshifted ray hits a
given observer at a given time, an axial ray emitted a while later will miss the
observer, apart from a few exceptional cases (see below). After a short while,
the rays seen by the observer in the same direction will be coming from other
regions of the same source and will be less blueshifted, i.e., they will form
the afterglow. See Sec. <ref> for a more detailed description.



A ray that hits QSS2 along the symmetry axis (from either the dipole maximum or
minimum side) is not deflected <cit.>. But the directions in the dipole
equator plane, ϑ = π/2, are also special. If a ray hits QSS2 there
having proceeded along ϑ = π/2 in the Friedmann region, then a
deflection will occur, but there will be no drift.[This effect was
verified numerically; its reason is not evident in the geodesic equations.] For
rays that hit QSS2 from other directions, the time-dependent deflection is to be
expected.

In the following, we will use two coordinate charts of the class (<ref>) –
(<ref>), one centered at the origin of QSS1, and the other centered at the
origin of QSS2. Their t will be the same, but the r-coordinates will differ.
Whenever confusion may arise, we will denote the r of QSS1 by r_ S1,
and the other one by r_ S2.

We consider a ray that was emitted from the LSH at the symmetry axis of QSS1 on
the dipole maximum side, and proceeds along this axis. The initial point is the
one that gives the maximally strong blueshift. (It is “ray A” in Sec. 9 of
Ref. <cit.> and “the ray” shown in Fig. <ref>). The ray then
passes all through QSS1 and emerges into the Friedmann region (after passing the
origin of QSS1 it becomes “ray B” of Ref. <cit.>). This part of the
ray's path is not shown in Fig. <ref>; in the following it will be called
the Upray. On it, we choose the point just behind the QSS1 region, with











    (r_ S1, t) = (r, t)_1 = 
           (0.0150907023847052114, -0.13926900571845113), 
           z_1 = -0.39590497158301252

(these numbers are taken from the numerical tables for ray B <cit.>).
We take (t_1, z_1) as the initial data for Ray 0. (Note: z_1 is the upward redshift between r_ S1 = 0 and r_ S1 = r_1. The true
(downward) redshift between these points is 1 + z_1  true = 1 / (1 +
z_1).)

For the initial r_ S2 on Ray 0 we take

    r_ S2 = r_i = 0.02.

We assume that the dipole maximum of QSS2 lies on ϑ = π in Fig.
<ref>, while the angle between the direction of the ray at (t_1, r_i)
and the dipole equator of QSS2 is

    δ = arctan(1/2);

see the uppermost inset in Fig. <ref>. The values of r_i and δ
were chosen such that the drift is clearly visible in the illustrations, while
the two QSS regions are not too close to each other. The angle of deflection
becomes larger when r_i is smaller and vice versa. The angle δ is
related to the initial k^ϑ by

    sinδ = k^ϑ_i Φ(t_i, r_i)/1 + z_1

(see derivation in Appendix <ref>), with Φ(t_i, r_i) calculated
from (<ref>) using (<ref>) and (<ref>), and (<ref>) for
t_B.

The upward redshift to t = 0 on Ray 0 is[For reproducibility of the
results, the numbers are quoted up to 17 decimal digits. Such precision was
needed to capture time differences of ≈ 10 min at the observer – see
Secs. <ref> – <ref>.]








    (1 + z)_ 3 up = 7.19256480334251602 × 10^-3.

The z_1 in (<ref>) is the upward z between r_ S1 = 0 and r_
S1 = r_1. Consequently, (<ref>) is the upward z between r_ S1 = 0
and the present observer. The proper z between the LSH and r_ S1 = 0
was calculated in Ref. <cit.>:

    1 + z_ ols3 = 1.11939135405414447 × 10^-7.

So, the proper 1 + z between the last scattering in QSS1 and the present time
is









    1 + z_3 = 1 + z_ ols3/1 + z_ 3 up = 1.55631737031266357
    × 10^-5.

This is near to 1.553 × 10^-5 of Ref. <cit.>; the difference
arose because in Ref. <cit.> the intervening QSS2 region was absent.
The 1 + z_3 of (<ref>) is within the range 1 + z < 1.689 × 10^-5
needed to blueshift the emission frequencies of hydrogen into the gamma sector
<cit.>.

Ray 0 overshot the present time t = 0 in consequence of numerical errors. The
calculation stopped at













    t_ fin = 1.29738987343870121 × 10^-10  NTU,

and the other coordinates of the endpoint were

    (r_ S2, ϑ) = (r, ϑ)_ fin = 
           (0.85217686701400219, 4.2386753628314251)

(recall: ϑ is counted from the X < 0 direction clockwise).

Rays 1 and 2 have their initial points at the same r_ S2 = r_i = 0.02 as
Ray 0, and the initial direction still determined by δ of (<ref>),
but the initial t later by 0.00002 NTU (Ray 1) and 0.00004 NTU (Ray 2). The
angle δ being the same means that the initial directions are
parallel-transported along the world line of the emitter; see Appendix
<ref> for more on this.

As seen in Fig. <ref>, each ray is deflected by a different angle when
passing through QSS2. This is the cosmic drift (non-repeatability of light
paths) of Refs. <cit.>. On comparing the two upper
insets in Fig. <ref> one sees that the projections of Rays 0 and 2
intersect at the edge of QSS2: Ray 0 goes further than Ray 2 around the BB hump.
Similarly, the projections of Ray 0 and Ray 1 intersect at point P of
coordinates




(X, Y) ≈ (0.019987, -0.02143).





Although the projections of Rays 0 and 1 on the (X, Y)-surface intersect at P,
the rays do not intersect in spacetime, as is seen from their projections on the
(X, t) and (Y, t) coordinate surfaces in Figs. <ref> and <ref>,
respectively. The intersecting (X, Y)-projections imply that an observer at P
would register both Ray 0 and (suitably later, from a different direction) Ray
1.

Figures <ref> and <ref> show [1 + z(r)] × 10^5 along Ray 0.
In short segments near the borders of QSS2, 1 + z decreases when followed
along the ray to the future. However, the times of flight of the rays through
QSS2 are made longer by the deflections. The net result is that the ray spends
more time in the redshift-generating region, and the final 1 + z is larger
than it would be in the absence of QSS2 – see (<ref>) and the remark under
it.








§ RATE OF ANGULAR DRIFT OF THE RAYS









Now we calculate the average rate of the angular drift of the rays applying a
few approximations. As the first approximation, we take the angles between the
images of the future-end segments of the rays in Fig. <ref> for the real
angles in spacetime. Since these angles are small, the departure from the exact
result should be tolerable. As shown in Appendix <ref>, the angle
between the end segments of Rays 2 and 0 is

    Δβ = 1.79491638375138503^∘.

The time difference between rays 0 and 2 at their initial points at r_ S2
= r_i is, from (<ref>) and (<ref>),

    Δ t = 0.00004  NTU = 39.2 × 10^5  y
           = 2.06192 × 10^12  min.

This implies the average rate of angular drift, with respect to the source time

    β t    =    4.57886832589639038 × 10^-7 ^∘/ y
       =    8.70507286292089426 × 10^-13 ^∘/ min
       =    0.15193218308397861 × 10^-13  rad/ min.

In (<ref>), the time is counted between the initial points of Rays 0 and 2,
at (r, ϑ)_ S2 = (r_i, π/2), with r_i given by (<ref>).
The upward redshift on Ray 0 from r_ S1 = 0 to the present time is the
z_ 3 up given by (<ref>). Hence, the proper redshift on Ray 0
between r_ S2 = r_i and the present time is

    1 + z_ ip = 1 + z_1/1 + z_ 3 up = 83.98881969561866230,

where z_1 is given by (<ref>). So, the time difference Δ t at
(r_i, t_1) results in the time difference (1 + z_ ip) ×Δ t
at the observer. The rate of angular drift at the observer is thus

    .α t|_ observer = .β t/
    (1 + z_ ip) 
           = 1.80895723543432844 × 10^-16  rad/min.

As a curiosity, in Ref. <cit.> (an off-center observer in a spherical
void, observing galaxies around her) the cosmic drift was ≤ 10^-6
arcsec/y ≈ 9.22 × 10^-18 rad/min, i.e. ≈ 19.6 times
smaller than (<ref>).

The observed GRBs typically last from less than a second to a few minutes
<cit.>. Let 10 minutes be the reference time. Imagine a ray in our model
reaching the observer 10 minutes later than Ray 0 from the same direction. Will
it still be a gamma ray?

A direct answer would require emitting a ray from the observer position back in
time later by

    τ = 10  min = 1.93993947388841468 × 10^-16  NTU

(from (<ref>) – (<ref>)) than the t_ fin of (<ref>), in
the direction opposite to Ray 0. But a ray having its initial point so close to
the endpoint of Ray 0 could be numerically distinguishable from Ray 0 only at
very high precision that is inaccessible to this author.[Such a ray
would travel through nearly the whole lifetime of the Universe. To capture time
differences counted in seconds, the numerical time-step would have to be of the
order of 1 s. Using (<ref>) and taking 13.819 × 10^9 y for the age
of the Universe <cit.>, the calculation would require ≈ 43
× 10^16 steps.] Therefore, we shall again resort to an approximate
estimate.

Ray 2 was sent from its initial point later by 0.00004 NTU than Ray 0. This
translates, by (<ref>), to

    Δ t_2 ≈ 0.00335955278782475  NTU

at the present time, t = 0. So, we send a ray (call it Backray 0) from the
observer position given by (<ref>) – (<ref>) backward in time, with
the initial t = t_ fin + Δ t_2, and with the initial direction
opposite to the final direction of Ray 0. There were difficulties with assigning
the correct number to the initial direction, see Appendix <ref>. The
number finally chosen was

    μ k^ϑ/k^r = 0.0107412585641537794

With this value, at the first contact with the edge of QSS2 Backray 0 and Ray 0
coincide to better than 10^-8 in both X and Y; see Appendix
<ref> again (but the coincidence is not as good all the way).



The integration of Backray 0 was stopped at the 2-surface S_b in spacetime
given by the equations

    r    =   r_i/sinϑ + tanδcosϑ≡r_i cosδ/sin (δ + ϑ), 
    φ   =    constant.

This surface is orthogonal to Ray 0 at its initial point S. It intersects the
plane of Fig. <ref> along the straight line OL, which has the
equation

    Y = r_i + X tanδ,

where r_i and δ are given by (<ref>) – (<ref>) and (X, Y)
are the coordinates in Fig. <ref>. The point where Backray 0 hit
S_b will be denoted F. Stopping Backray 0 there was convenient for prolonging
it into the QSS1 region.[The discontinuity in the slope of Backray 0 at
the edge of QSS2 in Fig. <ref> resulted from interpolation – only 1
in 600 calculated points is used by the drawing program. ]

The surface S_b and the line OL are at the same time tangent at S to the
circle of radius r_ S1 = r_1, where Ray 0 took over from Upray; see Fig.
<ref>.

The following parameters of Backray 0 at F will be needed for further
calculations:

    r_b    =    0.0196597503597857326, 
    
    t_b    =    -0.13922666531999014, 
    ϑ_b    =    1.5348829837350006, 
    
    1 + z_b    =    74.080313337692317, 
     R(t_b)    =    0.00617900818769279159.


Imagine a family of past-directed rays, all emitted from (r, ϑ) = (r,
ϑ)_ fin of (<ref>), in the direction opposite to the final
direction of Ray 0, with their emission instants continuously increased from t
= t_ fin of (<ref>) to t_ fin + Δ t_2 of (<ref>). Call
this family G (for “geodesics”). The earliest of these rays is the
time-reverse of Ray 0, the last coincides with Backray 0. Their intersections
with the surface S_b form an arc connecting point S to point F. The arc SF
represents the spatial drift of the endpoint of Backray 0 with respect to the
initial point of Ray 0.

The angle γ between the tangent vector to the arc SF at F and the
direction vector of Backray 0 at F is

    γ   =    1.5137231069990575  rad
       =    86.729945382475933^∘  ⟹
    γ      π/2 - γ = 0.0570732197958390142  rad
       =    3.2700546175240777^∘.

see Appendix <ref>. The (X, Y) coordinates of F are

    X_b    =    - r_b cosϑ_b = -7.05895596023793498 × 10^-4, 
    
    Y_b    =    r sinϑ_b = 0.0196470734618825471,

where r_b and ϑ_b are given by (<ref>) and (<ref>). The
time-coordinates of points S and F in spacetime differ by

    Δ t_1 = t_b - t_1 = 4.234039846099 × 10^-5  NTU.


Now take a ray of family G (call it Backray N) that was emitted back in
time from (r, ϑ)_ fin by U Δ t_2 later than the t_
fin of (<ref>), where U is small and Δ t_2 is given by
(<ref>). The point of intersection of Backray N with the surface S_b,
denoted P_N, lies on the arc SF near to the point S. Quantities referrring to
P_N will carry the subscript N. We make two more approximations:

1. The arc SP_N coincides with the straight line segment SP_N,
tangent to SP_N at S, where P_N is the projection of P_N on the
plane of Fig. <ref>.

2. The X_N coordinate of P_N, the angle γ_N between
Backray N and the normal to the arc SP_N at P_N, the Δ t_N  t_N -
t_1 and the redshift difference z_ ip - z_N are related to X_b,
γ, Δ t_1 and z_ ip - z_b all by the factor U.

It is not obvious how precise these approximations are because, as was seen in
Fig. <ref>, the directions of the rays do not necessarily change
uniformly: the projections of some rays intersect. We leave the problem of
verifying the correctness of these assumptions (and possible improvements upon
them) for the future. The approximation should be better, the shorter
SP_N is.

We wish to calculate by how much the past endpoint P_N and the direction of
Backray N at P_N will differ from those of Ray 0 at S when U = 10 
minutes / Δ t_2, i.e.,

    U =  5.77439795236701905 × 10^-14

from (<ref>) and (<ref>). Then we obtain at P_N

    X_N        -U r_b cosϑ_b 
       =    -4.07612208426468965 × 10^-17, 
    
    Y_N        r_i + X_N / tanδ = r_i + 2 X_N 
       =    1.99999999999999185 × 10^-2, 
    γ_N    =    U γ
       =    3.29563483524085620 × 10^-15  rad,    
    Δ t_N    =    U Δ t_1  
       =    2.44490310175544341 × 10^-18  NTU, 
    
    1 + z_N    =    1 + z_ ip - U (z_ ip - z_b) 
       =    83.988819695618090.

From (<ref>) – (<ref>), the (r, ϑ)_ S2 coordinates
of P_N are

    ϑ_N    =   arctan[Y_N/(- X_N)] 
       =    1.5707963267948946  rad, 
    
    r_N    =    Y_N / sinϑ_N 
       =    1.99999999999999185 × 10^-2.

The program that calculated this does not see the difference between Y_N and
r_N. This is because π/2 - ϑ_N = 2.03806104213235315 ×
10^-15, and the Fortran program finds sinϑ_N = 1 up to 16th
decimal place. (The difference π/2 - ϑ_N is not to be confused with
γ_N. The ϑ_N is the value of ϑ_ S2
at the endpoint of Backray N, while γ_N is the angle between
the direction of Backray N at the endpoint and the normal to SF there.)

The point P_N, lying on OL very near to the point S, lies very
nearly on the circle of radius r_ S1 = r_1 surrounding the QSS1 region.
So, we make one more approximation and assume that P_N lies on that
circle. The error in r committed thereby is smaller than 10^-17.



§ ACCOUNTING FOR THE SHORT-LIVEDNESS OF THE GAMMA-RAY FLASH




We now send a new ray (call it Downray 1), backward in time, in prolongation of
Backray N. What is the blueshift or redshift on this ray when it crosses the
LSH?

The numbers in (<ref>) and (<ref>) allow us to calculate the initial
direction for Downray 1 in terms of the (r, ϑ)_ S1 coordinates;
see the sketch in Fig. <ref>. It has its initial point at r_ S1 =
r_1 given by (<ref>), shifted by

    H       |X_N| / cosδ = √(5)|X_N| / 2 
       =    4.55724303250198597 × 10^-17

toward Y_ S1 > 0 with respect to the symmetry axis. Its initial direction
is inclined at the angle γ_N = π / 2 - γ_N to the OL
line. The ϑ_ S1 of its initial point is ε given by

    tanε = H/r_1 = 3.01990120560648060 × 10^-15.

The Fortran 90 program sees no difference between ε and tanε up to the level of 10^-17.



The tangent vector to Downray 1 at the point (t, r, ϑ) = (t_1 + Δ
t_N, r_1, ε) is

    k^r_N = √(1 - k r_1^2)(1 + z_N)/ R_N 
    D
               ×(- tanεsinγ_N - √(1
    - k r_1^2)cosγ_N),       
           k^ϑ_N = 1 + z_N/r_1  R_N  D
               ×(- √(1 - k r_1^2)sinγ_N +
    tanεcosγ_N),       
           D√(1 - kr_1^2 + tan^2 ε),    R_N  R(t_N);

see calculation and discussion in Appendix <ref>. In the numerical
calculations the factor (1 + z_N) was omitted – because smaller
values of k^r and k^ϑ result in smaller numerical steps of r and
ϑ which improves the numerical precision. The actual value of 1 + z
at the endpoint of Downray 1 was then calculated using (<ref>).

The initial t on Downray 1 is

    t_N = t_1 + Δ t_N = -0.013926900571845113   NTU.

Note that this is the same as t_1 of (<ref>) – with Δ t_N ≈
10^-18, the Fortran program does not see the difference. This confirms the
remark made under (<ref>) – if determined by direct numerical calculation,
the path of Backray N would become, at some point along the way,
indistinguishable from that of Ray 0.

Summarizing, the initial point of Downray 1 is (t, r, ϑ)_ S1 =
(t_N, r_1, ε), with the numbers given in (<ref>), (<ref>) and
(<ref>), respectively, while its initial direction is determined by
(<ref>) – (<ref>).



If continued from these initial data without corrections, Downray 1 would have
the (t, X) profile shown in Fig. <ref>. The t-coordinate along it
disagrees with that of the Upray by much more than should be expected, see
Appendix <ref> for details. This happens in consequence of a numerical
instability at the outer branch of the ERS, which this author was not able to
overcome. Therefore, somewhere along the way the t-coordinate of Downray 1 had
to be hand-corrected. This was done as follows. At the point of closest approach
of Downray 1 to the r = 0 line (see Appendix <ref>), the tangent
vector to the ray was parallel-transported from the actual time calculated along
the ray (t_o = -0.13932447906821049 NTU) to the initial time of the Upray,

    t'_o = -0.13932991589546649,

which is by 5.43682725598348959 × 10^-6 NTU earlier. The parallel
transport was done as in Appendix <ref>: z was thereby unchanged,
while k^r and k^ϑ were multiplied by R (t_o)/ R
(t'_o) = 1.3849964189634454. From there on, the ray was continued as Downray
2, see Figs. <ref> – <ref>. Consequences of other possible
modifications of Downray 1 will be discussed in Sec. <ref>.





The operation described above was done under two assumptions: (1) that the the
point of closest approach to r = 0 and the direction of Downray 1 there were
determined correctly, and (2) that the parallel transport could be done by the
Friedmannian rule of Appendix <ref>.

The first assumption is justified by the fact that Downray 1 proceeds close to
the ϑ = 0 axis, with |Y| < 4.404 × 10^-17 all the way between
P_N and X = 0, see Fig. <ref>. This means that any departure from
this path is smaller than the accuracy of the numerical code.



The second assumption is justified by the fact that at the closest approach to
the origin r ≈ 2.23 × 10^-22, E,_r/ E≈ 2.23
× 10^-19 and Φ≈ 2.15 × 10^-25, so the term in
(<ref>) that is responsible for the difference between the QSS and Friedmann
geometries is Φ E,_r/ E≈ 4.79 × 10^-44. The
Friedmann term Φ,_r is at this location 9.62 × 10^-4. The
non-Friedmannian correction is thus too small to show up at the level of
accuracy assumed here.

Near the origin some quantities become indeterminate of the type 0/0
<cit.>, so the integration of the ray equations had to include certain
precautions: (1) The numerical step had to be gradually decreased on approaching
the origin. (2) Downray 1 aimed at r = 0 with high precision, and at a certain
step its r-coordinate shot to r < 0. This had to be taken care of. See
Appendix <ref> for details.

Figure <ref> shows the (t, X) profiles of Downrays 1 and 2 across the
full diameter of the BB hump. Only the left half of the ERS is shown. Downray 1
at first goes along the steep slope of the BB hump, but then turns away and
crosses the LSH in the Friedmann region at

    ([       t;       r;       ϑ; 1 + z_P; 1 + z_o;         ])_D1LSH = ([ -0.13945067585653098  NTU;     0.0522717572229775690;        3.1415926411632658;        3.9974477620254634;       335.740919327408696;                           ])

(ϑ_D1LSH = 179.999999288^∘≈π). The z_P is relative
to the point P_N, z_o is relative to the present observer calculated from 1
+ z_o = (1 + z_N) (1 + z_P) using (<ref>).

Downray 2 never leaves the QSS1 region and finally hits the BB at the steep
slope – tangentially to a surface of constant r, as predicted in
(<ref>). It crosses the LSH at the point marked with the small square in
Fig. <ref> and with the dot in Fig. <ref>, at

    ([       t;       r;       ϑ; 1 + z_P; 1 + z_o;         ])_D2LSH = ([ -0.13939573683235620  NTU;     0.0149934194406606183;        3.6918768830001061;     0.0268536489146646434;     2.2554062768632145161;                           ]).


Figure <ref> shows the projection of Downrays 1 and 2 on the (X, Y)
coordinate plane. At this scale, both rays seem to proceed across the QSS1
region along the symmetry axis. Figures <ref> and <ref> show
details of the rays' paths. Downray 1 proceeds along a nearly straight line at
the angle ε of (<ref>) to the Y = 0 line and intersects it at
X ≈ -0.0030588. At (X, Y) ≈ (-4.322714809677376 × 10^-5,
-1.02236722643458 × 10^-17) it turns upwards and aims at r = 0. After
reaching (X, Y) ≈ (-2.230943471939892 × 10^-22,
-1.07257712221974956 × 10^-25), the Fortran program finds the next r
to be negative, see Appendix <ref> to see how this problem was handled.
On the other side of the origin, the two rays proceed at different angles.
Downray 1 stays within the 0 ≤ Y < 2.5 × 10^-10 strip up to the
boundary of the QSS1 region; see the upper panel of Fig. <ref>.
Downray 2 is initially inclined to the Y = 0 line at an angle that is ≈
10^12 times greater than in the X < 0 sector. On approaching the boundary
of the QSS1 region, it bends around and becomes nearly tangent to a constant-r
surface. It continues this way up to the intersection with the BB; see the lower
panel of Fig. <ref> and Fig. <ref>.



Finally, Fig. <ref> shows the {X, z} relation along Downray 2, with
z relative to the present observer. The observer would see blueshift if the
light sources were placed between X ≈ 0.0135 and X ≈ 0.0149, but
will not see any blueshift in the light emitted at the LSH. For z at LSH given
by (<ref>), and with the parameters of the spatially flat ΛCDM model
(H_0 = 69.6 km/(s × Mpc), Ω_m = 0.286, Ω_Λ =
0.714 <cit.>) the Cosmology Calculator of Refs.
<cit.> says the source of Downray 2 would be 8.726 Gyr to the
past of the present observer. In the model used here, it is ≈ 0.139 NTU
= 13.622 Gyr to the past. With 1 + z_o given by (<ref>), the whole range
of visible light is moved into near ultraviolet.



Assuming that the approximations were all valid, Eqs. (<ref>) and
(<ref>), and Fig. <ref> prove that 10 minutes after seeing a
gamma-ray flash the present observer will see radiation coming from the same
direction with frequency (1 + z_o)_D2LSH / (1 + z_3) ≈ 1.449 ×
10^5 times smaller. This will be in the ultraviolet range.



§ PERTURBATIONS OF DOWNRAY 2




Since the numerical calculations involved approximations, and they proceeded at
the borderline of numerical accuracy, it is useful to investigate the redshift
along rays propagating in proximity of Downray 2. For each ray described below,
the numerical step on approaching the minimum r and the jump through r = 0
were handled in the same way as described in Appendix <ref>.

Four sets of rays were considered:

1. Downrays 3 and 4 had the same initial (t, r, ϑ = ε) and
the same initial direction as Downray 2, but the corrected t at closest
approach to r = 0 larger (ray 3) and smaller (ray 4) than the t'_o of
(<ref>).

In the next three sets, the corrected t at closest approach to r = 0 was
(<ref>).

2. Downrays 5 – 7 had the same initial (t, r, ϑ = ε) as
Downray 2, but γ_N zero, minus that of Downray 2 and twice
that of Downray 2, respectively (γ_N = 0 means that the
initial direction is parallel to the Y = 0 line).

3. Downrays 8 – 10 had the same initial (t, r) as Downray 2, but
ε twice as large, and the same values of γ_N as
Downrays 5, 6 and 2, respectively.

4. Downrays 11 – 13 had the same initial (t, r) as Downray 2, but half the
value of ε, and the same set of values of γ_N as
in point 3.

For Downray 3 trial time-intervals τ_n = 10^-n NTU were added to the
t'_o of (<ref>) to verify the effect. With n = 18 the effect was
invisible, as expected. With n = 17, 16, …, 8 the changes in (t, r,
ϑ, 1 + z_P) at the LSH compared to those of Downray 2 were small and
qualitatively insignificant. Only at n = 7 a substantial change occurred. The
ray hit the BB at the boundary of the QSS1 region at earlier t than Downray 2,
but the program could not detect its intersection with the LSH. The ray with
τ_7.2 = 2 × 10^-7 NTU (Downray 3 in the figures) flew over the BB
hump and intersected the LSH in the Friedmann region, with (r, 1 + z_o) ≈
(0.03724, 8.0629).

For Downray 4, the trial intervals added to t'_o were τ_n = -
10^-n NTU. Similarly to Downray 3, a noticeable difference occurred only at
n = 6: the ray hit the LSH well before reaching the edge of the QSS1 region,
with (r, 1 + z_P) ≈ (0.0145, 182.256).

The conclusion is that the value of 1 + z at the LSH is not sensitive to the
values of τ_n and τ_n as long as τ≤ 10^-8 NTU =
980 years, i.e., it is not sensitive to numerical errors smaller than 10^-8
NTU.

Details of the other results are described in Appendix <ref>. Here is
just a short summary:

On rays 5 to 9 and 12, 1 + z_o was sufficient to move the whole visible light
range into the ultraviolet; see Table <ref> in Appendix
<ref> and Ref. <cit.> for the numbers. On Downray 10, the
visible range was moved into the near ultraviolet just adjacent to the visible
range. On rays 11 and 13, the visible range was slightly shifted toward violet
(Downray 11) or toward red (Downray 13), but most of it remained visible.

The final conclusion is that the result calculated in Sec. <ref> remains
valid for past-directed rays that reach the neighbourhood of the point P_N
with directions close to that of Downray N: on all those rays, the present
observer would see redshift or only weak blueshift (1 + z ≈ 0.91 on
Downray 11). A strong blueshift 1 + z ≈ 10^-5 would only be seen on
Ray 0 that goes along the symmetry axis of the QSS1 region.

In fact, this result was to be expected on the basis of the comments to eq.
(<ref>), but it is instructive to have it confirmed by explicit
computations.



§ CONCLUSIONS AND PROSPECTS




The present paper is a continuation and development of Ref. <cit.>. In
that paper, it was shown that if a light ray is emitted at the last-scattering
hypersurface in an axially symmetric quasispherical Szekeres region (denoted
QSS1) along the symmetry axis, with the frequency in the visible range, then it
will follow this axis all through QSS1 and might be blueshifted by so much that
the present observer will see it as a gamma ray. The axial direction is
unstable, and rays proceeding in a close neighbourhood of the symmetry axis will
be now observed in the ultraviolet range.

In the present paper, it was investigated what happens to the maximally
blueshifted ray (MBR) if, on its way between the emission point and the present
observer, it passes through another QSS region (denoted QSS2). Both QSS regions
are matched into a Friedmann background of negative spatial curvature, and the
observer is in this background.

The QSS metrics and some of their properties were introduced in Secs. <ref>
– <ref>. In Sec. <ref>, it was shown that if the MBR passes through
QSS2 along a generic path,[Non-generic paths are along the symmetry
axis and in the plane of the Szekeres dipole equator.] then is deflected by a
continuously changing angle via the cosmic drift effect
<cit.>. Therefore, any observer can find herself in this ray's
path only at a single instant of time. At all other times, the rays coming from
the same direction will have lower frequencies (i.e. will be redshifted or only
slightly blueshifted), while the MBR will miss the observer. The rate of angular
drift of the MBR was calculated in Sec. <ref>. Using this result, in Sec.
<ref> it was demonstrated by a numerical calculation that an MBR that
reached a present observer at some instant will no longer be visible to her
after 10 minutes. Instead, the ray coming from the same direction will be in the
ultraviolet range. Finally, in Sec. <ref> (with numerical details in
Appendix <ref>) it was shown that the result of Sec. <ref> is
not sensitive to perturbations smaller than 10^-8 NTU = 980 years of the
instant when the rays cross the origin of the QSS1 region, and smaller than
≈ 3.3 × 10^-15 rad = 1.89 × 10^-13 ^∘ of their
direction on exit from the QSS1 region.

In the present model the gamma-ray flash is instantaneous. As with any model,
this is an idealisation. In reality, mechanisms should exist that will stretch
the instant of observation of the gamma ray to a short but finite interval.
Modelling them may require using more general metrics and other BB profiles than
those used here. Refs. <cit.> and the present paper
provided a proof of existence of a mechanism that can shift the range of visible
light frequencies to the gamma range, and that will ensure short-livedness of
the high-frequency flash.

The two QSS regions were assumed axially symmetric for simplicity. Then the
geodesic equations show that the MBR must follow the symmetry axis all through
the QSS1 region, so in the numerical calculations it can be kept exactly on this
path, even though it is unstable <cit.>. But strong
blueshifts are generated also by fully nonsymmetric QSS regions
<cit.>. In the nonsymmetric case, the MBR will be changing its
direction with time in consequence of the cosmic drift within the source
<cit.>, even in the absence of the deflecting region. This
process is a more difficult numerical challenge: the maximally blueshifted ray
must be identified by numerics alone <cit.>.

The model presented here does not solve the problem of too-long-lasting
afterglows. This will be the subject of a separate paper.





§ CALCULATING Δ IN (<REF>)




The unit vector tangent to the line (t, ϑ, φ) = (t_1, π/2, 
constant) in the Friedmann region is, using (<ref>),

    n^α_2 = (0, √(1 + 2E)/ N, 0, 0).

The unit direction vector of a null vector k^α is <cit.>

    n^α_1 = u^α - k^α/k^ρ u_ρ.

In our case u^α = δ^α_0. Consequently, at the initial point
of the ray, using (<ref>), eq. (<ref>) implies

    n^α_1 = (0, k^r_i/1 + z_1,  k^ϑ_i/1 +
    z_1, 0),

where z_1 is given by (<ref>) (account was taken of k^φ = 0).
From (<ref>) with J_0 = 0 we have, again using (<ref>),

    k^r_i = ±√(1 + 2E_i)/ N_i√((1 + z_1)^2 -
    (Φ_i k_i^ϑ/ F_i)^2).

So, the angle δ between n^α_1 and n^α_2 is given by

    cosδ = - g_αβ n^α_1 n^β_2 = √(1 -
    [k^ϑ_i Φ_i/ F_i (1 +
    z_1)]^2).

Since the initial point of Ray 0 is in the Friedmann region, where F≡ 1, the above is equivalent to (<ref>).



§ PARALLEL TRANSPORT ALONG U^Α IN A FRIEDMANN
REGION




In a Friedmann region a vector k^μ parallel-transported along u^α
= δ_0^α obeys k^t = constant, k^i ℛ = constant, i
= 1, 2, 3, i.e., the ratio k^ϑ/k^r stays constant. Then

    k^t_ℓ t   =    k^t_et, 
    
    (k^r, k^ϑ)_e t   =    R (t_ℓ t)/ R (t_e
    t) (k^r, k^ϑ)_ℓ t,

where ℓ t and et denote “at a later time” and “at an earlier time”,
respectively. If k^μ is a null vector and we wish to reset z from
z_ℓ t to z_e t, then (<ref>) is modified as follows:

    (k^r, k^ϑ)_e t =  R (t_ℓ t)/ R (t_e t)(k^r, k^ϑ)_ℓ t×1 + z_et/1 + z_ℓ t.




§ CALCULATING THE ANGLES FOR (<REF>) – (<REF>)




The direction of Ray 0 in the (X, Y) surface at the future endpoint can be
read off from the table of values of X and Y found by the same program that
calculated the path of the ray. The last two pairs in this table are

    (X, Y)_1 
           = (0.38875803859076119, -0.75833536622992936), 
           (X, Y)_2 
           = (0.38875804873013375, -0.75833538643402199).


Hence the final direction of Ray 0 in the (X, Y) surface is inclined to the Y
= 0 line at the angle β_0, where

    |tanβ_0|       | X Y|_r_ fin = X_2
    - X_1/|Y_2 - Y_1|
       =    0.50184745960551459, 
       ⟹   β_0 = 26.64966989421893587^∘.

In the same way we find for Ray 2:

    β_2 = 28.44458627797032091

So, Δβ = β_2 - β_0 is as given by (<ref>).



§ CORRECTING THE DIRECTION OF BACKRAY 0 IN SEC.
<REF>




The value of μ = k^ϑ/k^r at t = t_ fin of (<ref>), as
calculated by the program that determined Ray 0, was

    μ k^ϑ/k^r = 0.0100796385641537794.

But with this μ at the initial point, the path of Backray 0 at the entry
point to QSS2 visibly differed (at the scale of the insets in Fig. <ref>)
from that of Ray 0. This had to be caused by numerical errors because there is
no cosmic drift in the Friedmann region <cit.>. Consequently, the μ
in (<ref>) had to be hand-corrected by trial and error. With μ given by
(<ref>), the two rays coincide at the edge of QSS2 up to the level of
accuracy[See footnote <ref> – the irregularity of Ray 0 at the
edge of QSS2 is an interpolation by the drawing program.] shown in Fig.
<ref>.





§ CALCULATING THE ANGLE Γ IN FIG. <REF>




Let us use ϑ as the parameter along the arc SF. The unit tangent
vector to SF at (r, ϑ) = (r_b, ϑ_b) is

    n_ SF^α = 1/ R(t_b) ℓ (0, . r
    ϑ|_ϑ = ϑ_b, 1, 0),

where R(t_b) is given by (<ref>), and . r
ϑ|_ϑ = ϑ_b is found by taking (<ref>) at
(r, ϑ) = (r_b, ϑ_b):

    . r ϑ|_r = r_b = - r_b  (ϑ_b + δ).

The normalizing factor ℓ is

    ℓ = r_b √(^2 (ϑ_b + δ)/1 - k r_b^2 + 1).

The unit direction vector of Backray 0 at the same point, calculated by analogy
to (<ref>), is

    n^α_ ray = (0, k^r_b/1 + z_b, k^ϑ_b/1 + z_b, 0),

where z_b is given by (<ref>). Consequently, the angle γ between
n_ SF^α and n^α_ ray is given by

    cosγ = r_b  R(t_b)/ℓ(1 + z_b) [-
    k^r_b  (ϑ_b + δ)/1 - k r_b^2 + r_b
    k_b^ϑ].

The values given in (<ref>) were calculated from (<ref>).



§ CALCULATING (<REF>) – (<REF>)




In the coordinates of the QSS1 region, the line OL obeys the equation X ≡
- r cosϑ = constant, so the unit tangent vector to it, n_
OL^α, obeys n_ OL^r = r tanϑ n_ OL^ϑ,
n_ OL^t = n_ OL^φ = 0 and g_αβ n_
OL^α n_ OL^β = -1. The angle γ_N in Fig. <ref>
is between the initial direction of Downray 1 and the unit tangent vector to OL
pointing downward, so

    n_ OL^r    =    - tanϑ/ R(t) ℓ_ OL,     
    n_ OL^ϑ = - 1/r  R(t) ℓ_ OL, 
    ℓ_ OL   =   √(tan^2 ϑ/1 - k r^2 + 1).

The unit direction vector of Downray 1 at P_N is, similarly to (<ref>) and
(<ref>),

    n^α_N = (0, k^r_N/1 + z_N, k^ϑ_N/1 +
    z_N, 0).

At P_N, - g_αβ n_ OL^α n_N^β = cosγ_N =
sinγ_N, so

    tanε/1 - k r_1^2 k^r_N + r_1 k^ϑ_N = - ℓ_N/ R_N (1 + z_N) sinγ_N,

where D is given by (<ref>) and

    ℓ_N       ℓ_ OL(r_1, ε) =  D/√(1 - k
    r_1^2), 
    
    t_N        t_1 + Δ t_N,    R_N  R(t_N).


The vector k^α obeys g_αβ k^α k^β = 0. Since
Downray 1 prolongs Backray N, Eq. (<ref>) continues to hold along it.
Consequently, at (r, ϑ) = (r_1, ε) we have

    (1 + z_N)^2 -  R_N^2/1 - k r_1^2(k^r_N)^2 - (r_1  R_N)^2
    (k^ϑ_N)^2 = 0.

Solving the set (<ref>), (<ref>) for k^r_N and k^ϑ_N we
obtain

    k^r_N = √(1 - k r_1^2)(1 + z_N)/ R_N 
    D
               ×(- tanεsinγ_N ±√(1 - k r_1^2)cosγ_N),       
           k^ϑ_N = 1 + z_N/r_1  R_N  D
               ×(- √(1 - k r_1^2)sinγ_N
    ∓tanεcosγ_N).

The double sign in (<ref>) follows from the ambiguity in the definition of
k^r_N: its direction is given only by the angle it forms with OL, so it can
point away from QSS1 or into QSS1. We need to follow the second possibility,
which corresponds to the lower sign and implies k^r_N < 0. The upper sign
would correspond to k^r_N > 0 because √(1 - k r_1^2)cosγ_N ≫tanεsinγ_N. This is
because √(1 - k r_1^2) > 1, cosγ_N ≈ 1 while
sinγ_N and tanε are both tiny in consequence
of (<ref>) and (<ref>).

The situation with the sign of k^ϑ_N is less clearcut. The double
sign in (<ref>) is coupled with that in (<ref>), and having chosen
“-” in (<ref>) we have to choose “+” in (<ref>) – but then
k^ϑ_N might still be positive or negative, depending on the
numerical values of the quantities in it. Since π/2 - γ_N =
γ_N > ε as seen from (<ref>) and (<ref>),
and ϑ increases clockwise, the initial direction of Downray 1 is
toward decreasing ϑ, so k^ϑ_N < 0 follows. This is
qualitatively reflected in Fig. <ref>. The Fortran 90 program,using
(<ref>), found k^ϑ_N = -3.39090742613403823 × 10^-12.

Equations (<ref>) – (<ref>) are (<ref>) – (<ref>) with the
signs already chosen.



§ NUMERICAL CALCULATIONS OF DOWNRAYS 1 AND 2




Throughout this section, the r-coordinate is meant relative to the origin of
QSS1, so r = r_ S1 in all cases.

The difference of the t-coordinates of Upray and Downray 1 at the smallest r
(≈ 2.23 × 10^-22) is ≈ 5.437 × 10^-6. Since
Downray 1 is supposed to represent the ray that reaches the observer 10 min
later than Upray, their time-difference at r = 0 should be 10 min × (1 +
z)_ 3 up = 7.19 × 10^-2 min = 1.395 × 10^-18 NTU, where
(1 + z)_ 3 up is the upward 1 + z between r = 0 and the present
observer given by (<ref>). So clearly, the difference found by the Fortran
program is unrealistically large. If Downray 1 were continued uncorrected, it
would fly above the BB hump and intersect the LSH in the Friedmann region with
1 + z as given in (<ref>); part of its path is shown in Fig.
<ref>. Taken literally, this would mean that the frequency on it is no
longer in the gamma range at the observation event. However, this would not
prove the intended point of this paper. It is known that the high-frequency
flash will end if the observer simply waits long enough, but the time-scale of
this process is too long <cit.>. What we wish to show here is that 10
min after the observer saw the gamma flash, the ray coming from the same
direction will be emitted off the symmetry axis of QSS1 and will have a lower
frequency for this reason. Therefore, the t-coordinate on Downray 1 at closest
approach to r = 0 had to be decreased by hand. The ray on which this was done
is labelled as Downray 2.



Downray 2 runs near to the symmetry axis of QSS1 and is nearly parallel to it.
When the ray approaches r = 0, at a certain step the next r becomes
negative. In an attempt to prevent this, the step Δλ in the affine
parameter λ was gradually decreased on approaching the origin. Table
<ref> shows the reductions in Δλ and the values of r
at which they were done. But the jump to r < 0 persisted. So, when the Fortran
program reached the limit of its precision, the jump to r < 0 was interpreted
as the jump from the last (r, ϑ) where r was still positive to (r,
π + ϑ) (see Figs. <ref> and <ref>), and the program
continued from there. As Downray 2 was receding from the origin and passed
through the values of r from Table <ref> in a reverse order,
Δλ was multiplied by the reverse of the factor from the table at
each border value of r, except the one in the last line. In addition, at r =
0.01 and X > 0, Δλ was multiplied by 100. Figures <ref>
– <ref> show the resulting graphs.








In order to locate the intersection of Downray 2 with the LSH with a
satisfactory precision, more modifications of Δλ were made. It was
multiplied by 1/100 at step 114,685,150, and then again by 1/1000 at step
115,226,000. The LSH and BB were reached in steps 114,698,353 and 115,634,944,
respectively, but the last number is not reliable. Since the ray becomes tangent
to a surface of constant r when approaching the BB, it is difficult to locate
the intersection point – the result strongly depends on the value of Δλ. But the exact point where the ray hits the BB has no physical meaning
because the QSS metric does not apply at times before the LSH.



§ DETAILS FOR SEC. <REF>




Table <ref> gives the redshifts/blueshifts on Downrays 5 to 13
between the LSH and the point P_N (1 + z_P), and between the LSH and the
present observer (1 + z_o).



Figure <ref> shows where Downrays 3 to 13 intersect the LSH; the
intersection points are marked with large dots. Four typical examples of the
rays are displayed in full. Downray 3 escapes from the QSS1 region and
intersects the LSH in the Friedmann region, beyond Fig. <ref>.
Downray 4 intersects the LSH between r = 0 and the boundary of the QSS1
region, and hits the BB still inside QSS1. Downray 7 turns toward Y < 0 near
the boundary of the QSS1 region, and intersects the LSH near the steep slope of
the BB hump. Downray 12 turns toward Y > 0, and intersects the LSH near
another point of the steep slope of the BB hump. Downray 13 follows a part of
the path of Downray 7, the other rays follow parts of the path of Downray 12.
They intersect the LSH at various ϑ, but at r nearly on the boundary
of the BB hump.





Figures <ref> and <ref> show the rays between their initial
points and the neighbourhood of r = 0. In this segment, they stay within the
strip -2 × 10^-16 < Y < 6 × 10^-16 and pass through r = 0 in
ways similar to that shown in Figs. <ref> and <ref>.







Figures <ref> and <ref> show the rays between r = 0 and r
= r_F given by (<ref>). Here they depart from the Y = 0 line much
further than in the X < 0 segment: most stay in the strip -7 × 10^-5 <
Y < 2 × 10^-5, except for Downray 12 that at the farthest distance from
the Y = 0 axis has Y ≈ -1 × 10^-4. Close to the boundary of
QSS1 they bend around and become nearly tangent to surfaces of constant r.
Downray 3 is an exception: it follows the nearly-constant-r path for a short
while, then escapes to larger r, as seen in Fig. <ref>. Downray 4
hits the BB at r distinctly smaller than the r_F of (<ref>), the
remaining rays hit it at r very close to r_F.

Acknowledgement For some calculations, the computer algebra system
Ortocartan <cit.> was used.

99
Kras2016a A. Krasiński, Cosmological blueshifting may explain the
gamma ray bursts. Phys. Rev. D93, 043525 (2016).

Kras2016b A. Krasiński, Existence of blueshifts in quasi-spherical
Szekeres spacetimes. Phys. Rev. D94, 023515 (2016).

Kras2017 A. Krasiński, Properties of blueshifted light rays in
quasispherical Szekeres metrics. ArXiv 1704.08145, accepted for Phys. Rev.
D.

Perlwww D. Perley, Gamma-Ray Bursts. Enigmatic explosions from the
distant universe. 

http://w.astro.berkeley.edu/dperley/pub/grbinfo.html

Szek1975 P. Szekeres, A class of inhomogeneous cosmological models.
Commun. Math. Phys. 41, 55 (1975).

Szek1975b P. Szekeres, Quasispherical gravitational collapse. Phys. Rev. D12, 2941 (1975).

PlKr2006 J. Plebański and A. Krasiński, An Introduction to
General Relativity and Cosmology. Cambridge University Press 2006, 534 pp, ISBN
0-521-85623-X.

Szek1980 P. Szekeres, Naked singularities. In: Gravitational
Radiation, Collapsed Objects and Exact Solutions. Edited by C. Edwards.
Springer (Lecture Notes in Physics, vol. 124), New York, pp. 477 – 487 (1980).

HeLa1984 C. Hellaby and K. Lake, The redshift structure of the Big
Bang in inhomogeneous cosmological models. I. Spherical dust solutions. Astrophys. J. 282, 1 (1984) + erratum Astrophys. J. 294, 702
(1985).

KrBo2011 A. Krasiński and K. Bolejko, Redshift propagation
equations in the β' ≠ 0 Szekeres models. Phys. Rev. D83,
083503 (2011).

QABC2012 C. Quercellini, L. Amendola, A. Balbi, P. Cabella, M.
Quartin, Real-time cosmology. Phys. Rep. 521, 95 – 134 (2012).

KoKo2017 M. Korzyński and J. Kopiński, Optical drift effects in
general relativity. J. Cosm. Astropart. Phys. 03, 012 (2018).

Hell1996 C. Hellaby, The nonsimultaneous nature of the Schwarzschild
R = 0 singularity. J. Math. Phys. 37, 2892 (1996). publication.

BoST1977 W. B. Bonnor, A. H. Sulaiman and N. Tomimura, Szekeres's
space-times have no Killing vectors, Gen. Rel. Grav. 8, 549 (1977).

DeSo1985 M. M. de Souza, Hidden symmetries of Szekeres
quasi-spherical solutions, Revista Brasileira de Física 15, 379
(1985).

HeKr2002 C. Hellaby and A. Krasiński. You cannot get through
Szekeres wormholes: Regularity, topology and causality in quasi-spherical
Szekeres models. Phys. Rev. D66, 084011 (2002).

Lema1933 G. Lemaître, L'Univers en expansion [The expanding
Universe]. Ann. Soc. Sci. Bruxelles A53, 51 (1933); English
translation: Gen. Relativ. Gravit. 29, 641 (1997); with an editorial
note by A. Krasiński: Gen. Relativ. Gravit. 29, 637 (1997).

Tolm1934 R. C. Tolman, Effect of inhomogeneity on cosmological
models. Proc. Nat. Acad. Sci. USA 20, 169 (1934); reprinted: Gen. Relativ. Gravit. 29, 935 (1997); with an editorial note by A.
Krasiński, in: Gen. Relativ. Gravit. 29, 931 (1997).

Kras2014d A. Krasiński, Blueshifts in the Lemaître – Tolman
models. Phys. Rev. D90, 103525 (2014).

Kras2014a A. Krasiński, Accelerating expansion or inhomogeneity? A
comparison of the ΛCDM and Lemaître – Tolman models. Phys.
Rev. D89, 023520 (2014) + erratum Phys. Rev. D89, 089901(E)
(2014).

unitconver Energy and Work Units Conversion, 

http://www.asknumbers.com/EnergyWorkConversion.aspx

Plan2014 Planck collaboration, Planck 2013 results. XVI. Cosmological
parameters. Astron. Astrophys. 571, A16 (2014).

Plan2014b Planck collaboration, Planck 2013 results. XV. CMB power
spectra and likelihood. Astron. Astrophys. 571, A15 (2014).

Elli1971 G. F. R. Ellis, Relativistic cosmology. In Proceedings
of the International School of Physics “Enrico Fermi”, Course 47: General
relativity and cosmology. Edited by R.K. Sachs, Academic Press, 1971, pp.
104-182. Reprinted: Gen. Relativ. Gravit. 41, 581 (2009); with an
editorial note by W. Stoeger, in Gen. Relativ. Gravit. 41, 575
(2009).

wripa E. L. Wright, A Cosmology Calculator for the World Wide Web.
Publ. Astr. Soc. Pac. 118, 1711 (2006).

calculator E. L. Wright, Ned Wright's Javascript cosmology calculator,


http://www.astro.ucla.edu/wright/ACC.html

emspec The Electromagnetic Spectrum, 

https://www.miniphysics.com/electromagnetic-spectrum_25.html

Kras2001 A. Krasiński, The newest release of the Ortocartan set of
programs for algebraic calculations in relativity. Gen. Relativ. Gravit.
33, 145 (2001).

KrPe2000 A. Krasiński, M. Perkowski, The system ORTOCARTAN –
user's manual. Fifth edition, Warsaw 2000.


