# A gauge-symmetrization method for energy-momentum tensors in high-order
electromagnetic field theories

Peifeng Fan  Key Laboratory of Optoelectronic Devices and Systems, College of
Physics and Optoelectronic Engineering, Shenzhen University, Shenzhen 518060,
China  Advanced Energy Research Center, Shenzhen University, Shenzhen 518060,
China  Jianyuan Xiao  School of nuclear science and technology, University of
Science and Technology of China, Hefei, Anhui 230026, China  Hong Qin
Princeton Plasma Physics Laboratory, Princeton University, Princeton, NJ
08543, USA

###### Abstract

For electromagnetic field theories, canonical energy-momentum conservation
laws can be derived from the underpinning spacetime translation symmetry
according to the Noether procedure. However, the canonical Energy-Momentum
Tensors (EMTs) are neither symmetric nor gauge-symmetric (gauge invariant).
The Belinfante-Rosenfeld (BR) method is a well-known procedure to symmetrize
the EMTs, which also render them gauge symmetric for first-order field
theories. High-order electromagnetic field theories appear in the study of
gyrokinetic systems for magnetized plasmas and the Podolsky system for the
radiation reaction of classical charged particles. For these high-order field
theories, gauge-symmetric EMTs are not necessarily symmetric and vice versa.
In the present study, we develop a new gauge-symmetrization method for EMTs in
high-order electromagnetic field theories. The Noether procedure is carried
out using the Faraday tensor $F_{\mu\nu}$, instead of the 4-potential
$A_{\mu}$, to derive a canonical EMT $T_{\text{N}}^{\mu\nu}$. We show that the
gauge-dependent part of $T_{\text{N}}^{\mu\nu}$ can be removed using the
displacement-potential tensor
$\mathcal{F}^{\sigma\mu\nu}\equiv\mathcal{D}^{\sigma\mu}A^{\nu}/4\pi$, where
$\mathcal{D}^{\sigma\mu}$ is the anti-symmetric electric displacement tensor.
This method gauge-symmetrize the EMT without necessarily making it symmetric,
which is adequate for applications not involving general relativity. For
first-order electromagnetic field theories, such as the standard Maxwell
system, $\mathcal{F}^{\sigma\mu\nu}$ reduces to the familiar BR super-
potential $\mathcal{S}^{\sigma\mu\nu}$, and the method developed can be used
as a simpler procedure to calculate $\mathcal{S}^{\sigma\mu\nu}$ without
employing the angular momentum tensor in 4D spacetime. When the
electromagnetic system is coupled to classical charged particles, the gauge-
symmetrization method for EMTs is shown to be effective as well.

##  I introduction

In classical field theories, one can derive canonical Energy-Momentum Tensors
(EMTs) $T_{\text{N}}^{\mu\nu}$ from the underpinning spacetime translation
symmetry using the Noether procedure (Noether, 1918). However, for classical
systems of electromagnetic field, the canonical EMTs are neither symmetric
with respect to tensor indices nor electromagnetic gauge invariant. Gauge
dependence is un-physical, and non-symmetric EMT is not consistent with
general relativity. In the present study, we will call a EMT symmetric if it
is symmetric with respect to tensor indices, and gauge symmetric if it is
gauge invariant. To date, much effort has been focused on symmetrizing the
EMTs (with respect to tensor indices), while constructing gauge-symmetric EMTs
is oftentimes a challenging task for general systems (Babak and Grishchuk,
1999; Gratus  _et al._ , 2012; Arminjon, 2016; Inglis and Jarvis, 2016;
Jiménez  _et al._ , 2018; Ilin and Paston, 2020).

The first method for symmetrizing EMTs was discovered by Belinfante
(Belinfante, 1939, 1940) and Rosenfeld (Rosenfeld, 1940), who added a
divergence-free tensor $\partial_{\sigma}\mathcal{S}^{\sigma\mu\nu}$ to obtain
a symmetric EMT, i.e.,

$\displaystyle
T_{\mathrm{BR}}^{\mu\nu}=T_{\text{N}}^{\mu\nu}+\partial_{\sigma}\mathcal{S}^{\sigma\mu\nu},$
(1)
$\displaystyle\partial_{\mu}\partial_{\sigma}\mathcal{S}^{\sigma\mu\nu}=0.$
(2)

Here, $T_{\mathrm{BR}}^{\mu\nu}$ is Belinfante-Rosenfeld (BR) EMT, and
$\mathcal{S}^{\sigma\mu\nu}$ is known as BR super-potential that depends on
the angular momentum tensor and is anti-symmetric with respect to $\sigma$ and
$\mu$ [see Eq. (52)]. General relativity suggests another method to generate
symmetric EMTs (Hawking and Ellis, 1973; Landau and lifshitz, 1975), which was
modified by Gotay and Marsden, who employed constraints to define symmetric
EMTs (Gotay and Marsden, 1992; Lopez  _et al._ , 2007). The relations between
these three types of symmetric EMTs have been discussed in the literature
(Zhang, 2005; Ilin and Paston, 2019; Baker  _et al._ , 2021).

In many systems, including the standard Maxwell system (6), the symmetrization
of $T_{\text{N}}^{\mu\nu}$ also renders it gauge symmetric. But for general
electromagnetic field theories with high-order field derivations, symmetry
with respect to tensor indices in general does not imply gauge symmetry and
vice versa. High-order electromagnetic field theories appear in the study of
gyrokinetic systems (Qin, 2005; Qin  _et al._ , 2007; Fan  _et al._ , 2020)
for magnetized plasmas and the Podolsky system (Bopp, 1940; Podolsky, 1942)
for the radiation reaction of classical charged particles. In the present
study, we propose a new method to gauge-symmetrize the canonical EMTs
$T_{\text{N}}^{\mu\nu}$ in general electromagnetic field theories with high-
order field derivations. Our method removes the gauge dependence, but does not
necessarily symmetrize the EMTs. In applications that don’t involve general
relativity, gauge-symmetrized EMTs are adequate.

We first reformulate the equation of motion for the field by the variational
principle with respect to the Faraday tensor $F_{\mu\nu}$, instead of the
4-potential $A^{\mu}$ as in the standard field theory. The Euler-Lagrange (EL)
equation is cast into an explicitly gauge-symmetric form. The canonical EMT is
then separated into a gauge-invariant part and a gauge-dependent part, the
later of which contains the anti-symmetric electric displacement tensor
$\mathcal{D}^{\mu\nu}$. We define a super-potential
$\mathcal{F}^{\sigma\mu\nu}\equiv\mathcal{D}^{\sigma\mu}A^{\nu}/4\pi$, called
displacement-potential tensor, whose divergence with respect to the first
index is divergence free with respect to the second index and removes the
gauge dependence in the canonical EMT. It is simpler to calculate the
displacement-potential tensor $\mathcal{F}^{\sigma\mu\nu}$ than the BR super-
potential $\mathcal{S}^{\sigma\mu\nu}$, and the former only gauge-symmetrizes
the canonical EMT without render it symmetric (with respect to tensor
indices). For first-order electromagnetic field theories, such as the standard
Maxwell system, $\mathcal{F}^{\sigma\mu\nu}$ reduces to the familiar BR super-
potential $\mathcal{S}^{\sigma\mu\nu}$, and the method developed here can be
used as a simpler procedure to calculate $\mathcal{S}^{\sigma\mu\nu}$ without
employing the angular momentum tensor in 4D spacetime. In addition, when the
electromagnetic system is coupled with classical charged particles, we find
that the method is effective as well, even though the Lagrangian density is
not gauge symmetric in general.

The paper is organized as follows. In Sec. II, we describe the gauge-
symmetrization method for the EMT in a general high-order electromagnetic
field theory, and highlight the difference in comparison with the BR method
using the example of the Podolsky system (Bopp, 1940; Podolsky, 1942). Section
III shows how the gauge-symmetrization method for the EMT works when the
electromagnetic system is coupled with classical charged particles.

##  II explicitly gauge-symmetric conservation laws for high-order
electromagnetic systems

###  II.1 Explicitly gauge-symmetric Euler-Lagrange equation

The Lagrangian density of a general electromagnetic system is written as

$\mathcal{L}_{F}=\mathcal{L}_{F}\left(x^{\mu},DA_{\mu},\cdots,D^{\left(n+1\right)}A_{\mu}\right),$
(3)

where $A=\left(\varphi,-\bm{A}\right)$ is the 4-potential and
$D=\left(\left(1/c\right)\partial_{t},\bm{\nabla}\right)$ is the derivative
operator over spacetime. The EL equation of the Lagrangian density is

$E_{A}^{\mu}\left(\mathcal{L}_{F}\right)=0,$ (4)

where the Euler operator $E_{A}^{\mu}$ of $A_{\mu}$ is defined by

$E_{A}^{\mu}\equiv\sum_{i=1}^{n+1}\left(-1\right)^{i}D_{\mu_{1}}\cdots
D_{\mu_{i}}\frac{\partial}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{i}}A_{\mu}\right)}.$
(5)

Note that the Lagrangian density $\mathcal{L}_{F}$ depends on derivatives of
$A$ with respect to the spacetime coordinates up to the $(n+1)$-th order. It
includes the standard Maxwell system, i.e.,

$\mathcal{L}_{F}=-\frac{1}{16\pi}F_{\mu\nu}F^{\mu\nu},$ (6)

as a special case, where

$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}.$ (7)

is the Faraday tensor. In Eq. (6), $\mathcal{L}_{F}$ depends only on first-
order derivatives of $A$. High-order electromagnetic field theories appear in
the study of gyrokinetic systems (Qin, 2005; Qin  _et al._ , 2007; Fan  _et
al._ , 2020) for magnetized plasmas and radiation reaction for classical
charged particles (Bopp, 1940; Podolsky, 1942). Physics requires that the EL
equation (4) is gauge symmetric, i.e, invariant under the gauge transformation
$A_{\mu}\mapsto A_{\mu}+\partial_{\mu}f$. In the present study, we assume that
$\mathcal{L}_{F}$ is explicitly gauge symmetric in the form of

$\mathcal{L}_{F}=\mathcal{L}_{F}\left(x^{\mu},F_{\mu\nu},DF_{\mu\nu},\cdots,D^{\left(n\right)}F_{\mu\nu}\right).$
(8)

From the variational principle,
$\delta\mathcal{A}=\delta\int\mathcal{L}_{F}d^{4}x=0$, we have

$\displaystyle
0=\delta\int\mathcal{L}_{F}d^{4}x=\int\left\\{E_{F}^{\mu\nu}\left(\mathcal{L}_{F}\right)\delta
F_{\mu\nu}\right\\}d^{4}x=\int\left\\{E_{F}^{\mu\nu}\left(\mathcal{L}_{F}\right)\left(\partial_{\mu}\delta
A_{\nu}-\partial_{\nu}\delta A_{\mu}\right)\right\\}d^{4}x$
$\displaystyle=-\int\partial_{\mu}\left\\{2E_{F}^{\left[\mu\nu\right]}\left(\mathcal{L}_{F}\right)\right\\}\delta
A_{\nu}d^{4}x,$ (9)

where the boundary term has been dropped, and $E_{F}^{\mu\nu}$ denotes the
Euler operator for the Faraday tensor $F_{\mu\nu}$ defined by

$E_{F}^{\mu\nu}\left(\mathcal{L}_{F}\right)=\frac{\partial\mathcal{L}_{F}}{\partial
F_{\mu\nu}}+\sum_{i=1}^{n}\left(-1\right)^{i}D_{\mu_{1}}\cdots
D_{\mu_{i}}\frac{\partial\mathcal{L}_{F}}{\partial\partial_{\mu_{1}}\cdots\partial_{\mu_{i}}F_{\mu\nu}}\,.$
(10)

In Eq. (9), superscript $\left[\mu\nu\right]$ represents anti-symmetrization
with respect to $\mu$ and $\nu$, i.e.,

$E_{F}^{\left[\mu\nu\right]}\left(\mathcal{L}_{F}\right)\equiv\frac{1}{2}\left[E_{F}^{\mu\nu}\left(\mathcal{L}_{F}\right)-E_{F}^{\nu\mu}\left(\mathcal{L}_{F}\right)\right].$
(11)

Due to the arbitrariness of $\delta A_{\nu}$ in Eq. (9), the equation of
motion for the system is

$\partial_{\mu}\mathcal{D}^{\mu\nu}=0,$ (12)

where

$\mathcal{D}^{\mu\nu}\equiv-8\pi
E_{F}^{\left[\mu\nu\right]}\left(\mathcal{L}_{F}\right)$ (13)

is the electric displacement tensor.

In Sec. III, we will consider electromagnetic systems coupled with charged
particles, and the Lagrangian density $\mathcal{L}$ will depend on 4-potential
$A_{\mu}$ i.e.,

$\mathcal{L}=\mathcal{L}\left(x^{\mu},A_{\mu},F_{\mu\nu},\cdots,D^{\left(n\right)}F_{\mu\nu}\right).$

Equation (12) then becomes

$\partial_{\mu}\mathcal{D}^{\mu\nu}=\frac{4\pi}{c}J_{f}^{\nu},$ (14)

where

$J_{f}^{\nu}\equiv-c\frac{\partial\mathcal{L}_{F}}{\partial A_{\nu}}$ (15)

is the free 4-current.

For the standard Maxwell system (6) without free 4-current, the electric
displacement tensor is the Faraday tensor, i.e.,
$\mathcal{D}^{\mu\nu}=F^{\mu\nu}$, and Eq. (12) reduces to Maxwell’s equation

$\partial_{\mu}F^{\mu\nu}=0.$ (16)

###  II.2 Infinitesimal criterion of symmetry and conservation laws

A continuous symmetry of the action $\mathcal{A}$ is a group of transformation

$\left(x^{\mu},A^{\nu}\right)\mapsto\left(\tilde{x}^{\mu},\tilde{A}^{\nu}\right)=g_{\epsilon}\cdot\left(x^{\mu},A^{\nu}\right),$
(17)

such that

$\int\mathcal{L}_{F}\left(\tilde{x}^{\mu},\tilde{F}_{\mu\nu},\tilde{D}\tilde{F}_{\mu\nu},\cdots,\tilde{D}^{\left(n\right)}\tilde{F}_{\mu\nu}\right)d^{4}\tilde{x}=\int\mathcal{L}_{F}\left(x^{\mu},F_{\mu\nu},DF_{\mu\nu},\cdots,D^{\left(n\right)}F_{\mu\nu}\right)d^{4}x,$
(18)

where $g_{\epsilon}$ constitutes a continuous group of the transformations
parameterized by $\epsilon$ (Olver, 1993). The infinitesimal generator of the
transformation group is

$\bm{v}\coloneqq\frac{d}{d\epsilon}|_{0}g_{\epsilon}\cdot\left(x^{\mu},A^{\nu}\right)=\xi^{\mu}\frac{\partial}{\partial
x^{\mu}}+\phi_{\mu}\frac{\partial}{\partial A_{\mu}}.$ (19)

By rewriting the symmetry condition (18) as

$\frac{d}{d\epsilon}|_{0}\int\mathcal{L}_{F}\left(\tilde{x}^{\mu},\tilde{F}_{\mu\nu},\tilde{D}\tilde{F}_{\mu\nu},\cdots,\tilde{D}^{\left(n\right)}\tilde{F}_{\mu\nu}\right)d^{4}\tilde{x}=0,$
(20)

we can derive the following infinitesimal version of the symmetry condition,

$\text{pr}^{\left(n+1\right)}\bm{v}\left(\mathcal{L}_{F}\right)+\mathcal{L}_{F}D_{\mu}\xi^{\mu}=0,$
(21)

where $\text{pr}^{\left(n+1\right)}\bm{v}$ is the prolongation of $\bm{v}$.
The standard prolongation formula for $\text{pr}^{\left(n+1\right)}\bm{v}$ can
be found in Ref. (Olver, 1993). In the present study, we rewrite the
prolongation formula with respect to $F_{\mu\nu}$, instead of $A_{\mu}$, as

$\displaystyle\text{pr}^{\left(n+1\right)}\bm{v}=\bm{v}+\left[G_{\sigma\rho}+\xi^{\sigma}D_{\sigma}F_{\sigma\rho}\right]\frac{\partial\mathcal{L}_{F}}{\partial
F_{\sigma\rho}}$ $\displaystyle+\sum_{i=1}^{n}\left[D_{\mu_{1}}\cdots
D_{\mu_{i}}G_{\sigma\rho}+\xi^{\sigma}D_{\sigma}D_{\mu_{1}}\cdots
D_{\mu_{i}}F_{\sigma\rho}\right]\frac{\partial\mathcal{L}_{F}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{i}}F_{\sigma\rho}\right)},$
(22)

where

$G_{\sigma\rho}=\partial_{\sigma}Q_{\rho}-\partial_{\rho}Q_{\sigma}\equiv
2\partial_{[\sigma}Q_{\rho]}$ (23)

and

$Q_{\nu}=\phi_{\nu}-\xi^{\sigma}D_{\sigma}A_{\nu},$ (24)

is a characteristic of the Lie algebra.

Combining Eqs. (12) and (21) generates the conservation law corresponding to
the symmetry,

$D_{\mu}\left\\{\mathcal{L}_{F}\xi^{\mu}-\frac{1}{4\pi}\mathcal{D}^{\mu\nu}\left(\mathcal{L}_{F}\right)Q_{\nu}+\mathbb{P}^{\mu}\right\\}=0,$
(25)

where

$\mathbb{P}^{\mu}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j+1}\left(D_{\mu_{j+1}}\cdots
D_{\mu_{i}}G_{\sigma\rho}\right)\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}_{F}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\sigma\rho}\right)}\right].$
(26)

The conservation law given by Eq. (25) is not gauge-symmetric in general.

###  II.3 Gauge-symmetrization of the canonical EMT

Now we assume the high-order electromagnetic field theory admits the spacetime
translation symmetry, i.e.,

$\frac{\partial\mathcal{L}_{F}}{\partial x^{\mu}}=0,$ (27)

and derive the corresponding energy-momentum conservation law. Because of Eq.
(27), the action is invariant under the spacetime translation

$\left(x^{\mu},A_{\nu}\right)\mapsto\left(x^{\mu}+\epsilon
X_{0}^{\mu},A_{\nu}\right),$ (28)

where $X_{0}^{\mu}$ is 4D constant vector field. The infinitesimal generator
$\bm{v}$, characteristic $Q^{\nu}$, and $G_{\sigma\rho}$ in Eq. (23) are

$\displaystyle\bm{v}=X_{0}^{\mu}\frac{\partial}{\partial x^{\mu}},$ (29)
$\displaystyle Q^{\nu}=-X_{0}^{\nu}\partial_{\nu}A_{\sigma},$ (30)
$\displaystyle G_{\sigma\rho}=-X_{0}^{\nu}\partial_{\nu}F_{\sigma\rho}.$ (31)

The Lagrangian density satisfies the infinitesimal criterion because

$X_{0}^{\mu}\frac{\partial\mathcal{L}_{F}}{\partial x^{\mu}}=0,$ (32)

which implies a conservation law. Substituting Eqs. (29)-(31) into to Eq.
(25), we obtain the canonical energy-momentum conservation law according the
standard Noether procedure,

$\displaystyle D_{\mu}T_{\text{N}}^{\mu\nu}=0,$ (33) $\displaystyle
T_{\text{N}}^{\mu\nu}=\mathcal{L}_{F}\eta^{\mu\nu}+\frac{1}{4\pi}\mathcal{D}^{\mu\sigma}\partial^{\nu}A_{\sigma}-\Sigma^{\mu\nu},$
(34)
$\displaystyle\Sigma^{\mu\nu}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j+1}\left(D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\partial^{\nu}F_{\sigma\rho}\right)\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}_{F}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\sigma\rho}\right)}\right].$
(35)

In Eq. (33), $T_{\text{N}}^{\mu\nu}$ is the canonical EMT derived from the
standard Noether procedure.

Obviously, $T_{\text{N}}^{\mu\nu}$ depends on the gauge as expected. In the
expression of $T_{\text{N}}^{\mu\nu}$ given by Eq. (34), the gauge dependence
comes from the second term, and the first and third terms are gauge symmetric.
Now we show how to gauge-symmetrize $T_{\text{N}}^{\mu\nu}$. Note that because
electric displacement tensor $\mathcal{D}^{\sigma\mu}$ is anti-symmetric, the
following equations hold,

$\displaystyle D_{\mu}\left(D_{\sigma}\mathcal{F}^{\sigma\mu\nu}\right)=0,$
(36)
$\displaystyle\mathcal{F}^{\sigma\mu\nu}\equiv\frac{1}{4\pi}\mathcal{D}^{\sigma\mu}A^{\nu}.$
(37)

Here, $\mathcal{F}^{\sigma\mu\nu}$ is a super-potential that is anti-symmetric
with respect to the first two indices. For easy reference, we will call
$\mathcal{F}^{\sigma\mu\nu}$ displacement-potential tensor. The divergence of
$\mathcal{F}^{\sigma\mu\nu}$ defines a divergence-free tensor, i.e.,

$T_{0}^{\mu\nu}\equiv
D_{\sigma}\mathcal{F}^{\sigma\mu\nu}=-\frac{1}{4\pi}\mathcal{D}^{\mu\sigma}\partial_{\sigma}A^{\nu},$
(38)

where used is made of Eq. (12). When $T_{0}^{\mu\nu}$ is added to
$T_{\text{N}}^{\mu\nu},$ the gauge dependence is removed, i.e.,

$\displaystyle D_{\mu}T_{\text{GS}}^{\mu\nu}=0,$ (39) $\displaystyle
T_{\text{GS}}^{\mu\nu}\equiv
T_{\text{N}}^{\mu\nu}+T_{0}^{\mu\nu}=\mathcal{L}_{F}\eta^{\mu\nu}+\frac{1}{4\pi}\mathcal{D}^{\mu\sigma}F_{\;\sigma}^{\nu}-\Sigma^{\mu\nu},$
(40)

where $T_{\text{GS}}^{\mu\nu}$ is the gauge-symmetric EMT.

It is worthwhile to mention that we derived the gauge-symmetrized EMT
$T_{\text{GS}}^{\mu\nu}$ from the expression of $T_{\text{N}}^{\mu\nu}$ in Eq.
(34), which is calculated from the prolongation with respect to $F_{\mu\nu}$.
On the other hand, had we started from Eq. (3) and calculated the EMT from the
prolongation with respect to $A_{\mu}$, we would have obtained a canonical EMT
in the form of

$T_{\text{N}}^{\mu\nu}=\mathcal{L}_{F}\eta^{\mu\nu}-\hat{\Sigma}^{\mu\nu},$
(41)

where

$\hat{\Sigma}^{\mu\nu}=\sum_{i=1}^{n+1}\sum_{j=1}^{i}\left(-1\right)^{j+1}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left(\partial_{\nu}A_{\sigma}\right)\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}_{F}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}A_{\sigma}\right)}\right].$
(42)

However, different from the situation in Eq. (34), every term in Eq. (42) is
gauge dependent, making the gauge symmetrization difficult, if not impossible.

For the standard Maxwell electromagnetic system specified by Eq. (6), the
electric displacement tensor reduces to the Faraday tensor $F^{\sigma\mu}$,
and the displacement-potential tensor $\mathcal{F}^{\sigma\mu\nu}$ reduces to
$F^{\sigma\mu}A^{\nu}/4\pi$, coinciding with the tensor used by Blaschke et
al. for the U$\left(1\right)$ gauge theory (Blaschke  _et al._ , 2016).

###  II.4 Comparison with the BR method

As described above, the method proposed in the present study employs
displacement-potential tensor $\mathcal{F}^{\sigma\mu\nu}$ to gauge-symmetrize
the EMT, while the BR method use the super-potential
$\mathcal{S}^{\sigma\mu\nu}$ to symmetrize the EMT. In this subsection, we
discuss the difference between the displacement-potential tensor
$\mathcal{F}^{\sigma\mu\nu}$ in Eq. (36) and the BR super-potential
$\mathcal{S}^{\sigma\mu\nu}$. To calculate $\mathcal{S}^{\sigma\mu\nu}$, we
need to first derive the 4D angular momentum conservation laws generated by
the Lorentz symmetry. Assume that system is invariant under rotational
transformation in 4D spacetime

$\left(x^{\mu},A^{\nu}\right)\mapsto\left(\tilde{x}^{\mu},\tilde{A}^{\nu}\right)=\left(\Lambda_{\epsilon}^{\mu\sigma}x_{\sigma},\Lambda_{\epsilon}^{\nu\sigma}A_{\sigma}\right),$
(43)

where $\left\\{\Lambda_{\epsilon}^{\mu\sigma}\right\\}$ is one-parameter
subgroup of the Lorentz group. The infinitesimal generator $\bm{v}$, the
characteristic $Q_{\rho}$, and the term $G_{s\rho}$ are calculated
respectively by Eqs. (19), (24), and (23) as

$\displaystyle\bm{v}=\frac{d}{d\epsilon}|_{0}\left(\Lambda_{\epsilon}^{\mu\sigma}x_{\sigma},\Lambda_{\epsilon}^{\nu\sigma}A_{\sigma}\right)=\left(\Omega^{\mu\sigma}x_{\sigma},\Omega^{\mu\sigma}A_{\sigma}\right),$
(44) $\displaystyle
Q_{\rho}=\phi_{\rho}-\xi^{\alpha}D_{\alpha}A_{\rho}=\Omega_{\rho\alpha}A^{\alpha}-\Omega_{\alpha\beta}x^{\beta}D^{\alpha}A_{\rho},$
(45) $\displaystyle
G_{s\rho}=\Omega_{\rho\alpha}F_{s}^{\;\alpha}-\Omega_{s\alpha}F_{\rho}^{\;\alpha}-\Omega_{\alpha\beta}x^{\beta}\partial^{\alpha}F_{s\rho},$
(46)

where the anti-symmetric tensor
$\Omega^{\mu\sigma}=\left[d\Lambda_{\epsilon}^{\mu\sigma}/d\epsilon\right]_{0}$
is the Lie algebra element of the Lorentz group. Substituting Eqs. (44)-(46)
into Eq. (25), we obtain the angular momentum conservation law in 4D
spacetime,

$\Omega_{\nu\sigma}D_{\mu}\left\\{x^{\sigma}T_{\text{N}}^{\mu\nu}+2E_{F}^{\left[\mu\nu\right]}\left(\mathcal{L}_{F}\right)A^{\sigma}+L^{\mu\nu\sigma}\right\\}=0,$
(47)

which can be rewritten as

$D_{\mu}\left\\{\left[x^{\sigma}T_{\text{N}}^{\mu\nu}-x^{\nu}T_{\text{N}}^{\mu\sigma}\right]+S^{\mu\nu\sigma}\right\\}=0.$
(48)

In above equations,

$\displaystyle
L^{\mu\nu\sigma}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j+1}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left(F_{\rho}^{\;\sigma}\right)D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\times$
$\displaystyle\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\rho\nu}\right)}-\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\nu\rho}\right)}\right]$
$\displaystyle-\left[\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j+1}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left(x^{\sigma}\partial^{\nu}F_{s\rho}\right)D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{s\rho}\right)}-x^{\sigma}\Sigma^{\mu\nu}\right]$
(49)
$\displaystyle\frac{1}{2}S^{\mu\nu\sigma}=E_{F}^{\left[\mu\nu\right]}\left(\mathcal{L}_{F}\right)A^{\sigma}-E_{F}^{\left[\mu\sigma\right]}\left(\mathcal{L}_{F}\right)A^{\nu}+\Delta^{\mu\nu\sigma},$
(50) $\displaystyle\Delta^{\mu\nu\sigma}\equiv 2L^{\mu\left[\nu\sigma\right]}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j+1}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left\\{\left(F_{\rho}^{\;\sigma}\right)D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\rho\nu}\right)}\right.\right.$
$\displaystyle\left.\vphantom{\left(F_{\rho}^{\;\sigma}\right)}-\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\nu\rho}\right)}\right]-\left(F_{\rho}^{\;\nu}\right)D_{\mu_{1}}\cdots
D_{\mu_{j-1}}$
$\displaystyle\left.\vphantom{\left(F_{\rho}^{\;\sigma}\right)}\times\left[\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\rho\sigma}\right)}-\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{\sigma\rho}\right)}\right]\right\\}$
$\displaystyle-\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j+1}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left(x^{\sigma}\partial^{\nu}F_{s\rho}-x^{\nu}\partial^{\sigma}F_{s\rho}\right)D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\times$
$\displaystyle\frac{\partial\mathcal{L}}{\partial\left(\partial_{\mu_{1}}\cdots\partial_{\mu_{j-1}}\partial_{\mu}\partial_{\mu_{j+1}}\cdots\partial_{\mu_{i}}F_{s\rho}\right)}+x^{\sigma}\Sigma^{\mu\nu}-x^{\nu}\Sigma^{\mu\sigma},$
(51)

where the superscript $\left[\mu\nu\right]$ denotes anti-symmetrization with
respect to $\mu$ and $\nu$.

The anti-symmetric BR super-potentail $\mathcal{S}^{\sigma\mu\nu}$ is defined
from the tensor $S^{\sigma\mu\nu}$ in Eq. (48) as (Belinfante, 1939, 1940;
Rosenfeld, 1940)

$\displaystyle\mathcal{S}^{\sigma\mu\nu}$
$\displaystyle\equiv\frac{1}{2}\left[S^{\sigma\nu\mu}-S^{\mu\nu\sigma}-S^{\nu\mu\sigma}\right].$
(52)

It is clear from Eqs. (49)-(51) that $\mathcal{S}^{\sigma\mu\nu}$ and
$\mathcal{F}^{\sigma\mu\nu}$ are related as follows,

$\mathcal{S}^{\sigma\mu\nu}=\mathcal{F}^{\sigma\mu\nu}+\frac{1}{2}\left[\Delta^{\sigma\nu\mu}-\Delta^{\mu\nu\sigma}-\Delta^{\nu\mu\sigma}\right].$
(53)

Equation (52) shows that in general $\mathcal{F}^{\sigma\mu\nu}$ is different
from $\mathcal{S}^{\sigma\mu\nu}$ when $\Delta^{\sigma\nu\mu}$ is non-
vanishing. For a first-order field theory, such as the standard Maxwell system
(6), $n=1$ and the last three terms vanish such that
$\mathcal{S}^{\sigma\mu\nu}=\mathcal{F}^{\sigma\mu\nu}$. In this situation,
adding $T_{0}^{\mu\nu}\equiv D_{\sigma}\mathcal{F}^{\sigma\mu\nu}$ to
$T_{\text{N}}^{\mu\nu}$ will render it both symmetric and gauge-symmetric, and
the method developed here can be used as a simpler procedure to calculate the
BR super-potential $\mathcal{S}^{\sigma\mu\nu}$ without the necessity to
calculate the angular momentum tensor in 4D spacetime.

As an example of high-order electromagnetic field theory, we consider the
Podolsky system (Bopp, 1940; Podolsky, 1942), which was proposed to study the
radiation reaction of classical charged particles. The Podolsky Lagrangian
density is

$\mathcal{L}_{\text{Po}}=\frac{1}{8\pi}\left\\{\bm{E}^{2}-\bm{B}^{2}+a^{2}\left[\left(\bm{\nabla}\cdot\bm{E}\right)^{2}-\left(\bm{\nabla}\times\bm{B}-\frac{1}{c}\partial_{t}\bm{E}\right)^{2}\right]\right\\}$
(54)

or in a manifestly covariant form

$\mathcal{L}_{\text{Po}}=-\frac{1}{16\pi}F_{\sigma\rho}F^{\sigma\rho}-\frac{a^{2}}{8\pi}\partial_{\sigma}F^{\sigma\lambda}\partial^{\rho}F_{\rho\lambda}.$
(55)

We substitute the Lagrangian density (55) into Eq. (34) to obtain the
canonical EMT

$\displaystyle 4\pi
T_{\text{N}}^{\mu\nu}=\left(-\frac{1}{4}F_{\sigma\rho}F^{\sigma\rho}-\frac{a^{2}}{2}\partial_{\sigma}F^{\sigma\lambda}\partial^{\rho}F_{\rho\lambda}\right)\eta^{\mu\nu}$
$\displaystyle+\left[F^{\mu\sigma}-a^{2}\left(\partial^{\mu}\partial_{\lambda}F^{\lambda\sigma}-\partial^{\sigma}\partial_{\lambda}F^{\lambda\mu}\right)\right]\partial^{\nu}A_{\sigma}+a^{2}\left(\partial^{\nu}F_{\;\rho}^{\mu}\right)\left(\partial_{\sigma}F^{\sigma\rho}\right),$
(56)

where use is made of the following equations,

$\displaystyle\frac{\partial}{\partial\left(\partial_{\sigma}F_{\mu\nu}\right)}\left[\partial_{\alpha}F^{\alpha\lambda}\partial^{\rho}F_{\rho\lambda}\right]=2\eta^{\sigma\mu}\partial_{\lambda}F^{\lambda\nu},$
(57) $\displaystyle
D_{\sigma}\frac{\partial}{\partial\left(\partial_{\sigma}F_{\mu\nu}\right)}\left[\partial_{\alpha}F^{\alpha\lambda}\partial^{\rho}F_{\rho\lambda}\right]=2\partial^{\mu}\partial_{\sigma}F^{\sigma\nu},$
(58) $\displaystyle
E_{F}^{\mu\sigma}=\frac{\partial\mathcal{L}_{\text{Po}}}{\partial
F_{\mu\sigma}}-D_{\rho}\frac{\partial\mathcal{L}_{\text{Po}}}{\partial\partial_{\rho}F_{\mu\sigma}}=-\frac{1}{8\pi}F^{\mu\sigma}+\frac{a^{2}}{4\pi}\partial^{\mu}\partial_{\lambda}F^{\lambda\sigma},$
(59) $\displaystyle
2E_{F}^{\left[\mu\sigma\right]}=-\frac{1}{4\pi}F^{\mu\sigma}+\frac{a^{2}}{4\pi}\left(\partial^{\mu}\partial_{\lambda}F^{\lambda\sigma}-\partial^{\sigma}\partial_{\lambda}F^{\lambda\mu}\right),$
(60)
$\displaystyle\Sigma^{\mu\nu}=\left(\partial^{\nu}F_{\sigma\rho}\right)\frac{\partial\mathcal{L}_{\text{Po}}}{\partial\left(\partial_{\mu}F_{\sigma\rho}\right)}=-\frac{a^{2}}{4\pi}\partial^{\nu}F_{\sigma\rho}\left[\eta^{\mu\sigma}\partial_{\lambda}F^{\lambda\rho}\right]=-\frac{a^{2}}{4\pi}\left(\partial^{\nu}F_{\;\rho}^{\mu}\right)\left(\partial_{\sigma}F^{\sigma\rho}\right).$
(61)

The displacement-potential tensor is

$\mathcal{F}^{\mu\nu\sigma}\equiv\frac{1}{4\pi}\mathcal{D}^{\sigma\mu}A^{\nu}=\frac{1}{4\pi}\left[-F^{\mu\sigma}+a^{2}\left(\partial^{\mu}\partial_{\lambda}F^{\lambda\sigma}-\partial^{\sigma}\partial_{\lambda}F^{\lambda\mu}\right)\right]A^{\nu},$
(62)

and

$4\pi\partial_{\sigma}\mathcal{F}^{\mu\nu\sigma}=\left[-F^{\mu\sigma}+a^{2}\left(\partial^{\mu}\partial_{\lambda}F^{\lambda\sigma}-\partial^{\sigma}\partial_{\lambda}F^{\lambda\mu}\right)\right]\partial_{\sigma}A^{\nu}.$
(63)

Adding Eq. (63) to Eq. (56), we obtain the gauge-symmetric EMT,

$\displaystyle 4\pi
T_{\text{GS}}^{\mu\nu}=\left[F^{\mu\sigma}F_{\nu\sigma}-\frac{1}{4}\left(F_{\sigma\rho}F^{\sigma\rho}\right)\eta^{\mu\nu}\right]-\frac{a^{2}}{2}\left(\partial_{\sigma}F^{\sigma\lambda}\partial^{\rho}F_{\rho\lambda}\right)\eta^{\mu\nu}$
$\displaystyle-a^{2}F_{\;\sigma}^{\nu}\left(\partial^{\mu}\partial_{\rho}F^{\rho\sigma}\right)+a^{2}F_{\;\sigma}^{\nu}\left(\partial^{\sigma}\partial_{\rho}F^{\rho\mu}\right)+a^{2}\left(\partial^{\nu}F_{\;\rho}^{\mu}\right)\left(\partial_{\sigma}F^{\sigma\rho}\right).$
(64)

It is easy to see that $T_{\text{GS}}^{\mu\nu}$ for the Podolsky system is not
symmetric, i.e., $T_{\text{GS}}^{\mu\nu}\neq T_{\text{GS}}^{\nu\mu}$.

To calculate the BR EMT $T_{\text{BR}}^{\mu\nu}$ for the Podolsky system, we
evaluate the $\Delta^{\mu\nu\sigma}$ term in Eq. (52). Using Eq. (51), we have

$\displaystyle\Delta^{\mu\nu\sigma}$
$\displaystyle=F_{\rho}^{\;\sigma}\left[\frac{\partial\mathcal{L}_{\text{Po}}}{\partial\left(\partial_{\mu}F_{\rho\nu}\right)}-\frac{\partial\mathcal{L}_{\text{Po}}}{\partial\left(\partial_{\mu}F_{\nu\rho}\right)}\right]-F_{\rho}^{\;\nu}\left[\frac{\partial\mathcal{L}_{\text{Po}}}{\partial\left(\partial_{\mu}F_{\rho\sigma}\right)}-\frac{\partial\mathcal{L}_{\text{Po}}}{\partial\left(\partial_{\mu}F_{\sigma\rho}\right)}\right]$
$\displaystyle=-\frac{a^{2}}{4\pi}\left[F^{\mu\sigma}\partial_{\rho}F^{\rho\nu}-\eta^{\mu\nu}\left(F_{\rho}^{\;\sigma}\partial_{\lambda}F^{\lambda\rho}\right)-F^{\mu\nu}\partial_{\rho}F^{\rho\sigma}+\eta^{\mu\sigma}F_{\rho}^{\;\nu}\left(\partial_{\lambda}F^{\lambda\rho}\right)\right].$
(65)

Substituting Eqs. (52) and (65) into Eq. (1), we obtain the BR EMT as

$\displaystyle 4\pi T_{\text{BR}}^{\mu\nu}=4\pi
T_{\text{GS}}^{\mu\nu}+2\pi\left[\Delta^{\sigma\nu\mu}-\Delta^{\mu\nu\sigma}-\Delta^{\nu\mu\sigma}\right]$
$\displaystyle=\left[F^{\mu\sigma}F_{\nu\sigma}-\frac{1}{4}\left(F_{\sigma\rho}F^{\sigma\rho}\right)\eta^{\mu\nu}\right]+\frac{a^{2}}{2}\left[\left(\partial_{\sigma}F^{\sigma\rho}\right)\left(\partial^{\lambda}F_{\lambda\rho}\right)-2F_{\rho}^{\;\sigma}\left(\partial_{\sigma}\partial_{\lambda}F^{\lambda\rho}\right)\right]\eta^{\mu\nu}$
$\displaystyle+a^{2}\left[F^{\nu\sigma}\left(\partial_{\sigma}\partial_{\rho}F^{\rho\mu}\right)+F^{\mu\sigma}\left(\partial_{\sigma}\partial_{\rho}F^{\rho\nu}\right)-\left(\partial_{\sigma}F^{\sigma\mu}\right)\left(\partial_{\rho}F^{\rho\nu}\right)\right.$
$\displaystyle\left.\vphantom{\left(\partial_{\rho}F^{\rho\nu}\right)}-F_{\;\sigma}^{\nu}\left(\partial^{\mu}\partial_{\rho}F^{\rho\sigma}\right)-F_{\;\sigma}^{\mu}\left(\partial^{\nu}\partial_{\rho}F^{\rho\sigma}\right)\right].$
(66)

It is easy to verify that $T_{\text{BR}}^{\mu\nu}$ for the Podolsky system is
both symmetric and gauge-symmetric.

##  III Gauge-symmetric EMTs for electromagnetic systems coupled with
classical charged particles

For self-consistent electromagnetic systems with free currents, the
electromagnetic fields are coupled with charged particles. In this section, we
apply the theory established in Sec. II to derive gauge-symmetric EMTs for
electromagnetic systems coupled with classical charged particles.

For practical applications, such as in the gyrokinetic theory (Qin, 2005; Qin
_et al._ , 2007; Fan  _et al._ , 2020) for magnetized plasmas, reduced
theoretical models are often adopted due to the intrinsic complexity of the
systems. The equations of motion for the systems are usually gauge invariant,
but the Lagrangian densities are not always specified by manifestly covariant
forms. For these systems, energy and momentum conservation laws need to be
derived separately. We demonstrate how the energy and momentum conservation
laws can be transformed into gauge-symmetric forms using the “3+1” form of Eq.
(36), i.e.,

$\frac{D}{Dt}\left\\{\frac{D}{D\bm{x}}\cdot\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi\right]\right\\}+\frac{D}{D\bm{x}}\cdot\left\\{\frac{D}{Dt}\left[-\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi\right]\right\\}=0$
(67)

and

$\frac{D}{Dt}\left\\{\frac{D}{D\bm{x}}\cdot\left[-\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}\right]\right\\}+\frac{D}{D\bm{x}}\cdot\left\\{\frac{D}{Dt}\left[\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}\right]\right\\}=0.$
(68)

###  III.1 Weak Euler-Lagrange equation and conservation law

The Lagrangian density of a generic classical electromagnetic field-charge
particle system assumes the form of

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{a}\mathcal{L}_{a}+\mathcal{L}_{F},$ (69)
$\displaystyle\mathcal{L}_{a}$ $\displaystyle=L_{a}\delta_{a},\quad
L_{a}=L_{a}\left(x^{\mu},\bm{X}_{a},\dot{\bm{X}}_{a};\varphi,\bm{A},\bm{E},\bm{B},D\bm{E},D\bm{B},\cdots,D^{n}\bm{E},D^{n}\bm{B}\right)$
(70)

where the subscript $a$ labels particles, $\bm{X}_{a}$ is its trajectory,
$\mathcal{L}_{a}$ is its Lagrangian density, and
$\delta_{a}\equiv\delta\left(\bm{x}-\bm{X}_{a}\right)$. Here, $\delta(x)$ is
the Dirac $\delta$-function.

In the “3+1” form, the equations of motion for the electromagnetic field are

$\displaystyle\nabla\cdot\bm{E}_{\bm{E}}\left(\mathcal{L}\right)=-\frac{\partial\mathcal{L}}{\partial\varphi},$
(71) $\displaystyle-\frac{1}{c}\frac{\partial}{\partial
t}\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]-\bm{\nabla}\times\left[\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]=\frac{\partial\mathcal{L}}{\partial\bm{A}}.$
(72)

In this study, it is assumed that the Lagrangian density $\mathcal{L}$ is
linear in terms of $\varphi$ and $\bm{A}$, and Eqs. (71) and (72) are thus
gauge symmetric. Specifically, we assume that $\mathcal{L}$ depends on
$\varphi$ and $\bm{A}$ only through the term
$-q_{a}\delta_{a}(\varphi+\bm{A}\cdot\dot{\bm{X}}_{a}/c)$, i.e., the
Lagrangian density can be written as

$\mathcal{L}=\sum_{a}\left[-\varphi+\frac{1}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\right]q_{a}\delta_{a}+\text{GSP}\left(\mathcal{L}\right),$
(73)

where “$\text{GSP}\left(\mathcal{L}\right)$” denotes the gauge-symmetric parts
of the Lagrangian density $\mathcal{L}$. The right hand side of Eqs. (71) and
(72) are the “3+1” form of Eq. (15), the free charge density $\rho_{f}$ and
current density $\bm{j}_{f}$, respectively. Using Eq. (73), we have

$\rho_{f}=-\frac{\partial\mathcal{L}}{\partial\varphi}=\sum_{a}q_{a}\delta_{a},\quad\bm{j}_{f}=c\frac{\partial\mathcal{L}}{\partial\bm{A}}=\sum_{a}q_{a}\dot{\bm{X}}_{a}\delta_{a}.$
(74)

The equation of motion for particles is also derived from the variational
principle. However, because particles and field reside on different manifolds,
the equation of motion for particles will be the weak EL equation (Qin  _et
al._ , 2014; Fan  _et al._ , 2018, 2019, 2020)

$\bm{E}_{\bm{X}_{a}}\left(\mathcal{L}\right)=\frac{D}{D\bm{x}}\cdot\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right),$
(75)

where $\bm{E}_{\bm{X}_{a}}$ is the Euler operator for the trajectory of the
$a$-th particle,

$\bm{E}_{\bm{X}_{a}}=\frac{\partial}{\partial\bm{X}_{a}}-\frac{d}{dt}\frac{\partial}{\partial\dot{\bm{X}_{a}}}.$
(76)

To derive a local conservation law from a symmetry, we need the infinitesimal
symmetry criterion for the Lagrangian density. A symmetry of the action
$\mathcal{A}\equiv\int\mathcal{L}dtd^{3}\bm{x}$ is defined by group
transformations

$\left(x^{\mu},\bm{X}_{a};\varphi,\bm{A}\right)\mapsto\left(\tilde{x}^{\mu},\tilde{\bm{X}}_{a};\tilde{\varphi},\tilde{\bm{A}}\right)=g_{\epsilon}\cdot\left(x^{\mu},\bm{X}_{a};\varphi,\bm{A}\right),$
(77)

such that

$\int\mathcal{L}\left(\tilde{x}^{\mu},\tilde{\bm{X}}_{a},\tilde{\bm{E}},\tilde{\bm{B}},\cdots,\tilde{D}^{n}\tilde{\bm{E}},\tilde{D}^{n}\tilde{\bm{B}}\right)d\tilde{t}d^{3}\tilde{\bm{x}}=\int\mathcal{L}\left(x^{\mu},\bm{X}_{a},\bm{E},\bm{B},\cdots,D^{\left(n\right)}\bm{E},D^{\left(n\right)}\bm{B}\right)dtd^{3}\bm{x}.$
(78)

The corresponding infinitesimal generator of (77) is

$\bm{v}=\xi^{\mu}\frac{\partial}{\partial
x^{\mu}}+\sum_{a}\bm{\theta}_{a}\cdot\frac{\partial}{\partial\bm{X}_{a}}+\phi_{0}\frac{\partial}{\partial\varphi}+\bm{\phi}_{\bm{A}}\cdot\frac{\partial}{\partial\bm{A}}.$
(79)

The infinitesimal criterion of the symmetry condition can be derived using the
same procedure in Sec. II.2,

$\text{pr}^{\left(n+1\right)}\bm{v}\left(\mathcal{L}\right)+\mathcal{L}D_{\mu}\xi^{\mu}=0.$
(80)

The prolongation of $\bm{v}$ now reads

$\displaystyle\text{pr}^{\left(n+1\right)}\bm{v}\left(\mathcal{L}\right)$
$\displaystyle=\bm{v}+\sum_{a}\left[\left(\dot{\bm{q}}_{a}+\xi^{t}\ddot{\bm{X}}_{a}\right)\cdot\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right]-\left[\bm{\nabla}Q_{0}+\xi^{\mu}D_{\mu}\left(\bm{\nabla}\varphi\right)\right]\cdot\frac{\partial\mathcal{L}}{\partial\bm{E}}$
$\displaystyle-\cdots-\sum_{i=1}^{n}\left[D_{\mu_{1}}\cdots
D_{\mu_{i}}\bm{\nabla}Q_{\alpha 0}+\xi^{\mu}D_{\mu}D_{\mu_{1}}\cdots
D_{\mu_{i}}\left(\bm{\nabla}\varphi_{\alpha}\right)\right]\cdot\frac{\partial\mathcal{L}}{\partial
D_{\mu_{1}}\cdots D_{\mu_{i}}\bm{E}}$
$\displaystyle-\left[D_{t}\bm{Q}_{\bm{A}}+\xi^{\mu}D_{\mu}\bm{A}_{,t}\right]\cdot\left(\frac{1}{c}\frac{\partial\mathcal{L}}{\partial\bm{E}}\right)$
$\displaystyle-\cdots-\sum_{i=1}^{n}\left[D_{\mu_{1}}\cdots
D_{\mu_{i}}D_{t}\bm{Q}_{\bm{A}}+\xi^{\mu}D_{\mu}D_{\mu_{1}}\cdots
D_{\mu_{i}}\bm{A}_{,t}\right]\cdot\frac{\partial\mathcal{L}}{\partial
D_{\mu_{1}}\cdots D_{\mu_{i}}\bm{E}}$
$\displaystyle+\left[\bm{\nabla}\bm{Q}_{\bm{A}}+\xi^{\mu}D_{\mu}\bm{\nabla}\bm{A}\right]:\left(\bm{\varepsilon}\cdot\frac{\partial\mathcal{L}}{\partial\bm{B}}\right)$
$\displaystyle+\cdots+\sum_{i=1}^{n}\left[D_{\mu_{1}}\cdots
D_{\mu_{i}}\bm{\nabla}\bm{Q}_{\bm{A}}+\xi^{\mu}D_{\mu}D_{\mu_{1}}\cdots
D_{\mu_{i}}\bm{\nabla}\bm{A}\right]:\left(\bm{\varepsilon}\cdot\frac{\partial\mathcal{L}}{\partial
D_{\mu_{1}}\cdots D_{\mu_{i}}\bm{B}}\right),$ (81)

where

$\bm{q}_{a}=\bm{\theta}_{a}-\xi^{t}\dot{\bm{X}}_{a}$ (82)

is another characteristic of $\bm{v}$ induced by particle’s trajectory. To
obtain the corresponding conservation law, we transform the infinitesimal
criterion into

$\displaystyle\partial_{t}\left[\mathcal{L}\xi^{t}-\frac{1}{c}\bm{Q}_{\bm{A}}\cdot\bm{E}_{\bm{E}}\left(\mathcal{L}\right)+\sum_{a}\left(\bm{q}_{a}\cdot\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)\right]+\bm{\nabla}\cdot\left[\mathcal{L}\bm{\kappa}-Q_{0}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)+\bm{Q}_{\bm{A}}\times\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]$
$\displaystyle+D_{\mu}\left[\mathcal{\mathbb{P}}_{1}^{\mu}+\mathcal{\mathbb{P}}_{2}^{\mu}\right]+\sum_{a}\left[\bm{q}_{a}\cdot\bm{E}_{\bm{X}_{a}}\left(\mathcal{L}\right)\right]+\left\\{\frac{\partial\mathcal{L}}{\partial\varphi}+\bm{\nabla}\cdot\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]\right\\}Q_{0}$
$\displaystyle+\left\\{\frac{\partial\mathcal{L}}{\partial\bm{A}}+\frac{1}{c}\partial_{t}\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]+\left[\bm{\nabla}\times\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]\right\\}\cdot\bm{Q}_{\bm{A}}=0,$
(83)

where

$\displaystyle\mathcal{\mathbb{P}}_{1}^{\mu}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left(\bm{\nabla}Q_{0}+\frac{1}{c}D_{t}\bm{Q}_{\bm{A}}\right)$
$\displaystyle\cdot\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial D_{\mu_{1}}\cdots
D_{\mu_{j-1}}D_{\mu}D_{\mu_{j+1}}\cdots D_{\mu_{i}}\bm{E}}\right],$ (84)
$\displaystyle\mathbb{P}_{2}^{\mu}$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\left(-\bm{\nabla}\times\bm{Q}_{\bm{A}}\right)$
$\displaystyle\cdot\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial D_{\mu_{1}}\cdots
D_{\mu_{j-1}}D_{\mu}D_{\mu_{j+1}}\cdots D_{\mu_{i}}\bm{B}}\right].$ (85)

The last two terms on the left-hand-side of Eq. (83) vanish due to Eqs. (71)
and (72), but the fourth term does not because of the weak EL equation (75).
If the characteristic $\bm{q}_{a}$ is independent of $\bm{x}$, $\bm{E}$, and
$\bm{B}$, the conservation law of the symmetry is established as

$\displaystyle\partial_{t}\left[\mathcal{L}\xi^{t}-\frac{1}{c}\bm{Q}_{\bm{A}}\cdot\bm{E}_{\bm{E}}\left(\mathcal{L}\right)+\sum_{a}\left(\bm{q}_{a}\cdot\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)\right]+\bm{\nabla}\cdot\left[\mathcal{L}\bm{\kappa}-Q_{0}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right.$
$\displaystyle\left.\vphantom{\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}}+\bm{Q}_{\bm{A}}\times\bm{E}_{\bm{B}}\left(\mathcal{L}\right)+\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right)\cdot\bm{q}_{a}\right]+D_{\mu}\left[\mathcal{\mathbb{P}}_{1}^{\mu}+\mathcal{\mathbb{P}}_{2}^{\mu}\right]=0.$
(86)

###  III.2 Gauge-symmetric energy conservation law

We first derive the gauge-symmetric energy conservation law, assuming that the
action $\mathcal{A}\equiv\int\mathcal{L}dtd^{3}\bm{x}$ is unchanged under the
time translation

$\left(t,\bm{x},\bm{X}_{a},\varphi,\bm{A}\right)\mapsto\left(t+\epsilon,\bm{x},\bm{X}_{a},\varphi,\bm{A}\right),\epsilon\in\mathbb{R}.$
(87)

The infinitesimal generator and characteristic are calculated as

$\displaystyle\bm{v}=\frac{\partial}{\partial
t},\medspace\xi^{t}=1,\bm{\kappa}=0,\>\bm{\theta}_{a}=0,\medspace\phi_{0}=\bm{\phi}_{\bm{A}}=0,$
(88)
$\displaystyle\bm{q}_{a}=-\dot{\bm{X}}_{a},\text{$\;$}Q_{0}=-\varphi_{,t},\medspace\bm{Q}_{\bm{A}}=-\bm{A}_{,t}.$
(89)

And the infinitesimal criterion (80) of the symmetry is

$\frac{\partial\mathcal{L}}{\partial t}=0.$ (90)

The corresponding energy conservation law is thus

$\displaystyle\partial_{t}\left[\mathcal{L}+\frac{1}{c}\bm{A}_{,t}\cdot\bm{E}_{\bm{E}}\left(\mathcal{L}\right)-\sum_{a}\left(\dot{\bm{X}}_{a}\cdot\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)\right]+\bm{\nabla}\cdot\left\\{\left[\varphi_{,t}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)-\bm{A}_{,t}\times\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]\right.$
$\displaystyle\left.\vphantom{\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}}-\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right)\cdot\dot{\bm{X}}_{a}\right\\}+D_{\mu}\left[\mathcal{\mathbb{P}}_{1}^{\mu}+\mathcal{\mathbb{P}}_{2}^{\mu}\right]=0,$
(91)

where

$\displaystyle\mathcal{\mathbb{P}}_{1}^{\mu}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\partial_{t}\bm{E}\cdot\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial D_{\mu_{1}}\cdots
D_{\mu_{j-1}}D_{\mu}D_{\mu_{j+1}}\cdots D_{\mu_{i}}\bm{E}}\right],$ (92)
$\displaystyle\mathcal{\mathbb{P}}_{2}^{\mu}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\partial_{t}\bm{B}\cdot\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial D_{\mu_{1}}\cdots
D_{\mu_{j-1}}D_{\mu}D_{\mu_{j+1}}\cdots D_{\mu_{i}}\bm{B}}\right].$ (93)

The energy density and flux in Eq. (91) are obviously gauge dependent. To
gauge-symmetrize the conservation law, we add Eq. (67) to Eq. (91) and obtain,

$\displaystyle\partial_{t}\left[\mathcal{L}-\frac{\partial\mathcal{L}}{\partial\varphi}\varphi-\sum_{a}\left(\dot{\bm{X}}_{a}\cdot\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)-\bm{E}\cdot\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]+\bm{\nabla}\cdot\left\\{c\bm{E}\times\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right.$
$\displaystyle\left.\vphantom{\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}}+c\frac{\partial\mathcal{L}}{\partial\bm{A}}\varphi-\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right)\cdot\dot{\bm{X}}_{a}\right\\}+D_{\mu}\left[\mathcal{\mathbb{P}}_{1}^{\mu}+\mathcal{\mathbb{P}}_{2}^{\mu}\right]=0.$
(94)

In deriving Eq. (94), we have rewritten the first and second terms of Eq. (67)
as

$\displaystyle\frac{D}{Dt}\left\\{\frac{D}{D\bm{x}}\cdot\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi\right]\right\\}=\frac{D}{Dt}\left\\{\bm{\nabla}\cdot\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]\varphi+\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\cdot\bm{\nabla}\varphi\right\\}$
$\displaystyle=\frac{D}{Dt}\left\\{-\frac{\partial\mathcal{L}}{\partial\varphi}\varphi+\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\cdot\bm{\nabla}\varphi\right\\},$
(95)
$\displaystyle\frac{D}{D\bm{x}}\cdot\left\\{\frac{D}{Dt}\left[-\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi\right]\right\\}=\frac{D}{D\bm{x}}\cdot\left\\{-\frac{\partial}{\partial
t}\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]\varphi-\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi_{,t}\right\\}$
$\displaystyle=\frac{D}{D\bm{x}}\cdot\left\\{c\bm{\nabla}\times\left[\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]\varphi+c\frac{\partial\mathcal{L}}{\partial\bm{A}}\varphi-\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi_{,t}\right\\}$
$\displaystyle=\frac{D}{D\bm{x}}\cdot\left\\{-\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\varphi_{,t}-c\bm{\nabla}\varphi\times\bm{E}_{\bm{B}}\left(\mathcal{L}\right)+c\frac{\partial\mathcal{L}}{\partial\bm{A}}\varphi\right\\},$
(96)

where use has been made of Eqs. (71) and (72). Adding Eqs. (95) and (96) to
Eq. (91) leads to Eq. (94).

We now prove that the energy density and flux in Eq. (94) are gauge symmetric.
It suffices to show that the following terms

$\displaystyle s_{1}$
$\displaystyle\equiv\mathcal{L}-\frac{\partial\mathcal{L}}{\partial\varphi}\varphi-\sum_{a}\left(\dot{\bm{X}}_{a}\cdot\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right),$
(97) $\displaystyle\bm{s}_{2}$ $\displaystyle\equiv
c\frac{\partial\mathcal{L}}{\partial\bm{A}}\varphi-\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right)\cdot\dot{\bm{X}}_{a}$

are gauge symmetric. Substituting Eq. (73) into the expression of $s_{1}$, we
have

$\displaystyle
s_{1}=\sum_{a}\left[-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right]+\text{GIP}\left(\mathcal{L}\right)$
$\displaystyle-\frac{\partial}{\partial\varphi}\left[\sum_{a}\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)\right]\varphi-\sum_{a}\left[\dot{\bm{X}}_{a}\cdot\frac{\partial}{\partial\dot{\bm{X}}_{a}}\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)\right]$
$\displaystyle=\text{GSP}\left(\mathcal{L}\right).$ (98)

Similarly, $\bm{s}_{2}$ is also gauge symmetric,

$\displaystyle\bm{s}_{2}=c\frac{\partial\mathcal{L}}{\partial\bm{A}}\varphi-\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right)\cdot\dot{\bm{X}}_{a}=c\frac{\partial}{\partial\bm{A}}\left[\sum_{a}\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)\right]\varphi$
$\displaystyle-\sum_{a}\left\\{\dot{\bm{X}}_{a}\frac{\partial}{\partial\dot{\bm{X}}_{a}}\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)-\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)\bm{I}+\text{GSP}\left(\mathcal{L}_{a}\right)\bm{I}\right\\}\cdot\dot{\bm{X}}_{a}$
$\displaystyle=\sum_{a}\text{GSP}\left(\mathcal{L}_{a}\right)\dot{\bm{X}}_{a}.$
(99)

###  III.3 Gauge-symmetric momentum conservation law

We now discuss how to derive a gauge-symmetric momentum conservation law,
assuming that the action $\mathcal{A}\equiv\int\mathcal{L}dtd^{3}\bm{x}$ of
the electromagnetic field-charged particle system is invariant under the space
translation

$\left(t,x,\bm{X}_{a},\varphi,\bm{A}\right)\mapsto\left(t,\bm{x}+\epsilon\bm{h},\bm{X}_{a}+\epsilon\bm{h},\varphi,\bm{A}\right),\epsilon\in\mathbb{R}.$
(100)

We emphasize that, different from the situation in standard field theories,
this symmetry group simultaneously translates both the spatial coordinate
$\bm{x}$ for the field and particle’s position $\bm{X}_{a}$ (Qin  _et al._ ,
2014; Fan  _et al._ , 2018, 2019). The infinitesimal criterion of this
symmetry is

$\frac{\partial\mathcal{L}}{\partial\bm{x}}+\sum_{a}\frac{\partial\mathcal{L}}{\partial\bm{X}_{a}}=0.$
(101)

From Eq. (100), the infinitesimal generator and its characteristic are

$\displaystyle\bm{v}=\bm{h}\cdot\sum_{a}\left(\frac{\partial}{\partial\bm{x}}+\frac{\partial}{\partial\bm{X}_{a}}\right),\medspace\xi^{t}=0,\bm{\kappa}=\bm{\theta}_{a}=\bm{h},,\medspace\phi_{0}=\bm{\phi}_{\bm{A}}=0,$
(102)
$\displaystyle\bm{q}_{a}=\bm{h},Q_{0}=-\bm{h}\cdot\bm{\nabla}\varphi,\medspace\bm{Q}_{\bm{A}}=-\bm{h}\cdot\bm{\nabla}\bm{A}.$
(103)

The corresponding momentum conservation law is obtained by substituting Eqs.
(102) and (103) into Eq. (86), i.e.,

$\displaystyle\partial_{t}\left[\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\cdot\left(\bm{\nabla}\bm{A}\right)^{T}+\sum_{a}\left(\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)\right]+\bm{\nabla}\cdot\left\\{\mathcal{L}\bm{I}+\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}-\mathcal{L}_{a}\bm{I}\right)\right.$
$\displaystyle\left.\vphantom{\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}}+\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{\nabla}\varphi+\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\times\left(\bm{\nabla}\bm{A}\right)^{T}\right]\right\\}+D_{\mu}\left[\mathcal{\bar{\mathbb{\bm{P}}}}_{1}^{\mu}+\mathcal{\bar{\mathbb{P}}}_{2}^{\mu}\right]=0,$
(104)

where

$\displaystyle\mathcal{\bar{\mathbb{\bm{P}}}}_{1}^{\mu}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\bm{\nabla}\bm{E}\cdot\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial D_{\mu_{1}}\cdots
D_{\mu_{j-1}}D_{\mu}D_{\mu_{j+1}}\cdots D_{\mu_{i}}\bm{E}}\right],$ (105)
$\displaystyle\mathcal{\bar{\mathbb{\bm{P}}}}_{2}^{\mu}=\sum_{i=1}^{n}\sum_{j=1}^{i}\left(-1\right)^{j}D_{\mu_{j+1}}\cdots
D_{\mu_{i}}\bm{\nabla}\bm{B}\cdot\left[D_{\mu_{1}}\cdots
D_{\mu_{j-1}}\frac{\partial\mathcal{L}}{\partial D_{\mu_{1}}\cdots
D_{\mu_{j-1}}D_{\mu}D_{\mu_{j+1}}\cdots D_{\mu_{i}}\bm{B}}\right].$ (106)

Again, the momentum density and flux in Eq. (104) are gauge dependent. We add
Eq. (68) to Eq. (104) to obtain a gauge-symmetric momentum conservation law,

$\displaystyle\partial_{t}\left[\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\times\bm{B}+\frac{1}{c}\frac{\partial\mathcal{L}}{\partial\varphi}\bm{A}+\sum_{a}\left(\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)\right]+\bm{\nabla}\cdot\left\\{-\frac{\partial\mathcal{L}}{\partial\bm{A}}\bm{A}+\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}\right)\right.$
$\displaystyle\left.\vphantom{\frac{\partial\mathcal{L}}{\partial\bm{A}_{\alpha}}}\left[\mathcal{L}_{F}-\bm{B}\cdot\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]\bm{I}+\left[-\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{E}+\bm{B}\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]\right\\}+D_{\mu}\left[\left(\mathcal{\bar{\mathbb{\bm{P}}}}_{1}^{\mu}+\mathcal{\bar{\mathbb{\bm{P}}}}_{2}^{\mu}\right)\right]=0.$
(107)

In the derivation of Eq. (107), we have rewritten the first and second terms
of Eq. (68) as

$\displaystyle\frac{D}{Dt}\left\\{\frac{D}{D\bm{x}}\cdot\left[-\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}\right]\right\\}=\frac{D}{Dt}\left\\{-\frac{1}{c}\bm{\nabla}\cdot\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]\bm{A}-\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\cdot\bm{\nabla}\bm{A}\right\\}$
$\displaystyle=\frac{D}{Dt}\left[\frac{1}{c}\frac{\partial\mathcal{L}}{\partial\varphi}\bm{A}-\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\cdot\bm{\nabla}\bm{A}\right],$
(108)
$\displaystyle\frac{D}{D\bm{x}}\cdot\left\\{\frac{D}{Dt}\left[\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}\right]\right\\}=\frac{D}{D\bm{x}}\cdot\left\\{\frac{1}{c}\frac{\partial}{\partial
t}\left[\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\right]\bm{A}+\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}_{,t}\right\\}$
$\displaystyle=\frac{D}{D\bm{x}}\cdot\left\\{-\bm{\nabla}\times\left[\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\right]\bm{A}-\frac{\partial\mathcal{L}}{\partial\bm{A}}\bm{A}+\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}_{,t}\right\\}$
$\displaystyle=\frac{D}{D\bm{x}}\cdot\left\\{-\bm{E}_{\bm{B}}\left(\mathcal{L}\right)\times\bm{\nabla}\bm{A}-\frac{\partial\mathcal{L}}{\partial\bm{A}}\bm{A}+\frac{1}{c}\bm{E}_{\bm{E}}\left(\mathcal{L}\right)\bm{A}_{,t}\right\\},$
(109)

where use has been made of Eqs. (71) and (72). Adding Eqs. (108) and (109)
into Eq. (104) gives Eq. (107).

To show that the momentum density and flux in Eq. (107) are gauge symmetric,
it suffices to show that the following terms are gauge symmetric,

$\displaystyle\bm{t}_{1}$
$\displaystyle\equiv\frac{1}{c}\frac{\partial\mathcal{L}}{\partial\varphi}\bm{A}+\sum_{a}\left(\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right),$
(110) $\displaystyle\bm{t}_{2}$
$\displaystyle\equiv-\frac{\partial\mathcal{L}}{\partial\bm{A}}\bm{A}+\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}\right).$

Substituting Eq. (73) into the expression of $\bm{t}_{1},$ we can see that it
is gauge symmetric, i.e.,

$\displaystyle\bm{t}_{1}=\frac{1}{c}\frac{\partial\mathcal{L}}{\partial\varphi}\bm{A}+\sum_{a}\left(\frac{\partial\mathcal{L}}{\partial\dot{\bm{X}}_{a}}\right)=\frac{1}{c}\sum_{a}\frac{\partial}{\partial\varphi}\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)\bm{A}+\sum_{a}\frac{\partial}{\partial\dot{\bm{X}}_{a}}\left[-q_{a}\varphi\delta_{a}\right.$
$\displaystyle\left.\vphantom{\left(\mathcal{L}_{a}\right)}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}+\text{GSP}\left(\mathcal{L}_{a}\right)\right]=\sum_{a}\frac{\partial}{\partial\dot{\bm{X}}_{a}}\left[\text{GSP}\left(\mathcal{L}_{a}\right)\right].$
(111)

Similarly, $\bm{t}_{2}$ is also gauge-symmetric,

$\displaystyle\bm{t}_{2}=-\frac{\partial\mathcal{L}}{\partial\bm{A}}\bm{A}+\sum_{a}\left(\dot{\bm{X}}_{a}\frac{\partial\mathcal{L}_{a}}{\partial\dot{\bm{X}}_{a}}\right)=-\frac{\partial}{\partial\bm{A}}\left(-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}\right)\bm{A}$
$\displaystyle+\sum_{a}\dot{\bm{X}}_{a}\frac{\partial}{\partial\dot{\bm{X}}_{a}}\left[-q_{a}\varphi\delta_{a}+\frac{q_{a}}{c}\bm{A}\cdot\dot{\bm{X}}_{a}\delta_{a}+\text{GSP}\left(\mathcal{L}_{a}\right)\right]=\sum_{a}\dot{\bm{X}}_{a}\frac{\partial}{\partial\dot{\bm{X}}_{a}}\left[\text{GSP}\left(\mathcal{L}_{a}\right)\right].$
(112)

##  IV conclusion

In this study, we developed a gauge-symmetrization method for the energy and
momentum conservation laws in general high-order classical electromagnetic
field theories, which appear in the study of gyrokinetic systems (Qin, 2005;
Qin  _et al._ , 2007; Fan  _et al._ , 2020) for magnetized plasmas and the
Podolsky system (Bopp, 1940; Podolsky, 1942) for the radiation reaction of
classical charged particles. The method only removes the electromagnetic gauge
dependence from the canonical EMT derived from the spacetime translation
symmetry, without necessarily symmetrizing the EMT with respect to the tensor
indices. This is adequate for applications not involving general relativity.

To achieve this goal, we reformulated the EL equation and infinitesimal
criterion in terms of the Faraday tensor $F_{\mu\nu}$. The canonical EMT
$T_{\text{N}}^{\mu\nu}$ is derived using this formalism, and it was found that
the gauge dependent part of $T_{\text{N}}^{\mu\nu}$ can be removed by adding
the divergence of the displacement-potential tensor, which is defined as

$\displaystyle\mathcal{F}^{\sigma\mu\nu}\equiv\frac{1}{4\pi}\mathcal{D}^{\sigma\mu}A^{\nu}.$
(113)

It was shown that the displacement-potential tensor
$\mathcal{F}^{\sigma\mu\nu}$ is related to the well-known BR super-potential
$\mathcal{S}^{\sigma\mu\nu}$ as

$\mathcal{S}^{\sigma\mu\nu}=\mathcal{F}^{\sigma\mu\nu}+\frac{1}{2}\left[\Delta^{\sigma\nu\mu}-\Delta^{\mu\nu\sigma}-\Delta^{\nu\mu\sigma}\right],$
(114)

where $\Delta^{\sigma\nu\mu}$ is defined in Eq. (51). Using the example of the
Podolsky system (Bopp, 1940; Podolsky, 1942), we show that
$\Delta^{\sigma\nu\mu}$ in general is non-vanishing for high-order field
theories. For a first-order field theory, such as the standard Maxwell system
(6), $\Delta^{\sigma\nu\mu}$ vanishes such that
$\mathcal{S}^{\sigma\mu\nu}=\mathcal{F}^{\sigma\mu\nu}$. In the case, the
method developed can be used as a simpler procedure to calculate the BR super-
potential $\mathcal{S}^{\sigma\mu\nu}$ without the necessity to calculate the
angular momentum tensor in 4D spacetime.

Lastly, we applied the method to derive gauge-symmetric EMTs for high-order
electromagnetic systems coupled with classical charged particles. Using the
“3+1” form of Eq. (36), we obtained the explicitly gauge-symmetric energy and
momentum conservation laws in a general setting [see Eqs. (94) and (107)].

###### Acknowledgements.

P. Fan was supported by Shenzhen Clean Energy Research Institute and National
Natural Science Foundation of China (NSFC-12005141). J. Xiao was supported by
the National MC Energy R&D Program (2018YFE0304100), National Key Research and
Development Program (2016YFA0400600, 2016YFA0400601 and 2016YFA0400602), and
the National Natural Science Foundation of China (NSFC-11905220 and 11805273).
H. Qin was supported by the U.S. Department of Energy (DE-AC02-09CH11466).

## References

  * Noether (1918) E. Noether, Nachr. König. Gesell. Wiss Göttingen, Math. -Phys. Kl. 235 (1918), also available in English at Transport Theory and Statistical Physics 1, 186-207 (1971).
  * Babak and Grishchuk (1999) S. V. Babak and L. P. Grishchuk, Physical Review D 62, 024038 (1999).
  * Gratus  _et al._ (2012) J. Gratus, Y. N. Obukhov, and R. W. Tucker, Annals of Physics 327, 2560 (2012).
  * Arminjon (2016) M. Arminjon, Advances in Mathematical Physics 2016, 1 (2016).
  * Inglis and Jarvis (2016) S. Inglis and P. Jarvis, Annals of Physics 366, 57 (2016).
  * Jiménez  _et al._ (2018) J. B. Jiménez, J. A. R. Cembranos, and J. M. Sánchez Velázquez, Journal of High Energy Physics 2018, 100 (2018).
  * Ilin and Paston (2020) R. Ilin and S. Paston, Universe 6, 173 (2020).
  * Belinfante (1939) F. Belinfante, Physica 6, 887 (1939).
  * Belinfante (1940) F. Belinfante, Physica 7, 449 (1940).
  * Rosenfeld (1940) L. Rosenfeld, Acad. Roy. Belg. Memoirs de Classes de Science 18, 1 (1940).
  * Hawking and Ellis (1973) S. W. Hawking and G. F. R. Ellis, “The large-scale structure of space-time,” (Cambridge University Press, New York, 1973) pp. 64–71.
  * Landau and lifshitz (1975) L. D. Landau and E. M. lifshitz, “The classical theory of fields,” (Butterworth-Heinemann, Oxford, 1975) pp. 46–89.
  * Gotay and Marsden (1992) M. J. Gotay and J. E. Marsden, Contemp. Math. 132, 367 (1992).
  * Lopez  _et al._ (2007) M. C. Lopez, M. J. Gotay, and J. E. Marsden,  (2007), arXiv:0711.4679 .
  * Zhang (2005) H. B. Zhang, Communications in Theoretical Physics 44, 1007 (2005).
  * Ilin and Paston (2019) R. V. Ilin and S. A. Paston, European Physical Journal Plus 134, 21 (2019).
  * Baker  _et al._ (2021) M. R. Baker, N. Kiriushcheva, and S. Kuzmin, Nuclear Physics B 962, 115240 (2021).
  * Qin (2005) H. Qin,  _A Short Introduction to General Gyrokinetic Theory_, Tech. Rep. (Princeton, NJ, 2005).
  * Qin  _et al._ (2007) H. Qin, R. H. Cohen, W. M. Nevins, and X. Q. Xu, Physics of Plasmas 14, 056110 (2007).
  * Fan  _et al._ (2020) P. Fan, H. Qin, and J. Xiao, “Discovering exact local energy-momentum conservation laws for electromagnetic gyrokinetic system by high-order field theory on heterogeneous manifolds,” (2020), arXiv:2006.11039 [physics.plasm-ph] .
  * Bopp (1940) F. Bopp, Annalen der Physik 430, 345 (1940).
  * Podolsky (1942) B. Podolsky, Physical Review 62, 68 (1942).
  * Olver (1993) P. J. Olver, “Applications of lie groups to differential equations,” (Springer-Verlag, New York, 1993) pp. 90–130, 242–283.
  * Blaschke  _et al._ (2016) D. N. Blaschke, F. Gieres, M. Reboud, and M. Schweda, Nuclear Physics B 912, 192 (2016).
  * Qin  _et al._ (2014) H. Qin, J. W. Burby, and R. C. Davidson, Physical Review E 90, 043102 (2014).
  * Fan  _et al._ (2018) P. Fan, H. Qin, J. Liu, N. Xiang, and Z. Yu, Frontiers of Physics 13, 135203 (2018).
  * Fan  _et al._ (2019) P. Fan, H. Qin, J. Xiao, and N. Xiang, Physics of Plasmas 26, 062115 (2019).