{'fraction_non_alphanumeric': 0.04035398230088495, 'fraction_numerical': 0.010226155358898722, 'mean_word_length': 4.478560655031243, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Simulation of Electromagnetic Solitons and their Interaction with Matter

A suitable correction of the Maxwell model brings to an enlargement of the space of solutions, allowing for the existence of solitons in vacuum. We review the basic achievements of the theory and discuss some approximation results based on an explicit finite-difference technique. The experiments in two dimensions simulate travelling solitary electromagnetic waves, and show their interaction with conductive walls. In particular, the classical dispersion, exhibited by the passage of a photon through a small aperture, is examined.


Introduction
It is known, although rarely stated in books, that Maxwell equations in empty space do not admit finite-energy solitary waves among their solutions. By this we mean smooth electromagnetic waves with compact support, travelling along straight-lines at the speed of light. One of the reasons to explain this fact can be attributed to the linearity of the Maxwell model, not allowing the due focussing of the signal on a constrained path. Such an impediment has stimulated the research of alternative nonlinear models, mostly based on modification of the Lagrangian (see for instance [5]), granting the existence of soliton-like solutions and justifying in this way electromagnetic phenomena such as photons.
A step ahead in the comprehension of electromagnetic solitary waves has been done in [11] and [12], where, in the framework of a self-consistent theory, explicit analytic solutions are carried out. The main argument is that the classical equations of the electromagnetism in vacuum are not capable to follow the evolution of finite-energy wave-fronts in the proper way (i.e., the one described for instance by the Huygens principle); this is due to the difficulty to impose both divE = 0 and divB = 0 at the same time on each point of the same front. As a matter of fact, if we define the wave-front as the surface envelope of the vector fields, one can easily show that requiring that the front evolves along its normal direction is incompatible with the two free-divergence conditions taken together. Therefore, the relation divE = 0 has been dropped, without this necessarily implying the existence of point-wise electric charges. Here, we circumscribe our exposition by outlining some of the peculiarity of the new approach, while we refer to [12] for a detailed explanation of the origin and the possible physical implications.
In particular, in this paper we devote our attention to some numerical simulations. It is soon evident that the removal of relation divE = 0 is important for the construction of numerical algorithms in general, because most of the difficulties encountered in simulations are indeed consequence of the imposition of such a constraint. These include the efforts made to build approximation spaces satisfying some divergence-free conditions (see for instance [2], [4], [7], [9], [14], [17]) or divergence corrections techniques (see for instance [15], [18], [19]). Some modifications of the Maxwell model, which may be in some way assimilated to the ones considered here, have been proposed at numerical level, in order to set up stable schemes (see [22], [24]), to handle boundary conditions (see [1], [3]), or for the treatment of wave propagation in linear non-dispersive lossy materials. These techniques are mainly adopted to overcome numerical troubles and are not intended to modifying the Maxwell model itself, as we are doing here. A survey of methods is given in [20], where the reader can also find an updated list of references.
For simplicity, we will use explicit finite-differences (in particular the Lax-Wendroff method). The aim is to validate the theory in [12] with a series of simple experiments. To this purpose we follow the evolution of solitary waves and see how they react when encountering matter. We will mainly concentrate to the physical appearance and significance of the solutions. Thus, no theoretical issues will be discussed. Moreover, the discretization parameters will be small enough to get sufficiently accurate approximations, finalized to recover qualitative information. In no way we shall claim that our numerical approach is competitive, being conscious that the one proposed here can be certainly improved upon, as far as per-formances versus costs are concerned.
The author would like to thank A. Ugolini (see [23]) for setting up the computational code.

The model equations
By denoting with c the speed of light, the classical Maxwell equations, in void three-dimensional space, in absence of electric charges, take the form: where the two fields E and cB have the same dimensions.
In [12], part of the analysis is devoted to electromagnetic free-waves. These are solutions of the following set of model equations: where ρ = divE, and V is a velocity vector field satisfying |V| = c. The field V is oriented as the vector field E × B. Note that relation (2.8) is certainly satisfied for all electromagnetic waves where E is orthogonal to B and |E| = |cB|. These requirements are standard. In addition, we expect the new set of equations to admit a large space of solutions. We will check later that, contrary to the Maxwell model, solitary waves with compact support, as well as perfect spherical waves, are now among the solutions. On the other hand, if divE is relatively small (as it actually happens in many practical circumstances), then (2.5) is as accurate as the corresponding standard Maxwell equation. Therefore, we expect the new model to be consistent with the existing ones, for a broad range of applications.
It should be clear that the added term ρV = cρJ, where |J| = 1, has the meaning of a current density flowing in the direction of the rays at speed c (Ampère term). We assume this to be true even if there is no presence of point-wise electric charges. As documented in [12], such a current is part of the wave itself. Moreover, equation (2.8) actually characterizes free-waves, since its says that there are no external 'forces' acting on the wave-fronts, as a result of external perturbations. It is a sort of Lorentz law (where the force turns out to be zero), in which moving charges are replaced by a balance of pure vector fields. We will see later, in section 4, how equation (2.8) should be modified in order to handle situations in which the wave interacts with matter.
It is important to observe that a conservation principle holds. In fact, by taking the divergence of equation (2.5), we come out with the following continuity equation: In addition, using that V is orthogonal to both E and B, one can easily prove the classical Poynting relation: Such properties confirm that what we are doing has physical meaning. It can be proven that the electromagnetic free-waves described by the new set of equations are perfectly compatible with the Huygens principle and the eikonal equation (see for instance [6]). In fact, if, for example, V is irrotational (i.e.: the rays are straight-lines), then V = ∇Ψ for some potential Ψ, so that the relation |V| = c is equivalent to |∇Ψ| = c, which is the eikonal equation. Moreover, (2.5) follows from minimizing the standard Lagrangian of classical electromagnetism, after imposing the constraint A = ΦJ to the electromagnetic potentials A and Φ (see [13] and [12], section 2.4). The covariant version of the equations and their invariance under Lorentz transformations are also discussed in [12], chapter 4.

The discretization method for 2-D solitons
Explicit expressions of wave-packets, moving at speed c and solving the system of equations (2.5)-(2.8), are available. The electromagnetic signals carried by these fronts can be very general. For the sake of simplicity, only examples in two dimensions will be examined here. This does not obscure, however, the general validity of the formulation.
We start by discussing a simple case. In the 3-D space with cartesian coordinates (x, y, z), we can define the following fields: Now, the continuity equation (2.9) becomes the simple transport equation: In practice, the example is two-dimensional since there is no dependency on the variable y. The fields describe the evolution of a solitary wavepacket shifting along the axis z. This is a group of parallel transversal fronts, travelling in the same direction V and modulated in intensity by the function g. It is important to remark that the fields E and B satisfy Maxwell equations (including (2.2)) only if the function f is constant (this check can be done by directly differentiating the expression in (3.1) and substituting in the model equations). This corresponds to a plane wave of infinite energy, which is the only Maxwellian front allowed by the above setting. If we want the energy to be finite, the only possible choice is f = 0, implying E = 0. The situation is different if we take into account equation (2.5). Now, as the reader can easily check by direct computation, the fields in (3.1) are solutions, and the two functions f and g can be completely arbitrary. In particular, f and g may have compact support. Note that, for simplicity, the wave is unbounded in the y direction, hence, in this case the energy of the wave-packet is not finite. Examples where the fronts are also bounded in the y direction can be easily constructed (recall to enforce the condition divB = 0), but the problem becomes three-dimensional and here we would prefer to avoid this complication.
Let us discuss a specific example. We start by fixing a point (x 0 , z 0 ). At time t = 0, we can assign the distribution: inside the rectangle centered at (x 0 , z 0 ) with sides of length 2τ and 2σ. Outside the same rectangle E is required to vanish. Accordingly, field cB is supposed to be orthogonal to the plane (x, z) and with the same intensity as E. When time passes, this signal-packet moves at speed c in the direction of the axis z: We do not need to use any approximation technique in order to follow the evolution of such a soliton, since the corresponding analytic expression exactly solves the set of equations (2.5)-(2.8). However, in view of more complex applications, let us introduce the Lax-Wendroff method. More precisely, we work with its 2-D version applied to a system of pure hyperbolic equations, as described in [21], section 4.9.
Let us suppose that, for some real constant d × d matrices M x and M z , our system can be written in the following form: where the d-components vector u depends on x, z and t.
We denote by (x i , z j ) the nodes of a uniform grid (of step-size h) on the plane (x, z). The time-discretization parameter is denoted by ∆t. Based on the stencil shown in figure 1, we pass from the approximated solution u k at time t k to the one at time t k+1 = t k + ∆t. This is done according to the scheme: where λ = ∆t/h and I is the identity matrix. Let us examine the case of a wave with the electric field laying on the plane (x, z). Therefore, the magnetic field, which is orthogonal to that plane, has only one component different from zero. The condition (2.4) is always satisfied, since our functions do not depend on the variable y. We organize the vector u in the following way (d = 3): For the classical Maxwell system (see (2.1) and (2.3)) it is enough to define the matrices: in order to recover (3.5). It is important to remark, however, that the condition (2.2) is not enforced (this should be done at the initial time, although, as we mentioned, fronts having divE = 0 and divB = 0 are extremely rare and out of interest for our applications). Note that the eigenvalues of both matrices are ±c, confirming that the system is of hyperbolic type. Regarding the modified system (see (2.5) and (2.6)), the above matrices have to be replaced by: Moreover, we require relation (2.8) to be verified. In other terms, V must be defined as: c(E × B)/|E × B|. Note that this time we do not need to satisfy (2.2). Finally, let us observe that, for both the matrices in (3.8), the modulus of the maximum eigenvalue turns out to be equal to c.
The matrices in (3.8) are not constant. Nevertheless,we will apply the Lax-Wendroff method to the new system also in this circumstance. Therefore, when implementing (3.6) at the k-th step, the entries of the matrices M x and M z will be evaluated at the central point of the stencil of figure 1. We are aware of the fact that this correction may deteriorate a bit the convergence properties of the method. In order to guarantee stability, we must impose the following constraint (CFL condition): In our experiments we always used the maximum ∆t allowed by (3.9), i.e.: ∆t = h/2c √ 2. With this choice we never encountered problems concerning stability.
We start by testing the discretization method on the soliton corresponding to the initial conditions in (3.3), with σ = τ . The initial condition for the discrete scheme is obtained by interpolation at the grid-points. It is interesting to note that, with this setting, the second component of E remains zero during the evolution. This is true both for the exact solution (see (3.4)) and the discrete one. As a matter of fact, the scheme (3.6) follows from repeated applications of the three matrices I, M x and M z . In this situation we have V 1 = 0 and V 3 = c. Every time we apply the matrix M x to a vector of the form ( * , * , 0), we get a vector of the same form, since E 1 = cB 2 . The same is true regarding the matrix M z . It is also interesting to note that relation (2.8) is preserved for any k. We can see in figure 2 the moving soliton at different times. The discretization parameters have been taken small enough to get rid of numerical disturbances. In other words, the parameter σ is large enough to allow the inclusion of a sufficiently high number of grid-points in the support of the wave (there are 13 × 13 nodes in this example). It is interesting to remark that the scalar function E 1 satisfies the wave equation: Note instead that the wave equation ∂ 2 ∂t 2 E = c 2 ∆E does not hold, since in (3.10) we are omitting the derivatives with respect to the variable x. On the other hand, we are not approximating a trivial transport problem for the unknown E 1 , but a coupled hyperbolic system involving both E and B. This is not a trivial remark, but a rather important achievement. Some longitudinal numerical dissipation is visible. It can be suitable reduced by using more appropriate numerical techniques. It is important however to remark that no numerical viscosity is introduced in the direction transverse to the motion, so that the solitary wave proceeds maintaining a well defined width. This is quite an important result. Suitable Sommerfeld type boundary conditions are assumed at the outflow boundary, in order to allow a smooth exit of the wave from the computational domain.
A completely different situation (see figure 3) is observed when using the same initial datum (3.3) and evolving the wave with the help of the two linear Maxwell equations (2.1) and (2.3), using the matrices in (3.7). In this case the condition (2.2) is not verified, since such a constraint is not satisfied by the initial datum and it is not enforced by the numerical scheme. As the reader can see, part of the solution is diffused all around, despite the fact that, in the initial dislocation, the vectors V are all oriented in the horizontal z direction. Surprisingly, the most of the energy does not move at all. Our guess is that the condition divE = 0 generates a kind of stationary charge. Of course, this behavior is not physically acceptable. Nevertheless, there is no way to modify the initial datum, maintaining f and g with compact support. In fact, as mentioned before, one can verify that in (3.1) the only possible initial condition compatible with (2.2) is E = 0. As far as solitons are concerned, this indicates that the poor performances of the classical equations of electromagnetism are mainly due to the total lack of initial compatible conditions. Figure 4: The first picture shows the level sets of |E|, for some initial datum of the form given by (3.11). The two successive enlargements show the initial distributions of E and V in the plane (x, z).
A less trivial example is obtained by taking an initial datum of the form: whose shape is given in figure 4. The magnetic field is along the direction of the y-axis and such that |E| = |cB|. The corresponding velocity field V is radially distributed with respect to a prescribed point. With this setting, during the evolution one gets a semi-cylindrical wave, simulated by the numerical scheme as shown in figure 5. There is a decay in intensity of the electromagnetic fields, due to the fact that the signal is spread out on front surfaces having increasing magnitude. This is in agreement with the equation of continuity and the energy preservation rules (see (2.9) and (2.10)), automatically imposed by the governing equations. The approximation is not excellent, since it is polluted by a secondary wave of small amplitude evolving backwards. We think that this is due to the poor performances of the Lax-Wendroff scheme. In fact, in the non-convex part of the support of the wave some oscillations are generated. They modify the sign of the electric field, but not the one of the magnetic field, so that the associated velocity field V (proportional to the Poynting vector E × B) locally changes orientation. This anomalous part is then correctly simulated by the numerical method, but, unfortunately it develops along a wrong path. This is a minor effect that should be however kept in mind when designing an alternative code.
After testing our algorithm on these simple examples, we are now ready to deal with a more serious situation.

Constrained waves
A generalization of the model is obtained by modifying equation (2.8) in the following way (see [12], chapter 5): where the constant µ > 0 is a charge divided by a mass, ρ = divE and p is a kind of pressure. As customary, the substantial derivative G = D Dt V is defined as D Dt V = ∂ ∂t V+(V·∇)V, so that G turns out to be an acceleration. Basically, if V is the normalized vector field tangent to the light rays, then the field G gives a measure of their curvature. Equation (4.1), even if there are no classical moving charges, can be assimilated to the Lorentz law. If G and p are zero, then one is dealing with a free electromagnetic wave, and the corresponding rays are straightlines. Actually, equation (2.8), where G = 0 and p = 0, says that the evolution of the wave is free from constraints, that is, there are no external factors ('forces') acting on it. When, for some reasons, (2.8) is not satisfied, then pressure develops. At the same time D Dt V is different from zero, so that V changes direction (the rays are curving) and the electromagnetic wave-fronts locally follow the evolution of the new normalized Poynting vector J = V/c. Thus, when G is different from zero, the wave is no longer free, and, following [12], it will be called constrained wave. Such a situation happens when the wave is subjected to external electromagnetic fields. This is true for instance during the interaction with matter at atomic level, as in reflection or diffraction. Therefore, our model equations are now able to provide the coupling between the curvature of the rays and the motion of the wave-fronts. Note that, during the change of trajectory of the rays, depending on the context, the polarization may also vary.
It is evident that equation (4.1) is inspired by the Euler equation for inviscid fluids (with an electromagnetic type forcing term). Put in other words, the vector field V evolves according to the laws of fluid dynamics. At the same time, it carries, on the transversal manifold, the two orthogonal fields E and B. In this peculiar way, the system of equations (2.5)-(2.6)-(2.7)-(4.1) describes the evolution of the triplet (E, B, V). We remark that our model is a system of vector equations. Therefore, although our analysis is based on simple examples, it has nothing to do with the equations modelling scalar solitons (see for instance [10]).
The exact physical meaning of (4.1) is carefully explained in [12], where an equation of state involving p, that uses the curvature of the space-time geometry, is also introduced. For the sake of simplicity, we do not discuss further the reasons of this construction and we refer to [12] for the details. Numerical experiments on constrained waves, concerning photons trapped in toroid shaped regions, are illustrated in [8].

Some experiments on diffraction
Our next step is to analyze the interaction of an incoming soliton with a perfectly conducting plane wall. We would like to discuss some experiments concerning diffraction, so that there is a hole in the wall and our soliton will be forced to pass through it. The barrier is simulated by imposing suitable boundary conditions, so that the impact is immediate and discontinuities are generated. According to figure 6, we expect the soliton to behave as a free-wave (hence G = 0 and p = 0 in (4.1)) before reaching the aperture. The reaction of the wall brings to an instantaneous diffusion of the signal, as also observed in real-life experiments with fluids or electromagnetic waves.
After the passage, the soliton continues his journey behaving as a free wave, following however non-parallel characteristic straight-lines.
Assuming that the path of the incoming wave is orthogonal to the boundary, there we impose V = 0 (no slip condition). Hence, for an instant, the last term on the right-hand side of equation (2.5) disappears. This way of imposing boundary conditions is suggested in [12], chapter 5, where the following Bernoulli type equation is deduced: Heuristically, when V goes to zero at the boundary, the pressure raises from zero to infinity, in order to maintain the energy balance. One can actually check that divV > 0 after the impact, bringing to a diffusive effect (with no dissipation of energy, however). As we shall see from the experiments, these assumptions seem to be correctly posed. From the numerical point of view the problem is stiff, therefore some unavoidable perturbations will be observed. A more realistic behavior could be obtained by smoothing the effects of the boundary. This could be simulated by slowing down the incoming soliton by varying V in a continuous way. We did not investigate however this possibility. In figures 7 and 8, we can see the evolution of the norm of the electric field. The width of the incoming soliton includes 32 grid points, while the size of the hole corresponds to 28 grid points in figure 7 and 26 grid points in figure 8. As expected, the no slip condition at the boundary produces a diffusive behavior. After hitting the obstacle, the development of the soliton is similar to that of figure 5 (with much less emphasis, the fronts assume a circular asset, and we know from the experiments of section 3 that the scheme can support this situation). There is no loss of energy or introduction of transversal dissipation. The spreading is a consequence of the new geometrical distribution of the rays, following non-parallel trajectories (in other words, as in figure 6, the stream-lines associated to the velocity vector V are diverging).
There are little wiggles scattered all around. These are mainly due to the exceeding part of the wave, reflected back by the wall. Therefore, these disturbances might be (in part) physically consistent. Meanwhile, they can also point out the limit of the algorithm and suggest better ways to handle boundary conditions (especially for the portion involved in the reflection process).
Finally, in figure 9, we tested the case when the incoming soliton is not symmetrically centered with respect to the aperture. Of course, we could show a multitude of other scattering experiments, regarding wave-packets of arbitrary shape interacting with a variety of obstacles. However, the purpose here was to provide a general qualitative insight. Figure 9: Evolution of a solitary wave-packet passing through the hole of a conductive wall. The aperture is now asymmetrically placed with respect to the trajectory of the wave-packet. Some minor disturbances are present. To some extent, they can be physically acceptable, but they are also due to the diffusive effects of the Lax-Wendroff method, contrasted by the sharp boundary constraints.
These experiments can certainly be ameliorated by using more appropriate approximations methods. Rather than accurate results, the interest here was to show that, with the help of the model equations (2.5)-(2.6)-(2.7)-(4.1), an electromagnetic wave can be treated exactly as a material fluid. We would like to remark once again that, although the pictures presented here are related to the evolution of the scalar quantity |E|, we are not solving a trivial scalar equation. Our solutions are electromagnetic localized emissions, presenting both the features of waves and particles (photons). Their internal structure is perfectly determined in terms of electric and magnetic fields. In the more complex and realistic 3-D case, they can even change polarization, reacting in this way like true electromagnetic entities, but with energies concentrated in finite regions of space. We believe that these peculiarities are important both for physical implications and technical applications.