  Simulating beam loading in radiofrequency accelerating structures is critical
for understanding higher-order mode effects on beam dynamics, such as beam
break-up instability in energy recovery linacs. Full wave simulations of beam
loading in radiofrequency structures are computationally expensive, while
reduced models can ignore essential physics and can be difficult to generalize.
We present a self-consistent algorithm derived from the least-action principle
which can model an arbitrary number of cavity eigenmodes and with a generic
beam distribution.
