{'fraction_non_alphanumeric': 0.04036052668692386, 'fraction_numerical': 0.011449158300960294, 'mean_word_length': 4.704944411936805, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relativistic Langevin Equation derived from a particle-bath Lagrangian

We show how a relativistic Langevin equation can be derived from a Lorentz-covariant version of the Caldeira-Leggett particle-bath Lagrangian. In one of its limits, we identify the obtained equation with the Langevin equation used in contemporary extensions of statistical mechanics to the near-light-speed motion of a Brownian particle in non-relativistic dissipative fluids. The proposed framework provides a more rigorous and first-principles form of the Langevin equation often quoted or postulated as ansatz in previous works. We then refine the aforementioned results by considering more terms in the particle-bath coupling, which improves the precision of the approximation for fully relativistic settings where not only the tagged particle but also the thermal bath motion is relativistic. We discuss the implications of the apparent breaking of space-time translation and parity invariance, showing that these effects are not necessarily in contradiction with the assumptions of statistical mechanics. The intrinsically non-Markovian character of the fully relativistic generalized Langevin equation derived here, and of the associated fluctuation-dissipation theorem, is also discussed.


I. INTRODUCTION
The Langevin equation is the cornerstone of modern nonequilibrium statistical mechanics [1,2], with profound implications for numerical computations [3,4]. It is the basis upon which one builds the diffusion equations [5][6][7][8], fluctuation-dissipation theorems [1,9], Brownian motion with state-dependent diffusion coefficient and its pathintegral formulation [10], and instrumental in areas overlapping with other fundamental theories of physics, such as the question of quantum decoherence [11,12]. Though it may appear to be a simple reformulation of Newton's second law that accounts for thermal noise and frictional damping, it encapsulates the modern view of the structure of matter and the dynamics of Brownian motion. For example, the Langevin equation summarizes our fundamental understanding of molecular motion in liquids as systems made up of stochastically moving molecules, on which short-range collisions and interactions with other surrounding molecules impart both a stochastic excitation and a viscous damping [13]. Extending the classical framework of Brownian motion, as embodied by the Langevin equation, to relativistic systems is vital to many areas of physics, which include relativistic fluids/plasmas [14][15][16], effective field theories of dissipative hydrodynamics [17], relativistic viscous electron flow in graphene [18,19] and other solids [20,21], quarkgluon plasma (QGP) [22] and nuclear matter [23], where thus far only the classical non-relativistic Langevin-type description is customarily used to describe fission and fusion-fission processes of hot nuclei [24][25][26][27][28]. Nonetheless, in all these cases it is clear that relativistic corrections would be vital. Furthermore, a Lorentz-covariant extension of the Langevin equation framework would be the first major step towards building a consistent relativistic version of nonequilibrium statistical mechanics [29] as well as of fluctuations theorems and thermodynamics in general [30,31].
One expects that the Langevin model generalizes to the limits of both fast-moving fluids (where the relativistic motion of constituent particles becomes important), and for relativistic Brownian particles. Several relativistic or quasi-relativistic extensions of the Langevin equation have been proposed, adopted and used [32][33][34][35]. However, rarely if ever any justification beyond that of an ansatz is provided, and no derivation from first principles is given for the case of relativistic environments. While there are many general approaches [29,36] to extending the Langevin equation to the weakly-or moderatelyrelativistic regime, several are either postulates of a dissipation model [34,36,37], or generalizations of causal fluid dynamics [38] to a case where they have limited applicability. Some works provide a justification for the case of the relativistic Brownian particle in a non-relativistic environment [33][34][35]. We, however, present bottom-up derivations of the relativistic extension of the Langevin equation that is valid for relativistic environments as well as relativistic particles, and which recovers the appropriate results when taking the respective limits.
Our approach is to extend a derivation of the Langevin equation from the Caldeira-Leggett model [39,40]. Section II is largely expository; we feel it is important to recall the particle-bath non-relativistic treatment in detail as the starting point of our derivation. It is useful arXiv:2107.07205v4 [cond-mat.stat-mech] 11 Jan 2022 in order to appreciate both the many parallels and the fundamental differences that arise with the relativistic treatment. We then (Sec. III) alter the (Galilean) Lagrangian to be explicitly Lorentz-covariant.
Due to the many difficulties associated with non-linear equations that result from the aforementioned Lorentzcovariant Lagrangian, the largest portion of the paper is dedicated to obtaining a compact approximation to the equations of motion of both the heat bath and the tagged particle. This procedure is described in Section IV, and a particularly detailed breakdown of the necessary steps is given in Section IV A. Broadly, we obtain a numerical and then a quasi-analytical closed-form approximation to the solution of the (Euler-Lagrange) equations of motion of the heat bath. We then use this solution to construct a coarser, but tractable, approximation that is both obviously a relativistic dual of the respective intermediate result of the Galilean derivation, and accounts for the most important relativistic corrections. By integrating the heat bath equations of motion, we obtain the final Lorentz-covariant dual of the Langevin equation, which is the central result of the present work.
Before we begin, we note that while some authors (cfr. [29,34,38]) use covariant notation, in our case, said notation would obscure the intent, and make the parallels to the non-relativistic derivation harder to see. For the same reason, we alter the straightforward derivation of the Galilean Langevin equation to be in the Lagrangian rather than Hamiltonian formulation.
We shall work with the coordinate time, of one frame referred to as the "lab frame" in the context of relativistic Brownian motion [29]. Inasmuch as the classical (generalised) Langevin equation derived from the particle-bath Caldeira-Leggett model is able to correctly describe Brownian motion in classical statistical mechanics, we shall assume that also our results may apply to general relativistic Brownian processes, at the level of special relativity. The final important distinction from previous attempts at extending classical Brownian motion to relativistic speeds, is that we aim to describe the behavior of a single Brownian particle (the Langevin picture), rather than of probability distributions (the diffusion or Fokker-Planck picture). Probabilities are rife with counter-intuitive phenomena that could be mistaken for causal propagation at superluminal speeds [41], thus we avoid these problems by excluding probabilistic considerations and working within the Langevin framework.
The main result of our paper is a fully-relativistic Langevin equation which describes fully-relativistic motion of both the tagged particle as well as of the environment (the bath oscillators). To our acknowledge this is the first time that such description is achieved, being all previous approaches either valid only for weaklyrelativistic systems or for a relativistic particle in a nonrelativistic environment.

II. GALILEAN LANGEVIN EQUATION
In this section we recall the derivation of the Langevin equation from a particle-bath Lagrangian under Galilean relativity. It is important to recall the derivation with some details as it will serve as the starting point for the new Lorentz-covariant derivation that we will present in the next Sections. Also, it is important to recall under which approximations the Langevin equation can be derived from a particle-bath model, which we shall elaborate on in this section.
Zwanzig [39], and later Caldeira and Leggett [40,42], considered a single (tagged) Brownian particle in linear dynamical coupling with a heat bath of harmonic oscillators. This choice is justified by three observations: (i) the coupling is an analytical function and usually sufficiently weak to be closely approximated by the linear term in its expansion, (ii) the stochastic motion of the environment is represent-able as a Fourier transform, and (iii) the behavior of the environment is only affected by the motion of the particle and no other external influence, while (iv) the behavior of the particle can often be described by a conservative force, and all dissipation is due to the interactions with the heat bath (i.e. no vector potentials). The weak coupling stems from the short-range nature of the interaction. Averaging the impulses to the time scales of the Brownian motion effectively reduces the force coming from the environment. This averaging also sufficiently smooths the function to be amenable to Fourier analysis.
It is important to note that while our method of using the Caldeira-Leggett Lagrangian is remarkably compact, it is neither the only, nor the most general derivation of the Langevin equation. The latter title belongs to the Zwanzig-Mori formalism [1,43], wherein only determinism of the underlying (microscopic and not necessarily classical) equations of motion is required. In Zwanzig-Mori formalism projection operators are used to decouple the motion of the entire coupled system into a relevant and irrelevant part, and obtain equations of motion for the relevant part. For the Caldeira-Leggett model, this step is not necessary, resulting in a compact derivation, which proves to be sufficient for our purposes.
These assumptions lead to the Caldeira-Leggett Lagrangian [42]: where the q i are the generalised coordinates of the reservoir modes, ω i -the natural frequency of the i-th mode and g i the specific coupling between the particle and the i-th mode. Finally x is the position of the tagged particle. Vectors are indicated with bold, and time derivatives with over-dots. Also recall, that we have chosen the Lagrangian over the Hamiltonian formulation [1], in order to take advantage of the Lagrangian formulation's explicit Lorentz invariance in later sections.
Ordinarily, the potential energy for the system of harmonic oscillators depends on m i ω 2 i q 2 i , plus we've stated that there is a coupling between the modes and the tagged particle, but one does not see why the tagged particle necessarily contributes to the potential energy. We impose this contribution, called the counter-terms [44] so that, the Lagrangian has manifest space and time translation invariance, as well as both space and time inversion symmetry. It is the reason why one can assume the leading order elements of a generalised coupling function to be linear in the difference in the generalised coordinates and x. This ad-hoc induction of symmetry is explored more in Sec. VIII D.
The Euler-Lagrange equations of motion for the particle and for the bath oscillators, respectively, are given by: where, without loss of generality, m i = 1, ∀i ∈ N . As we shall see later, there is an equivalent choice, that is more appropriate to the relativistic derivation, but for the sake of remaining comparable to the standard derivation, we shall continue with re-scaling the masses to unity. We should also highlight, that we treat the dependence of q on x as a parametric dependence, rather than functional. By contrast, in the following section, we shall make this dependence a functional one.
We assume that the trajectory x(t) is known, and solve for the motion of the modes q i (t; x), as a parametric function depending on x(t). Using the Laplace transformation of Eq. (2b), one obtains: which when integrated by parts yields The term, is known as the stochastic force [1]. Substituting Eq. (4) into Eq. (2a) we obtain the (Galilean) generalised Langevin equation [1,39] where is known as the friction kernel or the memory function [1]. It represents the ability of the medium to respond to temporally distributed perturbations, such as those encoded inẋ. The equation (6) is called generalised, because we have not yet imposed any restrictions on the behavior of the contributing terms: the stochastic force and the friction kernel. The following is an example of obtaining a special case known as the Markovian Langevin equation from Eq. (6). We start by noticing that Eq. (7) is a Fourier series decomposition of K(t) in terms of ω i , as is Eq. (5) a decomposition of F p . This series decomposition, can be replaced with a Fourier transform, for a continuous vibration spectrum, and a given density of states ρ(ω). The specific choice of ρ(ω) dictates the behavior of K(t) and F p , via: For the Markovian equation, ρ(ω) ∝ ω 2 , as per Debye's law and/or the Rayleigh-Jeans law for bosonic baths; and g(ω) = Const. Thus K(t) ∝ δ(t) as shown in [1]. The friction term in Eq. (6) becomes the classical viscous force ηẋ and Eq. (6) reduces to the special case of the Markovian Langevin equation [1].
This derivation not only obtains the Langevin equation Eq. (6) from general assumptions: weak coupling, spacetranslation invariance; but also provides a clear physical model for the origins of physical effects such as viscosity and the stochastic excitation. The most natural model for the environment -the harmonic heat bath, is also the same model used for describing quantum decoherence in the theory of open quantum systems [12, p. 656].
Importantly, this derivation demonstrates the classical deterministic worldview of statistical mechanics. Even in its final form the Galilean Langevin equation would be fully deterministic provided knowledge about initial conditions for all the degrees of freedom is at hand [1]. The memory function and the stochastic force can be fully predicted both forwards and backwards in time, provided the initial conditions at an arbitrary t = 0 are known. The reason why the memory function and the stochastic force appear to break time-inversion symmetry, is because (i) The number of Fourier modes (i.e. bath oscillators) is very large (in fact, merely recording O(10 9 ) initial conditions already requires at least 16Gib of RAM) and their initial conditions cannot be specified deterministically, hence they have to drawn from statistical distributions such as Boltzmann or Jüttner for the non-relativistic and relativistic limits respectively. (ii) The system is possibly chaotic due to the intrinsic non-linearity of the relativistic oscillators: minor perturbations to the initial conditions can greatly impact long-term behavior. As a result, while Eq. (6) is technically deterministic, it is effectively stochastic [1,39].
Finally, one must note the fluctuation-dissipation relations [1, see p. 23]: where . . . denotes ensemble or time averaging (which are equivalent for ergodic fluids). The Lagrangian (1) is manifestly Galilean covariant and not Lorentz covariant. What this means is that the Eq. (6) is not accurate in cases where either the environment and/or the tagged particle, are highly relativistic: e.g. stellar cores, the early universe, as well as the plasma in collider experiments.

III. RELATIVISTIC CALDEIRA-LEGGETT LAGRANGIAN
In this section we update the Caldeira-Leggett model and the derivation of the generalised Langevin equation to work in special relativity. While there are relatively many descriptions of a relativistic Brownian particle coupled to a non-relativistic heat bath [29,36], there are far fewer attempts at tackling the relativistic behaviour of the heat bath. Recall that in the Galilean case, the main difficulty lies in obtaining the equations of motion for the heat bath modes. This task is non-trivial, as there are many relativistic Lagrangians that reduce to Eq. (1) in the relevant limit.
Moreover, in the previous consideration we have encapsulated all external forces in V ext (x) = φ(x), which has causality implications under special relativity. By analogy with electromagnetism, we instead consider adding an extra degree of freedom in A and re-absorbing the "charge" equivalent into the definitions of A and φ.
Since electromagnetism is Lorentz covariant, we get familiar equations of motion and sidestep problems with action-at-a-distance.
In the original Galilean derivation by Zwanzig [39], the choice of harmonic equations of motion is a reflection of the efficacy of Fourier methods. It might be tempting to carry over the harmonic trajectory, rather than modifying the interaction. However this is mathematically equivalent to non-relativistic heat baths already considered in [29,36]. Thus, instead, the heat bath dynamics are relativistic equivalents of harmonic oscillations, rather than harmonic oscillations themselves. There are many models which reduce to simple harmonic oscillators, many of which are not solvable [45][46][47]. We chose to parameterise the relativistic equivalent of the harmonic potential, for the following reasons: (i) it is obviously a norm of a 4-vector, (ii) it obviously corresponds to the conventional coupled harmonic potential in the limit c → ∞ (iii) the constituent "events" can straightforwardly be interpreted as an interaction mediated by a (massless) particle, (iv) requiring that the "events" be null-separated is the relativistic (local) equivalent of Galilean (global) simultaneity, (v) when the interaction "events" are null-separated, many of the mathematical complications vanish. Bringing all of these considerations together we arrive at where and overdots signal differentiation with respect to corrdinate time. We must note the counter-intuitive presence of γ −1 rather than γ in the kinetic terms. This term in the Lagrangian has been discussed at length in standard literature [48, p. 323]. It pertains to the difference This leads to a renormalised frequency in Eq. (19).
between proper time and coordinate time, and the fact that the Euler-Lagrange equations are formulated using proper, rather than coordinate time (see Goldstein and Poole [48, p. 320]). The coupling term as discussed previously is formulated in terms of the Lorentz-invariant 4-interval between discrete interaction events (Eq. (13)), following the same logic as the introduction of counterterms in the original non-relativistic Caldeira-Leggett Lagrangian [40,44].
The equations of motion associated to the Lagrangian (13) are where the functions are evaluated at coordinate time t and F ext = ∇φ(t, x) + ∇( A·ẋ c ) −Ȧ c . In the equations of motion, this is the only term due to differentiation of the potential with respect to the first derivative that survives. The terms resulting from differentiating γ(ẋ) contribute a full interval, which is null in our model. The factors of γ −1 signal that the Lorentz contraction cannot and does not contribute to the force.

A. Problem setting
We proceed in complete analogy with the Galilean case, assuming that x(t) is a known trajectory, and reverse-engineer the equation satisfied by the heat bath. By expanding the left hand side of Eq. (15b) and using the product rule, we obtain To proceed, one projectsq along the directionq, decoupling the trajectories into longitudinal and transverse components.
By combining similar projections we obtain for the longitudinal, and for the transverse component. Here, notice that we have used q(t, x), while in reality q is only a function of time, and the x dependence comes indirectly from assuming knowledge of the trajectory x(t). Viewing this implicit dependence as an explicit depdendence on a free variable x is what allows us to solve Eqs. (17a) and (17b). Our focus for this section shall be developing a method for solving the two equations. At present, mathematics does not support an exact closed-form solution in terms of standard functions [46,47,49] to either of the known relativistic harmonic oscillators. Since Eq. (17b) and Eq. (17a) break the isotropy of the problem, one would expect that the components of the solution need to be tracked separately. However, looking at the McLaurin series expansion of Eq. (17a) in the limitq → 0 and comparing it to the expansion of Eq. (17b), we can see that the only difference is the numerical factor of powers ofq/c. This suggests that a functional form approximating the solution to Eq. (17b) can also approximate Eq. (17a), albeit with different numerical parameters. We shall be working under this assumption, and justify it in Appendix VI.
Our plan for solving the relativistic heat bath's equations of motion is as follows. Firstly, we shall produce a numerical solution to a simplified 1 + 1 dimensional problem, using the longitudinal component -Eq. (17a), as a base (Sec. IV B). We shall then use this numerical solution to construct a more compact closed-form solution that retains cardinal properties of the numerical solution, while remaining sufficiently similar to Eq. (4) (Sec. IV C). As it happens, Eq. (17a) can be solved in terms of elliptical functions, and while the solution is too cumbersome to even quote (much less manipulate), we shall use some of its properties to validate our previous step (Sections IV C and IV D). In Sec. IV D we find (and in Appendix VI verify), that the closed form approximation requires parameter fitting in order to be a good approximation to the solution of Eq. (17a). We carry out said fitting in Sec. IV E, and explain the necessary modifications to apply the same process to Eq. (17b) in Sec. IV F.

B. The numerical solution
By comparing the numerical solutions of Eq. (17a) and Eq. (3) in Fig. 1, we can see that relativistic corrections manifest as (i) a shift in eigenfrequency of plane wavelets (see Fig. 1) of the bath oscillators. (ii) a slight dependence of the wavelet phase of the bath oscillators on the position of the tagged particle x (see Figs. 2 to 4). (iii) smaller sub-oscillations along contours of constant phase, (Fig. 3). This effect is weak, and ignoring it only slightly reduces the χ 2 of the numerical fit, but greatly reduces the complexity of the following manipulations. So in order to build a Lorentz-covariant analogue of Eq. (3), we will need to incorporate these features.
We do this by introducing an adjusted or renormalised frequencyω, and a relativistic-correction phase ξ(x(t)) into the sinusoidal forms of the bath dynamics. These two corrections play a crucial role in building the solution to the dynamical problem, and they have a clear physical origin as summarised below.
The renormalization of the oscillator frequencies ω i is entirely expected since the dynamics of the relativistic oscillator is nonlinear (see Eq. (15b) versus Eq. (3)). As we know from classical mechanics, the eigenfrequency of an anharmonic oscillator is given by the eigenfrequeny of the corresponding harmonic oscillator plus a correction which depends on the anharmonic coefficients (as well as the square of the amplitude), see e.g. [50, pp.87-88].
Clearly, in our case the harmonic eigenfrequency is the one of the non-relativistic oscillator (γ(q) = 1), i.e. ω i , whereas the renormalised eigenfrequencyω i = ω i + δω i (γ(q)) is given by ω i plus a correction δω i that depends on the degree of anharmonicity (nonlinearity) of the system [50], i.e. on γ(q) in our case. We indeed found from numerical simulations for the initial conditions that we investigated thatω ∼ γ −3/4 (q).
The relativistic-correction phase ξ, instead, is required in order to bring the oscillatory/sinusoidal terms in the solution into a plane wave form. In turn, the so-obtained plane wave form is Lorentz-covariant thanks to Lorentzinvariance of the wave equation under the Lorentz transformations [51, p. 383]. We shall provide a detailed clarification and analysis of these terms in the next section.
We want our solution to be a manifestly Lorentzcovariant version of Eq. tioned corrections, we obtain: whereω i is the renormalised frequency that we introduced above, [52]. This analytical form provides a reasonably accurate, Lorentz-covariant match of the numerical solutions to Eq. (17a), thus eliminating the mismatch shown in Fig. 1. Firstly, unlike a naïve extension of the Galilean Eq. (3), in Eq. (19), we have locality and causality encapsulated in the phaseξ(t, x(t)), which is an integration constant. To understand the origin and physical meaning of this parameter, one must delve deeper into the solution's structure.
Secondly, we have two pairs of functions:ξ(x) and ξ(t, x(t)), which we shorten toξ(t); andω andω. The first pair represents the phase for two distinct situations, one is dependent only on x, while the other also explicitly depends on t. Similarly, the second pair are the frequencies for the same situations. The two formulations are equivalent as shown in Appendix C. We've added the overbar toξ for consistency.

C. Computationally-informed structure of the solution
In this section, we provide a justification of the aforementioned correction coefficients, i.e. the renormalised frequencyω and phaseξ, in reference to a quasianalytical solution of Eq. (17a) obtained in terms of elliptical functions in 1 + 1 dimensions. For compactness the indices i are dropped and reinserted only in Eq. (26). The computationally-obtained quasi-analytical solution is too long to be reported here (it can be found in the accompanying Mathematica notebooks [53]). Instead we use it to study the generic features of the solution in order to understand and justify the compact, analytical approximation, as well as explain the origins and methods of obtaining the renormalised frequencyω and a time dependent phase ξ from initial conditions.
We consider x = x as an independent variable of the scalar function q(t, x), and attempt to solve Eq. (17a), as if it were an ordinary differential equation in t, allowing x to vary. This allows for one extra degree of freedom, which permits approximating the solution in terms of elliptical functions. This solution, though cumbersome, contains two integration constants, whose roles are closely related to the roles ofω andξ in Eq. (19).
Firstly, the generic form of the Mathematica-generated quasi-analytical solution reads as where F is an expression whose exact structure is not important [53], save for being quasi-periodic. It is only approximately periodic in t i.e. no longer independent of x, as in Eq. (3), (see Fig. 4), due to the presence of C 2 (x) in the solution (and C 1 (x) to a lesser extent). By comparing Eq. (19) and Eq. (20) and matching the arguments at the local maxima of both functions,ξ(x)/c = C 2 (x) must be the integration constant responsible for the curvature of the wavefronts in the solution as x increases (see Fig. 3).
To understand the relationship between the other integration constants, one must dig more deeply into the structure of F . The form of F is the inverse of a combination of elliptical functions of a compound argument, that involves C 1 (x), t − C 2 (x), ω and g. Neither the expressions that contain C 1 (x) nor ω can be factored out. Thus when formulating Eq. (19), we reflect this by introducinĝ ω =ω(gx, ω, t, . . .), which encapsulates the fact that the solution is anharmonic, and that the frequency is shifted by a quantity that depends on C 2 (x). Thus both the phase and the frequency are position-dependent. This introduces an unnecessary complication which we can avoid by defining our integration constants slightly differently by exploiting a degree of freedom (see Appendix IV D). For our purposes it is easier to work with a time independentω with all of its time (and position) dependence added toξ, rather than withω andξ.
Bothξ andω play the role of integration constants. In principle, if one had a set of initial conditions to fix both C 1 (x) and C 2 (x), one could solve the equation where q(t, x) is defined in Eq. (19), to obtain bothξ and ω in closed form as functions of C 1 and C 2 . Notice, that there are two degrees of freedom corresponding to C 1 (x) and C 2 (x) that transform into two degrees of freedom in ξ andω, that are linked with only one equation. This is what allows us the freedom to choose to work withξ and ω, rather thanξ andω, which are more directly related to integration constants of Eq. (20). We also remind the reader that x enters expressions that straightforwardly generalize to (1+3) dimensional dynamics, while expressions containing x are specific to (1+1) dimensional considerations. Case in point, the concept of wavefronts is a convenient representation of the Fourier decomposition of the motions of q and x that only applies to the case of scalar x and q.

D. Integration constants and initial conditions
At this point it might appear that the phase-type term ξ(x) can be chosen arbitrarily, but that is not the case.
We are bound by causality, which links C 1 and C 2 , and gives us a clear understanding of the functional form of q i (x, t). By looking briefly at the functional form of F (λ) returned by Mathematica [53], one shall notice a deluge of terms, many of which fall into either or where we used the shorthand x = |x|. While it needs to be shown conclusively and rigorously, it is rather evident that causality constrains all such terms to be real. The solution also contains terms that mix the argument, (thus C 2 (x)) with C 1 (x) under a square root, so we can conclude that the integration constants are not independent of each other. An imprint of causality should be present on C 2 (x) as well as C 1 (x).
Finally, we must reconcile the integration constants, C 1 and C 2 , with the constant boundary conditions that we have assumed for the Galilean case. The constant boundary conditions could violate causality, thus proving problematic. So instead of assuming that the boundary conditions are constant across the entire domain of the solution, we instead assume that the boundary conditions are constant in the domain of events which are separated from the co-ordinate origin in a time-like fashion, i.e. within the light-cone of the origin, and zero otherwise.
As it turns out (see Appendix B), the relevant separation scales at which this could occur, are far in excess of the distances at which other assumptions we've made would break down. So we are justified in constraining the domain of x ∈ (0, 1), (g = 2 for convenience) and in choosing arbitrary initial conditions. For example, within the light-cone of the current event q(0, 0). Here we chose both numbers and constant initial conditions out of convenience, and for illustration purposes. q(0, x) = 1 is an arbitrary choice, while in the units of c = 1, q(0, x) = 0.85 is both sufficiently large to show relativistic effects, but also the largest number that doesn't result in numerical errors. Now suppose that we have imposed the relevant initial conditions and that the conditions are valid for our particle-bath system. We can thus obtain C 1 (x) and C 2 (x) from the closed-form approximation to the numerical solution detailed above. Moreover, we can determine what theξ(x, t) is for the relevantω = Const. via algebraic manipulations of elliptic and trigonometric functions, at least in principle. However, it is more relevant to attempt a different approach, that we have used to generate the relevant plots illustrating the match between the numerical solution and Eq. (19).  (19) is a generalised sinusoidal function, with the caveats forξ andω discussed above. The true solution to Eq. (17a) (and to Eq. (17b)), however, is not a perfect sinusoid, but can be approximated arbitrarily well by a suitable choice ofξ andω, which one can show by substituting Eq. (19) into Eq. (17a) (and Eq. (17b)). Of course, we are interested in the opposite, findingξ andω that minimize the mismatch.
Our first intuition is that the solution is fully periodic, so the frequency at x = 0 is independent of time, which is certainly the case for the regular sinusoid. Indeed we may obtainω under these assumptions, which we call frequency matching. The result is shown in Fig. 2.
However, as we can see from Fig. 2, the frequencymatched solution (blue) and the numerical solution (orange), will slowly drift apart, even though in the neighbourhood of x = 0 and t = 0 the agreement is nearly perfect. This "drift" is due to one of many differences between the relativistic and Galilean simple harmonic oscillators. The presence of γ(q) means that q andq are not always exactly π/2 out of phase, which can be viewed as a change of effective mass of the mode leading to a different frequency. We shall remove this drift from ω, and reabsorb it intoξ, as discussed in the previous section (see also Appendix D). Frequency-matching alone is therefore not sufficient to produce an accurate analytical description.
Instead, we must have two independent hypotheses for the forms ofξ andω and fit them simultaneously. Whilē ω can itself be treated as a scalar, we need two more parameters forξ.
The shape in Fig. 4 can be approximated by an offset parabola:ξ where A and B are the aforementioned parameters. Note that we assume thatξ(t = 0, x = 0) = 0, which we need for consistency with the initial conditions, but that ξ(t = 0, x = 0) = 0, and scales proportionally to time, as we would expect in the first approximation to the drift we saw in Fig. 2. By a least-square fit of the numerical solution using Eq. (19), for the case of g = 1 we obtain that the best fit parameters are B = 1.0000 ± 0.0003, A = 0.0045 ± 0.0005 andω = 0.78634 ± 0.0008.
Note thatω is not of the order of magnitude of γ for the relevant velocity. This shows that the change in frequency is not solely attributable to Lorentz time dilation.
Consequently, the results of the matching are presented in Fig. 5. Boundary conditions will require a different hypothesis, but a similar approach, which leaves unchanged the general form of the solution. Hence, with this general prescription we can obtain an analytical form which is consistent with the numerical solution, as demonstrated in Fig. 5.

F. Transverse component of the trajectory
It is worth noting that the procedure can be similarly repeated for Eq. (17b), which gives the transverse component of the acceleration of the heat bath mode. The differential equation that defines the equations of motion is different, by a factor of γ 2 (q), which one may think should significantly impact the applicability of Eq. (19), as an approximation. Our calculations show that the ansatz Eq. (19) is indeed a good approximation for this equation as well, with differentω andξ. For the case of constant boundary conditions, that we have discussed previously, the parabolic hypothesis forξ with linear time scaling shows similar agreement, but with different values of the parameters A and B.

V. TAGGED PARTICLE DYNAMICS
Having justified the general structure of the solution for the bath dynamics, based on the numerical results, we now proceed to integrating by parts to obtain the dynamical trajectory of the tagged particle using Eq. (19) calibrated as in the above section. In general, the sinusoidal behavior in the integrand is not independent of the trajectory x(s) (recall thatξ(s) =ξ(x(s), s), where the implicit dependence is the result of the parametrisation, and the explicit dependence is the effect of absorbing the x dependence fromω), because of the dependence encoded inξ(x(t)). From a mathematical point of view, however, we can treatξ(x(s)) as ultimately a generic function of s, and writeξ(x(s)) ≡ξ(s), which is an equivalent notation. This leads to the following integral for the last term on the r.h.s. of (19): where, in the right hand side, we use the Newton-Leibniz formula for a definite integral, expressed as two antiderivatives. We have slightly abused the notation here, in that the anti-derivative defines a set of functions which differ up to a constant term, which one can in principle determine from the boundary conditions. In the Galilean case, the result would have been the same regardless of the boundary conditions, but in the relativistic problem, we must determine the integration constant, and only then substitute s with the proper value. In this first step the necessity to find the integration constant may not be immediately evident, but will become apparent as we separate the terms.
We perform integration by parts in Eq. (25), and regroup the resultant terms. Upon replacing in Eq. (19), and following the steps reported in Appendix D, we finally get This is the sought-after result, and the relativistic analogue of Eq. (4).
At this point it may be useful to reiterate the intro-duced terms. Regarding the presence ofξ i , for each heat bath mode i, this is connected to establishing initial conditions in a statistical manner [1]. In fact, we (b) Weak coupling of g = 0.2. Note the near absence of x dependence of the phase. This is because the relevant quantity is gx rather than x, as one would expect from Eq. (17a).

FIG. 4:
Contour plots of the numerical solution q(x, t), where g is given in the sub-caption, and the remaining parameters (c, m, ω) were set to 1 in their respective units. Color gradients represent the value of q(x, t) from small (blue) to large (orange). These plots illustrate the x and t dependence introduced into the phaseξ of the relativistic wave solution (Eq. (19)).
will later assume that the initial conditions are drawn from a Maxwell-Boltzmann-type distribution as in the Galilean case [1,39], which mutatis mutandis is known as the Jüttner distribution in relativistic statistical mechanics [54,55]. Upon moving to a continuum of eigenfrequencies, the functionsξ i are replaced by an extra parametric dependence onω. To understand why, first recall how one moves to the continuum in the Galilean case. The mapping between an index and the mode's frequency is bijective, by which we meanω i =ω j if, and only if i = j. So the index is interchangeable with the frequency. Instead of stating thatξ i is the phase function belonging to the i-th mode, we state that it is a functionξ(t, x(t);ω i ) that depends on the frequency of the i-th mode. The summations with respect to i are replaced with integrations with respect to the continuum-spanning frequency variableω. The same is true of all the parameters which vary mode-to-mode, picking up an additionalω dependence.
Of course, here again we use the shorthand ξ(ω, t, x(t)) =ξ(ω, t), to make the notation more compact. The dependence on x is implicit and should be understood from the context.

VI. BAYESIAN JUSTIFICATION OF THE PROPOSED ANALYTICAL APPROXIMATION
Formal verification of the results that form the basis of Section IV B is needed. One could have in principle merely checked that substitutting Eq. (26) into Eqs. (17a) and (17b) for the appropriate parameters over a certain domain of the free variables, e. g. x, produces identity. Unfortunately, this forms a circular dependency, because the values ofξ andω are established assuming that the equations are correct to within machine precision.
Bayesian inference is widely used in statistical cosmology, particle physics and gradually also in other branches of physics. It is widely used for model comparison, particularly in cases where more fine-grained analysis of multiple theories is necessary. In addition to being able to reproduce any and all results obtained with frequentist statistical methods, Bayesian inference is able to automatically impose Occam's principle [56] in model comparison, thus avoiding "overfitting" (the famous "Fermi elephant").
We have based on Eq. following function: where q i is the approximate solution from Eq. (26), and ξ is defined in Eq. (24). Thus we have a hypothesis with three degrees of freedom: A and B from Eq. (24) andω, and two nuisance parameters: x and t.
For the likelihood we have chosen and the priors on the parameters were A ∈ (0, 100), B ∈ (0, 2) andω ∈ (0.001, 100), all uniform. Inference with PolyChord [57] yielded an uncorrelated Gaussian posterior with B = 1.0000 ± 0.0003, A = 0.0045 ± 0.0005 andω = 0.78634 ± 0.0008. and a Bayesian log-evidence of ln Z = −18.5742 ± 0.0005. This result is to be interpreted as follows: firstly, that our least-squares fit had been accurate. Secondly, the evidence is within 12σ of the result we would expect if the parameters were truly normal-distributed -ln Z expected = −18.5684. It does not mean that our analysis is invalid, quite the opposite, it suggests that the parameters' joint distribution is not exactly an uncorrelated multivariate normal distribution, but that it is a very good approximation: the difference is only 0.03% of the evidence. One interprets the evidence by comparing it to a better, but much more complex quasi-analytical solution to Eq. (17b) that can be obtained using Mathematica and that we do not report here as it is an extremely long expression (in the following referred to as "quasianalytical" solution). We use the same likelihood, but redefine q i in Eq. (27) to now be the quasi-analytical solution provided by Mathematica [53]. In this case the likelihood is independent of A, B andω. As a result of being independent of two parameters, the log-evidence jumps to ln Z = −10.11 ± 0.05, which is both significantly larger than the one we obtained earlier (due to the Occam penalty associated with A, B andω), but also indicative of the fact that the quasi-analytical solution is not perfect. If ∆ = 0 for all values of t, x, A, B and ω then the log-evidence is the log-volume of the prior space: ln Z = −9.9, which is the evidence one would have obtained with the "ground truth" solution. This discrepancy of 3% in the log-evidence tells us both that the quasi-analytical solution is a reasonable approximation to the solution to Eq. (17b). At the same time, however, it is not an exact analytical solution, since (10.11 − 9.90) = 3σ. For comparison, if we fix A, B andω, in our previous consideration, to their best fit values and combine the errors in quadrature, we obtain ln Z = −10.21 ± 0.08.
In quantitative terms, this means that the quasianalytical approximation is, with a e −10.11−9.9 = 81% Bayes ratio, the solution to Eq. (17b), and the approximation in Eq. (26) after parameter fitting -is the solution with a Bayes ratio 74%. Larger Bayes ratios signal a better fit. Paraphrasing Ref. [58], we have found strong evidence for the quasi-analytical approximation that one obtains from Mathematica, and marginally-strong evidence for the much more compact and manageable approximation in Eq. (26), provided that we do parameter fitting for the boundary conditions specified in Eqs. (23a) and (23b). This analysis certainly justifies the tradeoff between precision and compactness of the analytical expression in favour of Eq. (26), which will be used in the following to obtain a fully relativistic form of the Langevin equation.

A. Term classification
To obtain the final relativistic Langevin equation, we substitute Eq. (26) into Eq. (15a). However, we shall first analyze and classify the terms entering into the final equation.
Analysing Eq. (26), we can clearly notice parallels to the original Galilean Langevin equation Eq. (6), as well as fundamental differences. Let us start from the parallels. There is a frictional force term consisting of an integral, with the time derivative of the coordinateẋ, which multiplies the following memory function: This differs from the memory function of the Galilean case because of the renormalised frequencyω and because of the presence of the space-like parameterξ(t) needed to make the trigonometric functions manifestly Lorentzcovariant.
Similarly to the Galilean case in Zwanzig's treatment [1,39], one can identify terms that show dependencies on the boundary conditions, with the stochastic force: Here, exactly like in the original Galilean derivation by Zwanzig [39], the stochasticity of the force follows from "ignorance" about the boundary conditions [1]. Had we not coarse-grained away the information about the full microstate of the system, the stochastic force would appear fully deterministic [1]. It can be shown (in Section VIII F) that the above force term F p given by Eq. (29b) has zero average, F p (t) = 0, as required for the stochastic force in the Langevin equation [1]. We shall discuss this in more detail in subsequent sections. However, there are new terms that are not present in the Galilean equation (ifξ i (t, x) were to be treated as a small parameter, then all such terms would be O(ξ)), such as: Interestingly, this looks like a restoring force where the "spring constant" is dependent on the trajectory through ξ i (x(s)) -contained in the r.h.s. of the above equation -, which affects the overall dynamics of the system. In other words, Eq. (29c) describes a restoring force which retains "memory" of the dynamics, in a similar way as the memory encoded in the stochastic force. To our knowledge, this new, emergent "restoring force" induced by the coupling to the bath in the fully-relativistic regime has never been derived or discussed in previous literature. Thus the final equation takes the form where we use the prime to indicate that the primed terms are different in the relativistic formulation compared to the Galilean case. This equation is the central result of this paper, and represents a full-fledged relativistic Langevin equation. To our knowledge, this is the most general form of relativistic Langevin equation proposed so far, and recovers the previously known forms of Langevin equations (non-relativistic Generalised Langevin equation, Langevin eq. for relativistic tagged particle in non-relativistic bath, and weakly relativistic Langevin eq.) in the relevant limits.
In the next section we further discuss the structure of this equation and the significance of its terms, as well as the various limits that can be recovered. We will also address aspects related to symmetry implications connected with Eq. (30).

B.ξ andω in practice
It must be noted, that so far we have not used any specific value or functional form for eitherω i orξ i . Indeed, the two are "integration constants" in the sense that their value is fully determined by initial conditions. We have shown a functional form that can be obtained for a very specific set of initial conditions in Sec. IV E, and while these specific initial conditions (Eq. (23)) are chosen at random, the algorithm of determiningξ andω should be sufficiently general.
One can impose the initial conditions based on a specific model of interaction. For example if the heat bath comprises of a mass on a spring that enters elastic contact interaction with the tagged particle, at equilibrium, then the initial conditions are straightforward: While this particular situation is highly unlikely to occur in real systems, the principle can be generalised.

A. Stochasticity
We start our discussion of the Lorentz-covariant Generalised Langevin equation by noting that much like its Newtonian counterpart Eq. (6), it too is fundamentally deterministic. While some variants of the generalized Langevin equation share this property, not all versions are indeed deterministic or time-reversible. While what we have found is indeed a rather general statement about the behaviour of many systems, the true (and ultimate) derivation of the generalized Langevin equation must allow for intrinsic non-determinism in the relativistic limit.
In the Newtonian case, the effective non-determinism stems from two main arguments: the first is that the heat bath is too complex for the difference between deterministic motion and chaotic behaviour to be measurable. The same is true in special relativity, more so because of the issues of relative simultaneity.
The second argument notes that, in the limit of Brownian motion, the tagged particle is often much heavier than the constituents of its environment and thus the time-scales for the measurable behaviour of the Brownian particle are much longer than the time-scales needed to measure changes in the environment. In effect, the heat bath moves too quickly for its behaviour in the time scale of the tagged particle to be anything but random noise. This argument is less obviously true in special relativity: the faster moving environmental particles experience time more slowly than the tagged particle. However these complications vanish once one affixes the lab frame to the instantaneous rest frame of the tagged particle. In this frame of reference, the particles' subjective experience of time is irrelevant to the statistics of their motion.
Thus we fully expect the relativistic generalized Langevin equation to behave as "effectively stochastic" in many circumstances. It is possible, that the circumstances differ from the ones in which the Newtonian equation would have been effectively stochastic, which would require further investigation in future work.

B. Lorentz covariance
We are now equipped to consider the Lorentz covariance of the Langevin equation. Our first task is to show that the starting particle-bath Lagrangian is indeed Lorentz-covariant. This is tantamount to show that the action corresponding to Eq. (13) is Lorentz invariant: a scalar, a vector, etc. We thus consider the action associated with our starting particle-bath Lagrangian: The integral decouples into a series of terms: the kinetic terms of the bath eigenmodes and the tagged particle, the external Lorentz force acting on the external particle, and the interaction between bath eigenmodes and the tagged particle.
The kinetic terms simplify to integrals with respect to proper time of Lorentz scalars, e.g.
and similarly for the bath eigenmodes. The external force acts precisely like the electromagnetic force (which historically was the first Lorentz covariant interaction discovered), meaning that the V (x,ẋ, t) term is also Lorentz covariant.
Finally, since we model our interaction as a propagation of null-separated events, the interval, regardless of pre-factors would be a Lorentz scalar -zero. However, had we not restricted ourselves to null-separated events, the integral would become where in this case the δτ i is the Lorentz invariant proper time interval separating the interaction events. This being the final term in the starting Lagrangian proves that we have constructed a Lorentz-covariant theory of the interactions. It is well-known that the equations of motion which correspond to a Lorentz covariant Lagrangian, are themselves automatically Lorentz-covariant [59]. However, the Langevin equation (19) from which we obtain (30) is an approximate solution to the equations of motion (15a). Thus we must also demonstrate that Eq. (19) and the resultant Eq. (30) is Lorentz covariant. In practice, since the Langevin equation is a covariant equation for non-relativistic 3-vectors, what we must demonstrate is that the components on each side transform as spatial components of a Lorentz 4-vector.
For Eq. (19), after a boost with velocity u in some direction from the initial inertial reference frame S into the reference frame S , only the components of q along the direction u are affected. All trigonometric functions transform in the same way, and consistently, so that the equation holds for q (t ) in S . Namely, theω transforms as the inverse of a time-scale, t on both sides transforms as a time-scale. If we left it at that, the right-hand-side phase would differ from the corresponding phase on the left-hand-side. Fortunately, we can absorb the difference between the two into the laws of transformation ofξ, as we already did when introducing this parameter. What remains is a length at the initial time q(0), a 4-velocity multiplying a time scale 1/ω, and an integral. The integrand transforms like a length multiplied by two timescales: s and 1/ω, which we "fix" by noting that g i must transform like the inverse of the square of a time-scale. Consequently, the left-hand side and the right hand side transform as three-components of Lorentz four-vectors.
For Eq. (30), much like in the previous case, let's proceed in a term-by-term fashion. Firstly, the term is a standard [36, Eq. 12] Lorentz-covariant inertial term.
Adding the fourth (temporal) component [60] -E/c this becomes the tagged particle's four momentum. Thus, in order to demonstrate that Eq. (30) is compatible with special relativity, we must show that the right-hand side also transforms as a 4-vector of the same kind as the left-hand side [61]. The three forces, F p , F ext and F r are spatial components of contravariant 4-vectors, which we imposed by construction, with the temporal components of the 4-forces being equal to zero. Thus, the only term whose Lorentz covariance must be demonstrated is the memory term. Specifically it must transform as a (contravariant) 4-vector. We shall adopt the same approach as [36]. Namely, we shall introduce the memory tensor K µν , and re-write the integral as with implicit sums over repeated indices (µ). This is a contravariant 4-vector, provided that the memory tensor transforms as a Lorentz covariant rank-2 tensor, which we again, impose by construction following [36]. Specifically it's a diagonal matrix, with null temporal components, and diagonal entries all equal to [ This allows us to rewrite our relativistic Langevin equation (30) in explicitly Lorentz-covariant form: where we still used coordinate time t since we work in the instantaneous rest frame of the tagged particle, hence we do not have a difference between proper and coordinate time. Hence, we could equivalently replace t with τ in the above equation.
C. Slow (non-relativistic) heat bath limit Before we continue, it's important to consider whether Eq. (30) is compatible with the Langevin dynamics as is already known from previous works, by studying the relevant limits.
The first step is to analyze how Eq. (30) differs from the fully-Galilean limit, i.e. Eq. (6). There are three key differences: (i) a γ factor in the acceleration term, (ii) differences in how the stochastic force and the friction kernel are defined, and (iii) an extra term F r that is proportional to x(t). The first difference is trivial, in the non-relativistic (Newtonian) physics, γ → 1. The primed version of the friction kernel, and the stochastic force, Eq. (29a) and Eq. (29b), differ from Eq. (7) and Eq. (5), in a peculiar way. If we were to treatξ i (t, x(t)) as a small parameter, and expand the relevant equations, we would see that all the differences between the non-relativistic terms of the Langevin equation and the relativistic counterparts derived above are of order O(ξ i ).
That is, we have verified that in the Galilean nonrelativistic regime,ξ i → 0 andω → ω, upon taking the limit of the slow heat bath, i.e. γ(q) → 1. Hence, upon further setting γ = 1 in the acceleration term, the Galilean Langevin equation is correctly recovered in the appropriate limits by our Eq. (30).
Furthermore, starting from the general relativistic Eq. (30), one can choose to reduce the equations of motion of the bath Eq. (15b) to their Galilean counterpart Eq. (2b), while not taking the limit γ(ẋ) → 1. Therefore the Langevin equation reduces to where we notice that all the force terms are almost exactly the same as in the Galilean Langevin equation, and the term F r → 0. However, γ can still be 1, implying that the tagged particle moves at relativistic speeds in an otherwise non-relativistic bath. This limit given by Eq. (38), is equivalent to the equations derived by Plyukhin [35], Debbasch et al. [37], which describe the motion of a relativistic tagged particle embedded in a heat bath of non-relativistic oscillators.

D. Symmetries
During our derivation we have pointed out the necessity of lettingξ i (t, x(t)) be a function of time and of the trajectory of the tagged particle. This manifests as breaking of the following symmetries. If ∂ξ ∂t = 0, we no longer deal with harmonic waves, and the propagation is no longer exactly periodic, thus violating timetranslation invariance. For anything but the weakest coupling, ∂ξ ∂x = 0, if we were to shift the origins of both q and x, we wouldn't necessarily recover the same equation. Similarly, the inversions of time t, or of space, for q and x, would not recover the same equation. Nonetheless, the simultaneous inversion of both, would leave the equation unchanged if we were to extend the definition x). In other words we have lost parity covariance and time translation invariance.
However, as we shall see further below in this section, the loss of inversion symmetries is a natural consequence of combining statistical mechanics with special relativity, and the translation invariance isn't lost in the strict sense, just it doesn't appear to affect the length-like parameters of the problem.

Time inversion and parity
Let's first address the less controversial of the two apparently "lost" symmetries. Time and space inversion symmetry are violated individually. However, since Eq. (15b) is invariant with respect to a simultaneous inversion of both, and the function which approximates the general solution (Eq. (19)), is too, we can expect Eq. (30) to be invariant with respect to inversion of both time and space, combined together. This is not uncommon: CPT [62,63] is widely considered to be a fundamental symmetry of the relativistic theories while symmetries with respect to time and space inversion are individually violated even in simple Newtonian cases: a disk spinning clockwise is spinning counter-clockwise in the mirror universe, and a counterclockwise spinning disk is spinning clockwise if we reverse the arrow of time.
On the larger scales, one also has the thermodynamic arrow of time, which reflects the fact that irreversible processes flow in the direction of increasing entropy. This is a consequence of the breaking of time inversion invariance which is already at play in the Langevin equation [1]. However, we also predict a loss of parity invariance, which is consistent with recent observations [64]. It is yet to be explored at this point if the scale or extent of parity violations is similar for the experimental observations and our theory, but if it is, this might suggest that the thermodynamic arrow of time induced by dissipation leads to a thermodynamic chirality of the structures in the universe, that may have imprinted early in the universe's development, when the relativistic Langevin equation more accurately reflected its behavior.

Space-time translation invariance
In this section we shall briefly address the apparent breaking of translational invariance in Eq. (30). Let us first recall what kind of translational invariance is relevant here. Let us rescale the variables Then we substitute into Eq. (1), to yield see also Ref. [65]. This Lagrangian is invariant for V (x) = 0 under global translations of both the tagged particle and the bath, as is well known for the Caldeira-Leggett model [44,65]. Its equations of motion and the respective solutions, are left unchanged if the origins of x and Q are both displaced by the same arbitrary vector A ∈ R 3 . The same cannot be said of Eq. (30), primarily becausē ξ i depends on x, and the displacement of the origin of x would then lead to a change in the equation.
It is important to recall that, in the original Caldeira-Leggett Lagrangian or in general, we do not have enough terms to produce a complete square as in Eq. (40). The missing counter-terms are proportional to x 2 , and a priori we have no physical interaction that would provide them, apart from the special choice of coupling constants that leads to Eq. (40) shown above. So, as is common practice [44,65], rather than assuming a physical model for the interaction, we have imposed the translational invariance on the Lagrangian, and effectively swept the difference between the true Lagrangian corresponding to our model and the one that is translationally invariant into V (x). Thus, translational invariance is not a natural symmetry of the Caldeira-Legett model although it can be recovered for a suitable choice of coupling constants [65].
While the trick of adding counter-terms gives the illusion of translational invariance, it does not in fact resolve a real problem. We could have, by virtue of including a vector potential in Eq. (13) renormalised the terms to not contain any apparent violations of translational invariance, but that would obscure the fact that the spatial coordinates are not, in fact, coordinates. It should be clear that in Eq. (30) there is actually no real breaking of translational invariance with respect to coordinate reparametrization: x really is a displacement from a physically significant location, hence both x and q i should individually be invariant with respect to changes of coordinate origins. This consideration becomes even more relevant if one considers that, as will be elaborated on below in Sec. VIII E, the coordinates x (and to some extent also q i ) are meant to signify "displacements" from an initial thermodynamic state.

E. Continuous spectrum of the bath eigenmodes
In this section we shall attempt to provide a foundation for the relativistic extension of a cornerstone of nonequilibrium statistical mechanics, the fluctuation-dissipation theorem [31,66]. We compare the predictions of the Galilean and special relativity principles in the context of statistical mechanics, for simple models.
We start by focusing on the definitions of the stochastic force and memory function: Eq. (7) and Eq. (5). The two sums can be regarded as Fourier transforms if one replaces i with an integral ∞ 0 ρ(ω)dω, by introducing a density of states ρ(ω), and promotingω i to a continuous variable. By analogy, we should expect there to be Fourier representations of the memory function, of the restoring force and of the stochastic force. Namely, which completes our derivation.
Of note is the following difference between Eq. (41) compared to their discrete counterparts Eq. (29): the indices i are all replaced with a parametric dependence onω by suitably introducing a density of states of the bosonic bath vibrations [1]. However, as in the discrete case, no residual dependence on frequencies is observed also in the continuum case: everything which depends on ω is integrated out.

Markovian limit
Of particular interest is the case wherein ρ(ω) = αω 2 , as for bosonic particles, and g(ω) = Const. If applied in the limit of low velocities, ignoring all relativistic O(ξ) effects, the memory function's integral reduces to an integral from zero to infinite of a simple cosω i t, thus K (t) = K(t) ∝ δ(t), as shown in [1]. Therefore, the entire Langevin equation becomes Markovian, in the sense that the viscous response at time t depends on the velocityẋ(t) evaluated at the same time t, and there are no memory effects.
In the relativistic case, we do not have a simple cosine function inside the integral in Eq. (41a), but an integral which also contains the dependence on the bath oscillator trajectory viaξ. Upon inserting the form (24) that we obtain from numerics-assisted parametrization of trajectories, the inner integral in Eq. (41a) can be eval-uated analytically, and gives an expression of the type Hence, upon assuming ρ(ω) = αω 2 , and g(ω) = Const, it appears possible to retrieve that K (t) is proportional to δ(t) also in the relativistic case, just like in the Galilean case discussed in [1]. Unlike in the Galilean case, however, the assumption g(ω) = Const appears particularly strong and untenable in the relativistic case, because assuming that the tagged particle can interact with the same coupling strength with all the oscillator baths (irrespective of their separation in space-time from the tagged particle) is at odds with the principle of locality (which states that an object is directly influenced only by its immediate surroundings).
Therefore, we expect that since g(ω) = Const is imposed by the principle of locality in special relativity for a given physical system, the fluctuation-dissipation relation associated with our relativistic generalised Langevin equation (38) must be non-Markovian, as discussed in the next subsection.

F. Fluctuation-dissipation relations
Here we shall demonstrate that the force F p indeed has characteristics of thermal noise, thus qualifying Eq. (30) as a Langevin equation.
Consider that the bath-particle system is subject to a large number of computer simulations. In each such simulation, the heat bath initial conditions are taken at random, such that the probability density function (PDF) is where E is the total energy of the system for the microstate (x,ẋ); i.e. a Maxwell-Boltzmann [1] or Jüttner distributed random variable [54,55]. In the Galilean case, the system would have been in thermal equilibrium with respect to a frozen or constrained system coordinate x(0). Under these assumptions one can prove the fluctuation-dissipation theorem by direct evaluation of time averages [1]. However, special relativity makes this route somewhat more complex: one has more terms, and evaluating any term withξ dependence requires knowing the trajectory. With reference to Eq. (29b), the first two terms in the sum on the r.h.s. are identical to the terms that one has in the non-relativistic case (c.f. Eq. (5)) and it was shown by Zwanzig that they give F p (t) = 0 for the non-relativistic stochastic force [9,39]. Here we have an additional term, the last term in Eq. (29b) proportional to x(0). In general, this term is non-trivial to evaluate since it depends on the bath oscillator's trajectory throughξ i (s). However, we can show, by analytical integration, that for the case where the trajectory is parameterised by Eq. (24) (where ξ = As(x − B) 2 ), which was obtained from the numerics in our simulations, the result of the integral over s evaluated at s = 0 is simply cosω i t. Hence in this case also the last term in Eq. (29b) averages to zero, and F p (t) = 0 also for the relativistic case. Also in the hypothetical case ξ, one obtains sinusoidal functions i.e. a combination C sinω i (t + 1 4 ) + D cosω i (t + 1 4 ), with C and D some t-independent constants, which also leads to F p (t) = 0. Therefore, based on the numerical data that we have forξ(s), we can conclude that the term F p (t) has zero average and qualifies as the stochastic term or the noise in our relativistic Langevin equation. For most situations we expect that the integral over s in Eq. (29b) evaluates to a sinusoidal function of ω i t upon integrating away the s dependence and evaluating at s = 0.
Proceeding in a similar way, since the last term of Eq. (29b) for the numerically -evaluated trajectory form given by Eq. (24) (and presumably also for other more complicated forms) leads to a function proportional to cosω i t, evaluation of the time-correlation function of F p gives a product of sinusoidal functions for all the terms present in Eq. (29b). Following the same strategy as in [1] at page 23, and of [9], the products of two sinusoidal functions of argumentω i t andω i t respectively, by applying the standard trigonometric identity for the product of trigonometric functions, lead to a generic function K of argument (t − t ), Since coordinate time is relative, the right hand side of the fluctuation dissipation theorem must depend on the velocity of the observer. If one could find a physically significant frame of reference, however, proper time can be used and t can be replaced with τ in the above relation. Analogy with the Galilean case would suggest that the product inside the averaging brackets is an inner product of two 4-vectors, while the right hand side of the Galilean fluctuation dissipation theorem is a rank 2 tensor. Thus the product of two fluctuations must indeed also be the outer product of the force 4-vectors.
Hypothetical deviations from Eq. (43) which may occur for more non-trivial s-dependencies of ξ are considered in Appendix E.
In more general and complex settings, for complex trajectory dependence ofξ(t, x), simple methods of integration may not be sufficient to arrive at a FDT in closed form. However, one could opt for a path-integral approach to the problem along the general lines of Ref. [10]. That task is, however, well-beyond the scope of the current study and is left for future work, along with ascertaining whether the second additional term in Eq. (E1) can be observed in certain conditions.
Another approach would be to consider the less extreme relativistic behavior, i.e.ξ → 0. As we have mentioned previously, this case has been studied extensively and relations analogous to the fluctuation dissipation theorem (also known as Einstein-Sutherland relations), were obtained e.g. in Ref. [29]. These completely ignore any effects that may be due to O(ξ) terms, so, while a useful point of reference, they are not totally useful for fully relativistic conditions.

IX. CONCLUSION
In this largely expository paper, we have provided a first-principles derivation, from a microscopic Caldeira-Leggett particle-bath Lagrangian, of a full-fledged relativistic Langevin equation, Eq. (30). In some of its limits, this more general relativistic equation recovers commonly accepted [34,35,37] extensions of the Langevin equation to relativistic media and relativistic weakly-interacting particles. By relaxing some of the commonly used approximations, led to a full-fledged and more general representation of stochastic processes in relativistic media valid for more strongly relativistic conditions.
The new fully-relativistic generalised Langevin equation Eq. (30), or in fully covariant form Eq. (37), contains a new force term F r which is trajectory-dependent and requires further investigation in future work. Also, the Fourier modes associated with the bath oscillators, are modified into a plane wave-like form which is naturally covariant, with a relativistic correction length-scale ξ with is both time and trajectory dependent. Based on our numerical data, this dependence is an off-set parabola of the form ξ(x, t) = At(x − B) 2 . The possible generality of this form has to be further ascertained in future numerical studies.
The new fully-relativistic Langevin equation may predict fundamentally new physics. Firstly, there is some evidence that the universe is chiral at the very largest scales [64]. Our equation suggests a mechanism for parity violation due to dissipation, so analysing the results of the aforementioned experiment may be of utmost interest. The relativistic correction terms of the same order, also predict the presence of a weak restoring force that tends to bring the particle to its initial state/position, that may be observable by the same methods by which Sakharov oscillations [67], and Baryon Acoustic Oscillations (BAO) [68,69], are observed in the CMB radiation. Also, this effective force vanishes upon moving to speeds that are small compared to the speed of light, and is therefore a genuine new effect due to the interplay of relativistic motion and dissipation, which otherwise vanishes in the non-relativistic limit.
We have also, in a way, reached the limit of what the Caldeira-Leggett particle-bath models can tell us about the underlying dynamics of relativistic media. One could account for more phenomena by considering: (i) suitable incorporation of proper and co-ordinate times; (ii) couplings that depend on the position of the tagged particle; (iii) harmonic modes of other kinds (e.g. magnetic confinement potentials that lead to pure harmonic oscillations, etc.) (iv) pair production; (v) proper field--theoretic extensions of the Caldeira-Leggett Lagrangian model by means of path-integral formalism for state-dependent diffusion processes [10].
Of particular interest are the questions of pair production and field-theoretic extensions. By definition of rest energy, at velocities where corrections that we have neglected may become important, the energies of the medium and the particle are sufficient to randomly produce pairs of charged particles. This extension may be of interest both for examining super-heated exotic objects, as well as the early universe. Further important applications are in the context of nuclear physics, where so far only the Galilean generalised Langevin equation has been used to describe fission processes [26], while clearly relativistic corrections may be important for the fission of hot nuclei. Finally, the results presented here could be a first step for a more extensive and detail analysis of thermodynamic aspects of relativistic systems [70], including heat distribution [30], and fluctuation theorems [31].
In future work, one could also attempt to produce the General-Relativistic (GR) extension of the Langevin equation by the same methodology presented in this paper. Finally, it will be interesting to compare results from the above derivations with certain limits of Langevin equations obtained using the Schwinger-Keldysh formalism within AdS/CFT approaches [71,72].

ACKNOWLEDGMENTS
Matteo Baggioli and Bingyu Cui are gratefully acknowledged for discussions and assistance in the numerical calculations.
W. J. Handley is acknowledged in assistance in identifying suitable cosmological applications of the results.
A. Zaccone gratefully acknowledges financial support from the US Army Research Office through grant № W911NF-19-2-0055.
Appendix A: The dynamical coupling model in detail To answer this question we must more deeply elaborate on the origin of the modes. The real underlying model is that of hard collision potentials between the tagged particle and the constituent particles of the heat bath. There is a different probability of each collision between each species of particle, and multiple possible modes of interaction that we all ignore. Every potential that would lead to both the tagged particle and the heat bath confined to the region that we discussed, entails describing both the generalised q i and the dynamical x as displacements from a common equilibrium. Thusẋ is always interpreted as a relative velocity with respect to the bulk velocity of the heat bath.
In the relativistic case, the collisions are highly local, so one might think that they can safely assume t i ≈ t. However, the interaction with the heat bath modes entails the interaction of the particle that was hit by the tagged particle with the rest of the heat bath, and generate a bosonic excitation. This interaction is not local in general. Instead, to avoid substantial complications in the Euler-Lagrange equations, we shall note that the difference between q i (t i ) and q i (t) depends on the extent of the heat bath that we consider: the smaller the quantity of particles participating, the smaller the effective distance between t and t i . Thus we restrict ourselves to interactions between local modes, where q i (t i ) ≈ q i (t). This has the added benefit of making virtually any of the initial conditions that were deemed acceptable in the Galilean case, to be useful in the Lorentz covariant case (cfr. IV D). versa, weaker coupling allows for broader variations in x, as seen in Fig. 4.
Firstly, we see that the constraints set by causality can be implemented by constraining the range of x for a fixed g. Secondly, the units we are using are that in which the characteristic oscillation frequency is unity. In the LHC, the relevant frequency would the plasma frequency of the electrons in the cloud. By taking the average electron number density [73] n e ≈ 10 7 cm −3 = 10 13 m −3 we get a plasma frequency ω pe = n e e 2 0 m e + 3k B v 2 th,e ≈ 1.7 kHz, thus, x = 1 translates to approximately 150 km, six times the circumference of the LHC.For a star like e.g. Sol, the plasma frequency is much higher, and the size of the core is much closer to 1.0 × 10 5 km [74]. Thus one does not expect many relativistic effects at the LHC, but does at the core of the sun.
In either case, the distance at which the boundary conditions can be considered constant are far beyond the distance at which gx > 2 could cause issues.
It must also be noted, that g = 2 is an unrealistic coupling strength. The Caldeira-Leggett Hamiltonian corresponding to Eq. (1) and Eq. (13), is responsible for the heat capacity of the medium. If for the case of ω = 1, g = 2 were the true value, the heat capacity of the medium would be strongly affected by the presence of the tagged particle. This is not what we observe in the Galilean case, thus ruling out such a possibility for the Lorentz covariant case. The reason why we us values of g = 2 is purely illustrative, as the true solution to the Euler-Lagrange equations can only depend on gx and the plots are easier to read for reasonable values of x.
Appendix C: Equivalence of time-dependent frequency and phase Here we demonstrate that for an arbitrary phase relationship, the time-dependence of the frequency can be without loss of generality be shifted onto a time dependent phase.
Let φ(t) and ω(t) be some time-dependent functions. We can always re-write the phase of the form as where Ω = Const. To do so, recognise that the change of variables has two degrees of freedom, Ω and Φ, but only one constraint -Eq. (C2). Thus, by setting Φ(t) ≡ −ω(t) (t − φ(t)) − Ωt (C3) one can easily verify that Eq. (C2) holds as an identity. Thus one can always choose to deal with a timedependent phase in lieu of a time-dependent frequency, where, as in the following appendix, the integration is greatly simplified by not considering time-dependent frequencies in addition to an already implicitly timedependent phase.