Using the typical expansion functions (1-forms/edge-elements for E, and 2-forms/facet-elements for B) the formulations are basically the same after spatial discretization and you'd expect more or less the same accuracy. I do think they express slightly different opinions about time integration.
The mixed E/B formulation nudges you in the direction of ...
Regarding performance, Python is definitely the bottleneck. I have experienced the same issue with a 2D Euler code I had developed, even with vectorised operations everywhere possible. It was actually even worse, as I was using solve_ivp time schemes which reallocated memory at every step... You can try and profile your code to see where the bottlenecks are. ...
I am by no means experienced with the wave equation, but I think the issue comes from the imposition of the periodic BCs.
The periodic boundary conditions can be imposed by using ghost points: you do as if you were considering an extended system which, in Python terms, would have the state vector:
u_extend=[u[-1], u, u, ..., u[M-1], u[M], u]
If you use Julia's DifferentialEquations.jl it can automatically fix the index of your equations. A tutorial of this is shown in ModelingToolkit.jl. For example, we can write down an index-3 DAE:
function pendulum!(du, u, p, t)
x, dx, y, dy, T = u
g, L = p
du = dx
du = T*x