Superradiance from a Relativistic Source

We construct a model of cooperative superradiant emission from a highly relativistic multi-particle source. We revise an existing model of the literature for a relativistic two-level particle, and construct from it a Hamiltonian describing relativistic velocity dependent multi-particle superradiance. We adapt the standard diagrammatic framework to compute time evolution and density operators from our Hamiltonian, and demonstrate during the process a departure from standard results and calculation methods. In particular, we demonstrate that the so-called vertical photon result of the literature is modified by the relativistic Lorentz factor of the sample; we also introduce a set of coupled differential equations describing certain propagators in the velocity-dependent small sample framework, which we solve numerically via a hybrid fourth order Runge-Kutta and convolution approach. We demonstrate our methods for the simple case of two highly relativistic particles travelling with slightly differing velocities simulated at varying relativistic mean sample $\beta$ factors, and evaluate velocity coherence requirements for a sample to demonstrate enhanced superradiant emission in the observer frame. We find these coherence requirements to become increasingly restrictive at higher $\beta$ factors, even in the context of standard results of relativistic velocity differential transformations.


I. INTRODUCTION
In the quantum optics phenomenon of superradiance (SR) a collection of excited particles evolves into highly entangled states which generate greatly enhanced radiative emission over that predicted by the spontaneous decay of its constituent particles. In his seminal theoretical treatment of SR, Robert Dicke [1] first described the enhanced emission process as a cascade through symmetric states isomorphic to the maximal total angular momentum states from the elementary theory of the addition of spin-1/2 particles. If |S (N, m) denotes the symmetric superposition of the possible quantum states of N particles having m excited and N − m in the ground state (for example, |S (2, 1) = (|e, g + |g, e ) / √ 2), then the Dicke model describes the SR emission process by the cascade |S (N, N ) → |S (N, N − 1) → · · · → |S (N, 0) . (1) The probability of transition P k,k−1 = |S (N, k) → |S (N, k − 1) in the Dicke model is computed from a perturbation analysis of the Hamiltonian coupling the particle to the quantized radiation field via the relevant transition operator (electric dipole, magnetic dipole, or other), from which the central state k = N/2 [2] is found to couple most strongly to the radiation field. The full cascade through all states is characterized by the two salient SR features of a delay and enhancement in peak emission intensity. These theoretical predictions have been well verified by laboratory demonstrations of SR [3] since their original prediction by Dicke. Additionally, more sophisticated theoretical models of SR have been * cwyenber@uwo.ca developed which generalize the Dicke model to extended volumes; which introduce realistic dephasing, relaxation, and pumping terms; and which lend themselves to more efficient numerical simulation. SR indeed possesses a rich experimental and theoretical history at the scale of laboratory quantum optics. SR has also been recently applied to astrophysical environments, where it has been used to model maser flares [4][5][6][7] and fast radio bursts [8,9]; however, such environments lead to regimes of SR not realized in the laboratory and thus as of yet untreated theoretically. Some initial progress has been very recently made towards the theory of wide inhomogeneous broadening of SR, and two O (n) complex algorithms for simulating transients of n velocity channels have been derived [10,11] and applied [11] to astrophysical environments; however, there presently exists in the literature no relativistically covariant formulation of SR. Such a model could prove essential to analyzing, for example, the proposed highly relativistic SR events of [12]. The objective of this work is to develop a simple model of SR emission from a small sample of particles travelling with differing relativistic velocities, and with it to determine the degree of velocity coherence required in order to measure enhanced SR emission in the observer's frame.
There exist various representations of SR which suggest differing approaches to this problem. In the socalled Maxwell-Bloch model of a large population in an extended medium, for example, the complicated quantum mechanical time evolution is shown [13][14][15] to very quickly transition into a stage described by an ensemble of population inversion, polarization, and electric field trajectories determined by a relatively simple set of partial differential equations. In this later stage a correspondence to the classical electric and magnetic fields may be made, which could potentially be transformed to the observer's frame by the usual relativistic transformation of the electromagnetic field tensor [12].
As a simple and fundamental first step towards a relativistic description of SR, however, this work does not seek to model a large population. Moreover, it is not a simple exercise to generalize the aforementioned result [13][14][15] regarding the emergence of an ensemble of trajectories to observers of a relativistic source. We seek rather to develop a relativistic model of the proper quantum mechanical evolution of a small number of particles capable of SR emission.
Constructing a simple relativistic model of the quantum evolution of even a small number of particles presents numerous challenges. First, observer frame measurements cannot be described by the classical relativistic transformation of the electromagnetic field tensor. Instead, the time evolution of the observer frame quantum state must be computed from a Hamiltonian built upon a covariant formulation of the matter-field interaction. For this purpose we first introduce in Section III an existing covariant model [16] of a two-level particle established in the literature.
The second challenge faced by a relativistic model is the interpretation of the cascade through symmetric states. As detailed momentarily in Section II, the small sample Dicke model [1] identifies the intensity transient with the average emission rate over a distribution of stochastic cascades, where each realization involves transitions between the symmetric states at definite times [17]. In the case of photons emitted from a relativistic source, however, Dicke's interpretation becomes nuanced when considerations of time dilation or other relativistic effects are introduced into the stochastic cascade analysis. We therefore use the full diagrammatic method of quantum optics to evolve the matter-field system from the covariantly-derived Hamiltonian, and depart from the concept of stochastic discrete-time transition events in favour of a purely mathematical description of the system's time evolution operator.
The third challenge involves the choice of basis for the Hilbert space describing each particle's center of mass coordinate. SR is highly dependent upon the relative position of the particles, which therefore suggests the use of a position basis; however, we are interested in velocity coherence effects which are more naturally realized in a momentum basis. We define a compromised representation of the problem in Section IV A, and proceed to apply the diagrammatic method to a relativistic twoparticle SR system in Section IV C.
After building these foundations we turn in Section IV D to apply our revised diagrammatic method to the Hamiltonian introduced earlier in Section III. We finally solve the simple case of two particles travelling with offset relativistic velocities in Section V, where we establish a metric for SR enhancement as a function of velocity separation in the observer frame. We evaluate this metric for samples travelling at various highly relativistic speeds and we analyze the resulting velocity coherence requirements in the observer frame in the context of the standard relativistic transformation of velocity differentials.

II. EXISTING MODELS OF SUPERRADIANCE
Dicke [1] first analyzed the so-called small sample limit wherein all particles are contained within a volume of dimensions much less than the emission wavelength. He considered the expectation values ρ m of populations of all states |S (N, m) as evolving according to the transition rates {P m,m−1 }. The resulting responses of the populations ρ m were used to determine the total intensity transient which peaked after some delay at a maximum emission rate greatly exceeding that predicted by spontaneous emission. Dicke continued in his foundational work to describe SR within an extended volume, to make realistic revisions to his model, and to analyze the effect of momentum back-reaction during photon emission.
Later models of SR approached the problem from either the master equation of quantum optics or the Heisenberg picture evolution of coarse grained inversion and polarization operators coupled to the quantized radiation field. These approaches lead to the Maxwell-Bloch equations briefly mentioned in Section I. The Maxwell-Bloch equations are extremely useful for large samples which are intractable with purely quantum mechanical methods; however, they do not describe the initial quantum mechanical stage of SR. Moreover, our objective is to establish the most simple model possible for asking fundamental questions of a relativistic SR sample. The Maxwell-Bloch model is overly complicated for our present purposes.
We choose instead to generalize the small sample Dicke model [1] to a relativistic sample; however, as briefly discussed in Section I, we deviate from his original methods. The Dicke model imagines the system evolving via discrete transitions through definite quantum states, essentially constructing a stochastic time sequence of events. Such an analysis intrinsically involves a non-relativistic observer and it is not at all clear how to translate this stochastic time sequence into a relativistic framework. Rather than attempt to analyze the role of relativistic effects in the cascade, we choose to apply the rote diagrammatic method of quantum optics to compute the fully quantum mechanical time evolution operator from our relativistically-derived Hamiltonian, and to introduce observer frame intensities only at the end of the analysis as measurements upon the formally evolved quantum state of the system.

III. THE BOUSSIAKOU, BENNETT, AND BABIKER RELATIVISTIC MODEL OF A TWO-LEVEL PARTICLE
There exists in the literature a quantum mechanical model of a two-level particle [16] derived from a relativistically covariant Lagrangian describing the matterfield interaction. Canonical quantization of the resulting observer frame Hamiltonian yields a Klein-Gordon-like description of the particle's center of mass, as well as relativistic γ factor corrections to certain matter-field interaction terms. We refer to this Hamiltonian as the Boussiakou, Bennett, and Babiker (BBB) model of a relativistic two-level particle. We state the Hamiltonian below in this Section but leave its detailed synopsis to Appendix A.
In a confirmation of the BBB model's relativistic foundations, it is shown in Section IV to properly describe the relativistic time dilation of the spontaneous decay of an initially excited particle, now in our diagrammatic framework (a recovery of the result already obtained in [16]). The BBB model is well suited to the task of constructing a relativistic SR model. We describe in Appendix A its derivation at a high level as well as features and slight revisions pertinent to our work; for its detailed derivation, see [16]. The BBB Hamiltonian in the observer frame reads, where numerous terms require definition. The center of mass motion is described by the Klein-Gordon term for P = K the center of mass three-momentum of a Klein-Gordon particle of energy The operatorsπ † andπ raise and lower the internal energy state of the particle, respectively. Within the free field HamiltonianĤ 0 f the operatorsâ † k andâ k create and annihilate (respectively) photons of mode k and linear polarization ǫ k , where ǫ ⊥ k = (k/k)×ǫ k . The single-photon electric field coefficient is for V Q a fiducial field quantization volume and ǫ 0 the permittivity of free space. We imply in all summations over field modes k a summation also over two orthogonal polarizations.
Within the interaction termV the operatorsL † k and L k raise and lower the Klein-Gordon center of mass momentum by k, respectively; i.e., The notationL k (L † k ) is unique to this work. The corresponding terms in [16] read exp ±ik ·Q for a center of mass position operatorQ; i.e., For simplicity we have suppressed magnetization terms and restrict our work to the case of a particle possessing a rest frame electric dipole moment d ′ . The dipole moment vector is composed of d ′ and d ′ ⊥ parallel and perpendicular to the particle's velocity, respectively.

IV. RELATIVISTIC TWO-PARTICLE SR: THEORY
We develop in this section the diagrammatic representation of a relativistic two-particle SR sample. We assume familiarity with the diagrammatic representation of the SR time evolution operator; however, we provide a summary of the method in Appendix B.
The velocity dependent small sample SR problem is complicated by the localization of particles to within the wavelength of emission. Our development of the diagrammatic method therefore deviates from conventional methods, as we introduce in Section IV A a compromise between the position basis and the momentum basis which departs from a standard convolution structure [14].

A. Center of mass representation and localization considerations
Seeking in this work a simple model of relativistic SR, we restrict our analysis to the small sample limit defined such that all particles are contained within a volume of dimensions much smaller than the wavelength of emission λ. The Heisenberg uncertainty relation then implies that ∆P i ≫ / (2λ) in the initial momentum of any particle. This uncertainty could be realized, for example, by an initial Gaussian superposition in the center of mass Hilbert subspace. A particle initially excited with no photon in the radiation field could be described by the (unnormalized) state (12) Alternatively, the center of mass could be considered well-localized at some position R i such that The literature treatment of diagrammatic SR [14] conventionally adopts the latter approach, working in the position basis and choosing to neglect the changes in momentum caused by the actions ofL k orL † k in the interaction termV . Alternatively, in a perfectly formal execution of such calculations, the initial state of the system would be defined by a construction similar to Eq. (12), and modifications byL k andL † k to all P ′ in the superposition would be tracked during the system's time evolution. This exercise would greatly complicate calculations.
In our study of velocity dependent relativistic SR the initial particle momenta are critical to the velocity coherence analysis. We choose the following compromise between the strictly formal use of a superposition (such as Eq. (12)) and the literature approach of working in the position basis without tracking the particle momenta. We work in the momentum basis but we implicitly assume always that a center of mass momentum state |P is in fact a coarse distribution about a central value P. The intuitive picture is that momentum states overlap with coarseness much greater than /(2λ).
In practice the coarse momentum assumption is important to transition matrix elements. Suppose two particles starting with momenta P 1 and P 2 evolve according to the BBB Hamiltonian. We omit their momenta in our description of the state of the system, but imply always that their values are consistent with conservation of momentum through transition matrix elements. For example, consider the matrix elementṼ g,k;1,0 (t) for a two-particle system, where the particle energy state |1 corresponds to only the first particle excited, |g to both particles in the ground state, and |0 /|k to the zero photon / single photon of mode k states. This matrix element is implied to denote a transition between any two coarse momentum states consistent with momentum conservation in their central values; i.e., V g,k;1,0 (t) = g, k, P 1 − k, P 2 | Ṽ (t) |1, 0, P 1 , P 2 , (14) which is related to the non-interaction picture matrix element bỹ where α k1 = 1−cos θ k1 V 1 /c for θ k1 the angle between P 1 and k (this result is obtained in the next section, viz. Eq. (24)). Whether momentum conservation is expressed by modifying the momentum in the bra or the ket of Eq. (15) does not affect the results of later calculations. The situation is further nuanced by considerations of the cascade history. It is argued in the literature [14] that the density operator formalism used in Section IV D below does not need to reference multi-photon states; however, when particle momenta are carefully tracked (which is not the case in [14]), prior photons emitted during the SR cascade affect later transition matrix elements. For example, suppose that the sample begins doubly excited and that particle 2 first emits a photon of mode k ′ , such that its momentum is equal to P 2 − k ′ prior to the transition process of Eq. (14). The resulting matrix element of Eq. (15) would then be modified. In [14] the particle momenta are suppressed and such complications are avoided; this suggests that their methods are deficient for our present (relativistic) velocity dependent purposes, which necessarily reference the particle momenta.
We take the following approach. Instead of complicating the diagrammatic methods of the literature by tracking the particle momenta and generalising existing techniques to multi-photon states, we choose to suppress the particle momenta everywhere, save for transition matrix elements such as shown in Eq. (15). This is admittedly ad-hoc, but we are reassured by the success of the approach demonstrated in later sections. The model will be found to properly describe the continuum between the totally coherent limit where the particle velocity difference ∆v vanishes and the independent spontaneous emission limit ∆v → ∞; additionally, it will be found to accurately recover time dilation in the totally coherent case.

B. Convolution structure and the self-energy term under the BBB Hamiltonian
We now establish matrix elements for the time evolution operator, beginning from its standard series expansion identities of Appendix B. As demonstrated in the literature [14], Eqs. (B6) and (B7) possess convolution structures. For example, the U e,0;e,0 matrix element in the non-interaction picture may be expressed via Eq. (B6) as where f k depends upon the initial velocity and dipole moment orientation of the particle according to F e,0 is the free theory propagator for a particle in the excited state and F g,k the free theory propagator for a particle in the ground state with a single photon of mode k in the radiation field. A Laplace transform replaces the convolution structure of Eq. (16) with an algebraic relation; denoting the propagator Laplace transforms as F e,0 (s), F g,k (s), and U e,0;e,0 (s), we have the relation In the interaction picture Eq. (18) reads The summation is conventionally referred to as the (Laplace domain) selfenergy term, and it now deviates from the literature due to the velocity dependence as well as our use of the BBB Hamiltonian.
We require an expression for Ξ for two different orientations of the dipole moment d ′ relative to the initial velocity of the particle. Recall that the mode identifier k is implied to extend over two polarization directions, which we now explicitly identify as with associated perpendiculars by our earlier definition ǫ ⊥ k = (k/|k|) × ǫ k . When k is parallel to V i (and therefore also to d ′ ) Eq. (21) becomes undefined, but by the symmetry of such a situation two orthogonal polarization directions may be arbitrarily assigned (in the plane perpendicular to k) without affecting what follows.
From the BBB Hamiltonian we have that ω e,0 = γM c 2 / + ω ′ 0 /γ and that for an initial center of mass velocity V i (with which the γ factor is identified). We thus have for an angle θ k between k and V i . Let us define for We consider first the case where V i is parallel to d ′ , in which case only the ǫ kθ polarizations yield non-vanishing interaction matrix elements such that the self-energy term of Eq. (20) becomes .
(27) Repeating the arguments of [14] (now for our relativistic Hamiltonian), in the absence of interaction equation (19) possesses a pole at s = 0. We therefore expect, in the perturbed system, that the neighbourhood of s = 0 gives the main contribution to the sum over k such that we may approximate the self-energy factor of expression (27) by evaluating it at s = ǫ and taking the limit ǫ → 0 [14]. Using the Sokhotski-Plemelj identity [18] (where Pr denotes the Cauchy principal value) and expressing the sum over k as an integral via for Γ ′ 0 the rest frame spontaneous decay rate. In the time domain Eq. (30) reads We have discarded the imaginary principal value part leading to the Lamb-like energy level shift, which we assume has been included in the definition of ω ′ 0 . The fact that the self-energy Ξ obtained by our diagrammatic approach corresponds to half the (time dilated) spontaneous emission rate is understandable if we now use Ξ to solve forŨ e,0;e,0 in Eq. (19) and obtaiñ U e,0;e,0 (s) = 1 which has inverse Laplace transform Since |Ũ e,0;e,0 (t, t 0 ) | 2 corresponds to the probability of finding the particle excited after a time ∆t = t − t 0 , Eq. (33) tells us that the particle decays with time constant Γ ′ 0 /γ, which is the same decay behaviour calculated in [16]. Our time evolution operator approach, necessary to our SR application, indeed recovers the results obtained by Boussiakou et al. [16]. Note that this result is their demonstration that the BBB Hamiltonian successfully describes time dilation of spontaneous emission, which was the driving motive for our use of the BBB Hamiltonian in the present work.
It remains to yet evaluate Ξ when V i is perpendicular to d ′ . The calculation is similar to the above work and ultimately leads again to Ξ = Γ ′ 0 /(2γ). This is not a new result, but rather a reproduction of the conclusion of [16], unique only in our use of the time evolution operator and Laplace domain methods.

C. Relativistic two-particle SR propagators
We now have at our disposal all of the tools necessary to compute the propagators of a relativistic two-particle SR sample. The methods of this section depart from the standard diagrammatic methodology of the literature, as we shortly demonstrate that the velocity dependent problem in our coarse momentum framework (necessary to the small sample SR problem) no longer possesses a convolution structure in its propagator expansions.
We denote the doubly excited particle state |e , the state of only the first particle excited |1 , the state of only the second particle excited |2 , and the doubly ground state |g . In this section we in fact compute only the zero photon propagators, which will prove sufficient for our density operator work in Section IV D; however, we do provide the relativistic propagators for single particle transitions into single photon states in Appendix D. The results of Appendix D can be used to compute multiparticle photon emission propagators where future research demands. Because no photons are present in this section, we will omit the radiation field state in our notation by writing, for example,Ũ e;e instead ofŨ e,0;e,0 .

Doubly excited propagator
Let us first find the propagator for the two particles to both remain in the excited state. We obtain from the recurrence relation of Eq. (B6) the relatioñ We would be repeating ourselves to calculate the selfenergy parts within these convolutions: upon expanding the interaction matrix elements of Eq. (34) the central summations become the time domain equivalents to the calculation which followed Eq. (27) and yielded for γ 1 and γ 2 the relativistic factors of particle 1 and 2, respectively. We can approximate for γ the Lorentz factor of the average of the velocities, where the above approximation is correct to first order in ∆v/c for a velocity difference ∆v between the two particles. We finally have that Note that squaring the above result confirms that the doubly excited state again decays as expected by time dilation.

Singly excited propagators
We now calculate the propagator for the system to commence in the singly excited state of either particle and remain in that state at a later time. This situation differs from existing propagator calculations and results in the literature, as we shortly demonstrate a loss of convolution structure rooted in our use of the coarse momentum framework necessary to velocity dependent SR. We compute the propagator for the state |1 where only the first particle is excited; the propagator for |2 may be easily obtained by symmetry.
The recurrence relation of Eq. (B6) reads, for this matrix element, The first integration term on the right side of Eq. (39) describes the process of propagating from singly excited to singly excited, transitioning to the ground state, and transitioning back to the singly excited state. This term is familiar to us, becoming in the non-interaction picture a simple convolution structure containing the usual selfenergy expression Ξ (t 2 − t 1 ). It simplifies to The second line of Eq. (39) departs from the convolution structure. This line describes a process by which (i) the system transitions from the first particle singly excited state to the second particle singly excited state; (ii) the second particle emits a photon; and (iii) the first particle absorbs the photon. This physical picture of transferring a photon between the two particles suggests that this term will play a central role in SR velocity coherence.
It will greatly simplify what follows if we now choose a relative orientation between the initial velocities and dipole moments of the particles. Let us suppose that the two particles are travelling in the same direction with mean relativistic velocity v (|v| = v = βc) and velocity difference ∆v = v 2 − v 1 (|∆v| = ∆v) as measured in the observer frame, and that their dipoles are oriented parallel to each other but perpendicular to their shared velocity direction. The two interaction terms within the central summation in the second line of Eq. (39) then become where α k1/2 are the α k factors for v 1/2 , respectively.
The central two exponentials of Eq. (41) form the convolution violating term which cannot be expressed as a function of t 2 − t 1 . The occurrence of this term may be traced back to our use of the coarse momentum framework, wherein we chose not to evolve the distribution of momentum eigenstates describing a localized particle. Had we evolved such a distribution, we would have indeed retained a convolution structure between every adjacent occurrence of V in the series representations of propagators; however, we would also have been required to track the center of mass momenta in our kets, which would have dramatically increased the number of matrix elementsṼ nm . Despite removing the convolution structure, then, our approach remains more efficient than describing the state of the particle's center of mass via a superposition such as shown, for example, in Eq. (12). We next evaluate the summation over k in the second line of Eq. (39). The details of this lengthy calculation are left to Appendix C, where we find the entire second line of Eq. (39) to be equal to where We refer to C SR, rel β,∆v as the relativistic SR velocity coherence kernel, which is an important tool in this work. A plot of C SR, rel β,∆v 2 is shown in Fig. 1 for three values of β; we highlight features of C SR, rel β,∆v at the end of this section after establishing its relationship to the two-particle propagators.
Eq. (39) now reads which may be differentiated with respect to t to obtain the first order differential equation Repeating all of the above procedure from Eq. (39) onward for each of the possible single excitation propagators, we ultimately arrive at the set of coupled differential After first numerically calculating the relativistic SR velocity coherence kernel, Eqs. (47)-(50) may be solved by any forward stepping numerical method from the initial conditionsŨ j;i (t 0 , t 0 ) = δ ij . We use a fourth-order Runge Kutta stepping scheme for this purpose in later sections. We point out the following reassuring properties of the relativistic SR velocity coherence kernel (see Fig. 1): • In the limit as ∆v → ∞ the kernel C SR, rel β,∆v → 0; in this case Eqs. (47)-(50) are decoupled across the particles' identities such that they evolve independently.
• For non-relativistic β, C SR, rel β,∆v has characteristic width ω ′ 0 (∆v/c)t ′ ≈ 1 (to order of magnitude); i.e., it yields SR enhancement between two particles over an SR timescale T d only if their Doppler separation is less than the characteristic frequency 1/T d established by said timescale.
• The width of C SR, rel β,∆v narrows with increasing β; i.e., velocity coherence requirements become stricter at increasingly relativistic velocities.

D. Diagrammatic density operators for relativistic two-particle SR
Having in our possession propagators for the relativistic two-particle SR sample, we are now able to describe the system's density operator evolution. In this section we parallel the steps followed by Benedict et al. [14], but modify the standard method as our velocity dependent relativistic SR model requires.
A sample commencing in the initial state |Ψ (t 0 ) has a density operatorρ (t) = |Ψ (t) Ψ (t) | which evolves according to We will be interested shortly in the diagonal elements ρ m,m (t). Consider for a moment the simpler case of the single photon diagonal matrix elements from a singleparticle system evolving according to the BBB Hamiltonian's time evolution operator. Using the non-interaction picture version of Eq. (B5) to obtain the propagator we have that ρ g,k;g,k (t) = |f k | 2 2 [F g,k ⋆ U e,0;e,0 ] (t − t 0 ) × ρ e,0;e,0 (t 0 ) F * g,k ⋆ U * e,0;e,0 (t − t 0 ) . (53) Eq. (53) is not a convolution between the first and second lines; however, we are interested not in the diagonal single photon elements themselves, but in the trace over them. By a standard result of the literature [14] the trace operation restores the convolution structure. We now recover this result in our present BBB Hamiltonian and see how it is modified. The trace over Eq. (53) becomes ) U e,0;e,0 (t 1 , t 0 ) U * e,0;e,0 (t ′ 1 , t 0 ) ρ e,0;e,0 (t 0 ) ) U e,0;e,0 (t 1 , t 0 ) U * e,0;e,0 (t ′ 1 , t 0 ) ρ e,0;e,0 (t 0 ) symmetry, as well as recognizing that the integrand con- ; and from the fourth line to the fifth by exploiting the delta distribution in Ξ. Note that U g,0;g,0 (t, t 0 ) = e −iEP i (t−t0)/ for E Pi the kinetic energy of the Klein-Gordon particle possessing momentum P i .
Eq. (54) tells us that the so-called vertical photon result is retained in the BBB Hamiltonian, but that it now introduces a factor of Γ ′ 0 /γ. This result is critical to our density operator work to follow and is depicted pictorially in Fig. 2 Returning now to two-particle relativistic SR, the vertical photon result with accompanying Γ ′ 0 /γ factor is ubiquitous. Consider, for example, the quantity which represents the probability of finding the system at time t in the state with only the first particle excited. Similar calculations as preceding the single particle emission result of Eq. (54) find that the four diagrams of Fig.  3 together describe the quantity ρ 1;1 (t). These four diagrams exhaust the possible permutations of propagators obeying the following constraints: first, the leftmost side of a diagram must commence in both its top and bottom propagators in the excited state, corresponding to the starting state of the system at t = t 0 ; second, transition across any vertical photon line introduces a factor of Γ ′ 0 /γ; third, in any region demarcated by vertical photon partitioning, propagators must couple states of like total number of excited particles, decreasing by one when crossing a vertical photon partition; and fourth, the singly excited propagators must terminate in the state |1 where the first particle in particular is excited, at least for this particular calculation of ρ 1;1 (t). Note when interpreting Fig. 3 that a propagator matrix element U b;a describes transition from |a to |b ; i.e., from right to left in its indices, even though time flows from left to right in all diagrams.
The diagrammatic rules outlined above generalize to more complicated diagrams describing the quantity ρ g;g (t). One such (randomly chosen) diagram is shown FIG. 3. Diagrams corresponding to the quantity ρ1;1 (t) for a two-particle relativistic SR sample starting in the doubly excited state |e . Note that time flows from left to right in each diagram, but that a propagator conventionally denotes transition from the state of its right index to the state of its left index.
in Fig. 4. Propagators within the second (central) region of Fig. 4 are no longer constrained to terminate in any particular state, which increases the number of total possible diagrams contributing to ρ g;g (t) to 16. Actually, there is a great deal of symmetry to the set of differential Eqs. (47)-(50) and their initial conditions. From inspection, we see that U 1;2 = U 2;1 and that U 1;1 = U 2;2 . Let us refer to the result of multiplying a top propagator with a bottom propagator within a region between vertical photons as a block propagator. In light of these symmetries and the constraints placed upon diagrams outlined above, the only possible block propagators for our system are B e = |U e;e | 2 (56) as well as B * b . We are finally prepared to assess velocity coherence in a relativistic two-particle SR sample possessing some mean β and observer frame ∆v, by the following steps: 1. Solve the set of differential Eqs. (47)-(50) for the desired β and ∆v, 2. Compute the diagram for each of ρ 1;1 (t), ρ 2;2 (t) and ρ g;g (t) as follows: 4. In order to interpret a total emission power from the particle excitation time dependence, compute the time derivative of the total particle energy in the observer frame; i.e., (2ρ e;e +ρ 1;1 +ρ 2;2 ) ω ′ 0 /γ.
In terms of our block propagators and symmetries, we have that ρ e;e (t) = B e (t) (59) The above results are differentiated in step 4 above, which allows us to skip one convolution from any of the above calculations; for example,

V. RELATIVISTIC TWO-PARTICLE SR: RESULTS
The first result of our theory is a confirmation of the time dilation of two-particle SR in the ∆v → 0 limit. In this limit the relativistic coherence kernel goes to 1, so that all factors of γ in Eqs. by transforming to the variable τ = t/γ. Whatever the result of performing all steps required to compute a total emission power transient, then, that result will be a function of t/γ and thus time-contracted in the observer frame. This is not so much a useful result as it is a reassurance that our model behaves as expected.
The quantity of interest to this work is the behaviour of SR coherent enhancement as a function of the velocity difference of highly relativistic particles. Let us work in the variable q = Γ ′ 0 t/ (2γ) such that the only non symmetry-redundant quantities of Eqs.
We have also at the outset thatŨ e;e (q) = e −2q . A Python script was written to compute C SR, rel β,∆v , to solve Eqs. (64) and (65), to compute the block propagators, to convolve the appropriate diagrams for the density operator traces, and finally to obtain the total power transient for various particle velocity differences and relativistic β values.
In the simplest case, we first verify that our model demonstrates a transition of SR enhancement between the far separated ∆v = ∞ limit and the coherent ∆v = 0 limit. The photon emission rate at β = 0 is shown in Fig. 5 for the coherent case ∆v = 0 and for the well separated case ∆v = 100cΓ ′ 0 /ω ′ 0 . The former demonstrates SR enhancement, while the latter recovers the spontaneous emission rate from two independent particles. Fig. 5 does not differ substantially between the coherent SR (∆v = 0) and the independent spontaneous emission (∆v → ∞) cases. In order to obtain a more dramatic SR pulse demonstrating enhanced and delayed peak intensity one would need to simulate a larger number of particles. Such an exercise becomes very complicated very quickly, as the number of diagrams requiring evaluation grows extremely fast with the number of particles. Combinatoric considerations or computer methods could speed the process up and make a slightly larger population manageable, but we leave such exercises to future work. Instead, we focus now on extracting meaningful information from the two-particle case.
For the purpose of quantifying the departure of an SR system from the independent limit as a function of observer frame velocity difference, let I ∆v (t) = 2ρ e;e +ρ 1;1 +ρ 2;2 (66) for a sample having particle velocity difference ∆v, and define the relativistic velocity coherence metric which measures total departure of a given power transient from that transient which would be generated by independently emitting particles. The 1/A out front is a normalisation factor such that G (0) = 1. Fig. 6 shows G (∆v) for three values of β. The full-width half maxima (FWHM) of the various G (∆v) provide us with a measure of the velocity coherence required for two relativistic particles to demonstrate enhanced SR emission. For β = 0 the FWHM is approximately 9.1cΓ ′ 0 /ω ′ 0 , while for β = 0.95 the FWHM is approximately 0.52cΓ ′ 0 /ω ′ 0 . There are numerous ways to interpret the various FWHM. In the simplest approach, let us compare their values in the context of the relativistic transformation of velocity differentials. The usual relativistic velocity composition rules state that small velocity differences transform as ∆v ′ ≈ γ 2 ∆v between the observer and rest frames. The β = 0.95 plot of Fig. 6 corresponds to γ 2 ≈ 10.3; however, that plot's FWHM differed from the β = 0 case by a factor of 17.5. This suggests that the velocity coherence required between two particles to demonstrate enhanced SR emission in the observer frame is more restrictive than accounted for by simple considerations of the relativistic transformation of velocity spreads.
The situation is further nuanced if we ask which velocity scale should be assigned meaning when quantifying velocity coherence requirements. We have thus far compared velocity differences in terms of the relativistically invariant unit cΓ ′ 0 /ω ′ 0 ; however, the researcher may be more interested in absolute velocity differences in a selected reference frame. Whether velocity coherence requirements become more or less restrictive from a relativistic source is thus not well-defined in a general sense, and will depend upon what other interactions are of interest to the researcher. Moreover, many other relativistic effects may enter into the picture. We discuss some of these effects in Section VI.

VI. SUMMARY AND FUTURE WORK
We constructed a model of velocity dependent relativistic two-particle SR based upon an existing model [16] of relativistic spontaneous emission from a two-level particle. We applied the diagrammatic method of quantum optics to our relativistic SR Hamiltonian and found that certain standard results and techniques of the literature required revision. The introduction of velocity dependence caused a loss of convolution structure in zerophoton propagator calculations. We introduced a set of coupled differential equations, as well as the so-called relativistic velocity coherence kernel, which together enabled the calculation of SR propagators. The standard vertical photon result of the literature was re-proven and was shown to introduce an overall Lorentz factor into the density operator diagrams of SR.
We solved our relativistic SR model in Section V for two particles and characterized velocity coherence requirements for SR emission enhancement from relativistic sources at various β factors. Initial results suggested that velocity coherence requirements for enhanced SR emission in the observer frame become more restrictive as a source moves with higher relativistic velocity. Some caveats were made in Section V regarding the interpretation of velocity coherence requirements, which in reality becomes a nuanced exercise.
The model constructed in this paper forms a foundation for relativistic SR work. Future research could proceed in numerous directions. First, it remains to yet solve the present model for (i) dipoles aligned parallel to the particles' velocities, (ii) dipoles not oriented parallel to each other, and (iii) non co-linear particle velocities. Second, the diagrammatic method may also be used to calculate the spectrum and angular photon distribution from an SR source. This procedure could be applied to the present model and radiation intensity could be characterized across different directions from a relativistic SR source. Some initial calculation steps for this purpose are provided in Appendix D. Third, this model could be expanded to a larger number of particles. Although the diagrammatic method becomes quickly cumbersome with increasing population, computer methods exploiting combinatoric considerations could certainly simulate a moderately sized sample. Despite remaining physically negligible, a sample of four or more particles begins to demonstrate the salient SR features of delayed and enhanced peak intensity. A relativistic model of these features would be especially helpful to astrophysicists. Fourth, this model could be generalized to an extended sample by paralleling existing derivations of the Maxwell-Bloch equations, now with our relativistic SR Hamiltonian. Considerations of relativistic length contraction or time dilation of relaxation processes in such an extended model could lead to richer relativistic SR behavior. the charge separation in a frame co-moving with the particle's center of mass. Boussiakou et al. [16] select the Power-Zienau-Woolley (PZW) gauge for the electromagnetic field, which permits them to express the Lagrangian in the particle rest frame in a manifestly covariant form by contracting the field tensor with itself (the usual free field contribution to the Lagrangian) and contracting the field tensor with the polarization-magnetization tensor (a covariant description of the matter-field interaction). The polarization and magnetization fields comprising the latter tensor in the PZW gauge are constructed from the charge separation coordinate and its velocity via a prescription provided by existing PZW gauge theory in the literature [19].
The covariant Lagrangian formulation enables Boussiakou et al. to determine, via careful Lorentz transformations, an expression for the observer frame Hamiltonian in terms of the observer frame center of mass coordinates, charge separation coordinates, and radiation field; these observer frame quantities are coupled to the particle's rest frame polarization and magnetization properties. This form of the Hamiltonian ultimately enables the BBB model to relate observer frame behaviour to transition operators intrinsic to the particle. The BBB Hamiltonian in the observer frame reads, where a few outstanding terms have not yet been defined within the body of the paper. The position operatorsq andq ⊥ describe the charge separation coordinates (unprimed in the observer frame, primed in the rest frame) parallel and perpendicular to the center of mass velocity, the operatorsp andp ⊥ are their conjugate momenta operators, and the center of mass has an observer frame momentum operatorP. The Lorentz factor γ corresponds to the initial center of mass momentum of the particle. The operatorq ′ is the charge separation distance operator in the rest frame coordinates, µ is the rest frame reduced mass of the particle, andP ′ is the particle's rest frame polarization operator. For simplicity we have suppressed magnetization terms, consistent with the paper's restriction to the case of an electric dipole which couples to the observer frame electric (Ê) and magnetic (B) fields. The second line of Eq. (A3) is the so-called Röntgen interaction term.

Eigenfunctions of the unperturbed BBB Hamiltonian
We require the unperturbed eigenfunctions of the atomic Hamiltonian term defined by Eq. (A2). We present here the main results of the solution provided in [16]; the eigenfunction work below is adapted directly from their work.
Projecting onto the observer frame center of mass and charge separation coordinate basis and factoring the wavefunction as Ψ (Q, q) = e iK·Q ψ (q), one obtains (A4) for ǫ the (observer frame) energy of the state. Eq. (A4) may be expressed in rest frame coordinates as which is the usual hydrogenic Schrödinger equation, but where the energy eigenvalues of the conventional solutions correspond now to γ times the energy eigenvalues in the observer frame. Additionally, reversion back to the observer frame will introduce a length contraction to the conventional hydrogenic wavefunction shapes. Both of these results [16] are consistent with familiar results from the special theory of relativity [20].
It is important to keep in mind that the primed coordinates are simply a mathematical convenience for working with a wavefunction which is physically relevant to the observer frame. Although the unprimed and primed coordinates are related by the familiar Lorentz transformation between the observer and rest frame position coordinates, the primed coordinates serve only as a tool for solving Eq. (A4). The physics thus far developed has been based upon canonical quantization in the observer frame and meaningful conclusions from the present wavefunction can only be made with respect to that frame.
We now make an important point regarding normalization which affects our diagrammatic calculations in this work. The energy eigenfunction relation of Eq. (A4) becomes hydrogenic in the primed coordinates of (A5); we choose therefore to use the conventional hydrogenic wavefunctions in the primed coordinates, or their transformed version in the unprimed coordinates. Such functions are no longer normalized in the unprimed (observer) frame due to the factor of γ relating the coordinates parallel to the particle's motion. We might alternatively re-normalize the eigenfunctions in the observer frame; however, we will later be interested in expanding operators in basis states indexed by the primed (rest frame) coordinates. We therefore retain the conventional normalization of the hydrogenic wavefunctions relative to the primed coordinates, but remembering always that the completeness must then be asserted relative to the primed basis states; i.e., (A6)

The BBB Hamiltonian in its form useful to this work
In this section we summarize the work of [16] to express the BBB Hamiltonian in terms of operators and properties intrinsic to the particle, including internal raising and lowering operatorsπ † andπ, the internal dipole moment d ′ , and the rest frame transition energy ω ′ 0 . We also introduce a simplified form of the interaction Hamiltonian useful to this work.
The operatorP ′ appearing in the interaction operator V of Eq. (A3) refers to the polarization vector operator in the rest frame, and it couples to field operators in the observer frame. Boussiakou et al. [16] assert the dipole approximation to express the rest frame polarization operator asP ′ (r ′ ) =d ′ δ (r ′ − Q ′ ), whered ′ = eq ′ . With this approximation, we perform first the integration over d 3 r of Eq. (A3) in the observer frame; however, we must recognize that the arguments of the delta distribution are in the rest frame coordinates, and thus its support is length-contracted in the observer frame; i.e., the integration picks up a factor of 1/γ and we are left witĥ Let us prepare ourselves for the theory of dipole emission between two energy eigenstates |E 1,2 of the particle. We expand the dipole operator aŝ (A8) [14] based upon the original papers [21][22][23]. We work with the state of the system |Ψ (t) in the interaction picture defined as The state of the system at time t is related to its original state at time t 0 via the interaction picture time evolution operator Ũ (t, t 0 ); i.e., The Schrödinger equation forŨ (t, t 0 ) in the interaction picture reads We require in this work the exact first-order recurrence relation as well as the exact second-order recurrence relation We also make use of the infinite series expansion A matrix elementŨ nm = E n | Ũ |E m coupling two free theory energy eigenstates |E n and |E m may be interpreted by momentarily leaving the interaction picture.
The composition of adjacent occurrences of Ṽ within each integrand translates to the (k, l) th non-interaction picture matrix element where summation over p is implied and where ω q is the energy of a free theory energy eigenstate |E q . The term is the free theory propagator of the system in state |E p from time t = t n to time t = t n+1 . The sequence V kp F p (t n+1 , t n ) V pl within expression (B8) is interpreted from right to left as the physical process wherein the system (i) transitions from state |E l to state |E p at time t n , (ii) propagates freely in state |E p from time t n to time t n+1 , and (iii) transitions from state |E p to state |E k at time t n+1 . Each term in Eq. (B7) may be interpreted in this manner, where the order of a term in the expansion corresponds to the number of transition events. For each event order sequence, we are directed to integrate over every possible timing of the transition events which satisfies the ordering t > t n > t n−1 > · · · > t 0 . The summation over transition matrix elements during matrix multiplication of adjacent V kl terms directs us to sum over all possible intermediate transition states.
The particular form of the interaction of interest in this work, given by Eq. (A14), yields transitions from a state |e, 0 (excited particle with no photon in the radiation field) only to states of the form |g, k (ground state particle with one photon of mode k in the radiation field); and transitions from a state of the form |g, k only to the state |e, 0 . This interaction therefore greatly restricts the collection of possible transition event sequences; for example, a transition event of the form |e, k → |e, k ′ for k = k ′ is strictly forbidden [24].
Terms in the series expansion of a time evolution operator matrix element are conventionally brought into correspondence with diagrams depicting sequences of transition events. In light of our previous comments, all possible transitions of a single particle system evolving with the interaction defined by Eq. (A14) may be found in Fig. 7, which shows the first two non-vanishing terms in the diagrammatic representation of the matrix element corresponding to an initially excited particle being found in the excited state after time t − t 0 . dt 2 t2 t0 dt 1 k e i(ω ′ 0 /γ−α k1 ω k) t2 V 1,0;g,k (t 2 ) e −i(ω ′ 0 /γ−α k2 ω k) t1 V g,k;2 (t 1 )Ũ 2;1 (t 1 , t 0 ) .
by our prescription of Eq. (15). Let us assume that both dipoles are aligned with the x axis, that both velocities are in the z direction, and that the mode k has length k and direction Ω (with polar and azimuthal parts). Suppressing for a moment the t 1 , t 2 integrals andŨ 2;1 , the sum over k (including the leading [−i/ ] 2 multiplier) becomes whereω k = ω k α k with α k = 1 − β cos (θ) for θ the angle between the mode propagation direction and mean velocity v = (v 1 + v 2 ) /2 (or |v| = βc). The velocity difference is ∆v = v 2 − v 1 . We now assume that the integration over k contributes predominantly near where ω ′ 0 /γ −ω k = 0. We can therefore pull the slowly-varying part k 3 e ik∆vx(t1+t2)/2 out of the integral over k and evaluate it at k ′ 0 / (γα k ), make a change of integration variable, and formally extend our integration limits to obtain, Plugging this result back into our time integrals above, and recalling that t 1 only integrates up to t 2 (and thus only picks up half of the support of the Dirac delta distribution), we can immediately execute the integral over t 1 to obtain − Γ ′ which is our desired expression.

Appendix D: Photon Emission Propagator for the BBB Hamiltonian
We provide in this appendix the single particle propagators corresponding to evolution into non-zero photon states in the BBB Hamiltonian. These propagators are not necessary for the relativistic coherence study of the present work, but would be essential to potential future research into the emission spectrum from a relativistic SR system. We require for this purpose the elements V g,k;e,0 = g, k|V |e, 0 Let us consider first the case where V i is parallel to d ′ . Using the ǫ kθ and ǫ kφ polarizations of Section IV B, only the former yields a non-vanishing interaction matrix element V g,kθ;e,0 = ξ * k d ′ sin θ k /γ (D2) (where V denotes that this matrix element is computed for V i d ′ ) or, in the interaction picture, V g,kθ;e,0 = e iω g,kθ t (ξ * k d ′ sin θ k /γ) e −iωe,0t = e −i(ω ′ 0 /γ−α k ω k) t ξ * k d ′ sin θ k /γ (D3) which has Laplace transform V Q g,kθ;e,0 = ξ * k d ′ sin θ k /γ s + i (ω ′ 0 /γ − α k ω k ) .
The above was computed for the case where V i was parallel to d ′ ; if V i is now perpendicular to d ′ , we obtain the propagators (D11) where φ k is the angle between V i and the projection of k onto the plane perpendicular to d ′ .