# Dynamic Manipulation of a Deformable Linear Object: Simulation and Learning

Chen Qi Jing1, Timothy Bretl2, Nghia Vuong1, and Quang-Cuong Pham1,3 *Github
repo: https://github.com/qj25/der_mujoco.git (Videos:
https://youtu.be/Icb7Nl6voN0)1Nanyang Technological University, School of
Mechanical and Aerospace Engineering, 2University of Illinois at Urbana-
Champaign, 3Eureka Robotics, Singapore

###### Abstract

We show that it is possible to learn an open-loop policy in simulation for the
dynamic manipulation of a deformable linear object (DLO) – e.g., a rope, wire,
or cable – that can be executed by a real robot without additional training.
Our method is enabled by integrating an existing state-of-the-art DLO model
(Discrete Elastic Rods) with MuJoCo, a robot simulator. We describe how this
integration was done, check that validation results produced in simulation
match what we expect from analysis of the physics, and apply policy
optimization to train an open-loop policy from data collected only in
simulation that uses a robot arm to fling a wire precisely between two
obstacles. This policy achieves a success rate of 76.7% when executed by a
real robot in hardware experiments without additional training on the real
task.

##  I INTRODUCTION

The robotic manipulation of deformable linear objects (DLOs) is a popular area
of research which finds many applications in the manufacturing and medical
industry. Wire harness assembly is one such area which could benefit
significantly from an automated assembly process. The final assembly process
constitutes 70% of the production time while 90% of the process requires
manual work [1]. Automation could be a means of cheaper and more efficient
manufacturing in this area, especially when they are more widely used on
safety-critical systems [2] and require stricter manufacturing processes which
might not be suitable for less-experienced manual labour. An automated
manufacturing process experiences failures which are more systematic, which
are consequently easier to comprehend and anticipate. The objective of this
work is to improve on learning the dynamic manipulation of DLOs, which could
bring about future improvements in the automation of wire harness assembly.
The flinging task is chosen to demonstrate this.

The main difficulties being faced in the manipulation of DLOs are their
complex dynamics, the sensitivity of their movement to external and internal
parameters, and the underactuated nature of their system [3]. To tackle this
complex control task, forward sampling algorithms can be used. In recent
years, reinforcement learning (RL) has also gained popularity as a solution to
such problems. Both of these methods would benefit from an accurate simulator.

Dynamic manipulation of rigid objects could provide solutions to problems in
which quasi-static methods are guaranteed to fail. The same could be the case
when manipulating DLOs [4]. Other works on the dynamic manipulation of DLOs
focus on learning a task by training an agent with real data to ensure that
the training data obtained is accurate [5]. Works with training done in
simulation did not focus on the accuracy of the simulation but was successful
in bridging the sim-to-real gap by optimizing the actions of the system with a
learnt implicit policy and adaptive action sampling to achieve its goal [6].
For both these methods, data has to be collected through real experiments.

Figure 1: Simulation environment with a wire and two obstacles in MuJoCo for
the flinging task. The start (a) and a successful end (d) configuration are
shown.

The first contribution of this paper is the implementation of an existing
mathematical model of the Discrete Elastic Rods (DER) into MuJoCo*. The
processes carried out before each time step in MuJoCo, and the simulator’s
available functions, are taken into consideration when integrating the DER
model into it. The implementation will be discussed more under Section III.
This simulation will be useful for planning and learning accurate control
tasks for wires. With it, we have designed a novel dynamic wire manipulation
flinging task to demonstrate its use.

The second contribution is to show that the collection of real data for
training can be reduced with the use of an accurate wire model when carrying
out planning or learning. Real training data can be tedious to gather. With
stiffness accurately modelled in the simulation, our method is able to
successfully carry out the real task of flinging a wire111Video of real
experiment: https://youtu.be/Icb7Nl6voN0?si=wrrJlpGDIjeSxeZI&t=94, with a
76.7% success rate, using an agent that has only been trained in simulation.

##  II RELATED WORK

Simulation of DLOs have many uses apart from wire manipulation. DLOs have many
different representations in simulation. These representations have been used
in the study of polymers [7], musculoskeletal architectures [8], hair [9], and
DNA [10]. Finite element models (FEMs) are a common way to represent DLOs.
Many simulations make use of a mass-spring system. Despite their simplicity,
an interpretation of the system parameters of these models as physical
properties is lacking. To achieve the desired dynamic behavior could require
substantial effort in tuning. The mass-spring system in MuJoCo is applied
through ball joints between objects. The stiffness of the three rotational
axes of a ball joint is represented by a single parameter. Consequently, the
bending and twisting stiffness is represented by this stiffness parameter and
cannot be adjusted independently. This limits its ability to accurately model
DLOs. Position-based models [11] are widely used in computer animation for
their visual plausibility but are not known for their physical accuracy. For
the quasi-static manipulation of DLOs, their stable configurations are
required where both the bending (curvature) and twisting (torsion) energy are
at a local minimum. Minimal-energy curves have been successfully applied in
both 2D [12] and 3D [13] to manipulate DLOs. The discrete elastic rod model
[14] uses a similar concept of minimum energy to arrive at a stable
configuration of the wire. In addition to bending and twisting, the Cosserat
rod theory [15] [16] was able to include stretching and shearing to model
elastic rods which have a wider array of applications. As wire stretching and
shearing is relatively negligible for our application, the DER model will be
used in this paper to reduce computational cost. More details will be
discussed in Section III.

Manipulation of DLOs face many difficulties due to their complex dynamics and
under-actuated nature. Past works have explored learning manipulation of ropes
in a 2D space where the rope is placed flat against a table. In this way,
moving a rope to pose proved successful through numerous learning techniques
including learning from demonstrations [17] and learning from images [18].
Research has also been done to carry out more complex manipulation tasks on
DLOs such as knot tying [19] by representing the DLO in its topological rather
than geometrical state. In the 3D space, early works approached the
manipulation problem by first sampling gripper pose, then using numerical
simulation to approximate the wire configuration [13] [20]. Later, the set of
all local solutions of a Kirchhoff elastic rod was shown to be a smooth
manifold of finite dimension. As the smooth manifold can be parameterized by a
single chart, manipulation planning was able to be carried out using sampling-
based algorithms [21]. This was further brought to dual-arm manipulation of
the rod [22] which was carried out using the path-connected property of the
Kirchhoff elastic rod [23]. Another more specific manipulation task on DLOs
uses a tactile sensor to trace a wire then insert it [24].

Dynamic manipulation of DLOs proves a greater challenge due to its high
sensitivity to uncertainties. Iterative Residual Policy [6] was able
accomplish tasks with complex dynamics by learning an implicit policy using
delta dynamics and combining it with adaptive action sampling to optimize its
actions to reach a goal. This technique was able to generalize to other
scenarios better than alternative approaches. In Robots of the Lost Arc [5],
three novel dynamic DLO manipulation tasks were introduced – vaulting,
knocking, and weaving. A deep convolutional neural network along with
quadratic programming has been successfully deployed in a pipeline to learn
these tasks. When the learnt policy was deployed in real experiments, its
results showed substantial improvement over a fixed and human-specified
baseline. Their paper is what our work is largely inspired by. The simulation
designed in their paper was used as a platform to experiment with the tasks
and the learning algorithms before training is carried out with real
experiments. By including an accurately modelled wire in simulation, our work
will instead use the simulation to train a dynamic DLO manipulation policy.
The intention of this is to reduce real training data required for learning.
This will be explained further in Section IV.

##  III DYNAMIC WIRE SIMULATION

The method chosen in this paper for the simulation of a wire is using Discrete
Elastic Rods (DER) [14]. This section will briefly introduce the DER model,
along with some ideas behind it, followed by the structure of the simulation
program.

###  III-A Discrete Elastic Rod

Modelling a wire continuously in simulation is infeasible due to its infinite
dimensional property. DER theory splits a wire into discrete sections with
which the wire dynamics can be analysed. The more discrete sections there are,
the more closely the rod behaves like a real one. The DER paper differed from
other approaches of its in terms of the kinematic description, where the twist
of the material frame is represented by its angular deviation from the
parallel-transported natural Bishop frame, and the dynamic treatment, where
the material frame is treated as quasi-static while the centerline as dynamic.
The explicit representation of the centerline makes it suitable for our
simulation as it makes collision handling and rendering possible.

At each time step, the DER theory calculates the force on each node in the
discrete rod as negative of the derivative of the energy with respect to the
node position. This derivative can be obtained from the formulae of energy.
Generally, the direction the force acts in is the direction in which expected
node movement would most greatly decrease the sum of bending and twisting
energy. This value is dependent on both the nodes positions and the twist in
each discrete section. The formula is as follows. More details can be found in
the original paper.

$-\frac{dE(\Gamma)}{dx_{i}}=-\frac{\partial E(\Gamma)}{\partial
x_{i}}-\sum_{j=0}^{n}\frac{\partial
E(\Gamma)}{\partial\theta^{j}}\frac{\partial\theta^{j}}{\partial x_{i}}$
$None$

###  III-B DER in MuJoCo

After calculating the centerline forces, constraint handling and integration
of those forces will be carried out within MuJoCo [25]. Different physics
simulators operate differently and have their own set of built-in functions.
As such, they would require their own custom methods when integrating the DER
model within them. Below details specifically the framework used to integrate
the DER model into MuJoCo.

##### Modelling the Wire

Within the simulation, the wire is modelled as a continuous chain of discrete
capsules (one of the base object types that can be simulated in MuJoCo) where
the nodes are the joining points between these capsules. Each capsule is
joined with the ’connect’ equality constraint in MuJoCo, which allows for
rotational but not translational motions. These constraints will enforce the
overall wire inextensibility constraint. As the capsules are free to rotate
along the axis of their length without damping, an additional function is
implemented in the simulation which resets the rotational velocity of all
capsules at each time step. This is to ensure that the simulation does not
become unstable from an unrestricted increase in the capsules’ rotational
velocity.

##### Constraints-handling

In the original paper, an integration of the centerline using the forces
calculated is applied, followed by a manifold projection method that is used
to enforce constraint conditions, such as rod inextensibility. For this paper,
we will make use of the constraint solver within MuJoCo to enforce
constraints. The centerline forces calculated by the mathematical model is
applied at each of the corresponding joints which connect the discrete bodies.
A convex optimization problem is defined in the simulation where its global
solution is the constraint forces. The solved constraint forces are then
integrated to find the wire displacement.

###  III-C Implementation

Figure 2: Overall simulation program structure

The MuJoCo simulation environment is created with mujoco-py. The structure of
the overall program is reflected in Fig. 2. The program is first initialized
along with the MuJoCo environment. At each time step, the `update_force`
function is called which updates the program with the new twist angle,
$\theta_{n}$, and the positions of the nodes, $x_{i}$. From these updated
values, the program calculates the centerline forces on the wire. This force
is applied to the nodes of the wire with the MuJoCo function, `xfrc_applied`,
including these forces in the constraint equations of the simulator. The
resultant movement of the wire is then calculated from the global solution of
a formulated convex optimization problem solved by the MuJoCo simulator. The
time step of the simulator has to be small enough such that no instability is
caused when solving the optimization problem and to ensure that centerline
forces are updated. Delayed calculation of the centerline forces due to large
time steps might cause oscillations of the wire similar to a control system
with its control frequency set too low. A suitable time step which balances
the trade-off between the stability of the simulation and its computation time
is found to be around 0.002 seconds. Open-source code of the wire model in
MuJoCo is available on GitHub: https://github.com/qj25/der_mujoco.git.

##### Simulation Speeds

For both types of simulations, with and without DER calculations (`WD` and
`WOD`, respectively), the wire is modelled as a continuous chain of discrete
capsules, connected with the MuJoCo connect constraint, as described above.
The connect constraint does not have bending nor twisting stiffness and only
serves to keep the discrete capsules linked. `WD` includes the DER functions
that calculate and exert the centerline forces which simulate wire stiffness,
while `WOD` does not. These functions increase the computation time required.
The computation times are shown in Table I. The speed tests are a measure of
the real computation time taken (average of 300 environment runs) to simulate
one simulation time step (0.002s) with varying numbers of discrete rod
sections, with and without the DER force calculations. As can be seen, the
percentage increase in time taken decreases when more sections are simulated.
This means that the number of sections can be increased without concerns that
the time taken for stiffness calculations will increase exponentially.

Table I: Computation time to run one simulation time step (0.002s)

No. of Discrete Sections, $n$ |  Only Simulating Discrete Connected Sections in MuJoCo (s) |  With additional DER calculations (ms) |  Increase in time taken (%)
---|---|---|---
20 | 0.320 | 0.362 | 12.93
30 | 0.555 | 0.597 | 7.78
40 | 0.712 | 0.758 | 6.51
50 | 0.967 | 1.015 | 4.92
60 | 1.262 | 1.319 | 4.53

###  III-D Validation

To ensure that our simulation model of a straight isotropic wire is accurate,
some validation tests from [14] are carried out here. The referenced work
compares some testing results with available analytical solutions. Using the
localized helical buckling and Michell’s buckling instability tests,
consistency with those results will validate our simulation model.

##### Localized helical buckling

For a straight isotropic rod, localized helical buckling occurs when a twist
is introduced to one rod end and the rods ends are brought towards each other,
quasi-statically. The solution to this buckling can be described by an
analytical solution. Our model is shown to converge to this analytical
solution as the number of discrete section increases. A plot of the envelope
of the helix, $f(\varphi)$ against the dimensionless arc length, $s/s*$, is
shown in Fig. 3, where $\varphi=cos^{-1}(\mathbf{t}.\mathbf{e}_{x})$ is the
angular deviation of the tangent away from the axis passing through the rod
end points. $s$ is the rod arc length. The helix envelope can be simplified as
$f(\varphi)=\frac{cos(\varphi)-cos(\varphi_{0})}{1-cos(\varphi_{0})}$ [26],
where $s/s*=(\frac{\beta
m}{2\alpha}\sqrt{\frac{1-cos(\varphi_{0})}{1+cos(\varphi_{0})}})s$, and the
maximal angular deviation of the wire from the tangent can be calculated as
$\varphi_{0}=max_{s}\varphi(s)$. $\alpha$ and $\beta$ are the bending and
twisting stiffness, respectively. $m$ is the overall twist over the rod
length. The graph is shown to converge towards the analytical solution as more
discrete sections are used to represent the rod.

Figure 3: Localized Helical Buckling Test on a straight isotropic rod. A twist
of 27 turns is imposed on one end of the rod. The rod has properties as
follows: $L=9.29$, $\alpha=1.345$, $\beta=0.789$. The ends are brought closer
together by 0.3 units which causes helical buckling on the rod. The test is
carried out with different values of $n$ as shown in the legends. It can be
seen that the plot converges to the analytical values as $n$ increases. Table
II: Validation Results for Localized Helical Buckling Test

Number of Discrete Sections, $n$ | Average Error (unit per data point)
---|---
40 | 0.020736
60 | 0.008613
80 | 0.004083
110 | 0.001716
140 | 0.000782
180 | 0.000296

##### Michell’s Instability

Connecting the ends of a straight isotropic rod results in a loop being
formed. This loop lies in a single plane. Michell’s instability occurs when
one end is twisted along the tangent axis of the rod while the other is kept
still. As this twist angle increases buckling of the loop occurs, causing the
loop to fall out of the original plane. This phenomenon shows the coupling
between the rod twist and its equilibrium configuration. Buckling occurs at a
critical twist angle which can be determined analytically as
$\theta^{n}_{c}=2\pi\sqrt{3}/(\beta/\alpha)$ when the loop radius, $R=1$ [27]
[28]. Fig. 4 shows the graph of critical buckling angle against $\beta/\alpha$
that was obtained by testing our simulation model. It is shown that our
model’s results are similar to that of the results obtained analytically with
a maximum deviation of $1.68\%$ from the analytical plot when
$\beta/\alpha=0.5$.

Figure 4: Michell’s Instability Test on a straight isotropic rod with ends
connected to form an elastic ring. One end is twisted a certain amount until a
critical value is reached whereby the circular shape buckles into a non-planar
shape in the presence of small external disturbances. $\alpha=1$ and number of
discrete rods, $n=50$, while $\beta$ is varied to achieve different values of
$\beta/\alpha$ for which the critical twist, $\theta^{n}$ depends.

##  IV WIRE MANIPULATION

This section will focus on the dynamic manipulation of DLOs with robot arms.
More specifically, how flinging, a complex dynamic task, can be accomplished
with policy optimization in the environment created. In the environment, one
end of the wire is fixed to a point in the world space and the other is
attached to the end-effector of a 7R Panda robot arm.

###  IV-A Flinging Task

To demonstrate the usefulness of this simulation environment in robotics
planning and learning for the manipulation of DLOs, we will carry out a single
step policy optimization on the problem of flinging a 2-meter wire between two
obstacles on a table. The task can not be easily accomplished with quasi-
static manipulations due to the proximity of the wire ends and the restricted
movement of the robot arm. Because of the complex dynamics of a wire and the
under-actuated nature of manipulating one, classical methods have shown issues
of being resolution incomplete. Therefore, optimization on dynamic actions is
chosen as our solution. A render of the environment is shown in Fig. 1.
Flinging too softly will not allow the wire to go over either obstacles, while
flinging too hard will cause the wire to go over both obstacles. The movable
end of the wire is assumed to be firmly grasped by the robot arm end-effector.
The task is successful if, at the end of the robot arm motion, the wire comes
to rest between the two obstacles.

####  IV-A1 Policy Optimization

To learn the task, we make use of the stochastic policy optimization algorithm
within Proximal Policy Optimization (PPO) [29]. Our simulation environment is
designed with only a single step per episode. It provides its initial state as
an input to the algorithm. The algorithm then outputs an action vector back to
the environment to run. Depending on the results of that run, the policy is
updated. Instead of collecting a large random dataset for training, this
method updates its policy based on the results returned from taking the
stochastic actions decided by the policy itself. Essentially, this
implementation is a single step RL algorithm [30] which holds the benefit that
the resulting environment structure allows for other RL algorithms to be
easily applied when designing similar simulation experiments which require
closed-loop control.

The PPO algorithm used for policy optimization is provided by the Stable-
Baselines3 [31] package and trained on a gym [32] environment of our
simulation. The action space output, $U$, of the learning algorithm is a
custom parameterized dynamic trajectory. In [5], the dynamic manipulation
tasks were successfully carried out with a parameterized trajectory created,
using a quadratic programming solver, from an apex point and fixed start and
end points. Our paper uses multiple Cartesian sample points within the
restricted robot end-effector operational space, along with splines, to obtain
the final joint trajectory. This is done with the intention of exposing the
robot arm to a larger variety of movements while keeping the size of the
action space reasonable. These sample points are restricted to $y$ and $z$
axis translational and $y$ axis rotational pose differences from the initial
start position of the robot end-effector. A total of 3 such sample points are
taken, after which analytical inverse kinematic [33] is carried out to
determine the joint position at these sampled points. With a total of 4 sets
of joint positions (3 sampled + 1 start position), a B-spline interpolation
[34] is done on each of the 7 joints. Following that, the joint trajectory
obtained is fed into a joint position controller which commands the robot arm
to carry out the motion. The input to the algorithm is the bending and
twisting modulus of the wire, $\alpha$ and $\beta$, respectively. $d_{err}$ is
calculated as the closest distance between the position of the wire and the
center of the top of the two obstacles. If the wire is between the obstacles,
$d_{err}=0$. This error value is used to calculate the return or reward
function. The environment terminates after 75 controller command steps which
is enough for the wire to settle into a sufficiently stable state.

The reward function, $r$ is the negative of the minimum $d_{err}$ achieved
throughout the flinging. The smaller each $d_{err}$ is, the closer the reward
function is to 0. In addition, a reward of +10 is given when the task is
successfully carried out ($success_{true/false}=10$, else
$success_{true/false}=0$ if task fails).

The overall formula is as follows:

$r=max(-min(d_{err}),10*success_{true/false})$ $None$

Figure 5: Plot of mean episode reward against training time steps. Training is
done with a wire simulated using 30 discrete sections. The dashed line shows
the average rewards while the shaded region represents the mean standard
deviation during training. The solid line shows the agent which was able to
achieve the highest reward with the least time steps.

The environment is trained with 750,000 time steps total. Training is done on
a simulation environment which use 30 simulated rod sections. Results for the
optimization are shown in Fig. 5. The reward increases with training time
steps and approaches its maximum value (+10). Some of the policies are
observed to start improving later, creating a high variance for the overall
training.

##  V REAL EXPERIMENT

Real experiments222Video of real experiment:
https://youtu.be/Icb7Nl6voN0?si=wrrJlpGDIjeSxeZI&t=94 for the flinging task
are carried out to determine if the Panda robot arm is able to accomplish it
and reflect the results obtained from training in simulation. Following that,
wire positions for the real experiment will be obtained. This will be done
using a stereo vision camera to capture the experiment, specifically
identifying the wire’s features of interest and filtering out other features.
This information will be plotted and compared with the wire positions captured
in the MuJoCo environment to determine the accuracy of the model in simulating
the real flinging experiment.

Figure 6: The three wires of varying thickness and material used for the
flinging task.

For the real experiments, we will use the policy which was able to arrive at
the maximum reward with the least time steps. Three different wires (Fig. 6)
with varying thickness and material properties were used. The experiment was
repeated a total of 30 times for each wire using the policy trained. The
success rates are 76.6% (23/30 runs) for `bthin`, 10.0% (3/30 runs) for
`rthin`, and 83.3% (25/30 runs) for `bthick`. A possible reason `bthin` and
`bthick` be show more consistent results is that they are more uniformly
straight in their neutral position as compared to `rthin`. Another possible
reason is that the wires modelled in simulation have parameters which more
closely match the two black wires. All of the failures of `rthin` are caused
by the robot arm flinging the wire over both obstacles instead of just one.
The black wires experienced both under-flinging and over-flinging in their
failure cases.

###  V-A Experiment Setup

The experiment is set up as shown in Fig. 7. Placement of each object follows
the simulated environment. Only two meters of each wire will be used as per
the simulated environment. White tapes are placed at 20 cm intervals for ease
of feature detection. Experiments are reset to the same position each time as
shown in Fig. 7 (a). The stereo camera used is the Intel Realsense D435i,
which captures the left and right monochrome frames and a color frame. Each of
the cameras runs at a rate of 30 frames per second. The D435i cameras use a
global shutter which are suitable for capturing high speed objects without
distortion. The shutter speed has to be set relatively high to ensure that
there is minimal blur for ease of feature detection.

Figure 7: Flinging Experiment Setup along with an example of a successful
flinging motion carried out by the Panda Robot Arm (a to d).

###  V-B Stereo Vision

As the feature detection for the Intel Realsense software is designed for
general use, it is not able to reliably detect the relatively thin wires in
the images. Using the OpenCV Library [35], a stereo vision program was written
which detects the features in the left and right images, match them, filter
out unnecessary points, and connect the remaining points into a continuous
wire using splines. This is done on the experiment conducted with the `rthin`
wire. The results are as shown in Fig. 8.

Figure 8: Comparison of Real Experiment Results (right) using stereo vision
with the Simulation (left). The colors of the plot shows the time of the fling
at which the frame was captured, i.e. black (start) to red (end).

The stereo vision results might contain inaccuracies due to the limited depth
resolution available from the 1280x720p monochrome images. In addition,
inaccuracies might arise from the splines used to connect marker points as
they are just an approximation of the wire shape. The splines interpolate
between the identified points and could results in a loss of feature,
especially in areas where the curvature of the wire is high such as the end
position of the wire. Despite this, the general shape of the real wire
movement is similar to that obtained in simulation.

##  VI CONCLUSIONS

To accurately model a wire, the existing model of a DER has been integrated
into MuJoCo. Validation results show that the simulator performs similarly to
the model in the original DER paper [14]. Policy optimization was employed to
learn a dynamic task on the simulator. Its reliability was shown in the
results of the real experiments333Video of real experiment:
https://youtu.be/Icb7Nl6voN0?si=wrrJlpGDIjeSxeZI&t=94, where the task was
executed with a 76.7% success rate, without further training on real data. We
are able to compare the difference between the wire movement in the real and
simulated experiments. Although the general movement of the wire is observed
to be similar, more can be done to improve the accuracy of the simulation such
as accurately modelling friction of the wire and identifying the right
stiffness parameters to use.

###  VI-A Limitations and Future Work

One of the limitations is that the output trajectory from the policy is only a
guide for the robot arm to follow due to the nature of the joint trajectory
controller and how the expected trajectory is deployed. The actual trajectory
the robot takes in the simulator is slightly different. As such, the
controller of the real robot would have to be tuned precisely to ensure it
performs similarly to the one in the simulator.

Future work will be done on carrying out parameter identification of the wire.
Another option is to design the experiment in such a way that it is robust to
parameter uncertainties. The designed flinging task is an open-loop problem.
In future, the environment can be designed and employed to solve more complex
problems which require closed-loop control.

## References

  * [1] H. G. Nguyen, M. Kuhn, and J. Franke, “Manufacturing automation for automotive wiring harnesses,” _Procedia CIRP_ , vol. 97, pp. 379–384, 2021.
  * [2] M. Kuhn, H. G. Nguyen, H. Otten, and J. Franke, “Blockchain enabled traceability–securing process quality in manufacturing chains in the age of autonomous driving,” in _2018 IEEE International Conference on Technology Management, Operations and Decisions (ICTMOD)_. IEEE, 2018, pp. 131–136.
  * [3] J. Zhu, A. Cherubini, C. Dune, D. Navarro-Alarcon, F. Alambeigi, D. Berenson, F. Ficuciello, K. Harada, X. Li, J. Pan _et al._ , “Challenges and outlook in robotic manipulation of deformable objects,” _arXiv preprint arXiv:2105.01767_ , 2021.
  * [4] Q.-C. Pham, S. Caron, P. Lertkultanon, and Y. Nakamura, “Admissible velocity propagation: Beyond quasi-static path planning for high-dimensional robots,” _The International Journal of Robotics Research_ , vol. 36, no. 1, pp. 44–67, 2017.
  * [5] H. Zhang, J. Ichnowski, D. Seita, J. Wang, H. Huang, and K. Goldberg, “Robots of the lost arc: Self-supervised learning to dynamically manipulate fixed-endpoint cables,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 4560–4567.
  * [6] C. Chi, B. Burchfiel, E. Cousineau, S. Feng, and S. Song, “Iterative residual policy: for goal-conditioned dynamic manipulation of deformable objects,” _arXiv preprint arXiv:2203.00663_ , 2022.
  * [7] P. A. Wiggins, R. Phillips, and P. C. Nelson, “Exact theory of kinkable elastic polymers,” _Physical Review E_ , vol. 71, no. 2, p. 021909, 2005.
  * [8] X. Zhang, F. K. Chan, T. Parthasarathy, and M. Gazzola, “Modeling and simulation of complex dynamic musculoskeletal architectures,” _Nature communications_ , vol. 10, no. 1, pp. 1–12, 2019.
  * [9] F. Bertails, B. Audoly, M.-P. Cani, B. Querleux, F. Leroy, and J.-L. Lévêque, “Super-helices for predicting the dynamics of natural hair,” _ACM Transactions on Graphics (TOG)_ , vol. 25, no. 3, pp. 1180–1187, 2006.
  * [10] Y. V. Sereda and T. C. Bishop, “Evaluation of elastic rod models with long range interactions for predicting nucleosome stability,” _Journal of Biomolecular Structure and Dynamics_ , vol. 27, no. 6, pp. 867–887, 2010.
  * [11] N. Umetani, R. Schmidt, and J. Stam, “Position-based elastic rods,” in _ACM SIGGRAPH 2014 Talks_ , 2014, pp. 1–1.
  * [12] E. Jou and W. Han, “Minimal-energy splines with various end constraints,” in _Curve and surface design_. SIAM, 1992, pp. 23–40.
  * [13] M. Moll and L. E. Kavraki, “Path planning for deformable linear objects,” _IEEE Transactions on Robotics_ , vol. 22, no. 4, pp. 625–636, 2006.
  * [14] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, and E. Grinspun, “Discrete elastic rods,” in _ACM SIGGRAPH 2008 papers_ , 2008, pp. 1–12.
  * [15] E. Cosserat, _Théorie des corps déformables_. Librairie Scientifique A. Hermann et Fils, 1909.
  * [16] J. Spillmann and M. Teschner, “Corde: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects,” in _Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation_ , 2007, pp. 63–72.
  * [17] A. Nair, D. Chen, P. Agrawal, P. Isola, P. Abbeel, J. Malik, and S. Levine, “Combining self-supervised learning and imitation for vision-based rope manipulation,” in _2017 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2017, pp. 2146–2153.
  * [18] A. Wang, T. Kurutach, K. Liu, P. Abbeel, and A. Tamar, “Learning robotic manipulation through visual planning and acting,” _arXiv preprint arXiv:1905.04411_ , 2019.
  * [19] M. Saha and P. Isto, “Manipulation planning for deformable linear objects,” _IEEE Transactions on Robotics_ , vol. 23, no. 6, pp. 1141–1150, 2007.
  * [20] F. Lamiraux and L. E. Kavraki, “Planning paths for elastic objects under manipulation constraints,” _The International Journal of Robotics Research_ , vol. 20, no. 3, pp. 188–208, 2001.
  * [21] T. Bretl and Z. McCarthy, “Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations,” _The International Journal of Robotics Research_ , vol. 33, no. 1, pp. 48–68, 2014.
  * [22] A. Sintov, S. Macenski, A. Borum, and T. Bretl, “Motion planning for dual-arm manipulation of elastic rods,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6065–6072, 2020.
  * [23] A. Borum and T. Bretl, “The free configuration space of a kirchhoff elastic rod is path-connected,” in _2015 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2015, pp. 2958–2964.
  * [24] Y. She, S. Wang, S. Dong, N. Sunil, A. Rodriguez, and E. Adelson, “Cable manipulation with a tactile-reactive gripper,” _The International Journal of Robotics Research_ , vol. 40, no. 12-14, pp. 1385–1401, 2021.
  * [25] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ international conference on intelligent robots and systems_. IEEE, 2012, pp. 5026–5033.
  * [26] G. Van der Heijden and J. Thompson, “Helical and localised buckling in twisted rods: a unified analysis of the symmetric case,” _Nonlinear dynamics_ , vol. 21, no. 1, pp. 71–99, 2000.
  * [27] J. Michell, “On the stability of a bent and twisted wire,” _Messenger Math_ , vol. 11, pp. 181–184, 1889.
  * [28] A. Goriely, “Twisted elastic rings and the rediscoveries of michell’s instability,” _Journal of Elasticity_ , vol. 84, no. 3, pp. 281–299, 2006.
  * [29] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017.
  * [30] H. Ghraieb, J. Viquerat, A. Larcher, P. Meliga, and E. Hachem, “Single-step deep reinforcement learning for open-loop control of laminar and turbulent flows,” _Physical Review Fluids_ , vol. 6, no. 5, p. 053902, 2021.
  * [31] A. Raffin, A. Hill, M. Ernestus, A. Gleave, A. Kanervisto, and N. Dormann, “Stable baselines3,” 2019.
  * [32] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, “Openai gym,” _arXiv preprint arXiv:1606.01540_ , 2016.
  * [33] Y. He and S. Liu, “Analytical inverse kinematics for franka emika panda–a geometrical solver for 7-dof manipulators with unconventional design,” in _2021 9th International Conference on Control, Mechatronics and Automation (ICCMA)_. IEEE, 2021, pp. 194–199.
  * [34] C. De Boor and C. De Boor, _A practical guide to splines_. springer-verlag New York, 1978, vol. 27.
  * [35] G. Bradski and A. Kaehler, _Learning OpenCV: Computer vision with the OpenCV library_. ” O’Reilly Media, Inc.”, 2008.


#  How the perturbed vacuum might be completely decoupled from all physical
states, allowing for a well-defined Hamiltonian formulation of QED

Mads J. Damgaard111 B.Sc. at the Niels Bohr Institute, University of
Copenhagen. E-mail<EMAIL_ADDRESS>
###### Abstract

We carry out a Dirac sea reinterpretation of a discretized version of the
Hamiltonian of quantum electrodynamics (QED), and analyze the perturbed vacuum
in the continuum limit. We argue that if certain operators can be shown to be
the self-adjoint, the perturbed vacuum will have solutions that converge
nicely in this limit to states of the infinitesimal subspace of the Hilbert
space spanned by momentum eigenstates whose momenta sum to 0. This convergence
allows us to analyze what effect the perturbed vacuum has on the physical
states of the system. We then show that the interaction and the interference
between the perturbed vacuum and the physical states actually vanish in the
continuum limit, due to the fact that the states of the perturbed vacuum only
live in this vanishingly small part of the momentum space. We furthermore show
that this allows us to remove the terms that perturb the vacuum from the
Hamiltonian, without changing the dynamics of the physical states, thus
yielding us a well-defined Hamiltonian for the full theory of QED.

##  1 Motivation

It is conventional to use the path integral formulation of quantum mechanics
when working with quantum electrodynamics (QED), and other theories of quantum
field theory (QFT), rather than using the Hamiltonian formulation. Not much
attention is therefore typically given to the Hamiltonian of such theories.
However, since QFT assumes that the time evolution is unitary, a generator of
that time evolution must exist. Any mathematically consistent QFT must
therefore have an underlying Hamiltonian.

Before the negative-energy fermion solutions of QED are reinterpreted as
antiparticles, in accordance with the principles of the Dirac sea, the
Hamiltonian of the theory _can_ indeed be derived.222 See e.g. Weinberg [4],
or see Damgaard [1].  But after this Dirac sea reinterpretation, we seem to
get an infinite amount of vacuum fluctuations, as particles are created and
annihilated again from the bare vacuum. This makes it, at least at first
glance, seem like an almost impossible task to describe the theory using the
Hamiltonian formulation.

It is therefore quite natural to want to abandon this approach, especially
when the alternative one of the path integral formulation is available to us.
However, finding the full Hamiltonian of QED might lead to some new insights.
And in particular, it might give us more insight into how the perturbed vacuum
affects the dynamics of the _physical states_ of the theory, i.e. the states
which we might prepare and do experiments on in a laboratory.

In fact, the analysis of this paper into the Hamiltonian formulation of QED
suggests that the perturbed vacuum actually has no effect at all on the
physical states, and that the interaction terms that cause these vacuum
fluctuations can simply be removed for the theory, without changing its
dynamics.

##  2 The initial Hamiltonian

Before the Dirac sea reinterpretation, where the negative-energy fermion
eigenstates are reinterpreted as antiparticles, the Hamiltonian of QED is
formally given by2

$\displaystyle\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{D}}+\hat{H}_{\mathrm{C}}$
(1)

on a Fock space, $\mathbf{H}$, where $\hat{H}_{0}$ is the free energy
operator, and where $\hat{H}_{\mathrm{D}}$ and $\hat{H}_{\mathrm{C}}$ are the
operators of the Dirac interaction and the Coulomb interaction, respectively.
Note that while it is not typical to see the Coulomb interaction appear
explicitly in the path integrals of QFT, this is not because it is not there:
It is simply absorbed into the Feynman propagator of the photons (see e.g.
Weinberg [4]).

The exact definitions of $\hat{H}_{\mathrm{D}}$ and $\hat{H}_{\mathrm{C}}$ are
not important for the argument presented in this paper, but for the interested
reader, they are elaborated upon in Appendix A. The important thing for our
purposes is just to note that $\hat{H}_{\mathrm{D}}$ and
$\hat{H}_{\mathrm{C}}$ have the general forms given by

$\displaystyle\hat{H}_{\mathrm{D}}$
$\displaystyle=\sum_{\lambda=1}^{2}\sum_{s,s^{\prime}=1}^{4}\int\frac{d\mathbf{k}\,d\mathbf{p}}{(2\pi)^{6}}\,\frac{\alpha_{\lambda,s,s^{\prime}}(\mathbf{k},\mathbf{p})}{\sqrt{|\mathbf{k}|}}\big{(}\hat{a}_{\lambda}(\mathbf{k})+\hat{a}_{\lambda}^{\dagger}(\mathbf{-k})\big{)}\hat{b}_{s^{\prime}}^{\dagger}(\mathbf{p}+\mathbf{k})\hat{b}_{s}(\mathbf{p}),$
(2) $\displaystyle\hat{H}_{\mathrm{C}}$
$\displaystyle=\sum_{\lambda=1}^{2}\sum_{s,s^{\prime},t,t^{\prime}=1}^{4}\int\frac{d\mathbf{k}\,d\mathbf{p}_{1}\,d\mathbf{p}_{2}}{(2\pi)^{9}}\,\frac{\beta_{\lambda,s,s^{\prime},t,t^{\prime}}(\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2})}{\mathbf{k}^{2}}\hat{b}_{s^{\prime}}^{\dagger}(\mathbf{p}_{1}+\mathbf{k})\hat{b}_{t^{\prime}}^{\dagger}(\mathbf{p}_{2}-\mathbf{k})\hat{b}_{s}(\mathbf{p}_{1})\hat{b}_{t}(\mathbf{p}_{2}),$
(3)

with $\alpha_{\lambda,s,s^{\prime}}$ and
$\beta_{\lambda,s,s^{\prime},t,t^{\prime}}$ being some bounded complex
functions. Here $\hat{a}_{\lambda}^{\dagger}(\mathbf{k})$ and
$\hat{a}_{\lambda}(\mathbf{k})$ are of course the creation and annihilation
operators for the photons, and $\hat{b}_{s}^{\dagger}(\mathbf{p})$ and
$\hat{b}_{s}(\mathbf{p})$ are the creation and annihilation operators for some
fermions with 4 spin components.333 It is conventional to use the symbol $d$
instead of $b$ for the two spin indices of the negative-energy eigenstates,
$s\in\\{3,4\\}$, as is also seen in Appendix A. But here we use $b$ for all
spin indices in order to write these expressions more compactly.  And for the
free energy operator, we have

$\hat{H}_{0}=\sum_{\lambda=1}^{2}\int\frac{d\mathbf{k}}{(2\pi)^{3}}\,|\mathbf{k}|\hat{a}_{\lambda}^{\dagger}(\mathbf{k})\hat{a}_{\lambda}(\mathbf{k})+\int\frac{d\mathbf{p}}{(2\pi)^{3}}\,E_{\mathbf{p}}\Big{(}\sum_{s=1}^{2}\hat{b}_{s}^{\dagger}(\mathbf{p})\hat{b}_{s}(\mathbf{p})-\sum_{s=3}^{4}\hat{b}_{s}^{\dagger}(\mathbf{p})\hat{b}_{s}(\mathbf{p})\Big{)},$
(4)

where $E_{\mathbf{p}}=\sqrt{\mathbf{p}^{2}+m_{\mathrm{F}}^{2}}$, with
$m_{\mathrm{F}}$ being the mass of the fermions.

We will assume in this paper that $\hat{H}$ is self-adjoint on some dense
domain of $\mathbf{H}$, and we will furthermore assume that we are free to
make an ultraviolet cutoff, $|\mathbf{p}|,|\mathbf{k}|<\Lambda$, for
$\hat{H}$, as well as for all other Hamiltonians appearing in this paper, and
then expect the dynamics of these cut-off operators to converge to those of
the original ones when $\Lambda\to\infty$.

##  3 Discretizing the operators

Before we can make the Dirac sea reinterpretation, where we want to turn
$\hat{b}_{3}^{\dagger}(\mathbf{p}),\hat{b}_{4}^{\dagger}(\mathbf{p})$ into
$\hat{b}_{3}(-\mathbf{p}),\hat{b}_{4}(-\mathbf{p})$ and vice versa, we will
first need to discretize all the momenta. Let us therefore introduce the
discretized creation and annihilation operators,
$\hat{a}_{\lambda\mathbf{k}}^{\dagger}$, $\hat{a}_{\lambda\mathbf{k}}$,
$\hat{b}_{s\mathbf{p}}^{\dagger}$, $\hat{b}_{s\mathbf{p}}$,444 Such
discretized operators can be constructed using step functions (in momentum
space), which thus allows us to remain in the same Hilbert space.  where
$\mathbf{k}$ and $\mathbf{p}$ now range over a discrete and bounded set,
$\mathbb{K}$, which we can define by

$\mathbb{K}=\big{\\{}\mathbf{k}\in\mathbb{R}^{3}\,\big{|}\,\mathbf{k}=\delta
k\mathbf{z},\;\mathbf{z}\in\mathbb{Z}^{3}\setminus\\{(0,0,0)\\},\;|\mathbf{k}|<\Lambda\big{\\}}$
(5)

for some lattice spacing $\delta k$, and for the ultraviolet cutoff of
$\Lambda$.

With this discretization, we then need to replace the creation and
annihilation operators of Eqs. (2–4) such that all transition probabilities
converge when we let $\delta k\to 0$. Damgaard [1] (Section 8) shows a
heuristic derivation of the following rules to ensure this:

$\displaystyle\begin{gathered}\int\frac{d\mathbf{k}}{(2\pi)^{3/2}}\,\hat{a}_{\lambda}^{\dagger}(\mathbf{k})\to\sum_{\mathbf{k}\in\mathbb{K}}\delta
k^{3/2}\hat{a}_{\lambda\mathbf{k}}^{\dagger},\quad\int\frac{d\mathbf{k}}{(2\pi)^{3/2}}\,\hat{a}_{\lambda}(\mathbf{k})\to\sum_{\mathbf{k}\in\mathbb{K}}\delta
k^{3/2}\hat{a}_{\lambda\mathbf{k}},\\\
\int\frac{d\mathbf{k}}{(2\pi)^{3}}\,\hat{a}_{\lambda}^{\dagger}(\mathbf{k})\hat{a}_{\lambda}(\mathbf{k})\to\sum_{\mathbf{k}\in\mathbb{K}}\hat{a}_{\lambda\mathbf{k}}^{\dagger}\hat{a}_{\lambda\mathbf{k}}.\end{gathered}$
(8)

And these results can be generalized by observing that for any transitions
from an initial state,
$\ket{\mathbf{p}_{1},\ldots,\mathbf{p}_{n};\mathbf{k}_{1},\ldots,\mathbf{k}_{m}}$
to a final state,
$\ket{\mathbf{p}^{\prime}_{1},\ldots,\mathbf{p}^{\prime}_{n^{\prime}};\mathbf{k}^{\prime}_{1},\ldots,\mathbf{k}^{\prime}_{m^{\prime}}}$,
we need to have

$\braket{\mathbf{p}^{\prime}_{1},\ldots,\mathbf{p}^{\prime}_{n^{\prime}};\mathbf{k}^{\prime}_{1},\ldots,\mathbf{k}^{\prime}_{m^{\prime}}}{\hat{A}}{\mathbf{p}_{1},\ldots,\mathbf{p}_{n};\mathbf{k}_{1},\ldots,\mathbf{k}_{m}}\propto\delta
k^{\Delta/2+d}$ (9)

for the relevant matrix element of a given operator $\hat{A}$, where $\Delta$,
first of all, denotes the difference between the degrees of freedom of the
final and the initial state. This will typically be equal to
$n^{\prime}-n+m^{\prime}-m$, unless the initial or finial states have some
Dirac delta function(s) in them. And $d$, second of all, denotes the degrees
of freedom for the initial state in terms of how
$(\mathbf{p}_{1},\ldots,\mathbf{p}_{n},\mathbf{k}_{1},\ldots,\mathbf{k}_{m})$
is allowed to vary given that
$(\mathbf{p}^{\prime}_{1},\ldots,\mathbf{p}^{\prime}_{n^{\prime}},\mathbf{k}^{\prime}_{1},\ldots,\mathbf{k}^{\prime}_{m^{\prime}})$
is fixed for the final state.

The reasoning behind this general rule is that if the states in the
neighborhood of the initial point
$(\mathbf{p}_{1},\ldots,\mathbf{p}_{n},\mathbf{k}_{1},\ldots,\mathbf{k}_{m})$
have amplitudes that are proportional to $\delta k^{3a/2}$, then the states in
the neighborhood of the final point
$(\mathbf{p}^{\prime}_{1},\ldots,\mathbf{p}^{\prime}_{n^{\prime}},\mathbf{k}^{\prime}_{1},\ldots,\mathbf{k}^{\prime}_{m^{\prime}})$
need to get amplitudes that are proportional to $\delta k^{3(\Delta+a)/2}$ in
order for probabilities to be preserved in the continuum limit. So the initial
guess might be that the matrix elements should be proportional to simply
$\delta k^{\Delta/2}$. However, when there are $d$ degrees of freedom for the
initial states going to the same, fixed final state, the sum over these
degrees of freedom will yield a value that grows proportional to $\delta
k^{-d}$. This is why we also need to add $d$ to the exponent on the right-hand
side of Eq. (9) in order to get the correct normalization.

So for the discretized versions of these operators, we can see that we need to
have555 We use proportional-to signs here instead of keeping track of the
exact factors of $(2\pi)^{a}$.

$\displaystyle\hat{H}_{0}^{\prime}$
$\displaystyle=\sum_{\lambda=1}^{2}\sum_{\mathbf{k}\in\mathbb{K}}|\mathbf{k}|\hat{a}_{\lambda\mathbf{k}}^{\dagger}\hat{a}_{\lambda\mathbf{k}}+\sum_{\mathbf{p}\in\mathbb{K}}E_{\mathbf{p}}\Big{(}\sum_{s=1}^{2}\hat{b}_{s\mathbf{p}}^{\dagger}\hat{b}_{s\mathbf{p}}-\sum_{s=3}^{4}\hat{b}_{s\mathbf{p}}^{\dagger}\hat{b}_{s\mathbf{p}}\Big{)},$
(10) $\displaystyle\hat{H}_{\mathrm{D}}^{\prime}$
$\displaystyle\propto\sum_{\lambda=1}^{2}\sum_{s,s^{\prime}=1}^{4}\sum_{\mathbf{k},\mathbf{p}\in\mathbb{K}}\delta
k^{3/2}\frac{\alpha_{\lambda,s,s^{\prime}}(\mathbf{k},\mathbf{p})}{\sqrt{|\mathbf{k}|}}(\hat{a}_{\lambda\mathbf{k}}+\hat{a}_{\lambda(-\mathbf{k})}^{\dagger})\hat{b}_{s^{\prime}(\mathbf{p}+\mathbf{k})}^{\dagger}\hat{b}_{s\mathbf{p}},$
(11) $\displaystyle\hat{H}_{\mathrm{C}}^{\prime}$
$\displaystyle\propto\sum_{\lambda=1}^{2}\sum_{s,s^{\prime},t,t^{\prime}=1}^{4}\sum_{\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2}\in\mathbb{K}}\delta
k^{3}\frac{\beta_{\lambda,s,s^{\prime},t,t^{\prime}}(\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2})}{\mathbf{k}^{2}}\hat{b}_{s^{\prime}(\mathbf{p}_{1}+\mathbf{k})}^{\dagger}\hat{b}_{t^{\prime}(\mathbf{p}_{2}-\mathbf{k})}^{\dagger}\hat{b}_{s\mathbf{p}_{1}}\hat{b}_{t\mathbf{p}_{2}}.$
(12)

The factor of $\delta k^{3/2}$ here in $\hat{H}_{\mathrm{D}}^{\prime}$ thus
comes from having respectively $\Delta/2+d=-3/2+3$ and $\Delta/2+d=3/2+0$ for
the two terms. And the factor of $\delta k^{3}$ in
$\hat{H}_{\mathrm{C}}^{\prime}$ comes from having $\Delta/2+d=0+3$.

##  4 The Dirac sea reinterpretation

We will now try to turn
$\hat{b}_{3\mathbf{p}}^{\dagger},\hat{b}_{4\mathbf{p}}^{\dagger}$ into
$\hat{b}_{3(-\mathbf{p})},\hat{b}_{4(-\mathbf{p})}$ and vice versa in
accordance with the principles of the Dirac sea reinterpretation.

If we start by considering $\hat{H}_{0}^{\prime}$, we see that we thus need to
turn

$\sum_{s=3}^{4}\hat{b}_{s(-\mathbf{p})}^{\dagger}\hat{b}_{s(-\mathbf{p})}\to\sum_{s=3}^{4}\hat{b}_{s\mathbf{p}}\hat{b}_{s\mathbf{p}}^{\dagger}=\sum_{s=3}^{4}(1-\hat{b}_{s\mathbf{p}}^{\dagger}\hat{b}_{s\mathbf{p}}).$
(13)

Due to our ultraviolet cutoff, $|\mathbf{p}|<\Lambda$, the first term of $1$
within this sum on the right-hand side will only yield a constant energy, $C$.
And since we are free to remove this constant energy at every step when we let
$\delta k\to 0$, $\Lambda\to\infty$, without changing the dynamics of the
system, we can thus write

$\hat{H}_{0}^{\prime\prime}=\sum_{\lambda=1}^{2}\sum_{\mathbf{k}\in\mathbb{K}}|\mathbf{k}|\hat{a}_{\lambda\mathbf{k}}^{\dagger}\hat{a}_{\lambda\mathbf{k}}+\sum_{s=1}^{4}\sum_{\mathbf{p}\in\mathbb{K}}E_{\mathbf{p}}\hat{b}_{s\mathbf{p}}^{\dagger}\hat{b}_{s\mathbf{p}}$
(14)

for this discretized and reinterpreted $\hat{H}_{0}$.

If we then consider $\hat{H}_{\mathrm{D}}^{\prime}$ and
$\hat{H}_{\mathrm{C}}^{\prime}$, we see that the creation and annihilation
operators anti-commute almost everywhere, and we therefore do not need to add
or subtract any constant energies here. Thus, for these operators, we can
write similar formulas for the reinterpreted cases, namely

$\displaystyle\hat{H}_{\mathrm{D}}^{\prime\prime}$
$\displaystyle\propto\sum_{\lambda=1}^{2}\sum_{s,s^{\prime}=1}^{4}\sum_{\mathbf{k},\mathbf{p}\in\mathbb{K}}\delta
k^{3/2}\frac{\alpha_{\lambda,s,s^{\prime}}(\mathbf{k},\mathbf{p})}{\sqrt{|\mathbf{k}|}}(\hat{a}_{\lambda\mathbf{k}}+\hat{a}_{\lambda(-\mathbf{k})}^{\dagger})\hat{c}_{s^{\prime}(\mathbf{p}+\mathbf{k})}^{\dagger}\hat{c}_{s\mathbf{p}},$
(15) $\displaystyle\hat{H}_{\mathrm{C}}^{\prime\prime}$
$\displaystyle\propto\sum_{\lambda=1}^{2}\sum_{s,s^{\prime},t,t^{\prime}=1}^{4}\sum_{\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2}\in\mathbb{K}}\delta
k^{3}\frac{\beta_{\lambda,s,s^{\prime},t,t^{\prime}}(\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2})}{\mathbf{k}^{2}}\hat{c}_{s^{\prime}(\mathbf{p}_{1}+\mathbf{k})}^{\dagger}\hat{c}_{t^{\prime}(\mathbf{p}_{2}-\mathbf{k})}^{\dagger}\hat{c}_{s\mathbf{p}_{1}}\hat{c}_{t\mathbf{p}_{2}},$
(16)

but where $\hat{c}_{s\mathbf{p}}=\hat{b}_{s\mathbf{p}}$ when $s\in\\{1,2\\}$
and $\hat{c}_{s\mathbf{p}}=\hat{b}_{s(-\mathbf{p})}^{\dagger}$ when
$s\in\\{3,4\\}$.

##  5 Going back to the continuum limit

When we let $\delta k\to 0$, it would be nice if we could use the generalized
rule of Eq. (8–9) in reverse and get a nice operator back with no divergences.
However, this is not quite the case, as we will see in this section, since the
terms that consist either exclusively of creation operators or of annihilation
operators seem to cause us some trouble. These are in other words the terms
that include either

$\displaystyle\hat{a}_{\lambda\mathbf{k}}\hat{b}_{s^{\prime}(-\mathbf{p}-\mathbf{k})}\hat{b}_{s\mathbf{p}},\quad\hat{b}_{s^{\prime}(-\mathbf{p}_{1}-\mathbf{k})}\hat{b}_{t^{\prime}(-\mathbf{p}_{2}+\mathbf{k})}\hat{b}_{s\mathbf{p}_{1}}\hat{b}_{t\mathbf{p}_{2}},$
(17)

or their adjoint counterparts. Since these are the terms that perturb the
vacuum, let us refer to such terms as the _vacuum-perturbing_ terms in the
remainder of this paper.

Before we analyze these troubling vacuum-perturbing terms further, let us
discuss the other, less troubling terms first. We will refer to these as the
_physical_ terms of the operator for reasons that will be clear by the end of
this analysis.

If we disregard the vacuum-perturbing terms of
$\hat{H}_{\mathrm{D}}^{\prime\prime}$, we see that the rest of the terms
either have two incoming particles and one outgoing particle, or have one
incoming particle and two outgoing ones. We see that in all these cases, we
have $\Delta/2+d=3/2$, which fits the factor of $\delta k^{3/2}$. So the sum
of all these terms has a well-defined formal continuum limit, given by

$\displaystyle\hat{H}_{\mathrm{D,phys}}^{\prime\prime\prime}$
$\displaystyle=\sum_{\lambda=1}^{2}\sum_{s,s^{\prime}=1}^{4}\int\frac{d\mathbf{k}\,d\mathbf{p}}{(2\pi)^{6}}\,\frac{\alpha_{\lambda,s,s^{\prime}}(\mathbf{k},\mathbf{p})}{\sqrt{|\mathbf{k}|}}\big{(}\hat{a}_{\lambda}(\mathbf{k})+\hat{a}_{\lambda}^{\dagger}(\mathbf{-k})\big{)}\hat{c}_{s^{\prime}}^{\dagger}(\mathbf{p}+\mathbf{k})\hat{c}_{s}(\mathbf{p})$
(18) $\displaystyle\phantom{=}\quad\quad\quad-\,\text{vacuum-perturbing
terms}.$

(Here we take these continuous $\hat{c}_{s}$-operators to have the same
relation to the $\hat{b}_{s}$-operators as in the discrete case above.)

If we then look at $\hat{H}_{\mathrm{C}}^{\prime\prime}$, we see that, apart
from the vacuum-perturbing terms, we still have some terms with both two
incoming and two outgoing particles, for which $\Delta/2+d$ is thus still 3.
Additionally, there are some terms with either three incoming and one outgoing
particle, or with one incoming particle and three outgoing ones. Here we can
see that we respectively have $\Delta/2+d=-6/2+6$ and $\Delta/2+d=6/2+0$,
which also gives us $\Delta/2+d=3$ in both cases. This still fits the factor
of $\delta k^{3}$, and the sum all of these terms therefore has a well-defined
formal continuum limit as well, given by

$\displaystyle\hat{H}_{\mathrm{C,phys}}^{\prime\prime\prime}$
$\displaystyle=\sum_{\lambda=1}^{2}\sum_{s,s^{\prime},t,t^{\prime}=1}^{4}\int\frac{d\mathbf{k}\,d\mathbf{p}_{1}\,d\mathbf{p}_{2}}{(2\pi)^{9}}\,\frac{\beta_{\lambda,s,s^{\prime},t,t^{\prime}}(\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2})}{\mathbf{k}^{2}}\hat{c}_{s^{\prime}}^{\dagger}(\mathbf{p}_{1}+\mathbf{k})\hat{c}_{t^{\prime}}^{\dagger}(\mathbf{p}_{2}-\mathbf{k})\hat{c}_{s}(\mathbf{p}_{1})\hat{c}_{t}(\mathbf{p}_{2})$
(19) $\displaystyle\phantom{=}\quad\quad\quad-\,\text{vacuum-perturbing
terms}.$

Let us then look at the remaining, vacuum-perturbing terms, and see why we
unfortunately cannot simply use the general rule of Eq. (8–9) in reverse for
these. For the Dirac interaction, we here have either three incoming or
outgoing particles with nothing on the other side. In particular, we can see
that for the latter of these cases, due to momentum conservation, we have
$\Delta/2+d=(9-3)/2+0=3$. And this does not fit the factor of $\delta
k^{3/2}$. Similarly, for the Coulomb interaction, we can see that we have
$\Delta/2+d=(12-3)/2+0=9/2$ for the case when there are four outgoing
particles and no incoming ones, which also does not fit the factor of $\delta
k^{3}$. Thus, we cannot expect the vacuum-perturbing part of
$\hat{H}^{\prime\prime}$ to have converging dynamics in the continuum limit,
as its eigenvalues might very well diverge.

Before we abandon all hope of completing the Dirac sea reinterpretation,
however, let us note that such diverging energies for the perturbed vacuum
will only be a persisting problem if the states of the perturbed vacuum do not
decouple from what we here refer to as the “physical states” in the continuum
limit. If they do, on the other hand, we might be able to simply subtract this
vacuum energy at every step going towards the continuum limit, leading us to
converging dynamics for the physical states.

Motivated by this fact, let us consider the vacuum-perturbing terms
separately. And since these terms can only cause transitions from the bare
vacuum into other states with a combined momentum of 0, let us consider these
operators on the subspace of $\mathbf{H}$, call it
$\mathbf{H}_{\mathrm{vac}}$, spanned by momentum eigenstates where all the
momenta sum to 0.

On this subspace, we still have $\Delta/2+d=(9-3)/2+0=3$ for the terms where
there are three outgoing particles and no incoming particle, and
$\Delta/2+d=(12-3)/2+0=9/2$ for the terms where there are four outgoing
particles. And for the adjoint of these terms, we then get
$\Delta/2+d=(-9+3)/2+6=3$ and $\Delta/2+d=(-12+3)/2+9=9/2$, respectively. We
thus see that all these terms are exactly a factor of $\delta k^{3/2}$ shy of
fulfilling the rule of Eq. (9). So while
$\hat{H}_{\mathrm{D,vac}}^{\prime\prime}$ and
$\hat{H}_{\mathrm{C,vac}}^{\prime\prime}$ will not have converging formal
continuum limits on $\mathbf{H}_{\mathrm{vac}}$, we see that $\delta
k^{3/2}\hat{H}_{\mathrm{D,vac}}^{\prime\prime}$ and $\delta
k^{3/2}\hat{H}_{\mathrm{C,vac}}^{\prime\prime}$ _will_. These formal continuum
limits of $\delta k^{3/2}\hat{H}_{\mathrm{D,vac}}^{\prime\prime}$ and $\delta
k^{3/2}\hat{H}_{\mathrm{C,vac}}^{\prime\prime}$ on
$\mathbf{H}_{\mathrm{vac}}$, call them
$\hat{H}_{\mathrm{D,vac}}^{\prime\prime\prime}$ and
$\hat{H}_{\mathrm{C,vac}}^{\prime\prime\prime}$, will thus be given by

$\displaystyle\begin{aligned}
\hat{H}_{\mathrm{D,vac}}^{\prime\prime\prime}=\sum_{\lambda=1}^{2}\sum_{s,s^{\prime}=1}^{4}\int\frac{d\mathbf{k}\,d\mathbf{p}}{(2\pi)^{6}}\,\frac{\alpha^{\prime}_{\lambda,s,s^{\prime}}(\mathbf{k},\mathbf{p})}{\sqrt{|\mathbf{k}|}}\big{(}\hat{a}_{\lambda}(\mathbf{k})\hat{b}_{s^{\prime}}(-\mathbf{p}-\mathbf{k})\hat{b}_{s}(\mathbf{p})+\,\text{adjoint
term}\big{)},\end{aligned}$ (20) $\displaystyle\begin{aligned}
\hat{H}_{\mathrm{C,vac}}^{\prime\prime\prime}&=\sum_{\lambda=1}^{2}\sum_{s,s^{\prime},t,t^{\prime}=1}^{4}\int\frac{d\mathbf{k}\,d\mathbf{p}_{1}\,d\mathbf{p}_{2}}{(2\pi)^{9}}\,\frac{\beta^{\prime}_{\lambda,s,s^{\prime},t,t^{\prime}}(\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2})}{\mathbf{k}^{2}}\\\
&\phantom{=}\quad\quad\quad\big{(}\hat{b}_{s^{\prime}}(-\mathbf{p}_{1}-\mathbf{k})\hat{b}_{t^{\prime}}(-\mathbf{p}_{2}+\mathbf{k})\hat{b}_{s}(\mathbf{p}_{1})\hat{b}_{t}(\mathbf{p}_{2})+\,\text{adjoint
term}\big{)},\end{aligned}$ (21)

where all these creation and annihilation operators are then redefined as
operators on $\mathbf{H}_{\mathrm{vac}}$.666 One might note that the wave
functions of this reduced Hilbert space should not be normalized according to
$\int|\psi(\mathbf{p}_{1},\ldots,\mathbf{p}_{n};\mathbf{k}_{1},\ldots,\mathbf{k}_{m})|^{2}\,d\mathbf{p}_{1}\cdots
d\mathbf{p}_{n}\,d\mathbf{k}_{1}\cdots d\mathbf{k}_{m}=P$, where $P$ is some
probability, but should have the following Dirac delta function in the
integrand as well:
$\delta^{3(n+m)}(\mathbf{p}_{1}+\ldots+\mathbf{p}_{n}+\mathbf{k}_{1}+\ldots+\mathbf{k}_{m})$.
Note that $\alpha^{\prime}_{\lambda,s,s^{\prime}}(\mathbf{k},\mathbf{p})$ and
$\beta^{\prime}_{\lambda,s,s^{\prime},t,t^{\prime}}(\mathbf{k},\mathbf{p}_{1},\mathbf{p}_{2})$
here are not the same coefficients as before. (In particular, they are both 0
whenever the relevant terms would otherwise break the law that the difference
between the number of particles and antiparticles is conserved.)

Now, if the sum of these formal operators can be shown to be self-adjoint on
some domain in $\mathbf{H}_{\mathrm{vac}}$, the spectral theorem777See e.g.
Hall [3]. will tell us that
$\hat{H}_{\mathrm{D,vac}}^{\prime\prime\prime}+\hat{H}_{\mathrm{D,vac}}^{\prime\prime\prime}$
has $\varepsilon$-almost eigenstates888 In accordance with Hall [3]
(Definition 10.24), we take an _$\varepsilon$ -almost eigenstate_ of an
operator, $\hat{A}$, with $\varepsilon>0$, to refer to a state, $\psi$, for
which $\|\hat{A}\psi-\lambda\psi\|<\varepsilon\|\psi\|$ for some _eigenvalue_
$\lambda\in\mathbb{C}$.  with real-valued eigenvalues. This would then suggest
that, despite the divergences of
$\hat{H}_{\mathrm{D,vac}}^{\prime\prime}+\hat{H}_{\mathrm{C,vac}}^{\prime\prime}$,
the eigenstates of that operator, i.e. on the subspace of states whose momenta
sum to 0, will still converge in the continuum limit, and that it is only the
eigen _values_ of these states which will diverge.

Of course, having diverging energies in a theory is still a problem, and
therefore we still need for the perturbed vacuum to decouple from the physical
states in the continuum limit in order for this Dirac sea reinterpretation to
work. However, when we consider the perturbed vacuum, we see that it only
lives on a very small part of the momentum space, whose volume vanishes when
$\delta k\to 0$. Therefore it might indeed be the case that the interaction of
the perturbed vacuum with the physical states of the system vanishes in the
continuum limit. If this is true, and if the interference vanishes as well, we
get that the perturbed vacuum _does_ decouple from the physical states in the
continuum limit, in which case the diverging energies will have no effect on
the dynamics of the latter. This is what we want to show, given a few
mathematical assumptions, in the following section.

##  6 How the perturbed vacuum decouples from the physical states

Suppose that

$\displaystyle\hat{H}_{\mathrm{vac}}^{\prime\prime\prime}(\delta k)=\delta
k^{3/2}\big{(}\hat{H}_{0}^{\prime\prime\prime}+\hat{H}_{\mathrm{D,phys}}^{\prime\prime\prime}+\hat{H}_{\mathrm{C,phys}}^{\prime\prime\prime}+\delta
k^{-3/2}(\hat{H}_{\mathrm{D,vac}}^{\prime\prime\prime}+\hat{H}_{\mathrm{C,vac}}^{\prime\prime\prime})\big{)}$
(22)

is self-adjoint on some dense domain of $\mathbf{H}_{\mathrm{vac}}$ for all
$\delta k<0$ and for $\delta k\to 0$, where $\hat{H}_{0}^{\prime\prime\prime}$
is the formal continuum limit of $\hat{H}_{0}^{\prime\prime}$. Let us also
from now on refer to $\hat{H}^{\prime\prime}$ as
$\hat{H}^{\prime\prime}(\delta k)$ instead to make its dependency on $\delta
k$ explicit. Then if we consider the part of $\hat{H}^{\prime\prime}(\delta
k)$ that works only on states whose momenta sum to 0, we see that it
approximates $\delta
k^{-3/2}\hat{H}_{\mathrm{vac}}^{\prime\prime\prime}(\delta k)$ arbitrarily
well when $\delta k$ tends towards 0. It is therefore not unreasonable to
expect that for any generalized eigenstate, $\ket{\phi}$, of $\lim_{\delta
k\to 0}\hat{H}_{\mathrm{vac}}^{\prime\prime\prime}(\delta
k)=\hat{H}_{\mathrm{D,vac}}^{\prime\prime\prime}+\hat{H}_{\mathrm{C,vac}}^{\prime\prime\prime}$
with eigenvalue $E_{\phi}$, one can find a sequence of eigenstates999 Since
$\hat{H}^{\prime\prime}(\delta k)$ is discrete and bounded, its generalized
eigenstates will be actual eigenstates in $\mathbf{H}$.  of
$\hat{H}^{\prime\prime}(\delta k)$ such that the cross section of these
eigenstates on $\mathbf{H}_{\mathrm{vac}}$, when normalized at each step to
preserve the overall norm, will converge to $\ket{\phi}$, and such that their
eigenvalues, when divided by $\delta k^{-3/2}$ at each step, will converge to
$E_{\phi}$.

Let us assume this to be true for at least one such $\ket{\phi}$, and let
$(\ket{\phi_{\delta k}})_{\delta k^{-1}\in\mathbb{N}_{+}}$ be this sequence of
eigenstates of $\hat{H}^{\prime\prime}(\delta k)$ that approximates
$\ket{\phi}$. Let also $\delta k^{-3/2}E_{\phi_{\delta k}}$ be the eigenvalue
of each $\ket{\phi_{\delta k}}$, with $E_{\phi_{\delta k}}\to E_{\phi}$ when
$\delta k\to 0$.

Now, since $\ket{\phi}$ is a generalized state, it might contain Dirac delta
functions in its formula. Let us, however, first assume that it does not, and
analyze the interaction between $\ket{\phi_{\delta k}}$ and any physical
states with this assumption. This then implies that the amplitude of
$\ket{\phi_{\delta k}}$ over each point,
$(\mathbf{p}_{1},\ldots,\mathbf{p}_{n},\mathbf{k}_{1},\ldots,\mathbf{k}_{m})$,
in the discretized parameter space will go as $\delta k^{3(n+m-1)/2}$, where
the $-1$ in the exponent comes from the fact that
$\mathbf{p}_{1}+\ldots+\mathbf{p}_{n}+\mathbf{k}_{1}+\ldots+\mathbf{k}_{m}=0$
for the support of $\ket{\phi_{\delta k}}$.

Let us then consider the $\varepsilon$-almost eigenstates8 of the “physical”
part of $\hat{H}^{\prime\prime}(\delta k)$. Here we will further make the
assumption that

$\displaystyle\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}=\hat{H}_{0}^{\prime\prime\prime}+\hat{H}_{\mathrm{D,phys}}^{\prime\prime\prime}+\hat{H}_{\mathrm{C,phys}}^{\prime\prime\prime}$
(23)

is self-adjoint on some dense domain of $\mathbf{H}$. We can then expect that
for any $\varepsilon$-almost eigenstate, $\ket{\psi}$, of
$\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}$ with eigenvalue $E_{\psi}$, the
discretized version of this operator, i.e.

$\displaystyle\hat{H}_{\mathrm{phys}}^{\prime\prime}(\delta
k)=\hat{H}^{\prime\prime}(\delta k)-\,\text{vacuum-perturbing terms},$ (24)

will have a sequence of $\varepsilon_{\delta k}$-almost eigenstates, each with
the same eigenvalue, $E_{\psi}$, that converges to $\ket{\psi}$ when $\delta
k\to 0$, and where $\varepsilon_{\delta k}\to\varepsilon$ in the same limit as
well. Let $(\ket{\psi_{\delta k}})_{\delta k^{-1}\in\mathbb{N}}$ be such a
sequence, and note furthermore that we are also free to assume that each
$\ket{\psi_{\delta k}}$ includes no momentum eigenstates,
$\ket{\mathbf{p}_{1},\ldots,\mathbf{p}_{n};\mathbf{k}_{1},\ldots,\mathbf{k}_{m}}$,
where there are three or more momentum vectors that sum to exactly 0. Since
this only removes a combined state whose norm goes as $\delta k^{3/2}$, and
since $\hat{H}^{\prime\prime}_{\mathrm{phys}}(\delta k)$ approximates a nice,
smooth operator, this will only add an $O(\delta k^{3/2})$ contribution to
$\hat{H}^{\prime\prime}_{\mathrm{phys}}(\delta k)\ket{\psi_{\delta k}}$. And
therefore, it will not change the fact that $\varepsilon_{\delta k}$ can get
arbitrarily close to $\varepsilon$ when $\delta k\to 0$. With these
assumptions, we see that since $\ket{\psi}$ is an actual state in
$\mathbf{H}$, we get that the amplitude of $\ket{\psi_{\delta k}}$ over each
point
$(\mathbf{p}_{1},\ldots,\mathbf{p}_{n},\mathbf{k}_{1},\ldots,\mathbf{k}_{m})$
in the discretized parameter space will go as $\delta k^{3(n+m)/2}$.

Define then $\hat{A}^{\dagger}_{\phi_{\delta k}}$ and
$\hat{A}^{\dagger}_{\psi_{\delta k}}$ as the combined creation operators that
respectively create $\ket{\phi_{\delta k}}$ and $\ket{\psi_{\delta k}}$ from
the bare vacuum. In other words, let

$\hat{A}^{\dagger}_{\phi_{\delta k}}\ket{\,}=\ket{\phi_{\delta
k}},\quad\hat{A}^{\dagger}_{\psi_{\delta
k}^{\prime}}\ket{\,}=\ket{\psi_{\delta k}}.$ (25)

And define further $\ket{\chi_{\delta k}}$ by

$\ket{\chi_{\delta k}}=\hat{A}^{\dagger}_{\psi_{\delta
k}}\hat{A}^{\dagger}_{\phi_{\delta k}}\ket{\,}.$ (26)

We then hope to show that $\ket{\chi_{\delta k}}$ is an
$\varepsilon^{\prime}$-almost eigenstate of $\hat{H}^{\prime\prime}(\delta
k)$, where $\varepsilon^{\prime}$ can get arbitrarily small when $\delta k\to
0$, $\varepsilon\to 0$.

In order to do this, let us first note that each term of
$\hat{H}^{\prime\prime}(\delta k)$ includes an even number of
$\hat{b}_{s,\mathbf{p}}$ or $\hat{b}_{s,\mathbf{p}}^{\dagger}$-operators. This
means that if we want to commute the terms of $\hat{H}^{\prime\prime}(\delta
k)$ with the terms of $\hat{A}^{\dagger}_{\psi_{\delta k}}$, the
anticommutation relations of the $\hat{b}_{s,\mathbf{p}}$ and
$\hat{b}_{s,\mathbf{p}}^{\dagger}$-operators will work exactly as if they were
the corresponding _commutation_ relations instead. This tells us that the
difference between $\hat{H}^{\prime\prime}(\delta
k)\hat{A}^{\dagger}_{\psi_{\delta k}}$ and $\hat{A}^{\dagger}_{\psi_{\delta
k}}\hat{H}^{\prime\prime}(\delta k)$ is exactly equal to all the terms where
one or more annihilation operator of $\hat{H}^{\prime\prime}(\delta k)$ hits a
creation operator in $\hat{A}^{\dagger}_{\psi_{\delta k}}$ with the same
indices and turns into a 1. Now, if we look at the terms where _all_ (and at
least one) of the annihilation operators in $\hat{H}^{\prime\prime}(\delta k)$
hit some corresponding creation operators in $\hat{A}^{\dagger}_{\psi_{\delta
k}}$, we see that this is exactly the terms that make up

$\hat{H}_{\mathrm{phys}}^{\prime\prime}(\delta
k)\hat{A}^{\dagger}_{\psi_{\delta
k}}\ket{\,}=E_{\psi}\hat{A}^{\dagger}_{\psi_{\delta
k}}\ket{\,}+O(\varepsilon_{\delta k}).$ (27)

Here $O(\varepsilon_{\delta k})$ is a shorthand for some state with norm less
than $C\varepsilon_{\delta k}$ for some constant $C>0$. The reason why the
vacuum-perturbing terms do not contribute to these terms is that
$\hat{A}^{\dagger}_{\psi_{\delta k}}\ket{\,}$, by assumption, does not include
any momentum eigenstates where three or more momentum vectors sum to 0.
Therefore the terms that include only annihilation operators will not be able
to have all these operators hit corresponding creation operators. And of
course, since we need to have at least one annihilation operator hit a
creation operator, we can also exclude the terms that include only creation
operators. We therefore get that

$\displaystyle\hat{H}^{\prime\prime}(\delta k)\ket{\chi_{\delta k}}$
$\displaystyle=\hat{A}^{\dagger}_{\psi_{\delta
k}}\hat{H}^{\prime\prime}(\delta k)\hat{A}^{\dagger}_{\phi_{\delta
k}}\ket{\,}+E_{\psi}\hat{A}^{\dagger}_{\psi_{\delta
k}}\hat{A}^{\dagger}_{\phi_{\delta k}}\ket{\,}+\,\text{interaction
terms}\,+O(\varepsilon_{\delta k})$ (28) $\displaystyle=(\delta
k^{-3/2}E_{\phi_{\delta k}}+E_{\psi})\ket{\chi_{\delta k}}+\,\text{interaction
terms}\,+O(\varepsilon_{\delta k}),$

where we by “interaction terms” refer to all the terms where there is at least
one annihilation operator from $\hat{H}^{\prime\prime}(\delta k)$ that hits a
creation operator in $\hat{A}^{\dagger}_{\psi_{\delta k}}$ and also at least
one that hits a creation operator in $\hat{A}^{\dagger}_{\phi_{\delta k}}$. So
if we can only show that the norm of these combined “interaction terms” goes
to zero when $\delta k$ does, we obtain the result that we are looking for.

This is fortunately not too hard to do, as it turns out. As per the
assumptions above, a state,
$\ket{\mathbf{p}_{1},\ldots,\mathbf{p}_{n};\mathbf{k}_{1},\ldots,\mathbf{k}_{m}}$,
created by $\hat{A}^{\dagger}_{\phi_{\delta k}}$ will have amplitudes that go
as $\delta k^{3(n+m-1)/2}$. And a state state,
$\ket{\mathbf{p}_{1},\ldots,\mathbf{p}_{n^{\prime}};\mathbf{k}_{1},\ldots,\mathbf{k}_{m^{\prime}}}$,
created by $\hat{A}^{\dagger}_{\phi_{\delta k}}$ will have amplitudes that go
as $\delta k^{3(n^{\prime}+m^{\prime})/2}$. This means that the amplitudes of
any state,
$\ket{\mathbf{p}_{1},\ldots,\mathbf{p}_{n^{\prime\prime}};\mathbf{k}_{1},\ldots,\mathbf{k}_{m^{\prime\prime}}}$,
of $\ket{\chi_{\delta k}}$ will have amplitudes that go as $\delta
k^{3(n^{\prime\prime}+m^{\prime\prime}-1)/2}$, matching the local degrees of
freedom: $3(n^{\prime\prime}+m^{\prime\prime}-1)$. We can then analyze the
norms of the “interaction terms” using the same rule of Eq. (9) once again, as
long as we just remember that $\Delta$ will then always be plus 3 of what it
would normally be since the initial amplitude is proportional to $\delta
k^{3(n^{\prime\prime}+m^{\prime\prime}-1)/2}$ rather than to $\delta
k^{3(n^{\prime\prime}+m^{\prime\prime})/2}$.

For the Dirac interaction, the relevant terms here are those with two
annihilated particles and one created particle, as well as those with three
annihilated particles and zero created ones. In the first case, we see that we
have $\Delta/2+d=(-3+3)/2+0=0$. However, the matrix elements are proportional
to $\delta k^{3/2}$ rather than $\delta k^{0}$, which means that the combined
norm of these contributions to $\hat{H}^{\prime\prime}(\delta
k)\ket{\chi_{\delta k}}$ will vanish when $\delta k\to 0$. And if we consider
the terms where three particles annihilate each other with no particles
created, we see that we have $\Delta/2+d=(-9+3)/2+3=0$ as well, which means
that these contributions also vanish in the same limit. And for the Coulomb
interaction, we first of all have the terms with both two annihilated and two
created particles, for which we get that $\Delta/2+d=(0+3)/2+0=3/2$. Then
there are also the terms with three annihilated particles and a single
particle created, for which we have $\Delta/2+d=(-6+3)/2+3=3/2$. And lastly,
we have the terms where there are four annihilated particles and no particles
created, in which case we have $\Delta/2+d=(-12+3)/2+6=3/2$ yet again.
Therefore, since the matrix elements here are proportional to $\delta k^{3}$
rather than $\delta k^{3/2}$, the combined norm of these contributions will
also vanish when $\delta k\to 0$.

We thus obtain that

$\displaystyle\hat{H}^{\prime\prime}(\delta k)\ket{\chi_{\delta k}}=(\delta
k^{-3/2}E_{\phi_{\delta k}}+E_{\psi})\ket{\chi_{\delta k}}+O(\delta
k^{3/2})+O(\varepsilon_{\delta k}).$ (29)

We thus get that $\ket{\chi_{\delta k}}$ will indeed be an
$\varepsilon^{\prime}$-almost eigenstate of $\hat{H}^{\prime\prime}(\delta
k)$, where $\varepsilon^{\prime}$ can get arbitrarily small when $\delta k\to
0$, $\varepsilon\to 0$.

Before we move on to discuss what this means, however, we should also look at
what happens if $\ket{\phi}$ _does_ include some Dirac delta functions. Now,
we have just seen how the fact that the limit of $(\ket{\phi_{\delta
k}})_{\delta k^{-1}\in\mathbb{N}_{+}}$ includes a delta function when viewed
in $\mathbf{H}$, namely the one that makes the momenta sum to 0, results in an
interaction that vanishes in the continuum limit. And it is therefore natural
to expect that if this limit also includes additional delta functions, it will
only make the interaction vanish more rapidly, if anything. And indeed, it is
not hard to convince oneself of this. We can see that for any degree of
freedom that we remove from $\ket{\phi_{\delta k}}$ in the initial state, this
either removes a degree of freedom from the final state, in which case
$\Delta$ is invariant, or it will decrease $d$ by 1. And if $d$ decreases by
1, then even if $\Delta$ also increases by 1 at the same time, it will still
only ever reduce $\Delta/2+d$ in total. So having delta functions in
$\ket{\phi}$ will thus indeed only make the interaction terms vanish more
rapidly, if anything, and Eq. (29) will therefore still apply.

This is a great result, since it implies that the dynamics of the physical
states for $\hat{H}^{\prime\prime}(\delta k)$ converges to the dynamics of
$\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}$ when $\delta k\to 0$. The way
to see this is to note, first of all, that any state, $\ket{\Psi}$, of
$\mathbf{H}$ can be approximated arbitrarily well by a sum of a finite set of
$\varepsilon$-almost eigenstates,
$\\{\ket{\psi_{n}}\\}_{n\in\\{1,\ldots,N\\}}$, such that
$\ket{\Psi}\approx\sum_{n=1}^{N}\alpha_{n}\ket{\psi_{n}}$, where
$\alpha_{n}\in\mathbb{C}$. And for each $\ket{\psi_{n}}$ with eigenvalue
$E_{\psi_{n}}$, we can find a $\ket{\chi_{n,\delta k}}$, that is an
$\varepsilon^{\prime}$-almost eigenstate of the operator given by

$\displaystyle\hat{H}_{\phi_{\delta
k}}^{\prime\prime}=\hat{H}^{\prime\prime}(\delta k)-\delta
k^{-3/2}E_{\phi_{\delta k}},$ (30)

with an eigenvalue of $E_{\psi_{n}}$ as well. Thus, for any given final state,
$\ket{\Phi}$, approximated by
$\ket{\Phi}\approx\sum_{n=1}^{N}\beta_{n}\ket{\psi_{n}}$, we get that

$\displaystyle\bra{\Phi}e^{-i\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}}\ket{\Psi}$
$\displaystyle\approx\sum_{m=1}^{N}\sum_{n=1}^{N}\beta_{m}^{*}\alpha_{n}\bra{\psi_{m}}e^{-i\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}}\ket{\psi_{n}}\approx\sum_{n=1}^{N}\beta_{n}^{*}\alpha_{n}e^{-iE_{\psi_{n}}}$
(31)
$\displaystyle\approx\sum_{m=1}^{N}\sum_{n=1}^{N}\beta_{m}^{*}\alpha_{n}\bra{\chi_{m,\delta
k}}e^{-i\hat{H}_{\phi_{\delta k}}^{\prime\prime}}\ket{\chi_{n,\delta k}}$
$\displaystyle\approx\bra{\,}\hat{A}_{\phi_{\delta
k}}\hat{A}_{\Psi}e^{-i\hat{H}_{\phi_{\delta
k}}^{\prime\prime}}\hat{A}^{\dagger}_{\Psi}\hat{A}^{\dagger}_{\phi_{\delta
k}}\ket{\,},$

where $\hat{A}^{\dagger}_{\Psi}$ is the combined creation operator that
creates $\ket{\Psi}$. When we let $\delta k$ tend towards 0, and let
$N\to\infty$ and $\varepsilon,\varepsilon_{\delta k}\to 0$ as well, we
therefore see that any experiment in a system governed by
$\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}$ can simulated to arbitrary
precision using $\hat{H}_{\phi_{\delta k}}^{\prime\prime}$ instead, and vice
versa!

The dynamics of the physical states in a system governed by
$\hat{H}^{\prime\prime}(\delta k)$ thus converges to the dynamics of a system
governed by $\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}$. And since
$\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}$ includes no vacuum-perturbing
terms, we thus get the the perturbed vacuum decouples from all the physical
states in the continuum limit.

##  7 The Lorentz covariance of the resulting Hamiltonian

We have not yet argued why the procedure used in this paper for the Dirac sea
reinterpretation is a correct approach. This is, however, not hard to justify,
given our assumption above that in each of the steps of this procedure, going
from
$\hat{H}\to\hat{H}^{\prime}\to\hat{H}^{\prime\prime}\to\hat{H}^{\prime\prime\prime}_{\mathrm{phys}}$,
we can always impose an ultraviolet cutoff,
$|\mathbf{p}|,|\mathbf{k}|<\Lambda$, and expect the dynamics to converge when
we let $\Lambda\to\infty$. This assumption allows us to make Dyson expansions
of the time evolution operators $\exp(-i\hat{H}^{\prime})$ and
$\exp(-i\hat{H}^{\prime\prime})$, yielding us two scattering matrices that can
both be written out purely as a linear combination of creation and
annihilation operators. The point is then to note that all such terms of
$\exp(-i\hat{H}^{\prime\prime})$ will correspond to a similar term of
$\exp(-i\hat{H}^{\prime})$, with the same coefficients, since the Dirac sea
reinterpretation does not change the commutation and anticommutation relations
between the various creation and annihilation operators. This means that the
scattering matrix of $\hat{H}^{\prime\prime}$ can be derived directly from
that of $\hat{H}^{\prime}$. So if the dynamics of $\hat{H}^{\prime}$ become
Lorentz-covariant in the continuum limit, and when we also let
$\Lambda\to\infty$ at the same time, then the dynamics of
$\hat{H}^{\prime\prime}$ must also become Lorentz-covariant in the same
limit.101010 This argument can also be made even more compelling if one
changes the basis of the creation and annihilation operators to a position
basis.

And given that $\hat{H}$ of Eq. (1) is generally regarded as Lorentz-covariant
to begin with,111111An overall argument for this can also be seen in Damgaard
[1]. i.e. with formulas of $\hat{H}_{D}$ and $\hat{H}_{C}$ as shown in
Appendix A, we can therefore propose the following Hamiltonian of QED, call it
$\hat{H}_{QED}$, given formally on $\mathbf{H}$ by

$\hat{H}_{QED}=\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}.$ (32)

Here we recall that $\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}$ is the
operator obtained from $\hat{H}$ by first turning
$\hat{b}_{3}^{\dagger}(\mathbf{p}),\hat{b}_{4}^{\dagger}(\mathbf{p})$ into
$\hat{b}_{3}(-\mathbf{p}),\hat{b}_{4}(-\mathbf{p})$ and vice versa, and then
removing all the vacuum-perturbing terms afterwards.

##  8 On confirming the mathematical assumptions

In order to get to $\hat{H}_{QED}$, we had to make some mathematical
assumptions, in particular about the self-adjointness of
$\hat{H}_{\mathrm{vac}}^{\prime\prime\prime}(\delta k)$ and of
$\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}\equiv\hat{H}_{QED}$, as well as
their discretized versions. And we have also assumed that $\hat{H}$ is self-
adjoint on some domain to begin with. These propositions are by no means
trivial to show. However, Damgaard [2] does introduce a way that one might be
able to show the self-adjointness of such operators. That paper thus proves
the self-adjointness of a simplified Dirac interaction operator (with only one
fermion and only one spin component), and the technique used for this proof
might be extendable to more general operators as well. Thus, there is hope
that these propositions can be shown in a near future.

Additionally, there is also the matter of showing that we can always impose
ultraviolet cutoffs for these operators, and then expect the
$\varepsilon$-almost eigenstates, along with their corresponding eigenvalues,
to converge to those of the original operators. This is not a trivial task
either, but it might very well be possible. Appendix B briefly discusses a
potential approach for doing this.

##  9 Conclusion

We went through the process of making the Dirac sea reinterpretation for the
Hamiltonian of QED, and argued that while the energies of the perturbed vacuum
seem to diverge in the continuum limit, there likely exist sequences of
eigenstates of the perturbed vacuum which converge, at least given some
mathematical assumptions about the self-adjointness of certain operators. We
then showed that the interaction and interference between such perturbed
vacuum states and the “physical states” of the system, i.e. the states which
we might prepare and investigate in any realistic experiment, vanish in the
continuum limit, namely due to the fact that the perturbed vacuum states only
occupy a vanishingly small part of the momentum space. We discussed how this
means that the physical states of the discretized and reinterpreted
Hamiltonian, $\hat{H}^{\prime\prime}$, can then be simulated to arbitrary
precision by an operator,
$\hat{H}_{\mathrm{phys}}^{\prime\prime\prime}\equiv\hat{H}_{QED}$, with all of
the vacuum-perturbing terms removed. Finally, we argued that this operator
then ought to be Lorentz-covariant, at least if we also make the assumption
that the Hamiltonians in this paper can all be simulated to arbitrary
precision by ones with ultraviolet cutoffs imposed on them. We can therefore
propose this $\hat{H}_{QED}$ as a promising candidate for the Hamiltonian of
the full theory of QED.

##  Appendix A The Hamiltonian of QED before the Dirac sea reinterpretation

According to Damgaard [1], the continuous Hamiltonian of QED before the Dirac
sea reinterpretation, on a Fock space of photons and some Dirac spinor
fermions with a charge of $q_{\mathrm{F}}$, is formally given by

$\displaystyle\hat{H}=\hat{H}_{0}+\hat{H}_{\mathrm{D}}+\hat{H}_{\mathrm{C}},$
(33)

with

$\displaystyle\hat{H}_{0}$
$\displaystyle=\sum_{\lambda=1}^{2}\int\frac{d\mathbf{k}}{(2\pi)^{3}}\,|\mathbf{k}|\hat{a}_{\lambda}^{\dagger}(\mathbf{k})\hat{a}_{\lambda}(\mathbf{k})+\int\frac{d\mathbf{p}}{(2\pi)^{3}}\,E_{\mathbf{p}}\Big{(}\sum_{s=1}^{2}\hat{b}_{s}^{\dagger}(\mathbf{p})\hat{b}_{s}(\mathbf{p})-\sum_{s=3}^{4}\hat{d}_{s}^{\dagger}(\mathbf{p})\hat{d}_{s}(\mathbf{p})\Big{)}$
(34) $\displaystyle\hat{H}_{\mathrm{D}}$
$\displaystyle=-q_{\mathrm{F}}\sum_{\lambda=1}^{2}\int\frac{d\mathbf{k}\,d\mathbf{p}}{(2\pi)^{6}}\,\frac{1}{\sqrt{2|\mathbf{k}|}}\big{(}\hat{a}_{\lambda}(\mathbf{k})+\hat{a}_{\lambda}^{\dagger}(\mathbf{-k})\big{)}\hat{\bar{\psi}}(\mathbf{p}+\mathbf{k})\gamma^{\mu}\epsilon_{\mu\lambda}(\mathbf{k})\hat{\psi}(\mathbf{p}),$
(35) $\displaystyle\hat{H}_{\mathrm{C}}$
$\displaystyle=\int\frac{d\mathbf{k}\,d\mathbf{p}_{1}\,d\mathbf{p}_{2}}{(2\pi)^{9}}\,\frac{q_{\mathrm{F}}^{2}}{2\mathbf{k}^{2}}\big{(}\hat{\bar{\psi}}(\mathbf{p}_{1}+\mathbf{k})\otimes\hat{\bar{\psi}}(\mathbf{p}_{2}-\mathbf{k})\big{)}(\gamma^{0}\otimes\gamma^{0})\big{(}\hat{\psi}(\mathbf{p}_{1})\otimes\hat{\psi}(\mathbf{p}_{2})\big{)},$
(36)

and with

$\displaystyle\hat{\psi}(\mathbf{p})$
$\displaystyle=\frac{1}{\sqrt{2E_{\mathbf{p}}}}\Big{(}\sum_{s=1}^{2}u_{s}(\mathbf{p})\hat{b}_{s}(\mathbf{p})+\sum_{s=3}^{4}v_{s}(\mathbf{p})\hat{d}_{s}(\mathbf{p})\Big{)},$
(37) $\displaystyle\hat{\bar{\psi}}(\mathbf{p})$
$\displaystyle=\frac{1}{\sqrt{2E_{\mathbf{p}}}}\Big{(}\sum_{s=1}^{2}\hat{b}_{s}^{\dagger}(\mathbf{p})\bar{u}_{s}(\mathbf{p})+\sum_{s=3}^{4}\hat{d}_{s}^{\dagger}(\mathbf{p})\bar{v}_{s}(\mathbf{p})\Big{)}.$

Most of the symbols here are either explained above or are well-known to any
readers familiar with QFT. As mentioned in Footnote 3, we have simply used the
symbols $\hat{b}_{3}$ and $\hat{b}_{4}$ instead of $\hat{d}_{3}$ and
$\hat{d}_{4}$ in the body of this paper in order to write some of the
expressions more compactly. Apart from this, the most important thing to note
is just that the $\otimes$ symbols in Eq. (36) represent the matrix direct
product. Since $u_{s},v_{s}$ are 4-by-1 column vectors (namely the spinor
solutions to the Dirac equation) and $\bar{u}_{s},\bar{v}_{s}$ are 1-by-4 row
vectors, the integrand in this formula thus consists of a 1-by-16 row vector
multiplied by a 16-by-16 matrix, $\gamma^{0}\otimes\gamma^{0}$, multiplied by
a 16-by-1 column vector.

For further elaboration on these formulas, see Damgaard [1].

##  Appendix B On showing that the ultraviolet cutoffs are valid

In this appendix, we will briefly discuss a potential approach for showing
that the ultraviolet cutoffs used in this paper are valid, in the sense that
the $\varepsilon$-almost eigenstates and eigenvalues converge when the cutoffs
are lifted again.

Damgaard [2] uses an approach for constructing states that are mapped into
normalizable states by the Hamiltonian in which transitions that create a
particle (specifically a photon) are canceled by transitions that annihilate a
particle. And in that paper, the momenta of the annihilated particles in these
latter transitions take values in a set, $E$, whose radii grow very rapidly
with increasing momenta and increasing photon number. However, one might also
be able to use transitions where the momenta of the annihilated particles take
values from an $E$ set with constant radii instead.

Now, the trouble that we are facing for this task is that when we make an
ultraviolet cutoff, we might ruin an $\varepsilon$-almost eigenstate’s ability
to cancel its own “creation transitions,” as we might call them, since this
might remove the “annihilation transitions” that the state otherwise used for
this. So if we let the cut-off version of the Hamiltonian work on a given
$\varepsilon$-almost eigenstate, it might actually (somewhat
counterintuitively) result in state that is no longer normalizable.

But if the created particles can indeed also be canceled by using an $E$ set
with constant radii instead, the idea is then that for any
$\varepsilon$-almost eigenstate, $\ket{\psi}$, of the original operator, one
might be able to construct an $\varepsilon^{\prime}$-almost eigenstate,
$\ket{\psi_{\Lambda}}$, of the cut-off version by taking $\ket{\psi}$ and then
adding a small state (a small “tail,” so to speak) that conforms to this $E$.
The job of this small “tail” is then to cancel all the “creation transitions”
that are no longer canceled by the “annihilation transitions” due to the
cutoff. One might then be able to show that this “tail” can be made
arbitrarily small compared to the norm of $\ket{\psi}$ when
$\Lambda\to\infty$, and thus that the sequence of $\ket{\psi_{\Lambda}}$’s
will converge to $\ket{\psi}$, with $\varepsilon^{\prime}\to\varepsilon$ and
with eigenvalues of the sequence converging to that of $\ket{\psi}$ as well.

## References

  * [1] M. J. Damgaard, An alternative derivation of the Hamiltonian of quantum electrodynamics, viXra:2210.0044v2 [quant-ph].
  * [2] M. J. Damgaard, Self-adjointness of a simplified Dirac interaction operator without any cutoffs, arXiv:2311.12870 [quant-ph].
  * [3] B. C. Hall, Quantum Theory for Mathematicians (Springer, New York, 2013).
  * [4] S. Weinberg, The Quantum Theory of Fields, Vol. I: Foundations (Cambridge University Press, Cambridge, 1995).
  * [5] M. Srednicki, Quantum Field Theory (Cambridge University Press, Cambridge, 2007).
  * [6] T. Lancaster and S. J. Blundell, Quantum Field Theory for the Gifted Amateur (Oxford University Press, Oxford, 2014).