Relativity for games

We present how to implement special relativity in computer games. The resultant relativistic world shows the time dilation and Lorentz contraction exactly, not only for the player but also for all the nonplayer characters, who obey the correct relativistic equation of motion according to their own accelerations. Causality is explicitly maintained in our formulation by use of the covariant velocities, proper times, worldlines, and light cones. Faraway relativistic scenes can be accurately projected onto the skydome. We show how to approximate a rigid body consisting of polygons, which is ubiquitous in computer games but itself is not a relativistically invariant object. We also give a simple idea to mimic the Doppler effect within the RGB color scheme.


Introduction
The world is written in the language of quantum mechanics and relativity. They together provide the basis for the Standard Model of particle physics and for standard cosmology; see, e.g., Refs. [1] and [2] for reviews, respectively. For large length scales, 1 classical (i.e. nonquantum) mechanics works well. If all the objects have much lower speeds than the speed of light, nonrelativistic mechanics suffices. 2 When both conditions are met, the ordinary Newtonian mechanics becomes a good effective description of nature. Quantum mechanics is hard to visualize, as it is beyond our ordinary perception. On the other hand, there is no problem in visualizing relativistic effects, which can be significant even within the reach of classical (nonquantum) mechanics. The detection process of light by our eyes does not differ whether the light is emitted from the relativistic world or not. Therefore, there remains room to visualize a classical but relativistic world in computer games if we manage to formulate the implementation of a relativistic world. This is the goal of this work.
Up to now, virtually all games have been based on Newtonian mechanics, with only a few exceptions. In Refs. [6] and [7], a player can traverse a relativistic world with three and four spacetime dimensions, respectively, without any reaction from the nonplayer characters (NPCs). These two are the only games that have so far been found by the authors.
One may also find a good review on the visualization of relativistic effects in computers in Ref. [8]. Further developments after Ref. [8] can be found in Refs. [9,10,11,12,13,14]. In these works, what is implemented are the relativistic effects for a moving observer in a world that does not evolve in time in itself. In other words, one cannot affect the world in this kind of implementation.
In this paper, we formulate how to implement time evolution in a special-relativistic world, namely, interactivity between the NPCs, the player, and all the other objects. In particular, one will see that causality is maintained by heavy use of past light cones and worldlines in four spacetime dimensions. One may find concrete implementation of all the ideas in this paper in a primitive first-person shooter (FPS) game in Ref. [15]. Sample game screens are shown in Fig. 1, in which the player is accelerating in a fixed background without affecting the world. In Fig. 2, the player is fighting against enemies interactively. This paper is organized as follows. In Sec. 2, we review special relativity and spell out our notations. In Sec. 3, we show how to draw the fully special-relativistic world in computer games. Each observed object is drawn at the intersection between the observer's past light cone and the object's worldline. In Sec. 4, we present how to evolve the relativistic world, with an application to an FPS in mind. This is the first time that interactions between participants (namely, the player and the NPCs) are consistently taken into account without contradicting causality. Past and future light cones are heavily used for this. In Sec. 5, we give various ideas useful for a concrete implementation. In Sec. 6, we summarize our result. In Appendix A, we briefly summarize how the energy and momentum form a spacetime vector in relativity. In Appendix B, we show how the Lorentz contraction and dilation occur when a measuring rod is seen by an observer. In Appendix C, we present another formulation for the time evolution of the objects on the player's past light cone. This co-moving evolution is more The 3D world is gridded by the white lines, which are Cartesian coordinates in the reference frame introduced in Sec. 3. One can see that the coordinates appear bent as the player's velocity becomes sizable. The Doppler effect is visible such that the color of the Galactic Center is blue-shifted as the player accelerates. On the other hand, one can see that the backward scene, which is squeezed into the player's sight due to the so-called relativistic aberration, is red-shifted on the lower-right panel. Note that the implemented Doppler effect does not represent the real physical shift of the spectrum but is a mimicked one; see Sec. 5.4. We see that, as the player accelerates, the player's "proper time", corresponding to τ in Eq. (40), becomes delayed compared to the "world time" of the background at rest. ("FPS" in the second line on the screen denotes "frames per second" rather than "first-person shooter".) Figure 2: Game screens from the sample game in Ref. [15]. The player is shooting the yellow bullets. Both the enemy and bullet are assumed to emit light of given colors by themselves, rather than reflecting light from some source such as the Sun. The player's velocity is faster on the right so that the Cartesian coordinates, explained in the caption of Fig. 1, appear more bent. On the right, the enemy is shooting bullets towards the player according to the aiming algorithm in Sec. 4.3; we can also see the red shift on the enemy's bullets that come in the player's sight due to the relativistic aberration. We note again that the implemented Doppler effect does not represent the real physical shift of the spectrum but is a mimicked one; see Sec. 5.4. The time evolution of the player and enemy is determined as described in Sec. 4: The enemy determines its move according to the information on its own past light cone as in Sec. 3.1, and estimates others' move projected on its own future light cone as in Ref. 4.3. The larger and smaller square cursors, which are connected by the straight blue lines, indicate the enemy's position on the player's past light cone in Eq. (68) and the enemy's projected direction on the future light cone in Eq. (92), respectively. suitable for analytic treatment, and some of the results are shown too. In Appendix D, as a concrete example of the relativistic acceleration, we show how the rocket propulsion can be formulated. In Appendix E, we review the use of quaternions in order to handle the rotation in three spatial dimensions.
We have tried to make this paper more pedagogical and self-contained than ordinary physics papers, even for some elementary mathematics, since the subject is more interdisciplinary and we expect more computer-oriented readers.

Special relativity
We spell out our basic notations to describe a special-relativistic world.

Coordinates, vectors, worldlines
Usually computer games have d = 2 or 3 spatial dimensions. We write a position vector in the d-dimensional Cartesian (orthonormal) coordinate system as 3 where the d vector is written in the component and matrix notations in the second and third expressions, and we use curly and square brackets for them, respectively, unless otherwise stated. We generally write a d vector with a bold-style letter such as x. In the matrix notation, the transpose is given by The location of each particle in the world is represented by its position vector x at a given time t. In practice, in computer games, the "particle" may stand for a vertex in a polygon. The time dimension t is written as the 0th coordinate 4 where c = 299 792 458 m/s is the speed of light, which is a physical constant in special relativity. 5 We can always take the natural units c = 1, which we employ in this paper unless otherwise stated. In the natural units, time and length have the same dimension: 1 s = 299 792 458 m. When necessary, we can always come back from the natural units by recovering c by dimensional analysis. A unit length in computer programming is quite arbitrary, and we call it a "pixel" here just for concreteness. We may identify a pixel with an arbitrary physical length ; namely, is given in the units of m/pixel. Then the speed of light in the program is which has units of pixel/s. 6 In practice, it is most convenient to store all the time information in pixels too, by using the speed of light (4) in Eq. (3). Then the velocity and acceleration, which will be presented in Eqs. (41) and (53), are given without unit and with units of pixel −1 , respectively. 7 The spacetime coordinates are collectively written as a D-dimensional vector with D = d + 1: where we have again used the component and matrix notations in the second and third expressions in each line, respectively. We place the long arrow on top of any D vector, such as − → x , throughout this paper. The transpose in the matrix notation is In ordinary Newtonian mechanics, the trajectory of each particle from a time t ini to t fin can be written as This trajectory can be viewed as a worldline in D-dimensional spacetime: where − → x is the D vector given in Eq. (5) and the real number s parametrizes the worldline. 8 6 One may instead take whatever value of the speed of light in the unit of pixel/s, regardless of the choice of , if one simulates a different world with a different speed of light, as discussed in footnote 5.
7 When one only uses the acceleration to change the velocity, as in Secs. 4.1 and 4.2, it would be more practical to give the proper-time difference ∆τ in units of seconds. Then the time x 0 (as well as the position − → x ), the velocity − → u , and the acceleration − → a would be given in units of pixels, dimensionless, and in the units of s −1 , respectively. In that case, the time evolution in the nonnatural units becomes schematically where c is given in units of pixel/s as in Eq. (4).

Newtonian mechanics
In Newtonian mechanics, let us consider a particle at time t located at x(t), having a (nonrelativistic) velocity The particle's position after an infinitesimal time ∆t is given by or, to be more explicit, The velocity v at each moment can be chosen quite arbitrarily in nonrelativistic computer games. If we use Newtonian mechanics, v(t + ∆t) is determined from v(t) for a given input of the (nonrelativistic) acceleration a NR (t): where the acceleration at t is determined by Newton's equation of motion: The explicit form of the force F depends on the dynamics that you choose. In computer games, the trajectory and worldline become discrete sets of d and D vectors, respectively, due to discrete time steps:

Lorentz transformation
In spatial d dimensions, the ordinary inner product of a pair of d vectors, say x = x 1 , . . . , x d and y = y 1 , . . . , y d , is given by We may also write in the matrix notation. Note that We define the Lorentzian inner product of the D vectors − → x = x 0 , x 1 , . . . , x d and − → y = y 0 , y 1 , . . . , y d by 9 Note that − → x · − → y = − → y · − → x . We also use the notation interchangeably.
We also define the Lorentzian norm: 10 where is the ordinary d-dimensional (squared) norm; we may also write x 2 = x T x in the matrix notation; we also write Note that the Lorentzian norm (22) can be negative. Any D vector − → x is called spatial (spacelike), null (light-like), and temporal (time-like) when [ − → x ] 2 is positive, zero, and negative, respectively: It is straightforward to show that the Lorentzian norm behaves in the same way as the ordinary square for a sum of any D vectors: Let us consider a linear coordinate transformation Λ: 10 The extra square bracket is placed for a reader unfamiliar with relativistic notation: Do not confuse the Lorentzian norm [ − → x ] 2 with the 2-component of − → x which is denoted by x 2 . 11 More generally, we impose the invariance under the Poincaré transformation: The addition of a constant vector − → b is the translation, a constant shift of the coordinate origin, which is trivially realized and is not treated in this paper. A reader unfamiliar with the translation should only memorize the fact that any subtraction of two coordinates is invariant under any translation: where we have also shown more explicit component expressions in the second and third lines; Λ is a matrix, 12 Λ = Λ µ ν µ,ν=0,1,...,d .
Note that all the D vectors are transformed by Eq. (27) simultaneously: In particular − → x and − → y are transformed by the same Λ.
The transformation (27) is called the Lorentz transformation if it does not change the Lorentzian inner product (18): which is satisfied when and only when 13 where 14 This can be shown as follows: where the matrix notation is used in the second step as in Eq. (20). That is, Λ is a Lorentz transformation when and only when Eq. (30) is satisfied. It is important that any Lorentzian norm, say [ − → x ] 2 , is then Lorentz-invariant: Physically, we are looking for linear transformations that leave the speed of light invariant. If the speed of light is invariant, then a spherical wavefront of light that is emitted from the same point should be transformed to a spherical one. The wavefront of a spherical wave of 12 A reader unfamiliar with the distinction between the upper and lower indices does not have to be bothered by it. In this paper, we do not raise and lower the time index 0. We can freely raise and lower the spatial indices 1, . . . , d under the metric convention (21). That is, Λ i j = Λij = Λ ij for i, j = 1, . . . , d. 13 The condition (30) is analogous to the condition for the spatial rotation R T R = I, which is deduced from the rotational invariance x · y → x · y = x · y under the coordinate transformation x → x = R x. 14 In general, Λ has D (D − 1) /2 = d (d + 1) /2 degrees of freedom, and can be parametrized as where L is a boost matrix and in which R is an ordinary rotation matrix in d spatial dimensions and 0 T = 0 · · · 0 . We note that R obeys where I is the d-dimensional identity matrix; see also footnote 13.
= d degrees of freedom in Λ are in the boost matrix L. A concrete parametrization of L will be given in Sec. 2.5 after we introduce more physics.

Proper time and Lorentz-covariant velocity
Special relativity is defined as a theory that is invariant under Lorentz transformations. Therefore, it is important to write the theory in terms of Lorentz-invariant quantities. A Lorentz-invariant quantity can be formed as a Lorentzian inner product (20) or a norm (22). They are composed from D vectors, which are Lorentz covariant. Therefore, it is important that things, such as the velocity and acceleration, are written as Lorentz-covariant D vectors.
We start from a coordinate system { − → x }. Suppose that an infinitesimal shift of the worldline parameter s in Eq. (8) changes the position of the particle in D spacetime dimensions: We can define the displacement D vector: 15 The dilatation − → x → b − → x , with b = 0 being a constant, also leaves the condition [ − → x ] 2 = 0 unchanged. The invariance under dilations, called the scale invariance, is violated by the conformal anomaly in quantum field theories in general. Otherwise, the dilatation is excluded from the symmetry of the system of our consideration under the following two assumptions [18]: First, a boost along a direction with a given velocity v, combined with another boost along the same direction but with the negative velocity −v, becomes an identity transformation. Second, the length of a measuring rod perpendicular to the boost direction becomes identical to each other when its boost is with the velocity v and −v. i.e., There are three possibilities for the particle's move in classical physics: • An ordinary massive particle always goes time-like: − → ∆x 2 < 0, where the square denotes the Lorentzian norm (22).
• A massless particle, such as a photon of which light consists, always goes light-like: The existence of a tachyon indicates a pathology of the system, and violates causality in general. 16 Therefore, you must be careful in introducing a tachyon in your computer game.
Hereafter, we focus on massive and massless particles. The time difference in this step is ∆x 0 . We assume that none of the massive and massless particles go backward in time. 17 That is, we assume that ∆x 0 > 0. For a massive particle, we define its proper time τ such that it increases by an amount for an infinitesimal evolution − → ∆x. It is important to note that the proper time is, by definition, manifestly Lorentz invariant. 18 Note that the proper time is defined for each particle individually. Now we can define the D velocity for the massive particle: 16 In quantum field theory, the existence of a tachyon indicates that one is on a false vacuum, which will eventually roll down to a true vacuum that does not have a tachyon.
17 In quantum field theory, a particle going backward in time is identical to its antiparticle going forward in time. 18 For massive particles, we may e.g. use τ as the worldline parameter s. Instead, we may also choose an arbitrary monotonically increasing function of τ as s; see footnote 8. where ∆τ is given in Eq. (40). By definition, − → u is a Lorentz-covariant D vector. 19 That is, for where we have used the Lorentz invariance of ∆τ defined in Eq. (40). Note that, by definition, We hereafter choose the spatial component u as the d independent parameters, and then the time component of the D velocity is not independent: The evolution (37) for a small proper-time step ∆τ now reads We note that the ordinary Lorentz-noncovariant velocity v in this particular coordinate system is given by where ∆x and ∆x 0 are given in Eq. (39) and, in the last step, we have tentatively come back from the natural units for a reader unfamiliar with dimensional analysis; see Eq. (3). Using this Lorentz-noncovariant velocity, one can show that the spatial component of − → u becomes i.e., There is a one-to-one correspondence between u and v. We see that for a massive particle, the possible values of v are |v| < 1, or |v| < c when we recover c. On the other hand, |u| can be arbitrarily large. When the particle's velocity is much smaller than the speed of light |v| c = 1, we see that the covariant and noncovariant velocities match each other: where O( n ) is the Landau symbol representing the terms of order n and higher as → 0. From Eq. (40), we can write the infinitesimal time difference ∆x 0 in this particular coordinate system as In terms of u, the relation (50) can be written as where the Lorentz factor γ is defined by Note that γ(−u) = γ(u) and that γ(0) = 1. Hereafter, v does not play any role in the actual formulation of the relativity since v is a Lorentz-noncovariant quantity, and we simply call u (rather than v) the velocity unless otherwise stated. We also define the D acceleration, where ∆τ is given in Eq. (40). 20 Note that, by taking the and hence in any coordinate system. In particular, a 0 = 0 whenever u = 0.

Rest frame
Let us show that we can always take a coordinate system in which a particle appears to be at rest. Suppose that a particle has a velocity u in a particular coordinate system { − → x } at a given single moment when its proper time is τ . Then the following Lorentz transformation makes − → u → (1, 0): whereû := u/ |u|. 21 More explicitly, One can verify that L(u) indeed transforms − → u to rest, that L(−u) is the inverse matrix of L(u), and that L(u) satisfies the condition for the Lorentz transformation (30), respectively: In the new coordinate system namely the rest frame for the particle at the proper time τ , we have U = 0 at that moment. Hereafter, we write quantities in someone's rest frame in upper case in general, unless otherwise stated. During the subsequent infinitesimal time evolution, we get ∆τ = ∆X 0 because γ(0) = 1. That is, the proper-time flow of the particle is nothing but the time flow in each rest frame in each time step. Therefore, the proper time is the time felt by the particle itself. Now we can interpret from Eq. (51) that the time difference ∆x 0 in an arbitrary frame (other than the rest frame) is always larger than the proper time ∆τ felt by the moving particle itself. 22 This phenomenon is called time dilation.
Einstein's equivalence principle asserts that the physical law for a moving particle is the same as the one that is Lorentz-transformed from the law of the particle at rest. That is, one can write the relativistic equation of motion just by replacing dt by dτ : 23 where − → f is nothing but the D vector transformed from the rest-frame force felt by the particle: which reduces to the celebrated form appearing in Einstein's original paper [18], when written in terms of the noncovariant velocity v := u 1 / 1 + (u 1 ) 2 . 22 The time difference ∆x 0 is not an actual time difference seen by any observer. We will see that the actual time difference for an observer is the one deduced from the time foliation by the observer's past light cones rather than by an equal-time slice. See Sec. 3.1.
23 As both sides are covariant, the equation of motion becomes invariant under the Lorentz transformations: Note that the 0th (time) component of the D force in the rest frame must be zero in order for the equation of motion (59) to hold; see the discussion after Eq. (55). That is, the 0th component of the equation of motion does not have independent information.
Conversely, for a given acceleration A in the rest frame, we obtain the D acceleration in any frame in which the particle has the velocity u: That is, In particular, if there is no acceleration in the rest frame A = 0, the acceleration vanishes in any frame a = 0. We may also show that, when a = 0 in a frame, then the D acceleration vanishes in any frame. 24 We note that it is perfectly legitimate to work within the framework of special relativity in order to handle an acceleration unless it is due to gravitational interactions. When and only when the gravitational interactions are strong enough, we need general relativity for a full treatment.

Drawing the world
We show how to draw the relativistic world for a given set of worldlines of all the objects. A schematic figure for this section is presented in Fig. 3.
First we arbitrarily choose a reference frame { − → x } and store all the information as written in this coordinate system. In practice, it is convenient to choose a reference frame such that background objects are at rest in the frame. 25 Hereafter, we write the reference-frame quantities in lower case, unless otherwise stated.

Drawing the world on the observer's past light cone
For an observer at a spacetime point − → x O , the hypersurface consisting of possible light rays that can come to − → x O is the past light cone (PLC): 24 Suppose that a = 0 in a frame: − → a = a 0 , 0 . We write the velocity in this frame u. Then the rest-frame acceleration reads Since γ(u) ≥ 1, a 0 = 0 results from the time component of this equation. This means − → a = 0. Any Lorentz transformation of the zero vector is zero: Λ − → 0 = − → 0 . Therefore the D acceleration is zero in any frame if a = 0 in a frame. 25 In our real Universe, we can define the absolute rest frame such that the dipole component of the cosmic microwave background vanishes; see Ref. [19] for latest observational results on the cosmic microwave background. For later use, we also define the future light cone (FLC): We say that a spacetime point − → x is on the past side of PLC( − → nor on its past side. The observer can see the world sliced by its PLC( − → x O ), i.e., the observer can know things on PLC( − → x O ) and its past side. It is important that the PLC is Lorentz invariant in the sense that, if a particle is on the PLC, it remains so in any Lorentz-transformed frame. 26 This is also the case for the notion of the past and future sides of the PLC, and similarly for the FLC and both its sides.
Let W n be the worldline of the nth particle: where τ ini n and τ fin n are the initial and final proper times of the worldline. The observer at − → x O sees this particle at the intersecting point between PLC( − → x O ) and W n , which we write as hereafter. How to determine the intersection in a concrete implementation will be explained in Sec. 3.3.
An observer having a velocity u O in the reference frame { − → x } sees the world in the observer's rest frame: Each of all the other particles, say the nth one, is seen by the player as located at where is seen by the player as These are all needed to draw the world fully relativistically.

Discrete worldline
The program stores the worldline of each particle in the reference frame { − → x }. Due to the iterations, each worldline becomes a discrete set of its past spacetime points, just as in the Newtonian case (15): where s t (t = 1, . . . , N W ) is the parameter of the worldline, which can be identified as the proper time of the particle-we have abbreviated them as − → x t := − → x (s t )-and we always order from past to future: where

Intersection between worldline and PLC
As discussed above, it is important to compute the intersection between a worldline and a light cone. Let us spell out the method to obtain the intersection between the worldline (70) and PLC( − → x O ): The first point − → x t that violates this condition is the point that is closest to PLC( − → x O ) in W on the future side of PLC( − → x O ).
27 If one wants to place the observer at the origin, one may take, so to say, a central frame • Between the adjacent points − −− → x t −1 and − → x t , the worldline is obtained by linear interpolation: where 0 ≤ σ ≤ 1.
i.e., 29 where 30 This value of σ gives the intersecting point − →

Short summary
At each moment, we draw all the objects as if each of them is placed at the spatial position . Once all the spatial positions are given, concrete implementation of the drawing of the world is the same as in ordinary nonrelativistic 3D games. An approximate treatment of a rigid body will be presented in Sec. 5.3.

Time evolution
We show how to take into account the relativistic time evolution of the system, with an application to a first-person shooter (FPS) in mind. As reviewed above, a game in Newtonian mechanics draws the world at each equal-time slice. One problem is that such a time slice is not Lorentz invariant and is not compatible with relativity. Instead we employ the foliation of spacetime by the PLCs of the player, which is Lorentz invariant, as reviewed above. A schematic diagram is shown in Fig. 4.

Player's time evolution
What should the time evolution be? At each time step, the future-most surface of the world that is perceivable by the player is PLC( − → x P ), on which the particle's latest position − → x is located. In the next iteration, the game program calls the real time passed, ∆t, which is 29 The other point σ = β + β 2 − αγ /α gives the intersection with the future light cone. 30 be careful with the abuse of notation: γ here has nothing to do with the Lorentz factor.
Player time space W n : Worldline of nth particle Figure 4: Schematic diagram in 1 + 1 spacetime dimensions for the time evolution of the player and the nth particle under the PLC foliation.
identified as the player's proper time passed, ∆τ P . In the iteration, the player's proper time, position, and velocity move to, respectively, where A P is the acceleration of the player given at its rest frame, determined by in which F P is the ordinary Newtonian force felt by the player in its rest frame, coming both from the player's own input and from influences from its environment and others. Concrete examples of how to give A P will be shown in Sec. 4.5.
We update the values of the player's proper time, position, and velocity as in Eq. (79). Now the player's new position − → x P is determined, and we know the new PLC − → x P which is used to draw the world in the next iteration.
We add the new point to the last of the player's worldline set, as described in Sec. 3.2. The player may have also done some action at − → x P before its move, e.g., have shot a beam from − → x P , which will be treated later in Sec. 4.4. nth NPC (= player in this case)

Others' time evolution
The nth particle at − → x n PLC( − → x P) determines its own move, and extends its worldline W n , according to its own acceleration and influences from others, until the extended W n hits the new PLC − → x P . The new intersection point is drawn in the next iteration as above. Now we show how to determine the moves of NPCs. Suppose that an NPC, say the nth one, is at − → x n PLC( − → x P) , and has a velocity − → u and a proper time τ , all in the reference frame.
The NPC can see the world sliced by its own PLC( − → x n ), or, more precisely, can know things on PLC( − → x n ) and its past side. According to the available information, we determine the NPC's acceleration − → a n in the reference frame; see Sec. 4.5 for details. The NPC's move is where ∆τ n can be chosen to be ∆τ P , or else whatever (fixed or variable) value according to the NPC's reflexes. In general, the resultant − → x n + − → u n ∆τ n is still on the past side of PLC − → x P .
We iterate until it goes beyond PLC − → x P , namely, until the following condition is violated: In each step, the NPC may also do some action, as will be described in Secs. 4.4. Note that the NPC at − → x n sees the world cut by its own PLC( − → x n ), and cannot see further moves of others beyond it: Any other NPC, say the mth one, is seen by the one at − →

Aiming
We show how an observer O can estimate the future position of other objects. A schematic diagram is shown in Fig. 5. The observer at − → x O sees the world sliced by its own PLC( − → x O ), on which the perceivable future-most position of another object, say the nth NPC, 31 is located . From the given information, O wants to predict the nth NPC's future move.
Let − → u n be its velocity at − → x n PLC( − → x O) . We want to know the intersection between the FLC of O, and the would-be worldline of the nth NPC extended in the direction of − → u n : The value of τ at the intersection is obtained from where we have used the PLC condition [ − → x n − − → x O ] 2 = 0 and Eq. (43) in the last step. That is, the projected location of the nth NPC onto O's FLC is given by namely, at This point is nothing but the intersection with the FLC mentioned in footnote 29. The observer will see this coordinate in its rest frame at

Shooting
We continue the discussion from the previous subsection. Now O wants to shoot a beam or bullet, B, at the nth NPC from the location − → x O , when its velocity is − → u O . In O's rest frame { X }, we define the "future beam/bullet cone" (FBC): x O is the location of O in its rest frame, V is the speed of B satisfying 0 < V ≤ 1, andD is a possible direction in which B can be shot.
In O's rest frame, the would-be worldline of nth NPC (88) becomes From Eq. (95), we get and hence − → where we have defined the following (linear and bilinear) operations: Note that this operation is not Lorentz invariant. Solving Eq. (98) as in Sec. 3.3, we get 32 It is possible that the would-be worldline (94) does not intersect with FLC(X n ), in which case the expression inside the square root in Eq. (101) becomes negative. In such a case, O may still shoot in the direction given in Sec. 4.3, which always exists. O can shoot a beam/bullet whose worldline is given in O's rest frame by in whichD In the reference frame, the beam/bullet worldline (102) becomes where Do not confuse it with the norm (99): the distinction is made by the comma inside.

Acceleration
Now we exhibit how O's acceleration, used in Sec. 4.2, is determined. There can be several sources for O's acceleration: • To match ordinary human common sense, one may suppose that the spacetime is filled with (fictitious) air or fluid that is static in the reference frame, and may assume a friction force that is a function of the velocity of each particle in the reference frame: where f (u) > 0 can be any function of u, which is typically ∝ |u| and u 2 for the laminar and turbulent drags, respectively. Note that the time component can always be • O can accelerate in a direction that O wants. That is, O may self-accelerate by a magnitude A s in a directionD in its rest frame. Then O's self-acceleration in its rest frame is given by and in the reference frame by See also Appendix D for a more elaborate rocket propulsion.
For example, O may want to accelerate in the projected direction − → d of the nth particle on O's FLC given by Eq. (91): The direction of acceleration,D = D/ |D|, in O's rest frame can be obtained from the spatial component D of the vector • We show how to implement a collision, namely, a repulsive force exerted by the nth particle onto O, without contradicting causality. Suppose that the nth particle is on PLC( − → x O ), and we write its and O's locations as − → X n and − → X O , respectively, in the nth particle's rest frame. We write the corresponding velocities as − → U n = 1 0 and − → U O = L(u n ) − → u O , respectively. We write the displacement vector from n to O in n's rest frame as: We assume that the nth particle exerts a strong repulsive force −−−→ F n→O when the distance in its rest frame is smaller than its typical radius R n : where in which F > 0 is a constant and is the Heaviside step function. If this repulsive force is the only force acting on O, it Then the acceleration from this repulsive force is, in the reference frame, We sum up − −− → a n→O coming from all the particles on PLC( − → x O ).
In the end, O's total acceleration is given by the vector sum which is then put into Eq. (85).

Miscellaneous ideas
We show various ideas to make more realistically relativistic scenes.

Faraway scenes on the skydome
One puts a sufficiently large sphere that surrounds the player in the player's rest frame, and draws a background texture of the faraway scene on it: At an angle (Θ, Φ) on the sphere, we draw the corresponding texture at (S, T ) in the texture coordinates. Let us see how (S, T ) is determined for a given direction (Θ, Φ). Let { − → x } be the reference frame in which faraway background objects are at rest. Suppose that the player and a background object O are at − → x P and − → ξ , respectively, in this frame. In this subsection, we parametrize O's position by Using ordinary polar coordinates, the intersection between O's worldline and the player's PLC can be parametrized as where we have assumed that O is so far away that player's position − → x P can be regarded as being at the spacetime origin: Conversely, we may write the zenith and azimuthal angles as We may use texture mapping of the background image using the equirectangular projection. For example, in OpenGL, the background scene at the angle (θ, φ) in the reference frame corresponds to the texture coordinates (S, T ): where 0 ≤ S ≤ 1 and 0 ≤ T ≤ 1. Let us consider light that comes into the player's eyes from the (Θ, Φ) direction in the player's rest frame. Let − → Ξ be the position of the faraway object seen by the player in the player's rest frame: We parametrize it as where R is very large. 33 How is (Θ, Φ) related to (θ, φ), used to pick up the drawn texture point (S, T ) in Eq. (124)? Equation (125) implies that For a given direction (Θ, Φ), we determine − → Ξ by Eq. (126) up to the overall constant R, and then obtain − → ξ by Eq. (127). Once − → ξ is known, (θ, φ) are given by Eq. (123), and then the corresponding texture coordinates (S, T ) to be drawn are obtained by Eq. (124). Note that the overall normalization R drops out of the final expression since the coordinates − → Ξ always appear as ratios in Eq. (123). In practice, we may put R = 1.

2D background
2D games can be handled similarly. We explain how to draw a background using orthogonal projection in the 3 = 2 + 1D spacetime. We assume for simplicity that the player's view does not rotate with respect to the background. Suppose that the background is at rest at the reference frame { x }, where x = x 0 , x = x 0 , x 1 , x 2 is a 3(2 + 1)D spacetime vector, whose spatial components consist of a spatial vector x = x 1 , x 2 . 34 Let (S, T ) be the texture coordinates of the background, with 0 ≤ S ≤ 1 and 0 ≤ T ≤ 1. Suppose that the background texture is mapped onto a region x 1 min ≤ x 1 ≤ x 1 max and x 2 min ≤ x 2 ≤ x 2 max : When the player's velocity is u P in the reference frame, we want to know which point in the texture coordinate (S, T ) is picked up for a given point in the player's rest frame { X }: where we have chosen the origin of X to be the point corresponding to x P in the reference frame.
The player sees the world sliced by its PLC. Therefore, a natural way to draw the 2D world is to project the player's PLC onto the X 1 -X 2 plane. Let the drawn region on screen be X 1 min ≤ X 1 ≤ X 1 max and X 2 min ≤ X 2 ≤ X 2 max . A spatial point X = X 1 , X 2 in this region has the time coordinate on the player's PLC. The corresponding point in the reference frame is From this, we can read off the values of the spatial components x 1 , x 2 , which can be put in Eq. (128) to get the corresponding texture coordinate (S, T ). By this PLC formalism, one may improve the drawing of the world by an equal-time slice in the player's rest frame, employed, e.g., in Ref. [6].

Approximately rigid body
In the relativistic world, it is impossible to have an exactly rigid body, as it violates the locality. However, a thing is usually introduced as a rigid body, described by a set of polygons, in computer games. We show how to treat such a thing approximately.
Suppose that a body B has N vertices that are specified by spatial d vectors where a = 1, . . . , N . Each d vector ξ a specifies the position of the corresponding vertex measured from B's reference point, somewhere near its center. Hereafter, we call this point, somewhat sloppily, B's center. Let R t be a d × d rotation matrix that represents the orientation of B, which differs at each spacetime point − → x t that consists of the worldline { − → x 1 , − → x 2 , . . . }. 35 At each worldline point − → x t , the object B has the velocity (71) that is uniform for all the polygon vertices (132): B's position in its rest frame is given by More precisely, − → X t represents the position of B's center in its rest frame, at the tth point in B's worldline.
We approximate that B is a rigid body, i.e., R t rotates all the vertices (132) simultaneously: The position of each vertex, say the ath one, in B's rest frame is given by when B is at the jth point on its worldline. How can one determine its time coordinate X 0 a,t ? We approximate that each worldline of the vertex is parallel in spacetime to that of B's center: Then we can obtain the intersection between the worldline of the ath vertex and the player's PLC from the condition i.e., where X P is the player's position in B's rest frame. In the original reference frame, the position of the ath vertex is where the time and spatial components of − − → X a,t are given by Eqs. (134) and (136), respectively. This x a,t can be transformed to that in the player's rest frame as in Sec. 3.1 when the world is drawn.
Rotation of the rigid body, or, more explicitly, the information on R t above, can be most easily handled by using quaternions; see Appendix E for a review. We specify the rotation of a rigid body from its basic position by an angle θ around a directionˆ using the quaternion Θ ˆ , θ := cos where 0 ≤ θ ≤ π and ˆ = 1. There are three independent degrees of freedom in total in θ andˆ , which coincide with the physical degrees of freedom, namely, a unit vector to specify the front (upward, or whatever) direction, and an angle to determine the rotation around that direction. 2D games can be handled similarly. Usually an object is drawn by a rectangular picture. We may divide the picture into N × N cells, where N and N are appropriate numbers, typically of order 10 to 100. Then we can treat the N N vertices as above, and map the texture on each cell.

Doppler effect
Light that comes into our eyes has a spectrum; namely, its intensity I is a function of the angular frequency ω: where ω is related to the frequency f and the wavelength λ by Here and hereafter, we recover c from the natural units in the expressions involving ω. When light ray points in a directionn, where |n| = 1 is a unit d vector, it is known that the wave D vector 36 is proportional to the photon's D momentum and that − → k transforms as a D vector under the Lorentz transformation; see Appendix A.
Suppose that a source S and an observer O have their velocities − → u S and − → u O , respectively, in the reference frame. Let us consider light of an angular frequency Ω pointing towardsN with N = 1 in S's rest frame: The directionN can be obtained as follows: Then we obtainN The light (142) is observed by O as the wave vector The angular frequency Ω observed by O is then the time component of − → K , namely K 0 . Note that Ω is proportional to Ω: 36 In an nonrelativistic context, the d vector k/c = (ω/c)n is usually called the wave vector.
where C is a constant that can be computed from the above procedure. Then the spectrum (139) is perceived by O as a new spectrum I O : This is all needed to determine the Doppler effect on light.
Most fundamentally, one should consult the above procedure. In practice, it does not work. Why? In the RGB color model used in computers, a color is specified by three numbers (R, G, B) that roughly represent the perception of three types of cone cells in our eyes. For a given spectrum I(ω), the values of R, G, B are obtained, up to an overall normalization of (R, G, B), by where the CIE standard observer color matching functions r, g, and b are given in Ref. [20] and λ min = 380 nm and λ max = 780 nm in the CIE 1931 color space. If I(ω) were given to specify a color, instead of (R, G, B), then we could calculate the new values of (R , G , B ) for O: The problem is that one needs the full spectrum I(ω) for the exact relativistic computation, whereas one cannot know the functional form of I(ω) from merely the three numbers R, G, and B.
Here we present a naive procedure to mimic the Doppler effect. For blackbody radiation, (R, G, B) is known as a function of the temperature T . That is, we know the functions for the blackbody radiation; see, e.g., Ref. given set of (R, G, B), we may then naively compute color-specific temperatures (T R , T G , T B ) by inverting Eq. (155). 38 Physically, a temperature roughly behaves as an energy, which is the 0th component of the energy-momentum D vector, proportional to the wave vector (163) for photons. Therefore, the new color temperature for O would be estimated by T /C for each color. Then one may get the naive Doppler-shifted colors as We have shown how to mimic the Doppler effect within the RGB color scheme. As said, the full treatment requires the inclusion of the spectrum of the light for each pixel, instead of the current approximation using the three RGB numbers. That will allow visualization of the currently invisible ultraviolet and infrared lights when one is boosted significantly by a Lorentz transformation.

Summary
We have shown how to implement special relativity in computer games. The past light cone is used for the foliation of the world without violating causality. This is the first realization of interacting nonplayer characters that perceive the world and react to others, both relativistically without violating causality, under fully relativistic time evolution. In the formulation, the notion of the Lorentz-covariant velocity is extensively used, instead of the more widely used noncovariant velocity. We have shown several ideas to approximate the relativistic world such as the relativistic generalization of the skydome, rigid body, and Doppler effect.
You may find an implementation of the ideas presented in this paper in Ref. [15].
which is a fixed constant for each species of particle. 39 The D velocity is then given as − → u = − → p /m for a massive particle.
The time component of − → p is identified as the energy of the particle: From the definition of mass (158), we obtain When the particle is at rest p = 0, its energy reduces to its mass, and we get Einstein's celebrated formula for the rest energy where we have tentatively recovered the speed of light c. Note that the ratio of momentum to energy gives the ordinary noncovariant velocity in any coordinate system. For a massless particle such as a photon, of which light consists, the velocity − → u cannot be defined, and rather its momentum − → p is a suitable physical quantity. In quantum mechanics, where 10 −34 Js is the reduced Planck constant and, tentatively recovering c from the natural units, in which ω and k are the angular velocity and the wave vector, respectively. Note that where f is the frequency and λ i is the wavelength for the ith direction. Note also that the photon momentum D vector is light-like: The relation (163) is used when taking into account the Doppler effect in Sec. 5.4. 39 We have assumed that − → p is either time-like [ − → p ] 2 < 0 or light-like [ − → p ] 2 = 0. A particle having space-like momentum [ − → p ] 2 > 0 moves faster than the speed of light, and is the tachyon discussed in footnote 16.

B Lorentz contraction and dilation
We show that a measuring rod moving towards the observer looks longer than when it is at rest, because of the Lorentz "contraction", due to the PLC foliation. For simplicity, we work in the D = 1 + 1-dimensional spacetime with only d = 1 spatial dimension. We consider a measuring rod that has a length in its rest frame. Suppose it stays at rest forever in its rest frame. Then its "worldsheet" can be written in the rest frame as We call X 1 = 0 and the left and right ends of the rod, respectively. Suppose that the rod has a velocity − → u = γ u 1 , u 1 in the reference frame, where the Lorentz factor is Here and hereafter in this section, we do not make the unnecessary distinction between u and u 1 as we work in the d = 1 spatial dimension. See Fig. 6 for a schematic plot. The reference frame − → x is obtained from the rod's rest frame − → X by where We see that the time coordinate of the rest frame is written as X 0 = γ u 1 x 0 − u 1 x 1 and hence where v := u 1 /γ u 1 is the noncovariant velocity. In the reference frame { − → x }, the rod appears to move with the (noncovariant) velocity v, and its length measured in x 1 is contracted to be /γ u 1 as X 1 varies from 0 to . That is, the worldsheet of the rod can be written in the reference frame as where δx 1 := X 1 /γ u 1 . This is the celebrated Lorentz contraction. However, the contracted length /γ u 1 is the length at an equal-time slice x 0 = const., which is not really a length measured on any observer's PLC. The equal-time slices are represented by the blue dashed horizontal lines in Fig 6. To see the real measured length, let us prepare an observer staying at rest at x 1 = 0 in the reference frame. We restrict ourselves to the case of a right-moving rod u 1 > 0 without loss of generality. The observer sees the rod from − → x O = x 0 O , 0 . When the rod is moving away from the observer, x 0 O > 0, the left and right ends of the Therefore, the observed length of the rod is Using the identity γ(u) − u = (γ(u) + u) −1 , we see that this is shorter than the Lorentz contraction /γ u 1 . When the rod is moving towards the observer, x 0 O < − /u 1 , the left and right ends of the rod intersect with PLC( − → x O ) at − → x L and − → x R : The observed length of the rod is where we have used the identity 1/ (γ(u) − u) = γ(u) + u. We see that the observed length is longer than when the rod is moving towards the observer.
To summarize, when we write the departing velocity of the rod u dep , which is u 1 and −u 1 in Eqs. (174) and (176), respectively, the apparent length of the rod is In the left panel of Fig. 7, we show this factor as a function of u dep . We also show it in the right as a function of the noncovariant departing velocity v dep := u dep /γ(u dep ).

C Evolution in co-moving frames
To cultivate physical intuition, it is instructive to formulate how the player, in its rest frame at each moment, sees the motion of the others that are always on the player's PLC. The formulation in Sec. 4.2 is better than the one in this section in the sense that it takes into account the time dilation more accurately in the discretized time steps, while the formulation in this section makes the time evolution more tractable in an analytic treatment.
Suppose that the player's position and velocity at its proper time τ P are − → x P and − → u P , respectively, in the reference frame. After an infinitesimal time ∆τ P felt by the player, we shift the player's proper time, position, and velocity as 40 where − → a P is the player's acceleration at the spacetime point − → x P in the reference frame. In the PLC foliation, all the other objects are located on PLC( − → x P ): The nth-particle's position − → x n satisfies x 0 n < x 0 P and In the iteration (178), the nth particle moves onto where − → a n is the acceleration of the nth particle at − → x n in the reference frame. As the player has moved to − → x P , the hypersurface seen by the player in the next step is where we have used [ − → u n ] 2 = −1 and Eq. (179), as well as [ − → u P ] 2 = −1, in the first and second steps, respectively. Note that − → x n − − → x P , − → u n > 0. In the limit ∆τ P → 0, we get where ≈ indicates that the quadratic and higher-order terms of the infinitesimal quantities are neglected. 41 As in footnote 27, let us introduce the central-frame coordinates − → X in which the player at its proper time τ P looks at rest at the spacetime origin: − → X P (τ P ) = 0 and − → U P (τ P ) = 0. The corresponding central-frame position vector − → X to the reference-frame one − → x is The central frame differs at each proper time of the player τ P , and we call it the C(τ P ) frame; we also call them collectively the co-moving frames. In particular, the player is always at rest at the origin in the co-moving frames: − → X P = (0, 0) and − → U P = (1, 0). Here and hereafter, we write the quantities in the co-moving frames in a curly upper-case letter. We continue to write the quantities in the reference frame and in each rest frame (at a particular moment) in lower and upper cases, respectively. We consider the foliation by the player's PLCs and formulate the time evolution as always pulled back to the player's central frame at each moment.
At τ P , the co-moving-frame coordinates − → X are identified with the C(τ P )-frame coordinates − → X : where − → X is given in Eq. (185). For a given acceleration − → A P = (0, A P ) in the C(τ P ) frame, the player's next location and the velocity after an infinitesimal time ∆τ P , felt by the player, are 42 Figure 8: Schematic plot to show the difference between the situations where the nth particle is coming closer to the player, X n ·Û n < 0 in the player's central frame, and going away from the player, X n ·Û n > 0.
At τ P , every object is located on PLC( − → x P ). In the C(τ P ) frame, we may parametrize the location of the nth particle as due to the PLC condition − → X n 2 = 0, derived from − → X n − − → X P 2 = 0 with − → X P = 0. 43 After the infinitesimal time evolution, the nth particle moves onto PLC − → x P . In the C(τ P ) frame, where we see from Eq. (184) that 44 in which we have used − → X P = − → 0 , − → U P = (1, 0), and − → X n = (− |X n | , X n ) = |X n | −1,X n . We see that a particle coming toowards the player,X n · U n < 0, spends a longer proper time compared to a particle going away from the player,X n · U n > 0; see Fig. 8. 41 In an actual implementation with finite ∆τP, Eq. (183) is more usable than Eq. (184). 42 The reference-frame acceleration reads aP = AP + (γ(uP) − 1) (ûP · AP)ûP, from which we may obtain the reference-frame velocity in the next step: u P = uP + aP∆τP. 43 One might think that the PLC condition − → Xn 2 = 0 leads to − → Xn, − → Un ? = 0 by differentiating it with respect to τn. This is not the case because the nth-particle's trajectory is not on the PLC(xP) but away from it. That is, − → Xn 2 = 0 does not hold along the path of the nth particle. 44 As Eq. (184) is manifestly Lorentz invariant, it is applicable in any frame.
The Lorentz transformation from the C(τ P ) frame to the C(τ P ) frame is and the co-moving coordinates − → X at τ P are obtained from − → X in the C(τ P ) frame by the following Poincaré transformation: Accordingly, any velocity − → U in the C(τ P ) frame is transformed as Putting − → X n and − → U n into − → X and − → U in the above expression, respectively, we get where ∆τ n is given in Eq. (192). To summarize, the basic evolution equation in the co-moving frames is where dX n dτ P τ P := lim with X n | τ P and X n | τ P being given in Eqs. (196) and (189), respectively, and similarly for U n ; we have used U n | τ P = U n and A n | τ P = A n ; and we have abbreviated the symbol | τ P in the left-and right-hand sides of Eqs. (198) and (199).
The second terms in the right-hand sides of Eqs. (198) and (199) are the extra contribution from the Lorentz transformation that pulls back to the player's central frame in every time step. When the player accelerates towards X n , the nth-particle's apparent velocity (198) receives the extra contribution |X n | A P that kicks the nth particle away from the player. If you accelerate towards an object, its apparent position becomes farther away from you! On the other hand, the Lorentz transformation affects the apparent velocity (199) by −γ(U n ) A P , so the nth particle's apparent velocity tends to be pulled back to the player. These two effects compete with each other, and should eventually bring the nth particle back to the player if the player continues to accelerate towards it.
Let us see this in the following way. For simplicity, we let the nth particle stay at rest in the reference frame: A n = 0; see footnote 24. We consider a constant acceleration of the player towards the nth particle in the player's rest frame: where a > 0 is a constant. Note thatÂ P =X P . Putting this into Eq. (199), we obtain We set the initial condition such that the player had also been at rest in the reference frame before its acceleration. Then the directionX n remains constant throughout the time evolution until the nth particle collides with the player, and the apparent velocity of the nth particle is proportional to its direction:Û n ∝X n . We parametrize the apparent position and velocity as where U(τ P ) is positive and negative when the nth particle appears to be moving away and coming back to the player, respectively. Then Eq. (202) reads which has the solution under the initial condition U(0) = 0. Putting this into Eq. (198), we obtain dX dτ P = −e aτ P sinh(aτ P ) + aX , which is solved as under the initial condition X (0) = X 0 , where X 0 > 0 is the initial apparent distance. The apparent velocity and acceleration read dX (τ P ) dτ P = e aτ P (1 + aX 0 − e aτ P ) , = a e aτ P (1 + aX 0 − 2e aτ P ) .
We see that the particle indeed appears to be moving away from the player at the beginning: The apparent velocity is positive while τ P < τ c P , with the critical value τ c P being given by where, in the last step, we have recovered the speed of light from natural units. The apparent distance at τ c P is, in nonnatural units, When the apparent distance is small, X 0 c 2 /a, we get τ c P → X 0 /c. Then if one wants to see this effect for one second, the apparent distance must be larger than 3 × 10 8 m, which is roughly the distance between the Earth and Moon 4 × 10 8 m. Even when the apparent distance is of that order, an acceleration of the order of the gravitational acceleration on Earth, a ∼ 10 m/s 2 , elongates X (τ c P ) from X 0 only by a fraction of aX 0 /2c 2 ∼ 10 −8 , i.e., by aX 2 0 /2c 2 ∼ 10 m. This is a tiny effect for the distance and acceleration of everyday life. In computer games, one may consider a very large distance and/or acceleration to see such a relativistic effect.
The final proper time of collision is That is, the proper-time difference between the turnover and collision is The corresponding apparent speed is In a formal expansion of Eqs. (216)-(218), assuming that c = 1 were large, we obtain which coincides with the result for the constant acceleration in ordinary Newtonian mechanics.
On the other hand, the correct limit for c aτ P reads Again we see that the relativistic correction is significant for X 0 c 2 /a, as well as for X 0 cτ P .

D Rocket propulsion
As a concrete example of the relativistic acceleration, we show how the rocket propulsion is formulated. Suppose the following: A rocket has a mass m(τ ), is located at − → x (τ ), and has a velocity − → u (τ ) in a reference frame { − → x } at its proper time τ . Note that a mass is a Lorentz-invariant notion and does not depend on frames. The rocket is in an inertial frame, i.e., under no outer force. Let us first switch to the rest frame of the rocket at τ : Suppose that the rocket's mass, including its fuel, has changed by an amount during an infinitesimal proper time ∆τ . That is, its fuel is consumed by an amount |∆m|. At the proper time τ + ∆τ the rocket's mass becomes m + ∆m, with ∆m < 0, and its velocity − → U (τ + ∆τ ) = − → U (τ ) + − − → ∆U .
This relation applies in any frame. In the rest frame (226), we further obtain ∆U 0 = 0.
Suppose that the fuel |∆m| is converted into a jet of mass ∆m jet > 0 with velocity − → W . The energy and d-momentum conservation reads mU 0 = (m + ∆m) U 0 + ∆U 0 + ∆m jet W 0 + ∆E dis , mU = (m + ∆m) (U + ∆U ) + ∆m jet W , where ∆E dis > 0 is a possible energy dissipation. 46 Putting U 0 = 1, U = 0, ∆U 0 = 0, and W 0 = 1 + W 2 and dropping the quadratic infinitesimal terms, we obtain where we have assumed that C jet := −∆m jet /∆m > 0 and C dis := −∆E dis /∆m > 0 are constants. The final expression for the rocket's D acceleration in its rest frame is where d ln m dτ = 1 m dm dτ . Since Eq. (235) is written Lorentz covariantly, we can easily obtain its form in the reference frame at τ : To summarize, when one switches on the rocket and determines its thrust direction to beŴ in its rest frame, the given acceleration is Eq. (236) in the general reference frame. In terms of the acceleration in Eq. (110), the rocket propulsion corresponds to setting The (noncovariant) exhaust velocity of the jet is (238) When the jet is composed of massless photons, we may take the limit ∆m jet → 0 i.e. C jet → 0, in which the exhaust velocity approaches the speed of light V → 1. In the limit, Eq. (236) reads d − → u (τ ) dτ C jet →0 = (1 − C dis ) d ln m(τ ) dτ L(−u(τ )) 0 W (τ ) .
As an illustration, let us examine a rocket propulsion with the constant directionŴ (τ ) = (0, 0, −1). We set an initial condition u(0) = 0 and m(0) = m ini > 0. For a velocity − → u = √ 1 + u 2 , 0, 0, u , the Lorentz transformation to the rest frame (56) is 47 and the D acceleration is Solving the spatial component, we obtain the relation between the final velocity u fin and mass m fin : We see that the ratio of the initial rocket mass, m ini , to the final mass after using up all the fuel, m fin , governs the final covariant velocity u fin . When we set C dis = 0, the exhaust velocity (238) becomes V = 1 − C 2 jet , and the small velocity limit V 1 of Eq. (243) recovers the well known result u fin → V ln m ini m fin . In the limit m fin m ini , we can increase the final covariant velocity u fin → 1 2 (m ini /m fin ) (1−C dis ) 2 −C 2 jet without limit, and the rocket's noncovariant velocity approaches the speed of light v fin = tanh (1 − C dis ) 2 − C 2 jet ln m ini m fin → 1.

E Handling rotation with quaternions
We review how to implement a spatial rotation in three dimensions using the quaternion, which avoids the gimbal lock. This can be used to store the information on the direction of each 3D object introduced in Sec. 5.3. Let q be a three-vector: We parametrize a quaternion Q by four real numbers q, q 1 , q 2 and q 3 as Q = q; q 1 , q 2 , q 3 = q + q 1 i 1 + q 2 i 2 + q 3 i 3 =: q + q · i, where the unit quaternions obey in which i 2 i := (i i ) 2 = i i i i . That is, where ijk is the totally antisymmetric tensor 48 for (i, j, k) being even permutation of (1, 2, 3), −1 for (i, j, k) being odd permutation of (1, 2, 3), 0 otherwise.
Hereafter, we employ the shorthand notation ij = 3 i=1 3 j=1 etc. From the relation (247), it is straightforward to show that where the vector product is defined by (p × q) i = jk ijk p j q k . In particular, For any Q = q; q 1 , q 2 , q 3 , we define its conjugate Q by Note that We express a three-vector, say, the Cartesian coordinates x = x 1 , x 2 , x 3 by the quaternion 49 We want to rotate this vector by an angle θ around an axis n = n 1 , n 2 , n 3 . The resultant three-vector can be obtained from 48 To be explicit, the nonzero components of the totally antisymmetric tensor are 123 = 231 = 312 = 1 and 132 = 321 = 213 = −1. All the other components are zero. 49 Here and hereafter in this section, X stands for the quaternion (253) rather than the rest frame coordinate.
One may use the quaternion (255) to specify the direction of a rigid body. It is tedious but straightforward to show that 50 N XN = c 2 − s 2 n 2 x + 2sc n × x + 2s 2 (n · x) n · i, where s := sin θ 2 and c := cos θ 2 . When we take n to be a unit vectorn with |n| = 1, we obtain further N XN = [(cos θ) x + (sin θ)n × x + (1 − cos θ) (n · x)n] · i. (257) Suppose that we decompose x into the components that are parallel and perpendicular to the unit vectorn: where x := (n · x)n, x ⊥ := −n × (n × x) .
In other words, the three-vector x is transformed by X → X = N XN as where R ij =n inj (1 − cos θ) + δ ij cos θ − k ijknk sin θ.