A standard form of master equations for general non-Markovian jump processes: the Laplace-space embedding framework and asymptotic solution

We present a standard form of master equations (ME) for general one-dimensional non-Markovian (history-dependent) jump processes, complemented by an asymptotic solution derived from an expanded system-size approach. The ME is obtained by developing a general Markovian embedding using a suitable set of auxiliary field variables. This Markovian embedding uses a Laplace-convolution operation applied to the velocity trajectory. We introduce an asymptotic method tailored for this ME standard, generalising the system-size expansion for these jump processes. Under specific stability conditions tied to a single noise source, upon coarse-graining, the Generalized Langevin Equation (GLE) emerges as a universal approximate model for point processes in the weak-coupling limit. This methodology offers a unified analytical toolset for general non-Markovian processes, reinforcing the universal applicability of the GLE founded in microdynamics and the principles of statistical physics.

• Physics: Within statistical physics, particle motion in water is described by the generalized Langevin equation (GLE) [1].The GLE represents a quintessential non-Markovian stochastic model, capturing the hydrodynamic memory effect.
• Econometrics: The autoregressive integrated moving average (ARIMA) model stands as a recognized discretetime non-Markovian model describing many stylized structures of financial returns [3].
Central to these models is their ability to encapsulate long-memory effects inherent to various systems.This is typified by power law decaying autocorrelation functions (ACFs), transcending the conventional boundaries set by Markovian stochastic processes.A well-established analytical toolkit has been developed for Markovian stochastic processes [22,23], which includes stochastic differential equations (SDEs), master equations (MEs), and their asymptotic solutions.For instance, the theory of standard forms has been instrumental in the systematic classification of both SDEs and MEs.Given that MEs represent linear time-evolution equations for probability density functions and functionals (PDFs), they can be solved within the framework of linear algebra, particularly through methods like the eigenfunction expansion [22,23].
There are also various asymptotic methods tailored to MEs.Prominent among these are the system-size expansion [2,[24][25][26] and the Wentzel-Kramers-Brillouin (WKB) approximation [23,27].Notably, the system-size expansion stands as a historic cornerstone in the realm of statistical physics, especially concerning the Langevin equations.This is largely due to its role in extrapolating various Langevin equations from underlying microscopic physical dynamics.Hence, Markovian process theory offers a robust and structured foundation for statistical physics, at least in a formal sense.
In contrast to the structured theories for Markovian processes, those for non-Markovian processes remain more fragmented.A universally accepted master equation (ME) theory for non-Markovian processes is absent.Current MEs pertain specifically to particular non-Markovian SDE classes, such as GLE with exponential memories [28,29], GLE with linear potential [30], and semi-Markovian point processes [31].Without a standardized form for these MEs,

C. Notation for functionals
If the argument of a map f is a function {z(s)} s , f is called a functional.A functional is indicated by the square brackets f [{z(s)} s ].The functional notation f [{z(s)} s ] is sometimes abbreviated as f [z] if its meaning is obvious from the context.For a stochastic field variable {ẑ t (s)} s , the corresponding PDF is written as P t [z] = P t [{z(s)} s ], characterising the probability P t [z]Dz that {ẑ(s)} s ∈ s [z(s), z(s) + dz(s)), where the functional volume element is Dz := s dz(s).The ensemble average of any functional f [ẑ] is written as the path-integral representation On the basis of this notation, the PDF is formally rewritten by P t [z] = ⟨δ[z − ẑt ]⟩, where the δ functional is defined by δ[z − ẑt ] := s∈R + δ(z(s) − ẑt (s)).The concept of derivative can be generalised to the functional derivative, which is denoted by δf [z]/δz(s) (see Appendix A for the detail).

III. LITERATURE REVIEW: MARKOVIAN STOCHASTIC PROCESSES
This section offers a concise overview of the foundational theory of Markovian processes, serving as an introduction for readers less acquainted with Markovian processes and statistical physics.Specifically, we touch upon the standard form of the master equation (ME) for these processes.Experts who are solely focused on our primary findings may bypass this section, as the main results are presented in a standalone, comprehensive format.

A. Markovian stochastic differential equations
Let us consider a one-dimensional stochastic process characterised by the trajectory {v s } s≤t , where t is the current time.If the statistics of the infinitesimal future state vt+dt is completely characterised only by the current state vt , the stochastic dynamics is said to obey a Markovian stochastic process.However, a more general class of stochastic models can be considered that cannot be characterised only by the current state vt .For example, a stochastic model can depend on the full history {v s } s≤t .Such stochastic dynamics obey a non-Markovian stochastic process.This subsection reviews the theory of Markovian stochastic processes, in particular their standard forms and the system-size expansion.

Standard form of the white noise
White noise is a noise that is independent of its history.It is formally defined as the derivative of the Lévy process Lt , such that ξW t := d Lt / dt.According to the Lévy-Itô decomposition, any white noise can be decomposed as the sum of the white Gaussian noise and of the white Poisson noise (see Appendix B for a review), such that where m is the constant drift, σ is the standard deviation of the white noise, ξG is the white noise, and ξP λ(y) is the white Poisson noise with intensity distribution function λ(y) of the jump size y.

Standard form of stochastic differential equations
Any one-dimensional stochastic process can be constructed from white noise by introducing the state-dependence into the drift term a, standard deviation σ, and the intensity distribution function (IDF) λ(y), such that A one-dimensional stochastic Markovian process {v s } s≤t obeys the state-dependent SDE with the Itô interpretation assumed.We refer to this representation as the standard form of one-dimensional Markovian SDEs.The Markovian property is expressed by the fact that the right-hand side of Eq. ( 6) depends only on the current state vt .
FIG. 1: Schematic of the Markovian jump process (9).The path is piecewisely continuous and occasionally has jumps.A jump occurs during [t, t + dt) with probability λ(y|vt)dtdy with jump size ŷ ∈ [y, y + dy).The remarkable character of the Markovian jump process is that the intensity density λ(y|vt) depends only on the current state vt and does not depend on the whole history {vτ }τ<t.In this sense, the Markovian jump process is an history-independent Poisson process, in contrast to the non-Markovian jump process (or the history-dependent Poisson process) defined as Eq. ( 17).

Standard form of master equations
The SDE (6) describes the dynamics of stochastic systems for a single path.While the SDE are intuitive tools, they are not easy to handle because of their general nonlinear structure.For analytical calculations, the ME approach provides more systematic methods based on linear algebra.The ME is the equation governing the time-evolution of the PDF P t (v) := ⟨δ(v − vt )⟩ as follows with a linear operator L. The ME corresponding to the SDE ( 6) is given by This ME is known to covers all possible one-dimensional Markovian stochastic processes [22], and, thus, is called the standard form of the master equation in this report.The ME is very useful because it is always a linear dynamical equation 1 of the PDF P t (v).In other words, a standard approach to a given Markovian process is to consider its ME and solve the corresponding eigenvalue problem with linear algebra techniques.

Markovian jump process (history-independent Poisson process)
Markovian jump processes constitute a large subclass of Markovian SDEs, such that where the drift term and the white Gaussian noise term are absent, and only the jump term is present.Markovian jump processes depend only on the current state vt and can be called history-independent Poisson processes, in contrast to the non-Markovian jump processes (or history-dependent Poisson processes) defined in Eq. ( 17) in the main section.Markovian jump processes are popular models.For instance, the detailed description of physical Brownian motion is often modelled as a Markovian jump process for the velocity of the Brownian particle, where the velocity discontinuously changes due to molecular collisions.
1 Indeed, the formal solution of the ME ( 7) is given by with initial condition constants {c i } i , where µ i and ϕ i (v) are the ith eigenvalue and the corresponding eigenfunction, respectively.Jump size: FIG. 2: Schematic of the system-size expansion.(a) A typical trajectory of a Markovian jump process (10) where the jump size is scaled with ε, such that ŷi = ε Ŷi with ε-independent jump size Ŷi.(b) For small ε ≪ 1, the path becomes approximately continuous due to the small jump sizes and obeys the Langevin equation ( 13) under an appropriate stability condition around vt ≃ 0. This picture essentially applies even to the non-Markovian jump process (17) as shown in Sec.V.

B. The system-size expansion
Solving the eigenvalue problem is in general difficult, in particular when the linear operator L leads to an integrodifferential equation of the form (8). One of the systematic methods to obtain asymptotic solutions was invented by van Kampen, which is called the system-size expansion.Mathematically, the system-size expansion can be regarded as a weak-noise asymptotic limit for Markovian jump processes.This assumption is very natural particularly in the context of Brownian motions.This method provides a solid mathematical derivation of the Langevin equations from microscopic physical dynamics.In this subsection, we briefly review this methodology based on Refs.[24][25][26].

Sketch of the system-size expansion
Let us consider a Markovian jump process described by where λ(y|v t ) is the conditional intensity density of the jump size y, N (t) is the total number of jumps during [0, t), and tk is the k-th jump time, ŷk is the k-th jump size.Let us assume that the jump size is proportional to a small positive parameter ε > 0 (see Fig. 2(a)), so that we can write ŷk = ε Ŷk .(11) This implies that the noise term can be rewritten as ξP λ(y|vt) = ε ξP W (Y |vt) with the conditional intensity density W (Y |v t ) for the rescaled jump size Y .We thus obtain the SDE with a small jump-noise term: This is the scaling assumption for the system-size expansion.In other words, the small parameter ϵ can be interpreted as the small constant quantifying the weak-coupling with the stochastic environment.
In the small-noise asymptotic limit ε → 0 and for a broad variety of setups, assuming a stability condition around vt ≃ 0 (see Appendix C for details), the Markovian jump process reduces to the Langevin equation (see Fig. 2(b) for a schematic) where γ takes the meaning of a frictional constant and T is the temperature.See Appendix C for the detailed derivation.Thus, the system-size expansion is a celebrated mathematical foundation for the derivations of the Langevin equations from microscopic dynamics.

Collision
Velocity jump size FIG.3: The scaling assumption of the system-size expansion is natural on physical grounds for the description of a massive Brownian-particle motion.Let us denote the mass and initial velocity of the Brownian particle (respectively of a surrounding small particle) by M (respectively m) and V (respectively v).After an elastic collision, the velocity of the Brownian particle changes according to formula (14) deriving from the conservation of momentum.The velocity jump size y is proportional to the mass ratio ε := m/M , which is the small parameter of the problem in the limit of massive Brownian particle limit m ≪ M .
2. Physical validity of the scaling assumption.
The physical validity of the scaling ( 11) and ( 12) can be intuitively understood by considering a one-dimensional collision problem (Fig. 3).Let us prepare a small particle of mass m, velocity v and a large particle of mass M and velocity V .In a one-dimensional elastic collision, the post-collisional velocity V ′ of the large particle is given by Since the typical thermal velocities are given by v ≃ T /m and V ≃ T /M where T is the gas temperature, we have |V /v| ∝ ϵ 1/2 ≪ 1 with ϵ := m/M .For ϵ ≪ 1, we obtain the velocity jump y of the large particle as Thus, the velocity-jump size is proportional to ε, and satisfies the system-size expansion scaling (12) exactly.This example highlights that the scaling assumption of the system-size expansion is physically reasonable 2 when the Brownian particle in a gas is much heavier than the surrounding gas particles.
3. Scaling assumption in the master equations.
The scaling assumption y = εY at a trajectory level is equivalent to the scaling assumption for the ME: which is derived from the conservation of probability (i.e., the Jacobian relation), such that λ(y|v)dy = W (Y |v)dY .

C. Goal of this report
On the basis of the above theory regarding the standard forms of SDEs and MEs, our goals in this report are the following.
1. We derive the ME for the general non-Markovian jump process analogous to the standard form of the Markovian ME (8).
2. We asymptotically solve the ME for non-Markovian jump processes by generalising the system-size expansion; we finally obtain the GLE via a physically-reasonable coarse-graining approach.

Jump probability
History FIG.4: Schematic of the non-Markovian jump process (or history-dependent Poisson process).The probability that a velocity jump occurs during [t, t+dt) is given by λ(y|{vτ } τ ≤t )dtdy with jump size ŷ ∈ [y, y +dy).Here the intensity density λ(y|{vτ } τ ≤t ) explicitly depends on the whole history {vτ } τ ≤t and the system is thus truly non-Markovian.

IV. NON-MARKOVIAN MODEL AND FORMULATION
In this section, we first present the stochastic model studied in this report.We then introduce the Laplaceconvolution Markovian embedding that converts the original low-dimensional non-Markovian dynamics onto a Markovian field dynamics.Finally, the corresponding field ME is formulated.

A. Non-Markovian jump process (history-dependent Poisson process)
Let us consider a non-Markovian stochastic model that can encompass a large class of non-Markovian stochastic processes: We study the history-dependent compound Poisson process (see Fig. 4 for a schematic) where the intensity λ(y|{v τ } τ ≤t ) with jump size y is conditional on the full history of the system {v τ } τ ≤t .
The non-Markovian nature of this process makes the intensity a functional of the whole history.More technically, Eq. ( 17) implies that for any given history {v τ } τ ≤t with infinitesimal time evolution dv t := vt+dt −v t .Our aim is to provide the full analytical toolset for this history-dependent Poisson process by developing the corresponding field ME and by analysing its asymptotic solutions.

B. Markovian embedding
In this subsection, we apply the Markovian-embedding scheme to the history-dependent Poisson process (17).We finally obtain the SPDE govering the Markovian field dynamics and derive the corresponding field ME.

Basic idea
The idea of Markovian embedding is very simple: a low-dimensional non-Markovian dynamics can be converted onto a higher-dimensional Markovian dynamics by adding a sufficient number of auxiliary variables.This approach Markovian embedding FIG.5: Schematic of the Markovian embedding of the original one-dimensional non-Markovian jump process (17) onto the Markovian field dynamics {ẑt(s)}s.The auxiliary field variables {ẑt(s)}s are defined by Eq. ( 19) on the wave-number axis s ∈ (0, ∞) (i.e., one-dimensional field) and obey the first-order Markovian SPDE (24a).
dates back to Mori [36] around the mid 1960s3 .Also, the theory of the Kac-Zwanzig model [28,29,[37][38][39] can be regarded a theory of Markovian embedding between the generalised Langevin equation and the Hamiltonian-particles model with harmonic interaction.For example, the generalised Langevin equation with the sum of K-exponential memories can be thought of as a K-dimensional Markovian dynamics [28,29,33].This idea can be even applied to the Hawkes processes [4,40,41] for memory kernels expressed as a sum of exponential functions.Remarkably, this idea of Markovian embedding has been also applied to non-Markovian stochastic processes in quantum systems [42][43][44][45][46][47][48] in the context of the pseudomode approach around the mid 1990s.
The dimension needed for the Markovian embedding depends on the model but can be infinite in general.In this case, the dynamics can be regarded as a Markovian field dynamics.For instance, the GLE and the Hawkes processes have been converted onto Markovian field dynamics [32][33][34][35], which can be analysed by the field ME (which is a functional-differential equation for the probability density functional).
Markovian embedding is nontrivial and technically tricky for continuous-time stochastic processes, while Markovian embedding is rather straightforward for discrete-time stochastic processes (see Appendix D for brief clarification).This report aims at formulating a general embedding theory of the non-Markovian jump process (17), even though it is based on continuous time.

Variable set
Before proceeding with the derivation of the Markovian field dynamics, let us introduce a complete set of system variables useful for Markovian embedding.In the previous section, we used {v τ } τ <t as a naive complete set of system variables.This set is equivalent in information content to another set (v t ; {v τ } τ ≤t ) with acceleration ât := dv t / dt.Note that the acceleration can include the impulses described by the Dirac δ functions associated with the jumps.Let us now introduce the Laplace-convolution Markovian-embedding representation of the velocity trajectory as which is defined for s > 0 (see Fig. 5 for a schematic).We then adopt the variable set as a useful complete variable set.
The introduction of the auxiliary field variable {ẑ t (s)} s>0 is the technical but crucial trick to convert the general nonlinear non-Markovian model onto a Markovian field model.In the following, the wave number s is always considered strictly positive (s > 0), and the set {f (s)} s>0 for any function f (s) is sometimes abbreviated by {f (s)} s if its meaning is clear from the context.

Phase space
Let us introduce the state variables Γt as points in the phase space S, such that where R is the space of real numbers and F is the function space.In the following, we simplify the notation of functionals such as the intensity as a functional in terms of the history in the following way: In our notation, the functional argument (e.g., Γt and {ẑ t (s)} s ) follows other variables (e.g., y and vt ) after the separation by the semi-colon.Note that the ordinary Markov compound Poisson process corresponds to the case where the intensity does not depend on the historical velocities, such that with a nonnegative function λ(y|v t ).

C. Markovian field dynamics
The history-dependent compound Poisson process (17) characterising the original variables {v τ } τ ≤t is equivalent to the set of the following SDE and SPDEs characterising the new variables Γt := (v t , {ẑ t (s)} s ): where the jump term ξCP simultaneously acts on both vt and ẑt (s) for all s > 0 (see Fig. 5 for a typical configuration of the auxiliary field).The initial condition is given by This set of SDEs characterizes the complete dynamics of the phase point Γt = (v t , {ẑ t (s)} s ) in a closed form.
The essence of the trick is to use the Laplace-convolution transform, which encodes the whole history of vt (or equivalently its acceleration ât ) into a function of t (now, thus Markovian) and of an additional variable s.The function ẑt (s) dependent on s serves as the key device to render the system Markovian, utilizing an infinite series of equations for all ẑt (s).It is remarkable that this system is Markovian in the extended phase space Γt ∈ S, while the original one-dimensional process is non-Markovian.This means that we have successfully transformed the original non-Markovian dynamics into a Markovian dynamics by adding a sufficient number of variables.Since the resulting dynamics is Markovian, we can derive the corresponding ME for the PDF for the phase point Γt in the extended phase space.

Derivation
By directly solving Eq. (24a), we obtain with the dummy-variable transformation τ ′ := t + τ and t ′′ := t − t ′ from the second to the third line.Thus, the set of the SDEs (24a) is consistent with the Markovian-embedding representation (19).

D. Field master equation for the history-dependent Poisson process
The functional ME of the field corresponding to the SPDEs (24a) are given by with jump size vector with indicator function 1(s) = 1 for any s.

Derivation
We derive the field ME as follows.For any functional f at leading order 4 .By taking the ensemble average of both sides, up to the order of dt, we obtain with integral volume element dΓ := dvDz.By applying a variable transformation Γ + ∆Γ y = Γ ′ , the first term in the right-hand side is given by where the dummy variable Γ ′ is finally replaced with Γ.By applying the functional partial integration (A17), the third term in the right-hand side of Eq. ( 29) is given by By considering that the left-hand side of Eq. ( 29) is given by we finally obtain the following integral identity regarding any functional in the limit dt → 0. Since this relation holds for an arbitrary f [Γ], we obtain the field ME (26).

E. Functional Kramers-Moyal expansion
By applying the identity (see Eq. (A11) for the functional Taylor expansion) to the field ME (26), we obtain the functional Kramers-Moyal (KM) expansion: with the KM coefficient F. Remark: systematic calculations based on linear algebra Since the field ME ( 26) is linear, it can be analysed with the tools of linear algebra.Starting from the standard form (7 where µ i is the i-th eigenvalue and ϕ i [Γ] is the corresponding eigenfunction 5 .The time-dependent solution is then given by superposition of the eigenfunctions where the coefficients {c i } i are determined by the initial condition.The steady-state PDF corresponds to the zeroth eigenfunction ϕ i [Γ] with µ 0 = 0: Various physical quantities can be systematically calculated by the time-dependent (37) or steady-state solutions (38).For example, the correlation function is formally given by the path integral where A t and B t ′ are expressed as functionals of Γ.
V. APPLICATION 1: SYSTEM-SIZE EXPANSION AND GENERALIZED LANGEVIN EQUATION In this section, we illustrate the utilization of our field ME framework in relation to an asymptotic theory, drawing upon the system-size expansion.By enforcing a stability condition upon the system-size expansion, we ultimately infer the Generalized Langevin Equation (GLE) as a plausible coarse-graining process rooted in physical reasoning.
A. Assumptions 1. Small noise assumption: Let us consider the non-Markovian process with small jumps: where W is the ε-independent conditional intensity of the jump size Y (see Fig. 2 for the schematic of the small-jump assumption).This assumption is equivalent to the following scaling relation of the conditional intensity density This assumption leads to the scaling relation for the KM coefficients with the ε-independent 2. Linear stability: Let us additionally assume that the first-order KM coefficient A 1 [Γ] has a single stable point: We also assume that A 1 [Γ] is linearly stable around Γ = 0, such that where the rescaled wave number u is defined by u := s ε .3. Existence of the noise term: The noise term is assumed to be present even for ε → 0 and thus the variance term is nonzero: As a technical assumption, we assume that all the considered integrals converge.This assumption implies that physically singular processes, such as with long time tail with decaying speed slower than t −1 , are out of the scope of this paper.Note that the above stability assumptions parallel the conventional stability assumption for the Markovian jump process (see Appendix C for their comparison).Also, we note that all even-order KM coefficients are positive under assumption 3, such that A k [0] > 0 for all even k due to Pawula's theorem [23].

B. Asymptotic derivation of the functional Fokker-Planck equation
Given the following rescaled variables we obtain the following functional Fokker-Planck equation for ε → 0

Derivation
With the above assumptions, let us formulate the system size expansion for this model.We have and We also consider the functional Maclaurin series (A12) for the KM coefficient around Γ = 0 (50) with the dummy-variable arguments χ and ζ.This relation implies that and From the KM expansion (35a), by introducing t := εt, we obtain γV We thus obtain the functional Fokker-Planck equation (47) in the weak coupling limit ε → 0.

Equivalent stochastic dynamics
Furthermore, the functional Fokker-Planck ( 47) is equivalent to the stochastic dynamics described by with standard white Gaussian noise ξG t that is common to the stochastic dynamics of Vt and Ẑt (u).

C. Asymptotic derivation of the generalized Langevin equation
The stochastic dynamics ( 54) is equivalent to the generalized Langevin equation (GLE): with memory kernel M(t) given by expression (60) and colored Gaussian noise This result implies that the GLE is a minimal model for the coarse-grained description of general non-Markovian jump processes in the weak coupling limit ε → 0 under the stability condition (V A).
The memory kernel M(t) and the noise statistics can be explicitly derived as follows.Let us define the "matrix" K(u, u ′ ) and the corresponding eigenvalue µ and eigenfunctions {e(µ; u)} µ satisfying The matrix K(u, u ′ ) has the following properties (see Appendix E): (i) All of its eigenvalues are real and positive µ > 0. (ii) K(u, u ′ ) is a positive symmetric matrix and thus has an inverse matrix K −1 (u, u ′ ).We assume that the eigenfunctions {e(µ; u)} µ constitute a complete set, and have inverse matrices e −1 (u; µ ′ ) such that 6 With these notations, the memory kernel and the noise statistical properties are respectively given by 6 In N -dimensional linear algebra, the set of all eigenvectors {e i } i with e i = (e i1 , . . ., e iN ) of a symmetric matrix constitute a complete set.In addition, the matrix A := (e ij ) has the inverse matrix A −1 = (e −1 ij ), such that j e ij e −1 jk = δ ik and j e −1 ij e jk = δ ik .This property is a straightforward generalization from finite-dimensional to infinite-dimensional linear algebra.
Let us rewrite expression (54b) as By introducing the representation based on the eigenvectors {e(µ; we formally obtain the explicit representation of Ẑt (u), where we used Eq. ( 59).We thus obtain whose solution is given by which leads to the explicit form of Ẑt as from Eq. (64).From Eqs. (54a) and (67), we obtain This equation can be written as with the memory kernel and the colored Gaussian noise FIG. 6: Schematic paths of the intensity λt and the price vt described by the nonlinear Hawkes model (72) for the financial price dynamics.

VI. APPLICATION 2: PRICE DYNAMICS BASED ON NONLINEAR HAWKES PROCESSES
In this section, we illustrate another application of our formalism.We focus on modelling financial price dynamics based on a nonlinear Hawkes process.Linear Hawkes processes have become popular in econophysics as well as in econometrics of market microstructure.

A. Model
Let us consider a stochastic financial model based on the nonlinear Hawkes processes, which has recently become popular to describe the price dynamics of financial assets [49,50].Let us denote vt the logarithm of the price of some stock at time t.The price dynamics is given by where vt is kth jump size of the log price occurring at time tk .The amplitude of the jumps are independently and identically distributed with mark distribution ρ(y).The sequence of jumps defines the jump size series {ŷ k } k and the jump time series { tk }.We denote by Nt the total number of jumps during [0, t).We assume that both excitatory and inhibitory effects are balanced, which is realised when the mark distribution is symmetric: The intensity λt of the jumps is assumed to obey the nonlinear Hawkes process with non-negative intensity function g > 0 and memory kernel h(t).Recall that the intensity λt gives the probability per unit time for the next jump to occur: λt dt is the probability for the next jump to occur during [t, t + dt).See Fig. 6 for the schematic paths of this model regarding the intensity λt and the price vt .This model is an example of a history-dependent Poisson processes.Indeed, the following specific history-dependent Poisson process is equivalent to the nonlinear Hawkes price model (72) [34,35].

B. Markovian embedding
Our Laplace-convolution Markovian-embedding scheme (24) fully converts the nonlinear non-Markovian Hawkes process (72) into a Markovian field process.Indeed, by decomposing the memory kernel as the sum of exponentials the conditional intensity can be rewritten as This is equivalent to the Markov-embedding representation introduced in our previous works [32][33][34][35].

C. Field master equation
The field ME for the nonlinear Hawkes price model ( 72) is with By integrating out both sides over vt , this field ME reduces to the field ME for a marginal PDF P t [z] := P t [Γ] dx that was introduced in our previous works [32][33][34][35] (see Appendix F for the explicit derivation).In addition, the reduced field ME has been analytically solved in Refs.[32][33][34][35] for the asymptotic intensity PDF in the steady state.

D. Diffusive approximation
Let us apply the diffusive approximation by using the KM series (35a) for the the field ME (76) and truncating it at the second order.This leads to the following approximate Fokker-Planck equation with This field Fokker-Planck equation is equivalent to This recovers the standard Geometric Brownman Motion model of price dynamics for constant λt .For non constant λt , equation (79) recovers the general class of stochastic volatility models [51].Here, we derived that the volatility is proportional to the intensity λt of the underlying point process.In other words, our nonlinear Hawkes (72) combined with our Markovian embedding and the diffusive approximation provide an interpretation of the source of stochastic volatility, which is here interpreted as resulting from the underlying jump intensity and its nonlinear memory structure.

VII. DISCUSSION
This section delves into the ramifications of our research and outlines our perspective on several outstanding technical challenges yet to be addressed.

A. Comparison with the projection-operator formalism
Our formulation bears similarities to the projection-operator formalism, as both theories pertain to the derivations of the Generalized Langevin Equations (GLEs).In this subsection, we juxtapose the two approaches, evaluating their respective advantages and disadvantages.
The projection-operator formalism originated in the 1950s and 1960s, crafted by pioneers like Nakajima, Mori, Zwanzig, and Kawasaki [28,[52][53][54][55][56].Particularly, Mori's approach focuses on establishing a microscopic foundation for the Generalized Langevin Equations (GLEs).In the projection-operator formalism, the selection of several slow variables is necessitated, guided by physical intuitions or empirical findings, as these variables cannot be determined theoretically.Subsequently, a projection operator is defined to dissect the phase-space dynamics between the function space, exclusive to slow variables, and the remainder.
Through the application of integral identities associated with projection operators, GLEs are derived.A notable merit of this approach is the formal derivation of GLEs from microscopic dynamics, providing a rigorous connection to underlying physical processes.However, a significant drawback lies in the inherent ambiguity of the approximation involved.While all calculations are theoretically exact, eliciting nontrivial predictions mandates the approximate computation of noise statistics and friction coefficients.This level of approximation is notably more intricate compared to conventional statistical-physics theories.In fact, the determination of theoretical key perturbation/control parameters for the conclusive deduction of the GLEs from microscopic dynamics remains unambiguous, making this process elusive.
Within the foundational framework of statistical physics pertaining to the GLEs, a drawback of our theory is the requisite assumption of the one-dimensional non-Markovian jump process (17) as an initial standpoint.This assumption is fundamentally heuristic, primarily rooted in phenomenological considerations.Conversely, a significant advantage of our approach is the explicit definition of the key perturbation parameter.Specifically, the small-jump scaling parameter, ϵ -generally anticipated to represent the mass ratio between the Brownian particle and surrounding entities -serves a crucial and explicit role in our asymptotic computations.This is particularly coherent for modeling dynamics of massive Brownian particles.In this sense, we have successfully established the GLEs through a physically plausible coarse-graining process, pinpointing the essential control parameter for mathematical derivation, a contrast to the methodologies embedded in the projection-operator formalism.

B. Future issue 1: physical validation of the non-Markovian jump model
Our theory is premised on the non-Markovian jump model (17).While intrinsic to models in seismic activity, finance, and social science -for instance, the Hawkes process is a subset of this model -its applicability in physics remains indeterminate.Addressing this uncertainty will necessitate further theoretical or data-driven analysis in the future.
From a theoretical standpoint, the Markovian ME formalism (8) has been substantiated in the dynamics of Brownian particles amid dilute gases [57][58][59].Indeed, the linearized Boltzmann equation, derivable from Newtonian microscopic dynamics through the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy [57,58] in the low-density limit, is an instance of a Markovian jump process, thus allowing systematic theoretical validation of the Markovian ME formalism (8) via kinetic theory.
Conversely, theoretical validation for the non-Markovian jump model ( 17) is yet to be achieved.Formulating a statistical physics theory analogous to the Markovian kinetic theory to derive the non-Markovian jump process (17) from microscopic Hamiltonian dynamics is imperative.
C. Future issue 2: time-reversal symmetry of our field master equation Exploring time-reversal symmetry is crucial when examining stochastic dynamics influenced by equilibrium fluctuations.Regrettably, this symmetry is not upheld for the field ME (26).The indispensable condition for general master equations, fully detailed in Gardiner's textbook [22] and Appendix G, is invariably breached in our field ME (26).
The absence of time-reversal symmetry in our Laplace-type embedding representation can be intuitively understood by considering a typical path of {ẑ t (s)} s .Indeed, the SPDE (24a) states that the dynamics of ẑt (s) is composed of the excitation due to the Poisson jump ξCP The preservation of time-reversal symmetry in a ME is significantly contingent upon the choice of state variables [60].To illustrate, a Markovian jump model may maintain time-reversal symmetry in a one-variable representation, but it invariably loses this symmetry in a two-variable depiction.Let us assume a Markov jump process with the timereversal symmetry (see Fig. 8(a)): where ŷi and ti are the ith jump size and ith jump time, respectively.The ME for this one-variable representation can satisfy the time-reversal symmetry by setting an appropriate intensity density λ(y|v).
On the other hand, this process is equivalent to a discrete-time Markovian jump process with two variables where i is an integer time incremented at every jump event, Vi := vt i is the velocity at the ith jump time, and ∆ ti is the waiting time obeying an exponential distribution.This technique is called the subordination in the context of the continuous-time random walk theory [31] (see Fig. 8(b)).The ME for the two-variable representation Γi := ( Vi , ti ) has no time-reversal symmetry due to the monotonically increasing nature of ti , where ∆P i (Γ) := P i+1 (Γ) − P i (Γ) is the discrete-time difference operator, and L is a linear operator for the discrete-time ME.This mathematical fact suggests that there are several different Markovian-embedding formulations of the ME even if the stochastic dynamics is uniquely defined, and the time-reversal symmetry might be formulated for a specific Markovian-embedding representation.Therefore, another Markovian-embedding representation might be suitable for formulating the time-reversal symmetry.Resolving this issue will be instrumental in shaping the formulation of stochastic thermodynamics and energetics [61,62] for non-Markovian jump processes.This is reserved for the future.
where the left-hand side is the cross-correlation between ηt1 and ηt2 , T is the temperature and we have taken units where the Boltzmann constant is unity.This fluctuation-dissipation relation of the second kind is equivalent to the time-reversal symmetry of the GLE.This was one of the most important issues for the statistical-physics foundation of the GLE particularly within the context of linear response theory [1] and the projection-operator formalism [28].
Since our Markovian-embedding formulation does not yet convert the time-reversal symmetry of the field ME, the necessary and sufficient condition for this fluctuation-dissipation relation of the second kind is not yet identified.The identification of these conditions is also an important future challenge.
E. Future issue 4: formal relations to quantum field theory Our field master equation (Fokker-Planck) is formally related to quantum field theory.Indeed, the field Fokker-Planck equation for the GLE with time-reversal symmetry is equivalent to a non-Hermitian quantum field theory with the Hermitian part of its Hamiltonian describing a field of harmonic oscillators [33].Indeed, a similar renormalisation issue appears regarding the infinite zero-point energy of the field harmonic oscillators.Since the first-order contribution of the system-size expansion for the non-Markovian jump process leads to the GLE, its next-order perturbation theory might require the use of methods developed in quantum field theory, such as the Feynman-diagram expansion.Establishing such field-theoretical techniques will be an interesting future topic.

VIII. CONCLUSION
In conclusion, we have introduced a comprehensive stochastic framework through a field master equation, encompassing all one-dimensional non-Markovian jump processes.Utilizing the Laplace-convolution embedding representation, we have demonstrated the transformation of any non-Markovian jump process into Markovian-field dynamics.We subsequently derived the corresponding field master equation and procured an asymptotic solution using a generalized system-size expansion.In essence, this framework can be applied to any jump processes, assuming one-dimensional dynamics are driven by collisions.We posit that this model's flexibility makes it adept at accommodating a wide array of point-process data, proving invaluable for data analyses.
with df (z) := f (z + dz) − f (z) for infinitesimal dz := (dz 1 , . . ., dz K ).In the continuous limit, we apply the replacement to obtain the first-order functional Taylor expansion b. Full-order functional Taylor expansion For a finite-dimensional vector function f (z), the full-order Taylor expansion is given by with ∆z := (∆z 1 , . . ., ∆z K ).In the continuous limit based on the formal replacement (A7), we obtain the full-order functional Taylor expansion We note that this calculation can be readily generalised for a two-argument functional f [v; z] with a real value v and a function {z(s)} s , such that with small ∆v and {∆z(s)} s .Particularly, the Maclaurin series is given by In the small-noise limit ε → 0, we obtain the FP equation Here we show a trivial approach of Markovian embedding available only for discrete-time stochastic processes.Similar ideas are used in Econometrics [3] regarding the lag operator.If the time is discrete, this formulation can be straigtforwardly generalised even for K → ∞, where the embedding dimension is infinite and thus the dynamics is truly non-Markovian.

Technical contribution of the Laplace-convolution representation
This fact implies that Markovian embedding is trivial for discrete-time stochastic processes.However, a straightforward generalisation of this specific embedding is difficult for continuous-time stochastic processes.Indeed, it is challenging to generalise the shifting operator S for continuous-time representations, even at a formal level.
Let us attempt to write the formal continuous representation from the naive discrete-time embedding equation (D3).By considering the continuous limit with K → ∞ and dt → 0 for the time interval, let us write the phase-space vector as { Γt (s)} s≥0 parametrised with s ≥ 0 defined by This equation does not make sense even at the theoretical physics level due to the apparent singularity of the δ function and its derivative.Thus, the naive embedding (D1) for the discrete-time processes cannot be straightforwardly generalised toward the continuous-time processes, even at the formal level.The Laplace-convolution representation technically solves this problem.The shifting operator S has an analytically tractable representation in the Laplace-convolution space, and, thus, the original non-Markovian dynamics is mapped onto a first-order Markovian SPDE.This equation is equivalent to the field ME in Refs.[34,35].
Appendix G: Time-reversal symmetry of the master equation We review the necessary and sufficient condition of the validity of time-reversal symmetry according to Ref. [22].For a finite-dimensional Markovian stochastic processe x := (x 1 , . . ., x K ) T , the general master equation is given by   P t (x) + dy [λ(x|y)P t (y) − λ(y|x)P t (x)] , (G1) where A k is the drift term, B kl ≥ 0 is the diffusion term, and λ(y|x) ≥ 0 is the jump-intensity density for jumps from x to y.
Let us define the time-reversal operator ϵ k such that ϵ k = 1 if x k is an even variable, and ϵ k = −1 if x k is an odd variable.Typically, the velocity (position) is an odd (even) variable because it has the odd (even) parity under time reversal.The necessary and sufficient condition for time-reversal symmetry to hold is given by λ(y|x)P ss (x) = λ(ϵx|ϵy)P ss (y) (G2a)

FIG. 7 :FIG. 8 :
FIG.7: Absence of time-reversal symmetry for the auxiliary field variables {ẑt(s)}.According to the SPDE (24a), the path of ẑt(s) responds to the excitations due to jumps and then relaxes toward zero.The relaxation is time-irreversible, and the auxiliary field variables have thus no time-reversal symmetry.

D. Future issue 3 :
the fluctuation-dissipation relation of the second kind When the environment is in thermal equilibrium, the thermal fluctuation of the GLE must satisfies the fluctuationdissipation relation (FDR) of the second kind: ) where the dummy argument variables χ and ζ are introduced to distinguish the arguments involved in the derivatives from the arguments v and z of the function f [v; z].c.Variable transformation formula Let us consider a simple variable transformation s = as (A13) with a positive constant a. Considering the definition (A5)leads to the invariant integral relationship, ds δ δz(s) = ds δ δz(s) .