📚 Physics guide

This section provides the theoretical foundations behind Wakis, offering a clear and concise explanation of the physics principles and numerical methods that power the simulation engine.

We begin by revisiting Maxwell’s equations in their integral form and explaining how they are discretized using the Finite Integration Technique (FIT) on a structured Cartesian grid. The formulation naturally leads to Maxwell Grid Equations (MGEs), which are solved in time using a leapfrog scheme. The fields and material properties are represented on a Yee-type staggered lattice, and anisotropic or spatially varying materials are handled via sparse metric tensors.

Boundary conditions (PEC, PMC, periodic, PML) are discussed along with the treatment of sources and initial conditions. We also highlight how the solver supports geometry import and sub-pixel smoothing for embedded CAD models.

Finally, we cover the implementation aspects, including GPU acceleration with CuPy and parallelization with mpi4py, enabling high-resolution 3D simulations across multiple devices.

1. Introduction

🎯 Motivation

In modern accelerators, precise knowledge of beam-coupling impedance and wakefields is essential to ensure beam quality, mitigate heating, and optimize component design. Analytical methods, while powerful, often fall short for realistic 3D geometries — this is where full electromagnetic solvers like Wakis become indispensable.

📖 Background

The evaluation of beam-coupling impedance and wakefields is fundamental to the design and operation of particle accelerators. As charged particle bunches traverse beamline components, they interact with material boundaries and geometric discontinuities, generating electromagnetic fields collectively known as wakefields.

These wakefields can:

  • Degrade beam quality and cause coherent instabilities

  • Induce energy spread and emittance growth

  • Lead to beam-induced heating and power loss in accelerator components

In the frequency domain, the response of a structure to the passing beam is quantified by its beam-coupling impedance — a complex-valued function that encapsulates how each device stores and dissipates electromagnetic energy. The inverse Fourier transform of this impedance yields the wake function, describing the time-domain interaction between successive charged particles.

Accurate computation of the impedance and wake function is essential for:

  • Predicting collective beam dynamics, such as instabilities and bunch deformation

  • Estimating power deposition and guiding the mechanical design of beamline components

See also

  • 🎯 Beam dynamic simulations

To simulate collective beam dynamics and stability using impedance models, you can use the CERN-developed Python package Xsuite, an open-source macroparticle tracking simulation code suite.

See also

  • 💡 Beam Induced Heating

The dissipated power due to beam-induced heating can be estimated using BIHC, a tool within the Wakis ecosystem. It takes as input the impedance profile of a device and the beam spectrum, returning power loss predictions that are critical for vacuum design, cooling, and material compliance.

Analytical vs Numerical Approaches

Analytical methods offer elegant and accurate solutions for the beam-coupling impedance of simplified geometries. These are invaluable for physical insight and quick estimates in idealized cases such as:

However, real-world accelerator components like:

  • Beam Position Monitors (BPMs)

  • Injection kickers

  • RF-shielded bellows

  • Complex collimators or cavities

…often have no analytical solution. These must be addressed through full 3D numerical simulations of Maxwell’s equations Maxwell, 1865, using finite differences, finite elements, or finite integration.

Tip

  • 🧑‍🏫 Why Use Time-Domain Solvers?

Wakis employs a time-domain approach using the Finite Integration Technique (FIT), which offers key benefits:

  • Broadband response in a single simulation

  • Natural support for transient excitation (e.g., Gaussian bunch)

  • Efficient use of explicit solvers with GPU and MPI support

This makes Wakis well-suited for impedance characterization across a wide frequency range, complementing frequency-domain solvers like CST or HFSS.

For a broader overview of impedance modeling, see Metral et al., 2020.

3. Electromagnetic Formulation

⚡🧲 Maxwell’s Equations (Integral Form)

Wakis numerically solves Maxwell’s equations in their integral form, which is fundamental to the Finite Integration Technique (FIT). This approach preserves the physical laws in their conservative form and naturally fits the discretization on structured grids.

The time-domain integral form of Maxwell’s equations is:

\[\begin{split} \begin{align} \oint_{\partial A} \mathbf{E}\cdot \mathrm{d}\mathbf{s} &= -\iint_{A}\frac{\partial \mathbf{B}}{\partial t}\cdot \mathrm{d}\mathbf{A} \tag{1a}\\[6pt] \oint_{\partial A} \mathbf{H}\cdot \mathrm{d}\mathbf{s} &= \iint_{A}\left(\frac{\partial \mathbf{D}}{\partial t} + \mathbf{J}\right)\cdot \mathrm{d}\mathbf{A} \tag{1b}\\[6pt] \iint_{\partial V} \mathbf{B}\cdot \mathrm{d}\mathbf{A} &= 0 \tag{1c}\\[6pt] \iint_{\partial V} \mathbf{D}\cdot \mathrm{d}\mathbf{A} &= \iiint_{V}\rho\, \mathrm{d}V \tag{1d}\\[6pt] \mathbf{D} = \varepsilon \mathbf{E},\quad \mathbf{B} &= \mu \mathbf{H},\quad \mathbf{J} = \sigma \mathbf{E} + \rho\mathbf{v} \tag{1e} \end{align} \end{split}\]

These laws describe:

  • The evolution of electric \(\mathbf{E}\) and magnetic \(\mathbf{H}\) fields over time via their circulation around surfaces (Eqs. 1a–1b) and fluxes (magnetic flux density \(\mathbf{B}\), electric displacement field \(\mathbf{D}\))

  • The coupling to sources through current density \(\mathbf{J}\) and charge density \(\rho\) (Eqs. 1b, 1d),

  • The absence of magnetic monopoles (Eq. 1c),

  • And the constitutive relations of the materials (Eq. 1e), which relate the physical fields to the medium’s electromagnetic properties: permittivity \(\varepsilon\), permeability \(\mu\), and conductivity \(\sigma\), with \(\mathbf{v}\) denoting the velocity of moving charges.

In these equations, \(\varepsilon\), \(\mu\), \(\sigma\) can be considered tensors and frequency independent. To account for frequency or time dependency, the multiplication should be exchanged for a convolution \(\ast\).

🧱 Discretization with the Finite Integration Technique (FIT)

Wakis discretizes the integral form of Maxwell’s equations using the Finite Integration Technique (FIT) on a structured three-dimensional Cartesian grid.

\[ N_\text{cells} = N_x \times N_y \times N_z \]

This approach maps:

  • Line integrals → to grid edges

  • Surface integrals → to grid faces

  • Volume integrals → to grid cells

The resulting discretization yields the Maxwell Grid Equations (MGE), which evolve the fields on a staggered Yee-like lattice. Specifically:

  • \(\vec{E}\) and \(\vec{H}\) components are stored on edges

  • \(\vec{D}\) and \(\vec{B}\) components are defined on faces

  • Scalar quantities such as charge density reside at cell centers

This structure ensures that discrete curl, divergence, and gradient operators obey their continuous counterparts’ conservation properties, which is critical for numerical stability and accuracy. By adopting the Yee staggered grid formulation and initially charge-free conditions, the divergence equations (1c–1d) are satisfied implicitly by construction.

Maxwell Grid Equations (MGE)

Following the FIT numerical method, the continuous Maxwell equations are converted into discrete update rules for the electric and magnetic fields:

\[\begin{split} \begin{align} \mathbf{C}\mathbf{D}_s \, \mathbf{e} &= -\mathbf{D}_A \, \frac{\partial (\mathbf{M}_{\mu} \mathbf{h})}{\partial t} \tag{2a} \\[6pt] \widetilde{\mathbf{C}}\widetilde{\mathbf{D}}_s \, \mathbf{h} &= \widetilde{\mathbf{D}}_A \left( \frac{\partial (\mathbf{M}_{\varepsilon} \mathbf{e})}{\partial t} + \mathbf{M}_{\sigma} \mathbf{e} + \mathbf{j}_{\text{src}} \right) \tag{2b} \\[6pt] \end{align} \end{split}\]

Where:

  • \(\mathbf{C}\) is the discrete curl matrix

  • \(\mathbf{C}^T\) is its transpose (used for magnetic curl)

  • \(\mathbf{D}_s\), \(\widetilde{\mathbf{D}}_s\), \(\mathbf{D}_A\), and \(\widetilde{\mathbf{D}}_A\) are diagonal matrices representing cell edge lengths and face areas in the primal and dual~(\(\sim\)) grids.

  • The electromagnetic fields \(\mathbf{e}, \mathbf{h}, \mathbf{j}\) are stored in memory as 1D vectors of length \(\{3N_\text{cells}\}\) stored in lexicographic order, encapsulated in the Field class.

The numerical method is implemented in SolverFIT3D class. The curl matrices are \(\{3N_{cells}\times3N_{cells}\}\) sparse matrices with bands of +1 and -1. They are implemented efficiently in Wakis using scipy.sparse CSR format. The diagonal matrices also benefit from the scipy.sparse.diags object for optimized storage.

The electromagnetic fields are stored in a numpy-based Field object in Wakis:

  • Supports .toarray() and .fromarray() for optimized modification during the time-stepping

  • .from_matrix(), .to_matrix() to go from 1d to 3d matrix,automatically reshaped to the simulation grid dimensions.

  • Interoperates with CuPy and MPI through magic methods and flags.

  • .inspect() and other handy plotting methods

  • Custom magic methods for multiplication, addition, division: __div__, __add__, __mul__

  • __getitem__, __setitem__ Getters and setters to access the 3D coordinates on-the-fly: apply intial conditions, sources, save states…

Tip

Thank’s to the Field class, fields are stored in memory-continuous arrays for optimized performance, but can be accessed as a 3D matrix.

Some examples on how to access and operate Fields:

# modify field slice in z-direction
solver.E[100, 20:30, :, 'z'] = 0.

# access cell value of the 123456th cell in lexico-graphic order
solver.H[123456]

# sum or multiply two field objects, keeping the Field
E_tot = solver_1.E + solver_2.E

# Calculate the energy (T_00)
T_00 = 0.5*(solver.E.get_abs()**2 + solver.H.get_abs()**2)

# Inspect the 3 components of the field in one line
solver.J.inspect()

This formulation enables stable, explicit time stepping and modular extensions to lossy media, materials, and sources.

Material tensors and grid information

Wakis distinguishes between primal and dual grid geometries as part of its Finite Integration Technique (FIT) formulation. The grid operations are implemented in Wakis’ GridFIT3D class. Each quantity is mapped to a geometric entity and stored as a sparse diagonal matrix to enable fast, memory-efficient computations:

Quantity

Description

\(\mathbf{M}_\varepsilon^{-1}\)

Diagonal matrix of inverse permittivities

\(\mathbf{M}_\mu^{-1}\)

Diagonal matrix of inverse permeabilities

\(\mathbf{M}_\sigma\)

Diagonal matrix of electrical conductivities

\(\mathbf{D}_s\), \(\widetilde{\mathbf{D}}_s\), \(\mathbf{D}_A\), and \(\widetilde{\mathbf{D}}_A\)

Edge lengths and face areas (primal/dual)

To support anisotropic materials and imported geometries, Wakis stores the raw material data in structured Field objects — similar to 3D tensors — where values can be specified independently along the x, y, and z directions for each cell.

Before time-stepping, these directional fields are assembled into the corresponding sparse diagonal matrices using scipy.sparse.diags. This preserves the locality of FIT updates while enabling efficient CPU and GPU execution.

Wakis also supports spatially varying media, embedded materials from CAD imports, and subpixel smoothing, ensuring accurate representation of complex geometries and composite structures.

🕒 Time-Stepping Routine

Wakis uses the Leapfrog scheme, a second-order accurate and explicit time integrator. This method updates the magnetic and electric fields in a staggered fashion:

\[\begin{split} \begin{align} \mathbf{h}^{n+1} &= \mathbf{h}^n - \Delta t \, \widetilde{\mathbf{D}}_s \, \mathbf{M}_\mu^{-1} \, \mathbf{D}_A^{-1} \, \mathbf{C} \, \mathbf{e}^{n+0.5} \tag{3a} \\[6pt] \mathbf{e}^{n+1.5} &= \mathbf{e}^{n+0.5} + \Delta t \, \mathbf{D}_s \, \widetilde{\mathbf{M}}_\varepsilon^{-1} \, \widetilde{\mathbf{D}}_A^{-1} \, \widetilde{\mathbf{C}} \, \mathbf{h}^n - \widetilde{\mathbf{M}}_\varepsilon^{-1} \, \mathbf{j}_{\text{src}}^n - \widetilde{\mathbf{M}}_\varepsilon^{-1} \, \widetilde{\mathbf{M}}_\sigma \, \mathbf{e}^{n+0.5} \tag{3b} \end{align} \end{split}\]

Where:

  • \(\mathbf{C}\) and \(\widetilde{\mathbf{C}}\) are discrete curl matrices on the primal and dual grids

  • \(\mathbf{D}_s\), \(\mathbf{D}_A\) are diagonal matrices containing metric coefficients (edge lengths and face areas)

  • Material properties enter through \(\mathbf{M}_\mu\), \(\mathbf{M}_\varepsilon\), and \(\mathbf{M}_\sigma\)

The timestep \(\Delta t\) is constrained by:

  • The Courant–Friedrichs–Lewy (CFL) condition

  • Material relaxation times (for dispersive or lossy media)

Most matrix operations are precomputed and cached, enabling large-scale simulations with modest memory usage (~20M cells in <8 GB RAM/GPU).

🔌 Sources and Initial Conditions

Wakis allows arbitrary initial conditions on \(\vec{E}\), \(\vec{H}\), and \(\vec{J}\). Sources can be defined in multiple ways:

  • User-defined time-dependent callbacks placed after each step

  • Predefined source types: Gaussian beams, dipoles, plane waves, laser pulses, available in sources.py

Tip

⚙️ Source callbacks

A source callback can be easily created as:

  • A function like def update(solver, time) placed in a for loop after each step solver.one_step().

  • a class with the method Source.update(solver, time), passed to the solver.emsolve() routine.

Example: Gaussian Beam Current \(J_z\)

A rigid Gaussian bunch current is modeled as a line distribution:

\[ \mathbf{J}_z(x_{\text{src}}, y_{\text{src}}, \vec{z}) = \frac{q \beta c}{\sqrt{2\pi} \sigma_z} \, \exp\left( -\frac{(\vec{s} - s_0)^2}{2\sigma_z^2} \right) \]

with:

  • \(\vec{s} = \vec{z} - \beta c t\): beam-frame coordinate

  • \(s_0 = z_{\min} - \beta c t_{\text{inj}}\): center of bunch

  • \(q\) the charge in \(\text{nC}\)

  • \(\sigma_z\) the bunch length in \(\text{m}\)

This supports both ultra-relativistic (\(\beta \approx 1\)) and low-beta scenarios.

🧊🔚 Boundary Conditions

Wakis supports several boundary condition (BC) types:

  • PEC (Perfect Electric Conductor): enforces \(\vec{E}_{\parallel} = 0\)

  • PMC (Perfect Magnetic Conductor): enforces \(\vec{H}_{\parallel} = 0\)

  • Periodic BCs: implemented with synchronized ghost cells

  • PML (Perfectly Matched Layers): often referred to open or absorbing BC, are made of matched (\(\varepsilon = \mu\)) layers using graded conductivity \(\sigma\) profiles for reflection-free truncation of the computational domain

PML implementation

PMLs follow the formulation by Berenger [1994] and are ramped using smooth profiles [Oskooi et al., 2008] to reach adiabatic reflection.

Warning

The PML implementation description is under development!

📥🗿 Geometry Importing & Embedded Boundaries

Wakis integrates with PyVista to import CAD geometries in .STL, .STEP, or .OBJ formats. The mesh is overlaid onto the simulation domain and mapped onto the Cartesian grid using:

  • pyvista’s surface collision algorithm, based on VTK optimized ray-tracing, allows to detect where the input geometry intersects the primal and dual grids.

  • Assignment of material properties (\(\varepsilon_r\), \(\mu_r\), \(\sigma\)) in \(x\), \(y\), and \(z\) to the intersected cells using a first-order subpixel smoothing, inspired by the open-source solver MEEP (MIT).

Future versions aim to include a more advanced meshing algorithm for improved fidelity near corners and edges.

🚀 GPU and MPI Parallelization

Wakis supports heterogeneous architecture computing thanks to open-source packages like:

  • GPU acceleration using CuPy and cupyx.scipy.sparse

  • Drop-in replacement of NumPy/SciPy operations when use_gpu=True

  • MPI parallelization using mpi4py

  • Efficient longitudinal domain decomposition with ghost-cell synchronization

  • Seamless integration with multi-GPU setups using both cupy and mpi4py memory passing protocols

Note

👩‍💻 Developer Notes about Wakis

  • Fully open-source and available on GitHub

  • Packaged on PyPI

  • Documented with Sphinx and hosted on ReadTheDocs: https://wakis.readthedocs.io/

  • Includes CI/CD, with end-to-end tests running nightly on GitHub actions, tagged versioned releases, and numerous ready-to-run examples in both Python scripts and notebooks, inluding a dedicated playground repository.

3. Wake Potential and Impedance calculation

Wakis computes beam coupling impedance from time-domain electromagnetic field simulations by evaluating the wakefields generated by a moving charged particle (or bunch) as it traverses an accelerator structure.

📚 Physical Definition: Wake function and Impedance

The longitudinal wake function \(w(\vec{r_s}, \vec{r_t}, s)\) of an accelerator component can be defined as a Green function in the time domain (i.e., the component electromagnetic response to a pulse excitation):

\[ w(\vec{r_s}, \vec{r_t}, s) = \frac{1}{q_s q_t} \int_{-\infty}^{\infty} \vec{F}_{\text{Lorentz}} \ d\vec{z} = \frac{1}{q_s q_t} \int_{-\infty}^{\infty} E_z(z, t) + \beta c \, \vec{e}_z \times \vec{B}(z, t) \ d\vec{z} \]

With:

  • \(\vec{r_s} = (x_s, y_s, z_s)\) be the position of a source particle -or pulse excitation.

  • \(\vec{r_t} = (x_t, y_t, z_t)\) the test particle position -or integration point.

  • \(s = z_{min} - \beta c t\) the catch-up distance, with \(s_{max}\) beign the desired wakefield’s length.

The wake function is the input to beam-dynamics simulations. Its Fourier transform yields the longitudinal impedance, in frequency domain:

\[ Z_{\parallel}(\omega) = \int_{-\infty}^{\infty} w_{\parallel}(s) \, e^{-i \omega s / c} \frac{ds}{c} \]

In practice, since a beam is a distributed source, the wake function is not directly accessible through wakefield simulations, where the excitation is a gaussian-shaped current. Instead, in wakefield simulations we compute the wake potential \(W(s)\) generated by the full bunch distribution.

📈 Wake Potential from 3D electromganetic simulations

The wake potential, expressed in \(\text{V/pC}\), is calculated by integrating the electric and magnetic fields seen by a test particle as it follows behind the source:

\[ W(s) = \frac{1}{q_s} \int_{-\infty}^{\infty} \left[ E_z(z, t) + \beta c \, \vec{e}_z \times \vec{B}(z, t) \right]_{t = (s + z)/c} \, dz \]

For ultra-relativistic beams, the transverse component vanishes, and the expression simplifies to:

\[ W_\parallel(s) = \frac{1}{q_s} \int_{-\infty}^{\infty} E_z(z, t = (s + z)/c) \, dz \]

The transverse wake potential is recovered via the Panofsky-Wenzel theorem:

\[ W_{\perp,\alpha}(s) = \frac{\partial}{\partial \alpha} \int_{-\infty}^{s} W_\parallel(s') \, ds', \quad \alpha = x, y \]

Wakis implements this gradient using second-order finite differences.

Transverse Decomposition

Wakis supports transverse wake analysis:

\[ W_{\perp,x}(x, y, s) = W_C(s) + W_D(s) \Delta x_s + W_Q(s) \Delta x_t + \mathcal{O}(x^2) \]

  • \(W_D\): dipolar wake, linear in source offset

  • \(W_Q\): quadrupolar wake, linear in test offset

  • \(W_C\): coherent term, for asymmetric geometries

These are extracted by sampling field responses at multiple \((x_s, y_s, x_t, y_t)\) combinations, either by displacing the beam source (\(x_s, y_s\)) or the integration path (\(x_t, y_t\)).

🔁 From Wake to Impedance

Given the bunch profile \(\lambda(s)\) and the wake potential \(W(s)\), the beam coupling impedance is computed in Fourier space ia a deconvolution:

  • Longitudinal impedance in \(\Omega\):

\[ Z_\parallel(\omega) = \beta c \cdot \frac{\mathcal{F}[W_\parallel(s)]}{\mathcal{F}[\lambda(s)]} \]

  • Transverse impedance in \(\Omega/m\):

\[ Z_\perp(\omega) = -i \beta c \cdot \frac{\mathcal{F}[W_\perp(s)]}{\mathcal{F}[\lambda(s)]} \]

where \(\mathcal{F}\) denotes the Fourier transform. Wakis uses numpy.fft with a Hanning window and zero-padding for smooth frequency analysis.

Modularity

🧪 Open-Source Compatibility

The wake and impedance analysis is performed within the WakeSolver class in Wakesolver.py. Even if this module is tailored for Wakis, it is completely modular and can be used with other EM solvers besides Wakis’ SolverFIT3D output:

  • It has been tested with WarpX EM fields. WarpX is a powerful open-source PIC solver for full 3D EM fields

  • It has been benchmarked with CST Wakefield solver, using both EM field ouput (time-domain field monitors) and calculated wake potential and impedance. > See the full wake analysis benchmark in IPAC’23 proceedings, E. de la Fuente.

  • Interoperability: Wakis can read field maps in HDF5, CSV, or NumPy format

  • Subvolume extraction and interpolation are supported for field-based post-processing