A BBGKY-like Hierarchy for Quantum Field Theories

Current models of interacting charged particles rely on approximation methods derived from classical electrodynamics. These schemes enjoy great success in relatively weakly interacting environments; however, their predictive power is wanting at higher energies. Large radiative corrections occur when particles begin emitting photons of energy comparable to their mass blackburn . This is not to mention more exotic corrections from effects such as pair creation paircreation . With pulsed laser facilities currently achieving radiation intensities of over $10^{23}W/cm^{2}$ plasmachallenges , and next-generation burning plasma experiments set to probe comparable energy scales, more sophisticated models must be developed to achieve predictive capabilities. These models must be semi-classical and therefore derive from quantum electrodynamics.

A common starting part for many classical reduced models is the Bogoliubov–Born–Green–
Kirkwood–Yvon (BBGKY) hierarchy. The BBGKY hierarchy is a way of organizing the equations of motion for the $N$ particle distribution function as a hierarchical system of equations for the $k\leq N$ point functions. Through various means, this hierarchy can be reduced to a small number of equations, providing workable approximations to the exact time evolution of a system.

In MMW , the BBGKY hierarchy was shown to have a Hamiltonian structure. In addition to being an important theoretical revelation, this structure offered insight into how BBGKY hierarchies could be constructed for other Hamiltonian theories. Indeed, in Nparticle , and separately in otherpaper , a BBGKY-like hierarchy was obtained for $N$ boson quantum mechanics. Using a filtration on quantum operators, as opposed to phase space functions, both papers obtained equations for the $k$th order reduced density matrices very similar to those of the classical BBGKY hierarchy.

In this paper, we generalize the constructions of MMW , Nparticle , and otherpaper to arbitrary quantum field theories. Using a filtration of Hermitian polynomials in creation/annihilation, we are able to derive a Hamiltonian BBGKY-like hierarchy for the field-theoretic reduced-density matrices. Unlike both the classical and $N$-boson hierarchies, our hierarchy has infinitely many equations. Further, compared to the simple coupling of the other hierarchies, our hierarchy generically couples the $k$th order variable to the $k+1$ through $k+4$th order variables. Despite these differences, methods to truncate the classical BBGKY hierarchy can be applied to our hierarchy without issue, rendering our equations small in number and thus computable.

To the best of our knowledge, the hierarchy we present is novel. We note, however, that a hierarchy of n-particle Wigner functions was constructed for $\lambda\phi^{4}$ theory in scalar . While quite different from the approach we employ here, it remains open whether a connection between our hierarchies exists.

With our work, we hope to kickstart the field of using geometric mechanics and quantum field theory to study high-energy-density, many-body systems. Such a viewpoint seems not only natural but also powerful, as witnessed by the relative brevity of our constructions. With further investigation along these lines, we trust that a much greater understanding of high-energy-density environments can be achieved.

In III, we generalize the works of MMW , Nparticle , and otherpaper to provide a framework from which a wide class of quantum field theoretic BBGKY-like hierarchies may be constructed. These hierarchies preserve the Hamiltonian structure of the field theory and take a form almost identical to the classical BBGKY hierarchy. In IV, we specialize our discussion to the hierarchy formed from the filtration of ladder operator polynomials, which we dub the canonical hierarchy. We show that the canonical hierarchy has many desirable properties, making it particularly well-suited to approximation schemes. In particular, relevant observables such as the energy and spectral densities depend only on the lowest-level hierarchy variables and are thus easily computed without recourse to further approximations. We demonstrate this in V. Finally, in VI, we offer several avenues of further work including a method of obtaining Hamiltonian closures to our hierarchy.

For an arbitrary quantum field theory, we let $\mathfrak{g}$ denote the space of hermitian operators. This space forms a lie-algebra under the usual commutator bracket $-i[\cdot,\cdot]$. A distinguished subset of the dual space $\mathfrak{g}^{*}$ is the space of trace-class, positive-semidefinite, hermitian operators $\mathfrak{d}$ dualized under the trace map $\rho\in\mathfrak{d}:A\mapsto\left<\rho,A\right>\coloneqq\text{Tr}(A\rho)$. By normalizing the trace, we may identify every nonzero element of $\mathfrak{d}$ with a density matrix. Provided we restrict to unitary equations of motion, $\mathfrak{d}$ can therefore be thought of as the phase space of the quantum theory.

We equip $\mathfrak{d}$ with a Poisson structure by defining the bracket $\{\cdot,\cdot\}:C^{\infty}(\mathfrak{d})\times C^{\infty}(\mathfrak{d})\to C^{\infty}(\mathfrak{d})$

$$\{F,G\}(\rho)=\left<\rho,-i\left[\frac{\delta F}{\delta\rho},\frac{\delta G}{\delta\rho}\right]\right>,$$

(1)

where $\frac{\delta F}{\delta\rho}$ is the unique element of $\mathfrak{g}$ satifying

$$\frac{d}{d\epsilon}|_{\epsilon=0}F(\rho+\epsilon\delta\rho)=\left<\delta\rho,\frac{\delta F}{\delta\rho}\right>.$$

(2)

Given a Hamiltonian $H\in\mathfrak{g}$, the Hamiltonian flow of $H$ (viewed as a function on $\mathfrak{d}$) provides the usual quantum mechanical equations of motion. Indeed, if $\rho(t):J\subset\mathbb{R}\to\mathfrak{d}$ integrates the Hamiltonian flow, then for any $F\in C^{\infty}(\mathfrak{d})$,

$$\frac{d}{dt}F(\rho(t))=\left<\frac{\delta F}{\delta\rho},\frac{d}{dt}\rho(t)\right>=\{F,H\}(\rho(t))=-i\left<\rho(t),\left[\frac{\delta F}{\delta\rho},H\right]\right>=iTr\left(\frac{\delta F}{\delta\rho}\left[\rho(t),H\right]\right),$$

(3)

where the first equality follows from the chain rule and the last from the cyclic property of the trace. It is easy to see that this equality holds for all $F\in C^{\infty}(\mathfrak{d})$ if and only if

$$\frac{d}{dt}\rho(t)=i[\rho(t),H],$$

(4)

which is the desired relation.

Let $\mathfrak{g}_{1}\subset\mathfrak{g}_{2}\subset\mathfrak{g}_{3}\subset\ldots\subset\mathfrak{g}$ be a sequence of operator subspaces such that $-i[\mathfrak{g}_{n},\mathfrak{g}_{m}]\subset\mathfrak{g}_{n+m-1}$ and $\bigcup_{i}\mathfrak{g}_{i}$ is dense in $\mathfrak{g}$. We call such a sequence a filtration. For any filtration, we define

$$\mathfrak{G}\coloneqq\bigoplus_{i=1}^{\infty}\mathfrak{g}_{i}.$$

(5)

We equip $\mathfrak{G}$ with the lie-bracket

$$[(A_{i}),(B_{i})]_{\mathfrak{G}}=-i\left([A_{1},B_{1}],[A_{2},B_{1}]+[A_{1},B_{2}],\ldots,\sum_{n+m=k}[A_{n},B_{m}],\ldots\right).$$

(6)

As is the case for any direct sum, $\mathfrak{G}^{*}=\prod_{i=1}^{\infty}\mathfrak{g}^{*}_{i}$. Defining the subspaces
$\mathfrak{d}_{i}=\{\rho|_{\mathfrak{g}_{i}}:\rho\in\mathfrak{d}\}\subset\mathfrak{g}_{i}^{*}$, we have the distinguished subspace of $\mathfrak{G}^{*}$

$$\mathfrak{D}\coloneqq\prod_{i=1}^{\infty}\mathfrak{d}_{i}.$$

(7)

This space has the canonical Poisson bracket

$$\{F,G\}_{\mathfrak{D}}((\rho_{i}))=\left<(\rho_{i}),\left[\frac{\delta F}{\delta(\rho_{i})},\frac{\delta G}{\delta(\rho_{i})}\right]_{\mathfrak{G}}\right>.$$

(8)

Writing $\frac{\delta F}{\delta(\rho_{i})}=(F_{i})\in\mathfrak{G}$, this is equivalent to the formula

	$\displaystyle\left\{F,G\right\}((\rho_{i}))$	$\displaystyle=-iTr\left(\rho_{1}[F_{1},G_{1}]\right)-iTr(\rho_{2}[F_{2},G_{1}])-iTr(\rho_{2}[F_{1},G_{2}])-iTr(\rho_{3}[F_{3},G_{1}])+\ldots.$		(9)
		$\displaystyle=iTr\left(F_{1}[\rho_{1},G_{1}]\right)+iTr(F_{2}[\rho_{2},G_{1}])+iTr(F_{1}[\rho_{2},G_{2}])+iTr(F_{3}[\rho_{3},G_{1}])+\ldots.$

Let $\alpha:\mathfrak{G}\to\mathfrak{g}$ denote the natural inclusion $\alpha:(A_{i})\mapsto\sum_{i}A_{i}$. One can easily verify that $\alpha$ is a lie-algebra homomorphism, and consequently that the dual function $\alpha^{*}:\mathfrak{d}\to\mathfrak{D},\,\alpha^{*}:\rho\mapsto(\rho|_{\mathfrak{g}_{i}})$ is a Poisson map. If $(H_{i})\in\mathfrak{G}$ is any sequence of operators such that $\alpha((H_{i}))=H$, and $\rho(t)$ evolves according to the flow on $\mathfrak{d}$, then $\alpha^{*}\rho(t)$ integrates the flow of $(H_{i})$. This can be verified directly by checking that, for any $F\in C^{\infty}(\mathfrak{D})$,

$$\frac{d}{dt}F(\alpha^{*}\rho(t))=\left<(F_{i}),\alpha^{*}(\rho(t))\right>=\left\{F,(H_{i})\right\}(\alpha^{*}\rho(t)).$$

(10)

We may view the restricted dual operators $\rho|_{\mathfrak{g}_{i}}$ as equivalence classes $[\rho]_{i}$ of operators under the quotient $\rho_{a}\sim_{i}\rho_{b}$ if and only if $Tr(\rho_{a}A_{i})=\rho(\rho_{b}A_{i})$ for all $A_{i}\in\mathfrak{g}_{i}$. Letting $\rho_{i}\in[\rho]_{i}$ and $\pi_{k}:A\mapsto[A]_{k}$ be the quotient map, we define $[[\rho]_{i},A_{i}]\coloneqq\pi_{i}[\rho_{i},A_{i}]$ for any $A_{i}\in\mathfrak{g}_{i}$. Some simple algebra reveals that (10) holds if and only if

$$\frac{\partial}{\partial t}[\rho]_{i}=i[[\rho]_{i},H_{1}]+i[[\rho]_{i+1},H_{2}]+i[[\rho]_{i+2},H_{3}]+\ldots.$$

(11)

Provided that $H\in B_{N}$ for some small $N$, (11) reveals that the equations for $[\rho]_{i}$ form a hierarchical structure. Furthermore, these equations are Hamiltonian. Truncating or otherwise closing the hierarchy, we obtain approximate equations for $[\rho]_{1},\ldots,[\rho]_{N}$. Approximate expectation values are then obtained for every operator $O$ such that $Tr(O\rho)$ depends only on the equivalence class $[\rho]_{N}$.

Given some set of operators $\{X_{\alpha}\}$, the sequence of subspaces

$$B_{k+1}=\{\text{kth order Hermitian polynomials in }X_{\alpha},X_{\alpha}^{\dagger}\}$$

(12)

will always form a filtration of $\mathfrak{g}$. In particular, we will see it is desirable to set $\{X_{\alpha}\}=\{a_{\mathbf{p},I}\}$, the annihilation operators of the quantum field theory. Here, $\mathbf{p}$ represents the particle momentum index, and $I$ represents any discrete indices (e.g. species, helicity, etc.). This filtration is particularly suited to approximations since, for any renormalizable theory, we may write that $H=\sum_{i=2}^{5}H_{i}$ with $H_{i}\in B_{i}-B_{i-1}$. Defining the reduced density matrices

$$\Gamma^{(m,n)}_{(I^{\prime}_{1},\ldots,I^{\prime}_{n},I_{1},\ldots,I_{m})}[\rho](\mathbf{p}_{1},...,\mathbf{p}_{m};\mathbf{p}^{\prime}_{1},...,\mathbf{p}^{\prime}_{n})\coloneqq\text{Tr}(\rho a_{\mathbf{p}^{\prime}_{1},I_{1}^{\prime}}^{\dagger}...a_{\mathbf{p}^{\prime}_{n},I_{n}^{\prime}}^{\dagger}a_{\mathbf{p}_{1},I_{1}}...a_{\mathbf{p}_{m},I_{m}}),$$

(13)

this choice of filtration further allows the equivalence classes $[\rho]_{k}$ to be linearly encoded into a finite collection $\Gamma_{k}[\rho]=(\Gamma^{(m,n)}_{\mathcal{I}}[\rho]:m+n<k)$, say by the map $\epsilon_{k}:[\rho]_{k}\mapsto\Gamma_{k}[\rho]$. One easily verifies that $\epsilon_{k}$ is injective and hence invertible on its image. Since applying $\epsilon_{k}$ commutes with time derivatives, the equations of motion for $\Gamma_{k}$ can be obtained from $\eqref{eq:10}$ as

$$\frac{\partial}{\partial t}\Gamma_{i}=i\epsilon_{i}[\epsilon^{-1}_{i+1}(\Gamma_{i+1}),H_{2}]+i\epsilon_{i}[\epsilon^{-1}_{i+2}(\Gamma_{i+2}),H_{3}]+i\epsilon_{i}[\epsilon^{-1}_{i+3}(\Gamma_{i+3}),H_{4}]+i\epsilon_{i}[\epsilon^{-1}_{i+4}(\Gamma_{i+4}),H_{5}].$$

(14)

The power of these equations lies in the fact that they are hierarchical and can thus be closed at finite order using methods originally developed for the classical BBGKY hierarchy. The simplest such method is to set $\Gamma_{i}\equiv 0$ for $i>N$. In this case, the approximate equations of motion read

		$\displaystyle\frac{\partial}{\partial t}\Gamma_{1}=i\epsilon_{i}[\epsilon^{-1}_{i+1}(\Gamma_{2}),H_{2}]+i\epsilon_{i}[\epsilon^{-1}_{i+2}(\Gamma_{3}),H_{3}]+i\epsilon_{i}[\epsilon^{-1}_{4}(\Gamma_{i+3}),H_{4}]+i\epsilon_{i}[\epsilon^{-1}_{5}(\Gamma_{i+4}),H_{5}],$		(15)
	$\displaystyle\;\;\;\;\;\;\;\vdots$
		$\displaystyle\frac{\partial}{\partial t}\Gamma_{N-3}=i\epsilon_{i}[\epsilon^{-1}_{N-2}(\Gamma_{2}),H_{N-2}]+i\epsilon_{i}[\epsilon^{-1}_{N-1}(\Gamma_{3}),H_{N-1}]+i\epsilon_{i}[\epsilon^{-1}_{4}(\Gamma_{N}),H_{4}],$
		$\displaystyle\frac{\partial}{\partial t}\Gamma_{N-3}=i\epsilon_{i}[\epsilon^{-1}_{N-2}(\Gamma_{2}),H_{2}]+i\epsilon_{i}[\epsilon^{-1}_{N-1}(\Gamma_{3}),H_{3}]+i\epsilon_{i}[\epsilon^{-1}_{4}(\Gamma_{N}),H_{4}],$
		$\displaystyle\frac{\partial}{\partial t}\Gamma_{N-2}=i\epsilon_{i}[\epsilon^{-1}_{N-1}(\Gamma_{2}),H_{2}]+i\epsilon_{i}[\epsilon^{-1}_{N}(\Gamma_{3}),H_{3}],$
		$\displaystyle\frac{\partial}{\partial t}\Gamma_{N-1}=i\epsilon_{i}[\epsilon^{-1}_{N}(\Gamma_{2}),H_{2}],$
		$\displaystyle\frac{\partial}{\partial t}\Gamma_{N}=0.$

A more sophisticated method involves cluster expanding the reduced density matrices into their correlations below the $N$th order. This scheme similarly reduces $\eqref{eq:13}$ to $N$ equations. However, this method represents a superior weak-interaction limit since any relativistic quantum field theory has asymptotically vanishing correlations. In the simplest case of a scalar theory and $N=1$, this amounts to approximating

$$\Gamma^{(m,n)}(\mathbf{p}_{1},...,\mathbf{p}_{m};\mathbf{p}^{\prime}_{1},...,\mathbf{p}^{\prime}_{n})\approx\overline{\Gamma^{(1,0)}}(\mathbf{p}_{1}^{\prime})\ldots\overline{\Gamma^{(1,0)}}(\mathbf{p}_{n}^{\prime})\ldots\Gamma^{(1,0)}(\mathbf{p}_{1})\ldots\Gamma^{(1,0)}(\mathbf{p}_{m}).$$

(16)

Substituting this expression into $\eqref{eq:13}$ for $i=1$, we obtain an equation of the form $\partial_{t}\Gamma_{1}=F(\Gamma_{1})$ which can be solved and substituted into $\eqref{above}$.

A major fault in both of these methods is that they do not produce Hamiltonian equations of motion, and hence break the underlying structure of the quantum theory. Ideally, approximations to (14) should close the hierarchy in a Hamiltonian manner. However, a general method for doing so remains open. In VI, we offer several possible routes toward this goal.

We note that once an approximation method has been chosen, the equations of motion are most effectively computed using the relation

$$i\epsilon_{i}[\epsilon^{-1}_{j}(\Gamma_{\mathcal{I}}^{(m,n)}),A]=-i\text{Tr}(\rho\,[a_{\mathbf{p}^{\prime}_{1},I_{1}^{\prime}}^{\dagger}...a_{\mathbf{p}^{\prime}_{n},I_{n}^{\prime}}^{\dagger}a_{\mathbf{p}_{1},I_{1}}...a_{\mathbf{p}_{m},I_{m}},A]).$$

(17)

For example, consider a free scalar theory with the Hamiltonian

$$H=\int\frac{d^{3}\mathbf{q}}{(2\pi)^{3}(2E_{\mathbf{q}})}E_{\mathbf{q}}a^{\dagger}_{\mathbf{q}}a_{\mathbf{q}}\in B_{3},$$

(18)

where $E_{\mathbf{p}}=\sqrt{\mathbf{p}^{2}+m^{2}}$ is the particle energy. We work in natural units and assume the commutation relations

$$[a_{\mathbf{p},\mathcal{I}},a^{\dagger}_{\mathbf{p^{\prime}},\mathcal{I}^{\prime}}]=(2\pi)^{3}(2E_{\mathbf{p}})\,\delta_{I,I^{\prime}}\,\delta^{(3)}(\mathbf{p}-\mathbf{p}^{\prime}).$$

(19)

A simple calculation reveals that

$$\frac{\partial}{\partial t}\Gamma^{(n,m)}(t;\mathbf{p}_{1},\ldots,\mathbf{p}_{n};\overline{\mathbf{p}}_{1}\ldots,\overline{\mathbf{p}}_{m})=\left(-i\sum_{j=1}^{n}\mathbf{p}_{j}t+i\sum_{j=1}^{m}\overline{\mathbf{p}}_{j}t\right)\Gamma^{(n,m)},$$

(20)

with similarly trivial behavior for other free theories.

For interacting theories, $\eqref{eq:13}$ is no longer exactly solvable. If the scalar Hamiltonian was instead

$$H=\frac{1}{2}\int\frac{d^{3}\mathbf{q}}{(2\pi)^{3}(2E_{\mathbf{q}})}\frac{d^{3}\mathbf{s}}{(2\pi)^{2}(2E_{\mathbf{s}})}h(\mathbf{q},\mathbf{s})a_{\mathbf{q}}^{\dagger}a_{\mathbf{s}}^{\dagger}a_{\mathbf{q}}a_{\mathbf{s}}$$

(21)

with $h(\mathbf{q},\mathbf{s})$ a symmetric, real-valued function, then a similar computation reveals that

	$\displaystyle\frac{\partial}{\partial t}\Gamma^{(1,0)}(t;\mathbf{p})$	$\displaystyle=-i\int\frac{d^{3}\mathbf{q}}{(2\pi)^{3}(2E_{\mathbf{q}})}\;\;\;\Gamma^{(2,1)}(t;\mathbf{p},\mathbf{q};\mathbf{q})h(\mathbf{p},\mathbf{q}),$		(22)
	$\displaystyle\frac{\partial}{\partial t}\Gamma^{(1,1)}(t;\mathbf{p};\mathbf{p}^{\prime})$	$\displaystyle=-i\int\frac{d^{3}\mathbf{q}}{(2\pi)^{3}(2E_{\mathbf{q}})}\;\;\Gamma^{(2,2)}(t;\mathbf{p},\mathbf{q};\mathbf{p}^{\prime},\mathbf{q})\left(h(\mathbf{p},\mathbf{q})-h(\mathbf{p}^{\prime},\mathbf{q})\right),$
	$\displaystyle\vdots$

which is no longer trivial, highlighting the general need to seek approximations to $\eqref{eq:13}$.

Perhaps the simplest observable is the expected density of particles with discrete indices $I$ and momentum $\mathbf{p}$,

$$D_{I}(\mathbf{p})\coloneqq\left<\frac{1}{(2\pi)^{3}2E_{\mathbf{p}}}a_{\mathbf{p},I}^{\dagger}a_{\mathbf{p},I}\right>=\frac{1}{(2\pi)^{3}2E_{\mathbf{p}}}\Gamma_{(I,I)}^{(1,1)}(\mathbf{p};\mathbf{p}).$$

(23)

Another observable of interest is the total energy stored in a region of space. For any relativistic theory, we may write the Hamiltonian as

$$H=\int d^{3}\mathbf{x}\;\mathcal{H}(\mathbf{x}).$$

(24)

where $\mathcal{H}(\mathbf{x})$ is the $(0,0)$ component of the physical stress-energy tensor. The energy distribution of a system is then

$$\mathcal{E}(\mathbf{x})\coloneqq\left<\mathcal{H}(\mathbf{x})\right>=F(\Gamma_{5}).$$

where $F$ is some function depending on the Hamiltonian. For interacting theories, $\left<\mathcal{H}(\mathbf{x})\right>$ is in general a very complicated function. In many cases, however, the energy stored in particle interactions is much smaller than the energy of the particles themselves. For such cases, we may approximate

$$\left<\mathcal{H}(\mathbf{x})\right>\approx\left<\mathcal{H}_{\text{free}}(\mathbf{x})\right>=F_{\text{free}}(\Gamma_{3}),$$

(25)

with $\mathcal{H}_{\text{free}}$ being the Hamiltonian density of the free theory. For example, for an arbitrary scalar theory, we may use $\eqref{eq:17}$ to approximately compute that

	$\displaystyle\mathcal{E}(\mathbf{x})$	$\displaystyle\approx\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}(2E_{\mathbf{k}})}\;\frac{d^{3}\mathbf{p}}{(2\pi)^{3}(2E_{\mathbf{p}})}\;$		(26)
		$\displaystyle(m^{2}-E_{\mathbf{p}}E_{\mathbf{k}}+\mathbf{p}\cdot\mathbf{k})e^{i\mathbf{x}\cdot(\mathbf{k}+\mathbf{p})}\Gamma^{(2,0)}(t;\mathbf{k},\mathbf{p})+e^{-i\mathbf{x}\cdot(\mathbf{k}-\mathbf{p})}(m^{2}+E_{\mathbf{p}}E_{\mathbf{k}}+\mathbf{p}\cdot\mathbf{k})\Gamma^{(1,1)}(t;\mathbf{p};\mathbf{k}).$

For a massive theory, we may further consider the limit whereby the particle mass is much larger than the particle momenta, and hence that the energy of the fields resides entirely in the mass energy of the particles. Letting $\left<A\right>_{\mathbf{p}\to 0}$ denote the expectation value of an operator $A$ in such a limit, the number density of particles can thus be computed as

$$N(\mathbf{x})=\frac{\left<\mathcal{H}_{\text{free}}(\mathbf{x})\right>_{\mathbf{p}\to 0}}{m}.$$

(27)

In the case of a scalar theory, for example, the expected particle density is

$\displaystyle N(\mathbf{x})=\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}\sqrt{2E_{\mathbf{k}}}}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}\sqrt{2E_{\mathbf{p}}}}e^{-i\mathbf{x}\cdot(\mathbf{k}-\mathbf{p})}\Gamma^{(1,1)}(t;\mathbf{p};\mathbf{k}),$

(28)

which is verified as the expectation value of the operator

$$\mathcal{N}(\mathbf{x})=\int\frac{d^{3}\mathbf{k}}{(2\pi)^{3}\sqrt{2E_{\mathbf{k}}}}\frac{d^{3}\mathbf{p}}{(2\pi)^{3}\sqrt{2E_{\mathbf{p}}}}e^{-i\mathbf{x}\cdot(\mathbf{k}-\mathbf{p})}a_{\mathbf{k}}^{\dagger}a_{\mathbf{p}}.$$

(29)

Indeed, the interpretation of $\mathcal{N}(\mathbf{x})$ as a number-density is confirmed upon checking that $\int d^{3}\mathbf{x}\;\mathcal{N}(\mathbf{x})$ is the usual number operator. Similar operators for other theories are obtained from (29) by adding the appropriate discrete indices and external leg factors.

While an important first step, several more steps must be taken before our hierarchy can be used in the semi-classical modeling of high-energy environments. For one, before any computational work is undertaken, $\eqref{eq:13}$ must be recast as a differential equation and the desired initial conditions rewritten as reduced density matrices. This is straightforward in the simplest cases but requires considerable effort for realistic theories such as quantum electrodynamics. A more subtle issue is how nontrivial boundary conditions may be worked into our theory. It is clear this involves working in position space, but we leave the details to future work.

A more theoretical line of future work lies in finding Hamiltonian closures to our hierarchy equations. While it is desirable that approximations to the quantum evolution equations preserve the Hamiltonian structure of theory, it is not clear how this can be done. One promising line of investigation is to look for Poisson maps from some smaller Poisson manifold $P$ to $\mathfrak{D}$. Solutions to the Hamiltonian flow on $P$ would then push forward to solutions on $\mathfrak{D}$, while presumably being easier to solve. Given an operator ideal $\mathcal{I}$, a natural choice for $P$ would be the subspace of $\mathfrak{D}$ consisting of sequences $(\rho_{1},\rho_{2},\ldots)$ such that $\rho_{i}(\mathfrak{g}_{i}\cap\mathcal{I})=0$, with the Poisson bracket inherited from the quotient algebra $\mathfrak{g}/\mathcal{I}$. However, such ideals remain elusive if they exist.

In a similar vein, one could look for Poisson structures on $\prod_{i=1}^{N}\mathfrak{d}_{i}$ close to the Poisson structure on $\prod_{i=1}^{\infty}\mathfrak{d}_{i}$. Although no longer an exact solution, solving for the flow on the smaller space would still serve as a Hamiltonian approximation to the full equations of motion.