Constraints for electric charge from Maxwell's equations and boundary conditions

Suppose that the system occupies an arbitrary physical state. To derive our first and second boundary constraints, we begin by considering Gauss's law ³ in the form

$$\begin{eqnarray}&&{\unicode{x0222F}}_{{\mathscr{S}}}\langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)\cdot {{\rm{d}}}^{2}{\bf{r}}=\displaystyle \frac{\langle {\hat{Q}}_{{\mathscr{V}}}\rangle }{{\epsilon }_{0}},\end{eqnarray} \tag{ 2 }$$

where the closed surface ${\mathscr{S}}$ is such that the enclosed volume ${\mathscr{V}}$ is the fundamental domain. The periodic boundary conditions and corresponding restriction to allowed wavevectors see the electric field matched on opposite walls of the fundamental domain:

$$\begin{eqnarray*}&&\langle \hat{{\bf{E}}}\rangle (x=0,y,z,t)=\langle \hat{{\bf{E}}}\rangle (x=L,y,z,t)\end{eqnarray*}$$

with analogous results for the walls of constant y and z. It follows that the electric flux through ${\mathscr{S}}$ vanishes:

$$\begin{eqnarray*}\begin{array}{rcl}{\unicode{x0222F}}_{{\mathscr{S}}}\langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)\cdot {{\rm{d}}}^{2}{\bf{r}} & = & -{\displaystyle \int }_{0}^{L}{\displaystyle \int }_{0}^{L}\langle {\hat{E}}_{x}\rangle (x=0,y,z,t){\rm{d}}y{\rm{d}}z+{\displaystyle \int }_{0}^{L}{\displaystyle \int }_{0}^{L}\langle {\hat{E}}_{x}\rangle (x=L,y,z,t){\rm{d}}y{\rm{d}}z+...\\ & = & 0.\end{array}\end{eqnarray*}$$

According to Gauss's law (2), the total electric charge in ${\mathscr{V}}$ should also vanish:

$$\begin{eqnarray}&&\langle {\hat{Q}}_{{\mathscr{V}}}\rangle =0.\end{eqnarray} \tag{ 3 }$$

Explicit calculation reveals, however, that the total electric charge in ${\mathscr{V}}$ does not vanish immediately; it is dependent on electric charges and occupation numbers:

$$\begin{eqnarray}&&\langle {\hat{Q}}_{{\mathscr{V}}}\rangle =\sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}\langle {\hat{N}}_{{ \mathcal F }}-{\hat{N}}_{\overline{{ \mathcal F }}}+\sum _{{\bf{k}}}2\rangle .\end{eqnarray} \tag{ 4 }$$

The first contribution on the right-hand side of equation (4) is the total electric charge of the matter particles in ${\mathscr{V}}$, the second is the total electric charge of the matter antiparticles in ${\mathscr{V}}$ and the third is the total zero-point electric charge in ${\mathscr{V}}$, which is the analogue for electric charge of the total zero-point energy in ${\mathscr{V}}$ and emerges as a consequence of anticommutation relations in a similar way. For more information on the origin of zero-point electric charge, see the relevant discussions in [38, 39]. Comparing equations (3) and (4), we see that, firstly, we must satisfy

$$\begin{eqnarray}&&\sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}=0\end{eqnarray} \tag{ 5 }$$

if we are to have any hope of obeying Gauss's law (2) for a state describing finite numbers of matter particles and antiparticles, as the summation over wavevectors in equation (4) diverges ⁴ ; the finite cannot neutralise the infinite. Equation (5) is an embodiment of our first boundary constraint; the zero-point electric charge must vanish. Secondly, we must satisfy

$$\begin{eqnarray}&&\sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}\langle {\hat{N}}_{{ \mathcal F }}-{\hat{N}}_{\overline{{ \mathcal F }}}\rangle =0.\end{eqnarray} \tag{ 6 }$$

With our first boundary constraint (5) understood, equation (6) is an embodiment of our second boundary constraint; the total electric charge in the fundamental domain must vanish.

Suppose again that the system occupies an arbitrary physical state. To derive our third boundary constraint, we begin by considering the Ampère-Maxwell law in the form

$$\begin{eqnarray}&&{\oint }_{{\mathscr{C}}}\langle \hat{{\bf{B}}}\rangle ({\bf{r}},t)\cdot {\rm{d}}{\bf{r}}={\mu }_{0}\langle {\hat{I}}_{{\mathscr{S}}^{\prime} }\rangle (t)+\displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\iint }_{{\mathscr{S}}^{\prime} }\langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)\cdot {{\rm{d}}}^{2}{\bf{r}},\end{eqnarray} \tag{ 7 }$$

where the closed curve ${\mathscr{C}}$ lies on the walls of the fundamental domain at z = L/2, circling the z axis according to the right-hand rule, and the enclosed surface ${\mathscr{S}}^{\prime}$ is a square of cross-sectional area A = L². The periodic boundary conditions and corresponding restriction to allowed wavevectors see the magnetic field matched on opposite walls of the fundamental domain:

$$\begin{eqnarray*}&&\langle \hat{{\bf{B}}}\rangle (x=0,y,z,t)=\langle \hat{{\bf{B}}}\rangle (x=L,y,z,t)\end{eqnarray*}$$

with analogous results for the walls of constant y and z. It follows that the line integral of the magnetic field around ${\mathscr{C}}$ vanishes:

$$\begin{eqnarray*}\begin{array}{rcl}{\oint }_{{\mathscr{C}}}\langle \hat{{\bf{B}}}\rangle ({\bf{r}},t)\cdot {\rm{d}}{\bf{r}} & = & -{\displaystyle \int }_{0}^{L}\langle {\hat{B}}_{y}\rangle (x=0,y,z=L/2,t){\rm{d}}y+{\displaystyle \int }_{0}^{L}\langle {\hat{B}}_{y}\rangle (x=L,y,z=L/2,t){\rm{d}}y+...\\ & = & 0.\end{array}\end{eqnarray*}$$

According to the Ampère-Maxwell law (7), the sum of the total electric and displacement currents through ${\mathscr{S}}^{\prime}$ should also vanish:

$$\begin{eqnarray}&&\langle {\hat{I}}_{{\mathscr{S}}^{\prime} }\rangle (t)+{\epsilon }_{0}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\iint }_{{\mathscr{S}}^{\prime} }\langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)\cdot {{\rm{d}}}^{2}{\bf{r}}\,=\,0.\end{eqnarray} \tag{ 8 }$$

Explicit calculation reveals, however, that the sum of the total electric and displacement currents through ${\mathscr{S}}^{\prime}$ does not vanish immediately; it is dependent on electric charges, occupation numbers and velocities:

$$\begin{eqnarray}&&\langle {\hat{I}}_{{\mathscr{S}}^{\prime} }\rangle (t)+{\epsilon }_{0}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\iint }_{{\mathscr{S}}^{\prime} }\langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)\cdot {{\rm{d}}}^{2}{\bf{r}}=\displaystyle \frac{A}{V}\sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}{q}_{{ \mathcal F }}{v}_{{ \mathcal F }{\bf{k}}z}\langle {\hat{N}}_{{ \mathcal F }{\bf{k}}}-{\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}}\rangle (t)+...,\end{eqnarray} \tag{ 9 }$$

where, for the sake of brevity, we have refrained from showing terms due to Zitterbewegung [40, 41] explicitly. The first contribution on the right-hand side of equation (9) is proportional to the mean piece [42] of the total electric current of the matter particles in ${\mathscr{V}}$ and the second is proportional to the mean piece of the total electric current of the matter antiparticles in ${\mathscr{V}}$. Comparing equations (8) and (9), we see that we must satisfy

$$\begin{eqnarray}&&\displaystyle \sum _{{ \mathcal F }=1}^{{N}_{q}}\displaystyle \sum _{{\bf{k}}}{q}_{{ \mathcal F }}{v}_{{ \mathcal F }{\bf{k}}z} \langle {\hat{N}}_{{ \mathcal F }{\bf{k}}}-{\hat{N}}_{\bar{{ \mathcal F }}{\bf{k}}} \rangle (t)+...\,=\,0\end{eqnarray} \tag{ 10 }$$

to obey the Ampère-Maxwell law (7). Our argument can be repeated for other curves, leading to the overarching conclusion that we must satisfy ⁵

$$\begin{eqnarray}&&\displaystyle \sum _{{ \mathcal F }=1}^{{N}_{q}}\displaystyle \sum _{{\bf{k}}}{q}_{{ \mathcal F }}{{\bf{v}}}_{{ \mathcal F }{\bf{k}}} \langle {\hat{N}}_{{ \mathcal F }{\bf{k}}}-{\hat{N}}_{\bar{{ \mathcal F }}{\bf{k}}} \rangle (t)+...\,=\,0.\end{eqnarray} \tag{ 11 }$$

Equation (11) is an embodiment of our third boundary constraint; the total electric current in the fundamental domain must vanish.

Together, our first boundary constraint (5) and second boundary constraint (6) restrict the possible electric charges ${q}_{{ \mathcal F }}$ and occupation-number differences $\langle {\hat{N}}_{{ \mathcal F }}-{\hat{N}}_{\overline{{ \mathcal F }}}\rangle$ as follows, where we make use of the fact that each of the ${q}_{{ \mathcal F }}$ must be non-zero by construction (otherwise we wouldn't have N_q species of electrically charged matter) and let ${\mathfrak{e}}$ be an arbitrary non-zero quantity of electric charge:

• N_q = 1 is forbidden, as the requirement that q₁ = 0 from our first boundary constraint (5) would be at odds with the requirement that q₁ ≠ 0 by construction; it is not possible to have only one species of electrically charged matter in the theory, as there would be nothing to neutralise the zero-point electric charge. Evidently, a theory (without normal ordering) in which electrons and positrons constitute the only species of electrically charged matter cannot obey Gauss's law.
• If N_q = 2, we must have ${q}_{2}=-{q}_{1}={\mathfrak{e}}$ from our first boundary constraint (5) and thus $\langle {\hat{N}}_{1}-{\hat{N}}_{\overline{1}}\rangle =\langle {\hat{N}}_{2}-{\hat{N}}_{\overline{2}}\rangle$ from our second boundary constraint (6). Electric charge is quantised in this case, albeit somewhat trivially.
• If N_q = 3, we must have one of two possible scenarios; either ${q}_{3}=-{q}_{2}-{q}_{1}={\mathfrak{e}}$ and $\langle {\hat{N}}_{1}-{\hat{N}}_{\overline{1}}\rangle =\langle {\hat{N}}_{2}-{\hat{N}}_{\overline{2}}\rangle =\langle {\hat{N}}_{3}-{\hat{N}}_{\overline{3}}\rangle$ or
  $$\begin{eqnarray}&&{q}_{1}={\mathfrak{e}}( \langle {\hat{N}}_{2}-{\hat{N}}_{\bar{2}} \rangle - \langle {\hat{N}}_{3}-{\hat{N}}_{\bar{3}} \rangle ),\end{eqnarray} \tag{ 12 }$$
  $$\begin{eqnarray}&&{q}_{2}={\mathfrak{e}}( \langle {\hat{N}}_{3}-{\hat{N}}_{\bar{3}} \rangle - \langle {\hat{N}}_{1}-{\hat{N}}_{\bar{1}} \rangle )\end{eqnarray} \tag{ 13 }$$
  $$\begin{eqnarray}&&{q}_{3}={\mathfrak{e}}( \langle {\hat{N}}_{1}-{\hat{N}}_{\bar{1}} \rangle - \langle {\hat{N}}_{2}-{\hat{N}}_{\bar{2}} \rangle )\end{eqnarray} \tag{ 14 }$$
  with $\langle {\hat{N}}_{1}-{\hat{N}}_{\overline{1}}\rangle$, $\langle {\hat{N}}_{2}-{\hat{N}}_{\overline{2}}\rangle$ and $\langle {\hat{N}}_{3}-{\hat{N}}_{\overline{3}}\rangle$ differing from each other, as can be deduced from our first boundary constraint (5) and second boundary constraint (6) by regarding them as planes in (q₁, q₂, q₃) space and considering their intersection. Electric charge is non-trivially quantised in the latter scenario, assuming that $\langle {\hat{N}}_{1}-{\hat{N}}_{\overline{1}}\rangle$, $\langle {\hat{N}}_{2}-{\hat{N}}_{\overline{2}}\rangle$ and $\langle {\hat{N}}_{3}-{\hat{N}}_{\overline{3}}\rangle$ are integers. See figure 1.
• If N_q ≥ 4, more intricate possibilities exist.

In general, the theory cannot sustain electrically charged matter with only one sign of electric charge and it is sufficient but not necessary to have

$$\begin{eqnarray*}&&\sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}=0\ \ \ \ \langle {\hat{N}}_{1}-{\hat{N}}_{\overline{1}}\rangle =\langle {\hat{N}}_{2}-{\hat{N}}_{\overline{2}}\rangle =...=\langle {\hat{N}}_{{N}_{q}}-{\hat{N}}_{\overline{{N}_{q}}}\rangle .\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

For given electric charges and total occupation-number differences, our third boundary constraint (11) restricts the possible velocities.

Evidently, novel forms of electric-charge quantisation can arise even in the basic theory of quantum electrodynamics, as a result of our boundary constraints. This does not appear to have been recognised before. Note that the mechanism here has nothing to do with extra dimensions, magnetic monopoles and/or dyons, grand unification, classical constraints or quantum constraints.

The derivations above are simple and it is natural, therefore, to ask why such observations are not made routinely in the basic theory of quantum electrodynamics.

When Dirac first described the Dirac sea [43], he did recognise a 'difficulty' in that the 'infinite density of electricity' should, according to Gauss's law, 'produce an electric field of infinite divergence'. He went on, however, to ignore the difficulty, arguing that it 'seems natural ... to interpret [the electric charge density] as the departure from the normal state of electrification of the world'. This approach persists to the present day. In Cohen-Tannoudji, Dupont-Roc and Grynberg's textbook on quantum electrodynamics [38], Milonni's textbook on the quantum vacuum [39] and many other seminal works, the zero-point electric charge is circumnavigated without significant comment by focussing on the normal-ordered form(s) of the electric charge density and/or total electric charge. Weinberg highlights Dirac's difficulty in one of his textbooks [44] and, by omission, seems to convey the view that it has not yet been resolved. To derive our first boundary constraint (5), we have refrained from using such heuristic normal-ordering procedures, asking instead that the theory be self consistent in its natural form. In this author's view, the zero-point electric charge is to electric charge what the zero-point energy is to energy and, like the zero-point energy, we should seek to understand it rather than simply removing it. In section 5, we offer a fundamental resolution to Dirac's difficulty.

When periodic boundary conditions are used, they are usually regarded as a mere computational aid without physical meaning; ultimately, the (singular) limit V → ∞ m³ is taken and replacements like

$$\begin{eqnarray}&&\sum _{{\bf{k}}}\to \displaystyle \frac{V}{{\left(2\pi \right)}^{3}}{\iiint }_{\infty }{{\rm{d}}}^{3}{\bf{k}}\end{eqnarray} \tag{ 15 }$$

are made. To derive our second boundary constraint (6) and third boundary constraint (11), we have instead stayed true to the periodic boundary constraints. There are contexts in which it is appropriate to consider a closed space like this. In sections 4 and 6, we consider systems such as an ideal crystal, for example, where three-dimensional periodic boundary conditions are physically meaningful. In sections 5 and 6, we consider the possibility that the entire Universe is closed.

In the appendix, we show that our boundary constraints (5), (6) and (11) emerge via Heisenberg's equation of motion only when we are careful enough to exclude electromagnetic modes with wavevector k = 0; there is a factor of $1/\sqrt{| {\bf{k}}| }$ in the mode expansion of the potential four-vector ${\hat{A}}^{\alpha }({\bf{r}})$, which is not defined for k = 0. This subtlety is easy to overlook and is obscured when replacements like (15) are made, as $1/\sqrt{| {\bf{k}}| }$ multiplied by the continuous reciprocal-space volume element d³ k appears to be well defined for k = 0 in the usual spherical-type coordinates ∣k∣, ϑ and φ; ${{\rm{d}}}^{3}{\bf{k}}/\sqrt{| {\bf{k}}| }=| {\bf{k}}{| }^{3/2}\sin \vartheta {\rm{d}}| {\bf{k}}| {\rm{d}}\vartheta {\rm{d}}\varphi$.

The arguments we have made above in support of our boundary constraints (5), (6) and (11) are essentially geometrical in nature. The fundamental domain has the topology of a three-dimensional torus due to the periodic boundary conditions; an attempt to exit through one of its walls results in a re-entrance through the opposite wall. It is impossible for a non-vanishing electric flux to simultaneously exit and enter the fundamental domain or for a non-vanishing line integral of the magnetic field around the walls of the fundamental domain to simultaneously circulate in opposite directions.

Our first boundary constraint (5) also applies if, instead of a fundamental domain with three-dimensional periodic boundary conditions, we consider a space of infinite extent with the boundary condition that the electric field vanish at spatial infinity; it applies whether the space is closed or open. In contrast, our second boundary constraint (6) and third boundary constraint (11) apply no matter how large the volume V of the fundamental domain is but do not obviously apply in a space of infinite extent; they are emphatically due to the closed nature of the space. The limit V → ∞ m³ is singular [45] in this regard.

As our arguments are essentially geometrical in nature and assume only the validity of Maxwell's equations, it should be clear that derivations analogous to those presented above apply in other relevant theories, including theories with curved rather than flat spacetimes; our boundary constraints have a generality that transcends the specific formulation of the basic theory of quantum electrodynamics considered in this section.

When treating electrostatic problems using three-dimensional periodic boundary conditions, our second boundary constraint tells us to ensure that the total electric charge in the fundamental domain vanishes. This is true whether the periodic boundary conditions are used merely as a computational aid or to model a system with real three-dimensional quasi-periodicity. Possible examples of the former include the simulation of an electrically charged colloid [46, 47], the simulation of an electrically charged quasiparticle or defect in the solid state [48] and the simulation of an electrolyte solution [49]. Examples of the latter include an ionic atomic crystal [50–52], a Wigner crystal [53–55], a noble metal crystal [56] and an ionic colloidal crystal [57]. The requirement that the total electric charge in the fundamental domain vanishes can be seen in Ewald's original description of his summation technique [50] as well as Evjen's [51] and it is usually respected. The requirement is usually justified, however, on the grounds that it helps prevent the electrostatic potential and thus the Coulomb energy from diverging [48, 49] or on the grounds that it helps to ensure that a Ewald summation is independent of the chosen screening factor [58]. Although such arguments are complementary to the arguments we have made in section 3, they are nevertheless weaker, as they invite attempts to 'renormalise' the divergent quantities, thus leading to erroneous results [49, 58]. The arguments we have made in section 3 show clearly that such attempts are doomed to failure; it is geometrically impossible to have a non-vanishing total electric charge in the fundamental domain. This precludes the use of three-dimensional periodic boundary conditions to simulate or model non-neutral plasmas [59–61], for example.

When treating more general electromagnetic problems using three-dimensional periodic boundary conditions, our third boundary constraint tells us to ensure that, in addition to the above, the total electric current in the fundamental domain vanishes. Again, this is true whether the periodic boundary conditions are used merely as a computational aid or to model a system with real three-dimensional quasi-periodicity. A possible example of the former is the simulation of a plasma [62]. An example of the latter is a thin-wire metamaterial [63, 64].

Our first boundary constraint applies whether the Universe is closed or open and reads

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{{ \mathcal F }}{q}_{{ \mathcal F }}{d}_{{ \mathcal F }} & = & \displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{ \mathcal F }}{({\rm{k}}{\rm{n}}{\rm{o}}{\rm{w}}{\rm{n}})}}{q}_{{ \mathcal F }}+\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{ \mathcal F }}{({\rm{y}}{\rm{e}}{\rm{t}}\,{\rm{t}}{\rm{o}}\,{\rm{b}}{\rm{e}}\,{\rm{d}}{\rm{i}}{\rm{s}}{\rm{c}}{\rm{o}}{\rm{v}}{\rm{e}}{\rm{r}}{\rm{e}}{\rm{d}})}}{q}_{{ \mathcal F }}{d}_{{ \mathcal F }}\\ & = & 0,\end{array}\end{eqnarray} \tag{ 16 }$$

where the ${ \mathcal F }$ summation runs over all elementary species of Fermion in the Universe; ${q}_{{ \mathcal F }}$ is the electric charge of the ${ \mathcal F }$th species; ${d}_{{ \mathcal F }}$ with ${d}_{{ \mathcal F }}=1$ for a Weyl-type field is a spin-multiplicity factor, included for generality ⁶ , and we have partitioned the summation into a summation over known species and a summation over any species yet to be discovered. In this author's view, our first boundary constraint (16) embodies the fundamental resolution of Dirac's difficulty [43]; the zero-point electric charge poses no problem, as it vanishes by necessity. In 1930, there was little prospect of Dirac recognising this, as only a handful of particle types were known.

Pleasingly, we find that our first boundary constraint (16) is satisfied by the known elementary species of Fermion in the Universe, as described by the minimal standard model of particle physics:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{ \mathcal F }}{({\rm{k}}{\rm{n}}{\rm{o}}{\rm{w}}{\rm{n}})}}{q}_{{ \mathcal F }} & = & \displaystyle \sum _{f=1,2,3}\left\{\displaystyle \sum _{{\sigma }=L,R}\left[\displaystyle \sum _{c=r,g,b}\left(\displaystyle \frac{2e}{3}-\displaystyle \frac{e}{3}\right)-e\right]+0\right\}+\displaystyle \frac{1}{2}{\epsilon }e\left[\displaystyle \sum _{{\sigma }=L,R}\left(1-1\right)+1-1\right]\\ & = & 0,\end{array}\end{eqnarray} \tag{ 17 }$$

where the ${ \mathcal F }$ summation runs over all elementary species of Fermion in the theory and we have used equation (1) and table 1, leaving the value of the parameter unspecified in the interests of generality. This leads us to make a prediction about physics beyond the minimal standard model of particle physics; as equation (17) holds, our first boundary constraint (16) dictates that the appropriately weighted sum over the electric charges of any elementary species of Fermion yet to be discovered must vanish:

$$\begin{eqnarray}&&\displaystyle \sum _{\displaystyle \genfrac{}{}{0em}{}{{ \mathcal F }}{\left({\rm{y}}{\rm{e}}{\rm{t}}\,{\rm{t}}{\rm{o}}\,{\rm{b}}{\rm{e}}\,{\rm{d}}{\rm{i}}{\rm{s}}{\rm{c}}{\rm{o}}{\rm{v}}{\rm{e}}{\rm{r}}{\rm{e}}{\rm{d}}\right)}}{q}_{{ \mathcal F }}{d}_{{ \mathcal F }}=0.\end{eqnarray} \tag{ 18 }$$

A simple way to satisfy equation (18) is to have each electric charge in the summation be zero.

The possibility that there do indeed exist as-yet-unknown elementary species of Fermion with zero electric charge such as sterile right-handed neutrinos, for example, is appealing as they could serve as components of dark matter [15, 65, 66].

Our second boundary constraint applies if the Universe is closed and reads

$$\begin{eqnarray}&&\langle {\hat{Q}}_{{\mathscr{U}}}\rangle =0,\end{eqnarray} \tag{ 19 }$$

where $\langle {\hat{Q}}_{{\mathscr{U}}}\rangle$ is the total electric charge in the Universe. Equation (19) has already been reported multiple times, albeit with various interpretations; Landau and Lifshitz seem to have held the view that the total of a conserved quantity in a closed Universe is meaningless as the corresponding conservation law is trivial ('0 = 0') [67] and Li recently described equation (19) as a 'paradox' that necessitates modification of Maxwell's equations [68], for example. The identification made in this paper of equation (19) as a boundary constraint does not appear to have been made explicitly before.

Consider the minimal standard model of particle physics augmented with three-dimensional periodic boundary conditions to model a closed Universe. Combining equation (1) with the appropriate statement of our second boundary constraint (19), we obtain

$$\begin{eqnarray}\begin{array}{rcl}\langle {\hat{Q}}_{{\mathscr{U}}}\rangle & = & \sum _{{ \mathcal E }}{q}_{{ \mathcal E }}\langle {\hat{N}}_{{ \mathcal E }}-{\hat{N}}_{\overline{{ \mathcal E }}}\rangle (t)\\ & = & \langle {\hat{Q}}_{{\mathscr{U}}}^{{\prime} }\rangle +\displaystyle \frac{1}{2}\epsilon e\langle {\rm{\Delta }}{\hat{L}}_{{\mathscr{U}}}\rangle \\ & = & 0\end{array}\end{eqnarray} \tag{ 20 }$$

with

$$\begin{eqnarray}&&\langle {\hat{Q}}_{{\mathscr{U}}}^{{\prime} }\rangle =\sum _{{ \mathcal E }}{q}_{{ \mathcal E }}^{{\prime} }\langle {\hat{N}}_{{ \mathcal E }}-{\hat{N}}_{\overline{{ \mathcal E }}}\rangle (t)\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&\langle {\rm{\Delta }}{\hat{L}}_{{\mathscr{U}}}\rangle =\sum _{{ \mathcal E }}({{\ell }}_{m{ \mathcal E }}-{{\ell }}_{n{ \mathcal E }})\langle {\hat{N}}_{{ \mathcal E }}-{\hat{N}}_{\overline{{ \mathcal E }}}\rangle (t),\end{eqnarray} \tag{ 22 }$$

where the ${ \mathcal E }$ summations run over all elementary species in the theory and $\langle {\hat{N}}_{{ \mathcal E }}-{\hat{N}}_{\overline{{ \mathcal E }}}\rangle (t)$ is the total number of particles minus the total number of antiparticles of the ${ \mathcal E }$th species. According to equations (20)–(22), the parameter is not free as is usually understood [24, 25] but is instead fixed as a function of occupation-number differences:

$$\begin{eqnarray}&&\epsilon =-\displaystyle \frac{\langle {\hat{Q}}_{{\mathscr{U}}}^{{\prime} }\rangle }{e\langle {\rm{\Delta }}{\hat{L}}_{{\mathscr{U}}}\rangle },\end{eqnarray} \tag{ 23 }$$

assuming that the difference $\langle {\rm{\Delta }}{\hat{L}}_{{\mathscr{U}}}\rangle$ in total lepton numbers under consideration does not vanish. If one considers a realistic state with occupation-number differences chosen such that $\langle {\hat{Q}}_{{\mathscr{U}}}^{{\prime} }\rangle =0$, equation (23) dictates that

$$\begin{eqnarray*}&&\epsilon =0.\end{eqnarray*}$$

If none of the differences in total lepton number vanish, electric charge must thus be quantised with the familiar values.

Evidently, the combination of classical, quantum and boundary constraints can give rise to electric-charge quantisation with the familiar values in the minimal standard model of particle physics augmented with three-dimensional periodic boundary conditions. This does not appear to have been recognised before.

In this author's view, the preceding argument is interesting but unlikely to be the explanation for electric-charge quantisation. The minimal standard model of particle physics is not our final theory of nature and it seems reasonable to expect that electric-charge quantisation will arise more naturally in a more complete theory. If one extends the minimal standard model of particle physics by adding sterile right-handed neutrinos with Majorana masses, for example, electric-charge quantisation with the familiar values can be explained as a necessary consequence of the relevant classical and quantum constraints only [24, 25]. With electric charge already fixed in such a theory, our second boundary constraint (19) simply restricts the possible occupation-number differences.

The possibility that there do indeed exist sterile right-handed neutrinos with Majorana masses is appealing as they could give rise to the empirically observed neutrino masses via a seesaw mechanism [66].

A key question is whether or not the Universe is closed such that our second and third boundary constraints apply. At present, we cannot answer this definitively, as the shape of the Universe is not known with certainty.

In the minimal ΛCDM model of cosmology, the Universe is taken to be open rather than closed, with zero curvature and a trivial topology [15, 69]. Anomalous features observed in the cosmic microwave background have led some to suggest, however, that the Universe might be closed by virtue of having positive curvature and perhaps even a non-trivial topology [70–73].

It is encouraging to note that the electric charge and current densities appear to vanish on astronomical and grander scales ⁷ ; the apparent dominance of gravitational interactions over electromagnetic interactions in sculpting the Universe at large suggests an electric-charge imbalance of no more than 1 part in 10¹⁸ [74, 75] and more stringent albeit speculative bounds can be claimed based on the apparent isotropy of the Universe [75, 76]. Our boundary constraints offer an immediate theoretical explanation for this empirical observation; if the Universe is closed, it must exhibit such electrical neutrality to obey Maxwell's equations and the boundary conditions. In the absence of other explanations, one might reverse this to interpret the apparent electric neutrality of the Universe as indirect evidence that the Universe is closed.

The possibility that the Universe is indeed electrically neutral is appealing, perhaps, as it aligns with the hypothesis that the Universe is literally an intricate embodiment of nothing [77–79].

We have tacitly assumed above that Maxwell's equations apply to the Universe, in particular that the photon rest mass is zero [1, 80]. There is currently no direct empirical evidence to suggest otherwise, the validity of Maxwell's equations having been extremely well tested [81–83]. Astronomical observations might soon reveal whether or not the photon rest mass is indeed zero [84, 85]. Some pertinent consequences of a non-zero photon rest mass are considered in [68, 74, 86, 87].

The fundamental domain is a cube of volume V = L³ described by time t and right-handed Cartesian coordinates 0 ≤ x ≤ L, 0 ≤ y ≤ L and 0 ≤ z ≤ L with associated unit vectors $\hat{{\bf{x}}}$, $\hat{{\bf{y}}}$ and $\hat{{\bf{z}}}$. The position vector, del operator and Laplacian are

$$\begin{eqnarray*}&&{\bf{r}}=x\hat{{\bf{x}}}+y\hat{{\bf{y}}}+z\hat{{\bf{z}}},\ \ \ {\boldsymbol{\nabla }}=\hat{{\bf{x}}}\displaystyle \frac{\partial }{\partial x}+\hat{{\bf{y}}}\displaystyle \frac{\partial }{\partial y}+\hat{{\bf{z}}}\displaystyle \frac{\partial }{\partial z}\ \ \ \ {{\rm{\nabla }}}^{2}=\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {y}^{2}}+\displaystyle \frac{{\partial }^{2}}{\partial {z}^{2}},\end{eqnarray*}$$

respectively.

The periodic boundary conditions are ensured by only allowing wavevectors of the form

$$\begin{eqnarray*}&&{\bf{k}}=\displaystyle \frac{2\pi }{L}\left(l\hat{{\bf{x}}}+m\hat{{\bf{y}}}+n\hat{{\bf{z}}}\right),\end{eqnarray*}$$

where l, m, n ∈ {0, ± 1,...} are integers.

We consider N_q ∈ {2, 3,...} species of massive and electrically charged matter, labelled ${ \mathcal F }\in \{1,...,{N}_{q}\}$. For the ${ \mathcal F }$th species, we identify Fermionic annihilation operators ${\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}$ and creation operators ${\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}^{\dagger }$, where τ = + for particles, τ = − for antiparticles, ℏ k is an eigenvalue of linear momentum and ℏ m_s ∈ {ℏ/2, − ℏ/2} is an eigenvalue of the z component of mean spin. The ${\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}$ and ${\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}^{\dagger }$ satisfy the usual Fermionic anticommutation relations, in particular

$$\begin{eqnarray*}&&{\left[{\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}},{\hat{c}}_{{ \mathcal F }\tau ^{\prime} {\bf{k}}^{\prime} {m}_{s}^{\prime} }^{\dagger }\right]}_{+}={\delta }_{\tau \tau ^{\prime} }{\delta }_{{\bf{kk}}^{\prime} }{\delta }_{{m}_{s}{m}_{s}^{\prime} }.\end{eqnarray*}$$

The corresponding matter occupation numbers are

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{N}}_{{ \mathcal F }{\bf{k}}{m}_{s}} & = & {\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}^{\dagger }{\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}},\ \ \ {\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}{m}_{s}}={\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}^{\dagger }{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}},\\ {\hat{N}}_{{ \mathcal F }{\bf{k}}} & = & \sum _{{m}_{s}=\pm 1/2}{\hat{N}}_{{ \mathcal F }{\bf{k}}{m}_{s}},\ \ \ {\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}}=\sum _{{m}_{s}=\pm 1/2}{\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}{m}_{s}},\\ {\hat{N}}_{{ \mathcal F }} & = & \sum _{{\bf{k}}}{\hat{N}}_{{ \mathcal F }{\bf{k}}}\ \ \ \ {\hat{N}}_{\overline{{ \mathcal F }}}=\sum _{{\bf{k}}}{\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}}.\end{array}\end{eqnarray*}$$

For the electromagnetic field, we identify Bosonic annihilation operators ${\hat{a}}_{{\bf{k}}\alpha }$ and creation operators ${\hat{a}}_{{\bf{k}}\alpha }^{\dagger }$, where k ≠ 0 is a non-zero allowed wavevector, α = t for scalar photons and α = i ∈ {x, y, z} for Cartesian photons. The ${\hat{a}}_{{\bf{k}}\alpha }$ and ${\hat{a}}_{{\bf{k}}\alpha }^{\dagger }$ satisfy the usual Bosonic commutation relations, in particular

$$\begin{eqnarray*}&&\left[{\hat{a}}_{{\bf{k}}\alpha },{\hat{a}}_{{\bf{k}}^{\prime} \alpha ^{\prime} }^{\dagger }\right]={\delta }_{{\bf{kk}}^{\prime} }{\delta }_{\alpha \alpha ^{\prime} }.\end{eqnarray*}$$

The corresponding electromagnetic occupation numbers are

$$\begin{eqnarray*}&&{\hat{N}}_{\gamma {\bf{k}}\alpha }={\hat{a}}_{{\bf{k}}\alpha }^{\dagger }{\hat{a}}_{{\bf{k}}\alpha }\ \ \ \ {\hat{N}}_{\gamma {\bf{k}}}=\sum _{\alpha =t,x,y,z}{\hat{N}}_{\gamma {\bf{k}}\alpha }\ \ \ \ {\hat{N}}_{\gamma }=\sum _{{\bf{k}}\ne 0}{\hat{N}}_{\gamma {\bf{k}}}.\end{eqnarray*}$$

The Hilbert space is spanned by Fock states of the form

$$\begin{eqnarray*}&&| {1}_{A},{1}_{B},...,{n}_{I},{n}_{J},...\rangle ={\hat{c}}_{A}^{\dagger }{\hat{c}}_{B}^{\dagger }...\displaystyle \frac{{\left({\hat{a}}_{I}^{\dagger }\right)}^{{n}_{I}}{\left({\hat{a}}_{J}^{\dagger }\right)}^{{n}_{J}}...}{\sqrt{{n}_{I}!{n}_{J}!...}}| 0\rangle ,\end{eqnarray*}$$

where A, B,... stand for distinct matter modes (${ \mathcal F }\tau {\bf{k}}{m}_{s}$); I, J,... stand for distinct electromagnetic modes (k α); n_I , n_J ,... ∈ {0, 1,...} are electromagnetic-occupation-number eigenvalues and ∣0〉 is the unperturbed vacuum state, for which

$$\begin{eqnarray*}&&{\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}| 0\rangle =0\ \ \ \ {\hat{a}}_{{\bf{k}}\alpha }| 0\rangle =0.\end{eqnarray*}$$

The Hilbert space is equipped with the usual inner product

$$\begin{eqnarray*}&&\langle \psi | \phi \rangle =\langle \phi | \psi {\rangle }^{* },\end{eqnarray*}$$

where ∣ψ〉 and ∣ϕ〉 are arbitrary kets and 〈ψ∣ and 〈ϕ∣ are the corresponding bras. The usual Hermitian conjugate ${\hat{X}}^{\dagger }$ of an arbitrary operator $\hat{X}$ is defined such that

$$\begin{eqnarray*}&&\langle \psi | {\hat{X}}^{\dagger }| \phi \rangle =\langle \phi | \hat{X}| \psi {\rangle }^{* }.\end{eqnarray*}$$

The new ket ∣ψ⟫ and corresponding bra ⟪ψ∣ associated with an arbitrary ket ∣ψ〉 and corresponding bra 〈ψ∣ are

$$\begin{eqnarray*}&&| \psi \rangle\rangle =| \psi \rangle \ \ \ \ \langle\langle \psi | =\langle \psi | \hat{M},\end{eqnarray*}$$

where $\hat{M}$ is the indefinite metric, acting on the subspace of the Hilbert space pertaining to scalar photons as

$$\begin{eqnarray*}&&\hat{M}| {1}_{A},{1}_{B},\ldots ,{n}_{I},{n}_{J},...\rangle ={\left(-1\right)}^{{\sum }_{{\bf{k}}\ne 0}{n}_{{\bf{k}}t}}| {1}_{A},{1}_{B},\ldots ,{n}_{I},{n}_{J},...\rangle .\end{eqnarray*}$$

Note that $\hat{M}={\hat{M}}^{\dagger }={\hat{M}}^{-1}$ is Hermitian in the usual sense and unitary in the usual sense.

The new Hermitian conjugate ${\hat{X}}^{\unicode{x02021}}$ of an arbitrary operator $\hat{X}$ is defined such that

$$\begin{eqnarray*} &&\langle\langle \psi | {\hat{X}}^{\unicode{x02021}}| \phi \rangle\rangle = \langle\langle \phi | \hat{X}| \psi { \rangle\rangle }^{* }.\end{eqnarray*}$$

Note that the new Hermitian conjugate coincides with the usual Hermitian conjugate for the matter creation operators ${\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}^{\dagger }={\hat{c}}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}^{\unicode{x02021}}$ and the Cartesian-photon creation operators ${\hat{a}}_{{\bf{k}}i}^{\dagger }={\hat{a}}_{{\bf{k}}i}^{\unicode{x02021}}$ but not the scalar-photon creation operators ${\hat{a}}_{{\bf{k}}t}^{\dagger }=-{\hat{a}}_{{\bf{k}}t}^{\unicode{x02021}}$.

If the system occupies a state ∣Ψ(t)〉, the new mean of an arbitrary operator $\hat{X}$ is

$$\begin{eqnarray*}&&\langle \hat{X}\rangle (t)=\displaystyle \frac{ \langle\langle {\rm{\Psi }}(t)| \hat{X}| {\rm{\Psi }}(t) \rangle\rangle }{ \langle\langle {\rm{\Psi }}(t)| {\rm{\Psi }}(t) \rangle\rangle },\end{eqnarray*}$$

assuming that the new norm ⟪Ψ(t)∣Ψ(t)⟫ of ∣Ψ(t)〉 is non-zero.

We consider the Dirac matrices

$$\begin{eqnarray*}\beta =\left(\begin{array}{cc}{1}_{2\times 2} & {0}_{2\times 2}\\ {0}_{2\times 2} & -{1}_{2\times 2}\end{array}\right),\ \ \ {\boldsymbol{\alpha }}=\left(\begin{array}{cc}{0}_{2\times 2} & {\boldsymbol{\sigma }}\\ {\boldsymbol{\sigma }} & {0}_{2\times 2}\end{array}\right)\ \ \ \ {\boldsymbol{\Sigma }}=\left(\begin{array}{cc}{\boldsymbol{\sigma }} & {0}_{2\times 2}\\ {0}_{2\times 2} & {\boldsymbol{\sigma }}\end{array}\right),\end{eqnarray*}$$

where

$$\begin{eqnarray*}{\boldsymbol{\sigma }}=\left(\begin{array}{cc}0 & 1\\ 1 & 0\end{array}\right)\hat{{\bf{x}}}+\left(\begin{array}{cc}0 & -{\rm{i}}\\ {\rm{i}} & 0\end{array}\right)\hat{{\bf{y}}}+\left(\begin{array}{cc}1 & 0\\ 0 & -1\end{array}\right)\hat{{\bf{z}}}\end{eqnarray*}$$

is a vector of Pauli matrices.

The ${ \mathcal F }$th species of matter is embodied by the Dirac-type field

$$\begin{eqnarray*}&&{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}})=\sum _{{\bf{k}}}\sum _{{m}_{s}=\pm 1/2}({\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{{ \mathcal F }+{\bf{k}}{m}_{s}}+{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{{ \mathcal F }-{\bf{k}}-{m}_{s}})\end{eqnarray*}$$

with rest mass ${m}_{{ \mathcal F }}$ and electric charge ${q}_{{ \mathcal F }}$, where the

$$\begin{eqnarray*}&&{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}=\displaystyle \frac{1}{\sqrt{V}}\exp \left\{-{\rm{i}}\left[-\displaystyle \frac{{\rm{i}}\beta {\boldsymbol{\alpha }}\cdot {\bf{k}}}{2| {\bf{k}}| }{\tan }^{-1}\left(\displaystyle \frac{{\hslash }| {\bf{k}}| }{{{cm}}_{{ \mathcal F }}}\right)\right]\right\}{\left[{\delta }_{\tau +}{\delta }_{{m}_{s}1/2},{\delta }_{\tau +}{\delta }_{{m}_{s}-1/2},{\delta }_{\tau -}{\delta }_{{m}_{s}1/2},{\delta }_{\tau -}{\delta }_{{m}_{s}-1/2}\right]}^{{\rm{T}}}\end{eqnarray*}$$

are spinors corresponding to modes of definite energy, linear momentum and mean-spin projection along the z axis in that

$$\begin{eqnarray*}\begin{array}{rcl}{h}_{{ \mathcal F }}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}} & = & \tau \sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}},\ \ \ {\bf{p}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}={\hslash }{\bf{k}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}\\ {s}_{{ \mathcal F }z}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}} & = & {\hslash }{m}_{s}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}},\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{h}_{{ \mathcal F }} & = & -{\rm{i}}{\hslash }c{\boldsymbol{\alpha }}\cdot {\boldsymbol{\nabla }}+{c}^{2}{m}_{{ \mathcal F }}\beta ,\ \ \ {\bf{p}}=-{\rm{i}}{\hslash }{\boldsymbol{\nabla }}\\ {{\bf{s}}}_{{ \mathcal F }} & = & \displaystyle \frac{{\hslash }}{2}{\boldsymbol{\Sigma }}-\displaystyle \frac{{{\hslash }}^{2}c\beta {\boldsymbol{\alpha }}\times {\boldsymbol{\nabla }}}{2\sqrt{-{{\hslash }}^{2}{c}^{2}{{\rm{\nabla }}}^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}}+\displaystyle \frac{{{\hslash }}^{3}{c}^{2}{\boldsymbol{\nabla }}\times ({\boldsymbol{\Sigma }}\times {\boldsymbol{\nabla }})}{2\sqrt{-{{\hslash }}^{2}{c}^{2}{{\rm{\nabla }}}^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}(\sqrt{-{{\hslash }}^{2}{c}^{2}{{\rm{\nabla }}}^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}+{c}^{2}{m}_{{ \mathcal F }})}\end{array}\end{eqnarray*}$$

are the first-quantised free-field Dirac Hamiltonian, linear momentum and mean spin, respectively. Note that the ${u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}$ satisfy the orthonormality condition

$$\begin{eqnarray*}&&{\iiint }_{{\mathscr{V}}}{u}_{{ \mathcal F }\tau {\bf{k}}{m}_{s}}^{\dagger }{u}_{{ \mathcal F }\tau ^{\prime} {\bf{k}}^{\prime} {m}_{s}^{\prime} }{{\rm{e}}}^{{\rm{i}}({\bf{k}}^{\prime} -{\bf{k}})\cdot {\bf{r}}}{{\rm{d}}}^{3}{\bf{r}}={\delta }_{\tau \tau ^{\prime} }{\delta }_{{\bf{kk}}^{\prime} }{\delta }_{{m}_{s}{m}_{s}^{\prime} }.\end{eqnarray*}$$

The electric charge and current densities are

$$\begin{eqnarray*}&&\hat{\rho }({\bf{r}})=\sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}^{\dagger }({\bf{r}}){\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}})\ \ \ \ \hat{{\bf{J}}}({\bf{r}})=\sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}^{\dagger }({\bf{r}}){\bf{v}}{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}}),\end{eqnarray*}$$

respectively, where

$$\begin{eqnarray*}&&{\bf{v}}=c{\boldsymbol{\alpha }}\end{eqnarray*}$$

is the first-quantised Dirac velocity. The total electric charge and total electric current in the fundamental domain ${\mathscr{V}}$ are

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{Q}}_{{\mathscr{V}}} & = & {\iiint }_{{\mathscr{V}}}\hat{\rho }({\bf{r}}){{\rm{d}}}^{3}{\bf{r}}\\ & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}\sum _{{m}_{s}=\pm 1/2}{q}_{{ \mathcal F }}({\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}^{\dagger }{\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}+{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}^{\dagger })\\ & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}{q}_{{ \mathcal F }}({\hat{N}}_{{ \mathcal F }}-{\hat{N}}_{\overline{{ \mathcal F }}}+\sum _{{\bf{k}}}2)\\ {\hat{{\bf{K}}}}_{{\mathscr{V}}} & = & {\iiint }_{{\mathscr{V}}}\hat{{\bf{J}}}({\bf{r}}){{\rm{d}}}^{3}{\bf{r}}\\ & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}\sum _{{m}_{s}=\pm 1/2}{q}_{{ \mathcal F }}{{\bf{v}}}_{{ \mathcal F }{\bf{k}}}({\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}^{\dagger }{\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}+{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}^{\dagger })+\hat{{\boldsymbol{\Delta }}},\\ & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}{q}_{{ \mathcal F }}{{\bf{v}}}_{{ \mathcal F }{\bf{k}}}({\hat{N}}_{{ \mathcal F }{\bf{k}}}-{\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}})+\hat{{\boldsymbol{\Delta }}},\end{array}\end{eqnarray*}$$

respectively, where

$$\begin{eqnarray*}&&{{\bf{v}}}_{{ \mathcal F }{\bf{k}}}=\displaystyle \frac{{\hslash }{c}^{2}{\bf{k}}}{\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}}\end{eqnarray*}$$

is a mean-velocity eigenvalue and

$$\begin{eqnarray*}&&\begin{array}{l}\hat{{\boldsymbol{\Delta }}}=\sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}\displaystyle \frac{{{cq}}_{{ \mathcal F }}}{\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}(\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}+{c}^{2}{m}_{{ \mathcal F }})}\\ \times \,\left\{\left[\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}(\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}+{c}^{2}{m}_{{ \mathcal F }})\hat{{\bf{x}}}-{c}^{2}{{\hslash }}^{2}{k}_{x}{\bf{k}}\right]\right.\\ \times \,({\hat{c}}_{{ \mathcal F }+{\bf{k}}1/2}^{\dagger }{\hat{c}}_{{ \mathcal F }--{\bf{k}}1/2}^{\dagger }+{\hat{c}}_{{ \mathcal F }+{\bf{k}}-1/2}^{\dagger }{\hat{c}}_{{ \mathcal F }--{\bf{k}}-1/2}^{\dagger })\\ +\left[\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}(\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}+{c}^{2}{m}_{{ \mathcal F }})\hat{{\bf{y}}}-{c}^{2}{{\hslash }}^{2}{k}_{y}{\bf{k}}\right]\\ \times \,(-{\rm{i}}{\hat{c}}_{{ \mathcal F }+{\bf{k}}1/2}^{\dagger }{\hat{c}}_{{ \mathcal F }--{\bf{k}}1/2}^{\dagger }+{\rm{i}}{\hat{c}}_{{ \mathcal F }+{\bf{k}}-1/2}^{\dagger }{\hat{c}}_{{ \mathcal F }--{\bf{k}}-1/2}^{\dagger })\\ +\left[\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}(\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}+{c}^{2}{m}_{{ \mathcal F }})\hat{{\bf{z}}}-{c}^{2}{{\hslash }}^{2}{k}_{z}{\bf{k}}\right]\\ \times \,\left.({\hat{c}}_{{ \mathcal F }+{\bf{k}}1/2}^{\dagger }{\hat{c}}_{{ \mathcal F }--{\bf{k}}-1/2}^{\dagger }-{\hat{c}}_{{ \mathcal F }+{\bf{k}}-1/2}^{\dagger }{\hat{c}}_{{ \mathcal F }--{\bf{k}}1/2}^{\dagger })+\mathrm{usualHermitianconjugate}\right\}\end{array}\end{eqnarray*}$$

is due to Zitterbewegung.

For each non-zero allowed wavevector k ≠ 0, we identify the rescaled reciprocal-space Fourier components

$$\begin{eqnarray*}&&{\hat{\lambda }}_{{\bf{k}}t}=\displaystyle \frac{1}{| {\bf{k}}| }\sqrt{\displaystyle \frac{1}{2{\hslash }{\epsilon }_{0}{cV}| {\bf{k}}| }}{\iiint }_{{\mathscr{V}}}\hat{\rho }({\bf{r}}){{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{{\rm{d}}}^{3}{\bf{r}}\ \ \ \ {\hat{{\boldsymbol{\lambda }}}}_{{\bf{k}}}=\displaystyle \frac{1}{c| {\bf{k}}| }\sqrt{\displaystyle \frac{1}{2{\hslash }{\epsilon }_{0}{cV}| {\bf{k}}| }}{\iiint }_{{\mathscr{V}}}\hat{{\bf{J}}}({\bf{r}}){{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}{{\rm{d}}}^{3}{\bf{r}}.\end{eqnarray*}$$

${\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}})$ obeys the Dirac equation

$$\begin{eqnarray*}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial {\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}})}{\partial t}=\left\{{h}_{{ \mathcal F }}+{q}_{{ \mathcal F }}\left[\hat{{\rm{\Phi }}}({\bf{r}})-{\bf{v}}\cdot \hat{{\bf{A}}}({\bf{r}})\right]\right\}{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}}).\end{eqnarray*}$$

The electromagnetic field is embodied by the potential four-vector

$$\begin{eqnarray*}&&{A}^{\alpha }({\bf{r}})=(\hat{{\rm{\Phi }}}({\bf{r}})/c,\hat{{\bf{A}}}({\bf{r}})),\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&\hat{{\rm{\Phi }}}({\bf{r}})=\sum _{{\bf{k}}\ne 0}\sqrt{\displaystyle \frac{{\hslash }c}{2{\epsilon }_{0}V| {\bf{k}}| }}({\hat{a}}_{{\bf{k}}t}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}-{\hat{a}}_{{\bf{k}}t}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}})\ \ \ \ \hat{{\bf{A}}}({\bf{r}})=\sum _{{\bf{k}}\ne 0}\sqrt{\displaystyle \frac{{\hslash }}{2{\epsilon }_{0}{cV}| {\bf{k}}| }}({\hat{{\bf{a}}}}_{{\bf{k}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}+{\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}})\end{eqnarray*}$$

are the scalar and vector potentials, respectively. Note that $\hat{{\rm{\Phi }}}({\bf{r}})={\hat{{\rm{\Phi }}}}^{\unicode{x02021}}({\bf{r}})$ is Hermitian in the new sense. The four-divergence of A^α (r) is

$$\begin{eqnarray*}\begin{array}{rcl}\hat{{\rm{\Lambda }}}({\bf{r}}) & = & \displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{\partial \hat{{\rm{\Phi }}}({\bf{r}})}{\partial t}+{\boldsymbol{\nabla }}\cdot \hat{{\bf{A}}}({\bf{r}})\\ & = & {\hat{{\rm{\Lambda }}}}^{(+)}({\bf{r}})+{\hat{{\rm{\Lambda }}}}^{(-)}({\bf{r}}),\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{{\rm{\Lambda }}}}^{(+)}({\bf{r}}) & = & \sum _{{\bf{k}}\ne 0}\sqrt{\displaystyle \frac{{\hslash }}{2{\epsilon }_{0}{cV}| {\bf{k}}| }}{\rm{i}}\left[| {\bf{k}}| ({\hat{\lambda }}_{{\bf{k}}t}-{\hat{a}}_{{\bf{k}}t})+{\bf{k}}\cdot {\hat{{\bf{a}}}}_{{\bf{k}}}\right]{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}\\ {\hat{{\rm{\Lambda }}}}^{(-)}({\bf{r}}) & = & -\sum _{{\bf{k}}\ne 0}\sqrt{\displaystyle \frac{{\hslash }}{2{\epsilon }_{0}{cV}| {\bf{k}}| }}{\rm{i}}\left[| {\bf{k}}| ({\hat{\lambda }}_{{\bf{k}}t}^{\dagger }+{\hat{a}}_{{\bf{k}}t}^{\dagger })+{\bf{k}}\cdot {\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger }\right]{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}\end{array}\end{eqnarray*}$$

are the positive- and negative-frequency parts, respectively. Note that $\hat{{\rm{\Lambda }}}({\bf{r}})={\hat{{\rm{\Lambda }}}}^{\unicode{x02021}}({\bf{r}})$ is Hermitian in the new sense and that the reciprocal-space Fourier components ${\hat{\lambda }}_{{\bf{k}}t}={\hat{\lambda }}_{-{\bf{k}}t}^{\dagger }$ make no overall contribution to $\hat{{\rm{\Lambda }}}({\bf{r}})$, as their contributions to ${\hat{{\rm{\Lambda }}}}^{(+)}({\bf{r}})$ and ${\hat{{\rm{\Lambda }}}}^{(-)}({\bf{r}})$ cancel.

The electric and magnetic fields are

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\bf{E}}}({\bf{r}}) & = & -{\boldsymbol{\nabla }}\hat{{\rm{\Phi }}}({\bf{r}})-\displaystyle \frac{\partial \hat{{\bf{A}}}({\bf{r}})}{\partial t}\\ & = & \sum _{{\bf{k}}\ne 0}\sqrt{\displaystyle \frac{{\hslash }c}{2{\epsilon }_{0}V| {\bf{k}}| }}{\rm{i}}\left[(-{\bf{k}}{\hat{a}}_{{\bf{k}}t}+| {\bf{k}}| {\hat{{\bf{a}}}}_{{\bf{k}}}){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}-({\bf{k}}{\hat{a}}_{{\bf{k}}t}^{\dagger }+| {\bf{k}}| {\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger }){{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}\right]\end{array}\end{eqnarray} \tag{ A1 }$$

$$\begin{eqnarray}\begin{array}{rcl}\hat{{\bf{B}}}({\bf{r}}) & = & {\boldsymbol{\nabla }}\times \hat{{\bf{A}}}({\bf{r}})\\ & = & \sum _{{\bf{k}}\ne 0}\sqrt{\displaystyle \frac{{\hslash }}{2{\epsilon }_{0}{cV}| {\bf{k}}| }}{\rm{i}}{\bf{k}}\times ({\hat{{\bf{a}}}}_{{\bf{k}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}-{\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger }{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}}),\end{array}\end{eqnarray} \tag{ A2 }$$

respectively. Note that $\hat{{\bf{E}}}({\bf{r}})={\hat{{\bf{E}}}}^{\unicode{x02021}}({\bf{r}})$ is Hermitian in the new sense.

Let us emphasise here that we are working in the Gupta-Bleuler (Lorenz-gauge) formalism [36, 37]. In this formalism, Maxwell's equations do not all hold in terms of operators but are expected to hold with respect to the new mean. Equations (A1) and (A2) ensure the validity of Gauss's law for magnetism and the Faraday-Lenz law

$$\begin{eqnarray*}&&{\boldsymbol{\nabla }}\cdot \hat{{\bf{B}}}({\bf{r}})=0\ \ \ \ {\boldsymbol{\nabla }}\times \hat{{\bf{E}}}({\bf{r}})=-\displaystyle \frac{\partial \hat{{\bf{B}}}({\bf{r}})}{\partial t},\end{eqnarray*}$$

respectively. Instead of Gauss's law and the Ampère-Maxwell law, however, we obtain

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \hat{{\bf{E}}}({\bf{r}})=\displaystyle \frac{\hat{\rho }({\bf{r}})}{{\epsilon }_{0}}-\displaystyle \frac{{\hat{Q}}_{{\mathscr{V}}}}{{\epsilon }_{0}V}-\displaystyle \frac{\partial \hat{{\rm{\Lambda }}}({\bf{r}})}{\partial t}\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\times \hat{{\bf{B}}}({\bf{r}})={\mu }_{0}\hat{{\bf{J}}}({\bf{r}})-\displaystyle \frac{{\mu }_{0}{\hat{{\bf{K}}}}_{{\mathscr{V}}}}{V}+{\boldsymbol{\nabla }}\hat{{\rm{\Lambda }}}({\bf{r}})+\displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{\partial \hat{{\bf{E}}}({\bf{r}})}{\partial t},\end{eqnarray} \tag{ A4 }$$

respectively. Note that the '${\hat{Q}}_{{\mathscr{V}}}$' and '${\hat{{\bf{K}}}}_{{\mathscr{V}}}$' terms in equations (A3) and (A4) emerge only when we are careful enough to exclude electromagnetic modes with wavevector k = 0, as

$$\begin{eqnarray*}&&\begin{array}{l}\displaystyle \frac{1}{{\epsilon }_{0}V}\sum _{{\bf{k}}\ne 0}{\iiint }_{{\mathscr{V}}}\hat{\rho }({\bf{r}}^{\prime} ){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot ({\bf{r}}-{\bf{r}}^{\prime} )}{{\rm{d}}}^{3}{\bf{r}}^{\prime} =\displaystyle \frac{1}{{\epsilon }_{0}\sqrt{V}}\sum _{{\bf{k}}}\left[\displaystyle \frac{1}{\sqrt{V}}{\iiint }_{{\mathscr{V}}}\hat{\rho }({\bf{r}}^{\prime} ){{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}^{\prime} }{{\rm{d}}}^{3}{\bf{r}}^{\prime} \right]{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}\\ -\displaystyle \frac{1}{{\epsilon }_{0}V}{\iiint }_{{\mathscr{V}}}\hat{\rho }({\bf{r}}^{\prime} ){{\rm{d}}}^{3}{\bf{r}}^{\prime} \\ =\displaystyle \frac{\hat{\rho }({\bf{r}})}{{\epsilon }_{0}}-\displaystyle \frac{{\hat{Q}}_{{\mathscr{V}}}}{{\epsilon }_{0}V}\\ \displaystyle \frac{{\mu }_{0}}{V}\sum _{{\bf{k}}\ne 0}{\iiint }_{V}\hat{{\bf{J}}}({\bf{r}}^{\prime} ){{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot ({\bf{r}}-{\bf{r}}^{\prime} )}{{\rm{d}}}^{3}{\bf{r}}^{\prime} =\displaystyle \frac{{\mu }_{0}}{\sqrt{V}}\sum _{{\bf{k}}}\left[\displaystyle \frac{1}{\sqrt{V}}{\iiint }_{{\mathscr{V}}}\hat{{\bf{J}}}({\bf{r}}^{\prime} ){{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{r}}^{\prime} }{{\rm{d}}}^{3}{\bf{r}}^{\prime} \right]{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{r}}}\\ -\displaystyle \frac{{\mu }_{0}}{V}{\iiint }_{{\mathscr{V}}}\hat{{\bf{J}}}({\bf{r}}^{\prime} ){{\rm{d}}}^{3}{\bf{r}}^{\prime} \\ ={\mu }_{0}\hat{{\bf{J}}}({\bf{r}})-\displaystyle \frac{{\mu }_{0}{\hat{{\bf{K}}}}_{{\mathscr{V}}}}{V}.\end{array}\end{eqnarray*}$$

It follows from equation (A4) that the total electric current through the surface ${\mathscr{S}}^{\prime}$ satisfies

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{I}}_{{\mathscr{S}}^{\prime} } & = & {\iint }_{{\mathscr{S}}^{\prime} }\hat{{\bf{J}}}({\bf{r}})\cdot {{\rm{d}}}^{2}{\bf{r}}\\ & = & \displaystyle \frac{1}{V}{\iint }_{{\mathscr{S}}^{\prime} }{\hat{{\bf{K}}}}_{{\mathscr{V}}}\cdot {{\rm{d}}}^{2}{\bf{r}}-\displaystyle \frac{1}{{\mu }_{0}}{\iint }_{{\mathscr{S}}^{\prime} }{\boldsymbol{\nabla }}\hat{{\rm{\Lambda }}}({\bf{r}})\cdot {{\rm{d}}}^{2}{\bf{r}}-{\epsilon }_{0}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\iint }_{{\mathscr{S}}^{\prime} }\hat{{\bf{E}}}({\bf{r}})\cdot {{\rm{d}}}^{2}{\bf{r}}.\end{array}\end{eqnarray*}$$

The system is governed by the Hamiltonian

$$\begin{eqnarray*}&&\hat{H}={\hat{H}}_{\mathrm{matter}}+{\hat{H}}_{\mathrm{EMfield}}+{\hat{H}}_{\mathrm{EMinteractions}},\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{H}}_{\mathrm{matter}} & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}{\iiint }_{{\mathscr{V}}}{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}^{\dagger }({\bf{r}}){h}_{{ \mathcal F }}{\hat{{\rm{\Psi }}}}_{{ \mathcal F }}({\bf{r}}){{\rm{d}}}^{3}{\bf{r}}\\ & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}\sum _{{m}_{s}=\pm 1/2}\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}({\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}^{\dagger }{\hat{c}}_{{ \mathcal F }+{\bf{k}}{m}_{s}}+{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}{\hat{c}}_{{ \mathcal F }-{\bf{k}}{m}_{s}}^{\dagger })\\ & = & \sum _{{ \mathcal F }=1}^{{N}_{q}}\sum _{{\bf{k}}}\sqrt{{{\hslash }}^{2}{c}^{2}| {\bf{k}}{| }^{2}+{c}^{4}{m}_{{ \mathcal F }}^{2}}({\hat{N}}_{{ \mathcal F }{\bf{k}}}-{\hat{N}}_{\overline{{ \mathcal F }}{\bf{k}}}-2),\\ {\hat{H}}_{\mathrm{EMfield}} & = & {\iiint }_{{\mathscr{V}}}\displaystyle \frac{1}{2}\left[-\displaystyle \frac{{\epsilon }_{0}}{{c}^{2}}\displaystyle \frac{\partial \hat{{\rm{\Phi }}}({\bf{r}})}{\partial t}\displaystyle \frac{\partial \hat{{\rm{\Phi }}}({\bf{r}})}{\partial t}+{\epsilon }_{0}\displaystyle \frac{\partial \hat{{\bf{A}}}({\bf{r}})}{\partial t}\cdot \displaystyle \frac{\partial \hat{{\bf{A}}}({\bf{r}})}{\partial t}-{\epsilon }_{0}{\boldsymbol{\nabla }}\hat{{\rm{\Phi }}}({\bf{r}})\cdot {\boldsymbol{\nabla }}\hat{{\rm{\Phi }}}({\bf{r}})\right.\\ & & \left.+\sum _{i=x,y,z}\sum _{j=x,y,z}\displaystyle \frac{1}{{\mu }_{0}}{\partial }_{i}{\hat{A}}_{j}({\bf{r}}){\partial }_{i}{\hat{A}}_{j}({\bf{r}})\right]{{\rm{d}}}^{3}{\bf{r}}\\ & = & \sum _{{\bf{k}}\ne 0}\displaystyle \frac{1}{2}{\hslash }c| {\bf{k}}| ({\hat{a}}_{{\bf{k}}t}^{\dagger }{\hat{a}}_{{\bf{k}}t}+{\hat{a}}_{{\bf{k}}t}{\hat{a}}_{{\bf{k}}t}^{\dagger }+{\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger }\cdot {\hat{{\bf{a}}}}_{{\bf{k}}}+{\hat{{\bf{a}}}}_{{\bf{k}}}\cdot {\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger })\\ & = & \sum _{{\bf{k}}\ne 0}{\hslash }c| {\bf{k}}| ({\hat{N}}_{\gamma {\bf{k}}}+2)\\ {\hat{H}}_{\mathrm{EMinteractions}} & = & {\iiint }_{{\mathscr{V}}}\left[\hat{\rho }({\bf{r}})\hat{{\rm{\Phi }}}({\bf{r}})-\hat{{\bf{J}}}({\bf{r}})\cdot \hat{{\bf{A}}}({\bf{r}})\right]{{\rm{d}}}^{3}{\bf{r}}\\ & = & \sum _{{\bf{k}}\ne 0}{\hslash }c| {\bf{k}}| ({\hat{\lambda }}_{{\bf{k}}t}^{\dagger }{\hat{a}}_{{\bf{k}}t}-{\hat{\lambda }}_{{\bf{k}}t}{\hat{a}}_{{\bf{k}}t}^{\dagger }-{\hat{{\boldsymbol{\lambda }}}}_{{\bf{k}}}^{\dagger }\cdot {\hat{{\bf{a}}}}_{{\bf{k}}}-{\hat{{\boldsymbol{\lambda }}}}_{{\bf{k}}}\cdot {\hat{{\bf{a}}}}_{{\bf{k}}}^{\dagger })\end{array}\end{eqnarray*}$$

describe the matter, the electromagnetic field and electromagnetic interactions, respectively. Note that $\hat{H}={\hat{H}}^{\unicode{x02021}}$ is Hermitian in the new sense.

The state ∣Ψ(t)〉 of the system evolves according to Schrödinger's equation

$$\begin{eqnarray*}&&{\rm{i}}{\hslash }\displaystyle \frac{{\rm{d}}| {\rm{\Psi }}(t)\rangle }{{\rm{d}}t}=\hat{H}| {\rm{\Psi }}(t)\rangle .\end{eqnarray*}$$

The time derivative ${\rm{d}}\hat{X}/{\rm{d}}t$ of an arbitrary operator $\hat{X}$ with no explicit time dependence is given by Heisenberg's equation of motion

$$\begin{eqnarray*}&&\displaystyle \frac{{\rm{d}}\hat{X}}{{\rm{d}}t}=\displaystyle \frac{{\rm{i}}}{{\hslash }}\left[\hat{H},\hat{X}\right].\end{eqnarray*}$$

To recover all of Maxwell's equations, we must focus on the 'physical' subspace of the Hilbert space, defined such that

$$\begin{eqnarray*}&&{\hat{{\rm{\Lambda }}}}^{(+)}({\bf{r}})| {\psi }_{\mathrm{physical}}\rangle =0,\end{eqnarray*}$$

where ∣ψ_physical〉 is an arbitary physical ket. If the system occupies an arbitrary state ∣Ψ_physical(t)〉 comprised solely of the ∣ψ_physical〉, the new means of equations (A3) and (A4) are

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)=\displaystyle \frac{\langle \hat{\rho }\rangle ({\bf{r}},t)}{{\epsilon }_{0}}-\displaystyle \frac{\langle {\hat{Q}}_{{\mathscr{V}}}\rangle }{{\epsilon }_{0}V}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\times \langle \hat{{\bf{B}}}\rangle ({\bf{r}},t)={\mu }_{0}\langle \hat{{\bf{J}}}\rangle ({\bf{r}},t)-\displaystyle \frac{{\mu }_{0}\langle {\hat{{\bf{K}}}}_{{\mathscr{V}}}\rangle (t)}{V}+\displaystyle \frac{1}{{c}^{2}}\displaystyle \frac{\partial \langle \hat{{\bf{E}}}\rangle ({\bf{r}},t)}{\partial t},\end{eqnarray} \tag{ A6 }$$

respectively, assuming that the new norm ⟪Ψ_physical(t)∣Ψ_physical(t)⟫ of ∣Ψ_physical(t)〉 is non-zero. Equations (A5) and (A6) reduce to Gauss's law and the Ampère-Maxwell law, respectively, if and only if

$$\begin{eqnarray*}&&\langle {\hat{Q}}_{{\mathscr{V}}}\rangle =0\ \ \ \ \langle {\hat{{\bf{K}}}}_{{\mathscr{V}}}\rangle =0,\end{eqnarray*}$$

from which our boundary constraints (5), (6) and (11) follow.