Wakefields

Resistive Wall Wakefield Classes

pmd_beamphysics.wakefields.ResistiveWallWakefield dataclass

ResistiveWallWakefield(radius, conductivity, relaxation_time, geometry=ROUND)

Bases: ResistiveWallWakefieldBase, ImpedanceWakefield

Accurate impedance-based resistive wall wakefield model.

This class computes the short-range resistive wall wakefield by numerically evaluating the full longitudinal impedance Z(k) and transforming to the wakefield W(z) via FFT. This approach captures all physical effects including the anomalous skin effect at high frequencies.

Physics Background

When a relativistic charged particle travels through a conducting beam pipe, electromagnetic fields penetrate into the conductor due to the finite conductivity. The skin depth δ decreases with frequency:

\[\delta = \sqrt{\frac{2}{\omega \mu_0 \sigma}}\]

At very high frequencies (relevant for short-range wakes), the electron mean free path becomes comparable to the skin depth, and the simple Ohmic model breaks down. The Drude model accounts for this through a frequency- dependent conductivity with relaxation time τ:

\[\sigma(k) = \frac{\sigma_0}{1 - ikc\tau}\]

The longitudinal impedance per unit length for a round pipe is:

\[Z(k) = \frac{Z_0}{2\pi a} \cdot \frac{\zeta(k)}{1 - i k a \zeta(k) / 2}\]

where ζ(k) is the surface impedance and a is the pipe radius.

The wakefield is obtained via cosine transform (for z ≤ 0):

\[W(z) = \frac{2c}{\pi} \int_0^{\infty} \text{Re}[Z(k)] \cos(kz) \, dk\]

Characteristic Scales

• s₀: Characteristic length scale:

$\(s_0 = \left( \frac{2a^2}{Z_0 \sigma_0} \right)^{1/3}\)$

For copper with a = 2.5 mm, s₀ ≈ 8 µm.

• W₀: Characteristic wake amplitude at z = 0:
• Round: \(W_0 = c Z_0 / (\pi a^2)\)
• Flat: \(W_0 = c Z_0 \pi / (16 a^2)\)

Geometry

• Round (circular pipe): Closed-form impedance expression.
• Flat (parallel plates): Requires numerical integration over transverse modes; more computationally expensive.

Parameters:

Name	Type	Description	Default
radius	float	Radius of the beam pipe [m]. For flat geometry, this is the half-gap.	required
conductivity	float	Electrical DC conductivity of the wall material [S/m].	required
relaxation_time	float	Drude-model relaxation time τ of the conductor [s]. Typical values: Cu ≈ 27 fs, Al ≈ 8 fs.	required
geometry	str	Geometry of the beam pipe: 'round' or 'flat'. Default is 'round'.	ROUND

Attributes:

Name	Type	Description
s0	float	Characteristic length scale [m]
W0	float	Characteristic wake amplitude [V/C/m]

Notes

• This model is ~10-20× slower than ResistiveWallPseudomode but more accurate, especially for the first oscillation peak.
• The FFT-based convolution uses O(N log N) complexity for density convolution, compared to O(N) for the pseudomode model.
• For flat geometry, impedance calculation requires numerical integration and is significantly slower than round geometry.

See Also

ResistiveWallPseudomode : Fast approximate model using damped sinusoid

References

.. [1] K. Bane and G. Stupakov, "Resistive wall wakefield in the LCLS undulator beam pipe," SLAC-PUB-10707 (2004). https://www.slac.stanford.edu/cgi-wrap/getdoc/slac-pub-10707.pdf

.. [2] A. Chao, "Physics of Collective Beam Instabilities in High Energy Accelerators," Wiley, 1993, Chapter 2.

Examples:

::

wake = ResistiveWallWakefield.from_material(
    "copper-slac-pub-10707", radius=2.5e-3, geometry="round"
)
wake.wake(-10e-6)  # Wake at 10 µm behind source
wake.impedance(1e5)  # Impedance at k = 100/mm

Functions

pmd_beamphysics.wakefields.ResistiveWallWakefield.impedance

impedance(k)

Evaluate the impedance at wavenumber k.

Parameters:

Name	Type	Description	Default
k	float or ndarray	Wavenumber [1/m].	required

Returns:

Name	Type	Description
Z	complex or ndarray	Impedance [Ohm/m].

Source code in pmd_beamphysics/wakefields/resistive_wall/impedance.py

def impedance(self, k):
    """
    Evaluate the impedance at wavenumber k.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    if self.geometry == Geometry.ROUND:
        return longitudinal_impedance_round(
            k, self.radius, self.conductivity, self._ctau
        )
    else:  # FLAT
        return longitudinal_impedance_flat(
            k, self.radius, self.conductivity, self._ctau
        )

pmd_beamphysics.wakefields.ResistiveWallWakefield.plot

plot(zmax=None, zmin=0, n=200, normalized=False, ax=None)

Plot the resistive wall wakefield W(z).

Parameters:

Name	Type	Description	Default
zmax	float	Maximum trailing distance [m]. Defaults to 100 * s0.	None
zmin	float	Minimum trailing distance [m]. Default is 0.	0
n	int	Number of points. Default is 200.	200
normalized	bool	If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀. Default is False.	False
ax	Axes	Axes to plot on. If None, creates a new figure.	None

Source code in pmd_beamphysics/wakefields/resistive_wall/impedance.py

def plot(self, zmax=None, zmin=0, n=200, normalized=False, ax=None):
    """
    Plot the resistive wall wakefield W(z).

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. Defaults to 100 * s0.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    normalized : bool, optional
        If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀.
        Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 20 * self.s0

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)

    if ax is None:
        _, ax = plt.subplots()

    if normalized:
        ax.plot(zlist / self.s0, Wz / self.W0)
        ax.set_xlabel(r"$|z|/s_0$")
        ax.set_ylabel(r"$W(z)/W_0$")
    else:
        ax.plot(zlist * 1e6, Wz * 1e-12)
        ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
        ax.set_ylabel(r"$W(z)$ (V/pC/m)")

pmd_beamphysics.wakefields.ResistiveWallWakefield.plot_impedance

plot_impedance(kmax=None, kmin=0, n=500, ax=None)

Plot the resistive wall impedance Z(k).

Parameters:

Name	Type	Description	Default
kmax	float	Maximum wavenumber [1/m]. Defaults to 10/s0.	None
kmin	float	Minimum wavenumber [1/m]. Default is 0.	0
n	int	Number of points. Default is 500.	500
ax	Axes	Axes to plot on. If None, creates a new figure.	None

Source code in pmd_beamphysics/wakefields/resistive_wall/impedance.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 500, ax=None
):
    """
    Plot the resistive wall impedance Z(k).

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. Defaults to 10/s0.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 500.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 10 / self.s0
    super().plot_impedance(kmax=kmax, kmin=kmin, n=n, ax=ax)

pmd_beamphysics.wakefields.ResistiveWallPseudomode dataclass

ResistiveWallPseudomode(radius, conductivity, relaxation_time, geometry=ROUND)

Bases: ResistiveWallWakefieldBase, PseudomodeWakefield

Fast pseudomode-based resistive wall wakefield model.

This class models the short-range resistive wall wakefield using a damped sinusoidal "pseudomode" approximation. When a relativistic charged particle travels through a conducting beam pipe, it induces image currents in the pipe walls. Due to the finite conductivity of the wall material, these currents penetrate into the conductor and dissipate energy, creating a wakefield that acts on trailing particles.

Physics Background

The resistive wall impedance arises from the skin effect in the conducting walls. At high frequencies (short distances), the AC conductivity of metals deviates from DC behavior due to the Drude relaxation time τ:

\[\sigma(\omega) = \frac{\sigma_0}{1 - i\omega\tau}\]

This frequency dependence causes the wakefield to oscillate and decay exponentially behind the source particle. The wakefield can be well approximated by a single damped sinusoid (pseudomode):

\[W(z) = A \cdot e^{d \cdot z} \cdot \sin(k_r z + \phi) \quad (z \le 0)\]

where the parameters \(k_r\) (oscillation wavenumber) and \(Q_r\) (quality factor, related to decay rate \(d = k_r / 2Q_r\)) depend on the dimensionless relaxation parameter \(\Gamma = c\tau / s_0\).

Characteristic Scales

• s₀: Characteristic length scale where the wake transitions from the \(1/\sqrt{z}\) DC behavior to oscillatory AC behavior:

$\(s_0 = \left( \frac{2a^2}{Z_0 \sigma_0} \right)^{1/3}\)$

• Γ: Dimensionless relaxation parameter \(\Gamma = c\tau / s_0\). For copper, Γ ≈ 0.8; for aluminum, Γ ≈ 0.2.

• W₀: Characteristic wake amplitude at z = 0:

• Round: \(W_0 = c Z_0 / (\pi a^2)\)
• Flat: \(W_0 = c Z_0 \pi / (16 a^2)\)

Algorithm

The pseudomode form enables an O(N) algorithm for computing particle kicks, compared to O(N²) or O(N log N) for general wakefields. This makes it suitable for multi-pass tracking simulations.

Parameters:

Name	Type	Description	Default
radius	float	Radius of the beam pipe [m]. For flat geometry, this is half the gap.	required
conductivity	float	Electrical DC conductivity of the wall material [S/m].	required
relaxation_time	float	Drude-model relaxation time τ of the conductor [s]. Typical values: Cu ≈ 27 fs, Al ≈ 8 fs.	required
geometry	str	Geometry of the beam pipe: 'round' or 'flat'. Default is 'round'.	ROUND

Attributes:

Name	Type	Description
s0	float	Characteristic length scale [m]
Gamma	float	Dimensionless relaxation parameter Γ = cτ/s₀
kr	float	Resonant wavenumber [1/m]
Qr	float	Quality factor of the pseudomode

Notes

• The pseudomode approximation has ~10-20% error compared to the full impedance model, primarily in the first oscillation peak.
• The polynomial fits for k_r and Q_r are valid for Γ ≲ 3. A warning is issued for larger values.
• Materials with known conductivity and τ values are available via from_material().

See Also

ResistiveWallWakefield : Accurate impedance-based model (slower)

References

.. [1] K. Bane and G. Stupakov, "Resistive wall wakefield in the LCLS undulator beam pipe," SLAC-PUB-10707 (2004). https://www.slac.stanford.edu/cgi-wrap/getdoc/slac-pub-10707.pdf

.. [2] D. Sagan, "Bmad Manual," Section 24.6: Short-Range Wakefields. https://www.classe.cornell.edu/bmad/manual.html

.. [3] A. Chao, "Physics of Collective Beam Instabilities in High Energy Accelerators," Wiley, 1993, Chapter 2.

Examples:

::

wake = ResistiveWallPseudomode.from_material(
    "copper-slac-pub-10707", radius=2.5e-3, geometry="round"
)
wake.wake(-10e-6)  # Wake at 10 µm behind source

Attributes

pmd_beamphysics.wakefields.ResistiveWallPseudomode.Gamma property

Gamma

Dimensionless relaxation time Γ = c * τ / s₀.

pmd_beamphysics.wakefields.ResistiveWallPseudomode.Qr property

Qr

Dimensionless quality factor Q_r of the wakefield pseudomode.

pmd_beamphysics.wakefields.ResistiveWallPseudomode.kr property

kr

Real-valued wave number k_r of the wakefield pseudomode [1/m].

Functions

pmd_beamphysics.wakefields.ResistiveWallPseudomode.plot

plot(zmax=None, zmin=0, n=200, normalized=False, ax=None)

Plot the resistive wall wakefield W(z).

Parameters:

Name	Type	Description	Default
zmax	float	Maximum trailing distance [m]. Defaults to 10 decay lengths.	None
zmin	float	Minimum trailing distance [m]. Default is 0.	0
n	int	Number of points. Default is 200.	200
normalized	bool	If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀. Default is False.	False
ax	Axes	Axes to plot on. If None, creates a new figure.	None

Source code in pmd_beamphysics/wakefields/resistive_wall/pseudomode.py

def plot(self, zmax=None, zmin=0, n=200, normalized=False, ax=None):
    """
    Plot the resistive wall wakefield W(z).

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. Defaults to 10 decay lengths.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    normalized : bool, optional
        If True, plot dimensionless Ŵ(z/s₀) = W(z)/W₀ vs z/s₀.
        Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 1 / (self.kr / (2 * self.Qr)) * 10

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)

    if ax is None:
        _, ax = plt.subplots()

    if normalized:
        ax.plot(zlist / self.s0, Wz / self.W0)
        ax.set_xlabel(r"$|z|/s_0$")
        ax.set_ylabel(r"$W(z)/W_0$")
    else:
        ax.plot(zlist * 1e6, Wz * 1e-12)
        ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
        ax.set_ylabel(r"$W(z)$ (V/pC/m)")

pmd_beamphysics.wakefields.ResistiveWallPseudomode.plot_impedance

plot_impedance(kmax=None, kmin=0, n=500, ax=None)

Plot the resistive wall impedance Z(k).

Parameters:

Name	Type	Description	Default
kmax	float	Maximum wavenumber [1/m]. Defaults to 10/s0.	None
kmin	float	Minimum wavenumber [1/m]. Default is 0.	0
n	int	Number of points. Default is 500.	500
ax	Axes	Axes to plot on. If None, creates a new figure.	None

Source code in pmd_beamphysics/wakefields/resistive_wall/pseudomode.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 500, ax=None
):
    """
    Plot the resistive wall impedance Z(k).

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. Defaults to 10/s0.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 500.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 10 / self.s0
    super().plot_impedance(kmax=kmax, kmin=kmin, n=n, ax=ax)

pmd_beamphysics.wakefields.ResistiveWallPseudomode.to_bmad

to_bmad(file=None, z_max=100, amp_scale=1, scale_with_length=True, z_scale=1)

Export wakefield in Bmad format.

Parameters:

Name	Type	Description	Default
file	str	Output file path. If None, returns string only.	None
z_max	float	Trailing z distance for Bmad.	100
amp_scale	float	Amplitude scaling factor.	1
scale_with_length	bool	Whether to scale with length.	True
z_scale	float	Z scaling factor.	1

Returns:

Type	Description
str	Bmad-formatted wakefield string.

Source code in pmd_beamphysics/wakefields/resistive_wall/pseudomode.py

    def to_bmad(
        self, file=None, z_max=100, amp_scale=1, scale_with_length=True, z_scale=1
    ):
        """
        Export wakefield in Bmad format.

        Parameters
        ----------
        file : str, optional
            Output file path. If None, returns string only.
        z_max : float
            Trailing z distance for Bmad.
        amp_scale : float
            Amplitude scaling factor.
        scale_with_length : bool
            Whether to scale with length.
        z_scale : float
            Z scaling factor.

        Returns
        -------
        str
            Bmad-formatted wakefield string.
        """
        s = f"""! AC Resistive wall wakefield
! Adapted from SLAC-PUB-10707
!    Material        : {self.material_from_properties()}
!    Conductivity    : {self.conductivity} S/m
!    Relaxation time : {self.relaxation_time} s
!    Geometry        : {self.geometry.value}
"""

        if self.geometry == Geometry.ROUND:
            s += f"!    Radius          : {self.radius} m\n"
        elif self.geometry == Geometry.FLAT:
            s += f"!    full gap        : {2*self.radius} m\n"

        s += f"!    s₀              : {self.s0}  m\n"
        s += f"!    Γ               : {self.Gamma} \n"
        s += "! sr_wake =  \n"

        s += f"{{{z_scale=}, {amp_scale=}, {scale_with_length=}, {z_max=},\n"
        s += PseudomodeWakefield.to_bmad(self) + "}\n"

        if file is not None:
            with open(file, "w") as f:
                f.write(s)

        return s

Base Classes

pmd_beamphysics.wakefields.WakefieldBase

Bases: ABC

Abstract base class for longitudinal wakefield models.

All wakefield implementations should inherit from this class and implement the required abstract methods.

Sign Convention

The wakefield W(z) is defined for a trailing particle at position z relative to the source particle at z=0:

• z < 0: Behind the source (trailing particle feels the wake)
• z > 0: Ahead of the source (causality requires W=0)

The wake is positive for an energy-losing interaction (the trailing particle loses energy).

Methods:

Name	Description
wake	Evaluate the wakefield at position z [V/C/m]
impedance	Evaluate the impedance at wavenumber k [Ohm/m]
convolve_density	Compute integrated wake from a charge density distribution [V/m]
particle_kicks	Compute momentum kicks for a distribution of particles [eV/m]
plot	Plot the wakefield over a range of z values
plot_impedance	Plot the impedance over a range of wavenumbers

See Also

ParticleGroup.apply_wakefield : Apply wakefield kicks to a ParticleGroup

Functions

pmd_beamphysics.wakefields.WakefieldBase.convolve_density

convolve_density(density, dz, offset=0, include_self_kick=True, plot=False, ax=None)

Compute integrated wakefield from a charge density distribution.

This computes the causal wake potential:

\[V(z_i) = \sum_{j>i} Q_j \cdot W(z_j - z_i) + \frac{1}{2} Q_i \cdot W(0)\]

where only particles ahead (larger z index) contribute to the wake felt by each particle. This is mathematically a correlation, not a convolution, and is computed efficiently using FFT.

Parameters:

Name	Type	Description	Default
density	ndarray	Charge density array [C/m]. Index 0 is the tail (back), index n-1 is the head (front).	required
dz	float	Grid spacing [m].	required
offset	float	Offset for the z coordinate [m]. Default is 0.	0
include_self_kick	bool	Whether to include the half self-kick term. Default is True.	True
plot	bool	If True, plot the density profile and wake potential. Default is False.	False
ax	Axes	Axes to plot on. If None and plot=True, creates a new figure.	None

Returns:

Name	Type	Description
integrated_wake	ndarray	Integrated longitudinal wakefield [V].

Source code in pmd_beamphysics/wakefields/base.py

def convolve_density(
    self,
    density: np.ndarray,
    dz: float,
    offset: float = 0,
    include_self_kick: bool = True,
    plot: bool = False,
    ax=None,
) -> np.ndarray:
    """
    Compute integrated wakefield from a charge density distribution.

    This computes the causal wake potential:

    $$V(z_i) = \\sum_{j>i} Q_j \\cdot W(z_j - z_i) + \\frac{1}{2} Q_i \\cdot W(0)$$

    where only particles ahead (larger z index) contribute to the wake
    felt by each particle. This is mathematically a correlation, not
    a convolution, and is computed efficiently using FFT.

    Parameters
    ----------
    density : np.ndarray
        Charge density array [C/m]. Index 0 is the tail (back),
        index n-1 is the head (front).
    dz : float
        Grid spacing [m].
    offset : float, optional
        Offset for the z coordinate [m]. Default is 0.
    include_self_kick : bool, optional
        Whether to include the half self-kick term. Default is True.
    plot : bool, optional
        If True, plot the density profile and wake potential. Default is False.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None and plot=True, creates a new figure.

    Returns
    -------
    integrated_wake : np.ndarray
        Integrated longitudinal wakefield [V].
    """
    from scipy.signal import correlate

    n = len(density)
    charge = density * dz  # Convert to charge per bin [C]

    # Build causal wake array: W[k] = wake at lag k*dz
    # W[0] = W(0), W[1] = W(-dz), W[2] = W(-2*dz), ...
    z_wake = -np.arange(n) * dz + offset
    W_causal = self.wake(z_wake)

    # Correlation: result[i] = sum_k charge[i+k] * W[k]
    # This gives the wake from particles ahead (larger indices)
    corr_result = correlate(charge, W_causal, mode="full")

    # Extract the portion corresponding to non-negative lags
    # correlate output has length 2*n-1, with zero lag at index n-1
    integrated_wake = corr_result[n - 1 : 2 * n - 1]

    # Correlation includes full W[0] contribution from self
    # Remove it, then add half if self-kick is requested
    W0 = W_causal[0]
    integrated_wake = integrated_wake - charge * W0
    if include_self_kick:
        integrated_wake = integrated_wake + 0.5 * W0 * charge

    if plot:
        self._plot_convolve_density(density, dz, integrated_wake, ax=ax)

    return integrated_wake

pmd_beamphysics.wakefields.WakefieldBase.impedance abstractmethod

impedance(k)

Evaluate the impedance at wavenumber k.

Parameters:

Name	Type	Description	Default
k	float or ndarray	Wavenumber [1/m].	required

Returns:

Name	Type	Description
Z	complex or ndarray	Impedance [Ohm/m].

Source code in pmd_beamphysics/wakefields/base.py

@abstractmethod
def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Evaluate the impedance at wavenumber k.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    pass

pmd_beamphysics.wakefields.WakefieldBase.particle_kicks

particle_kicks(z, weight, include_self_kick=True)

Compute wakefield-induced energy kicks per unit length.

This is a default O(N²) implementation. Subclasses may override with more efficient algorithms.

Parameters:

Name	Type	Description	Default
z	ndarray	Particle positions [m].	required
weight	ndarray	Particle charges [C].	required
include_self_kick	bool	Whether to include the self-kick term. Default is True.	True

Returns:

Name	Type	Description
kicks	ndarray	Energy kick per unit length at each particle [eV/m].

Source code in pmd_beamphysics/wakefields/base.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
) -> np.ndarray:
    """
    Compute wakefield-induced energy kicks per unit length.

    This is a default O(N²) implementation. Subclasses may override
    with more efficient algorithms.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        Whether to include the self-kick term. Default is True.

    Returns
    -------
    kicks : np.ndarray
        Energy kick per unit length at each particle [eV/m].
    """
    z = np.asarray(z)
    weight = np.asarray(weight)
    n = len(z)
    kicks = np.zeros(n)

    # O(N²) algorithm: sum contribution from all particles
    for i in range(n):
        for j in range(n):
            if i == j:
                if include_self_kick:
                    # Self-kick uses W(0⁻) / 2
                    kicks[i] -= 0.5 * weight[j] * self._self_kick_value()
            else:
                dz = z[i] - z[j]
                if dz < 0:  # j is ahead of i, so i feels wake from j
                    kicks[i] -= weight[j] * self.wake(dz)

    return kicks

pmd_beamphysics.wakefields.WakefieldBase.plot

plot(zmax=None, zmin=0, n=200, ax=None)

Plot the wakefield over a range of z values.

Parameters:

Name	Type	Description	Default
zmax	float	Maximum trailing distance [m]. If None, uses a sensible default.	None
zmin	float	Minimum trailing distance [m]. Default is 0.	0
n	int	Number of points. Default is 200.	200
ax	Axes	Axes to plot on. If None, creates a new figure.	None

Source code in pmd_beamphysics/wakefields/base.py

def plot(self, zmax: float = None, zmin: float = 0, n: int = 200, ax=None):
    """
    Plot the wakefield over a range of z values.

    Parameters
    ----------
    zmax : float, optional
        Maximum trailing distance [m]. If None, uses a sensible default.
    zmin : float, optional
        Minimum trailing distance [m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if zmax is None:
        zmax = 1e-3  # 1 mm default

    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)  # Negative z for trailing particles

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(zlist * 1e6, Wz * 1e-12)
    ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
    ax.set_ylabel(r"$W(z)$ (V/pC/m)")

pmd_beamphysics.wakefields.WakefieldBase.plot_impedance

plot_impedance(kmax=None, kmin=0, n=200, ax=None)

Plot the impedance over a range of wavenumbers.

Parameters:

Name	Type	Description	Default
kmax	float	Maximum wavenumber [1/m]. If None, uses a sensible default.	None
kmin	float	Minimum wavenumber [1/m]. Default is 0.	0
n	int	Number of points. Default is 200.	200
ax	Axes	Axes to plot on. If None, creates a new figure.	None

Source code in pmd_beamphysics/wakefields/base.py

def plot_impedance(
    self, kmax: float = None, kmin: float = 0, n: int = 200, ax=None
):
    """
    Plot the impedance over a range of wavenumbers.

    Parameters
    ----------
    kmax : float, optional
        Maximum wavenumber [1/m]. If None, uses a sensible default.
    kmin : float, optional
        Minimum wavenumber [1/m]. Default is 0.
    n : int, optional
        Number of points. Default is 200.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates a new figure.
    """
    if kmax is None:
        kmax = 1e6  # 1/µm default

    k = np.linspace(kmin, kmax, n)
    Z = self.impedance(k)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(k * 1e-3, np.real(Z), label=r"Re[$Z$]")
    ax.plot(k * 1e-3, np.imag(Z), label=r"Im[$Z$]")
    ax.set_xlabel(r"$k$ (1/mm)")
    ax.set_ylabel(r"$Z(k)$ (Ω/m)")
    ax.legend()

pmd_beamphysics.wakefields.WakefieldBase.wake abstractmethod

wake(z)

Evaluate the wakefield at longitudinal position z.

Parameters:

Name	Type	Description	Default
z	float or ndarray	Longitudinal position [m]. Negative z is behind the source.	required

Returns:

Name	Type	Description
W	float or ndarray	Wakefield value [V/C/m]. Returns 0 for z > 0 (causality).

Source code in pmd_beamphysics/wakefields/base.py

@abstractmethod
def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at longitudinal position z.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0 (causality).
    """
    pass

pmd_beamphysics.wakefields.PseudomodeWakefield

PseudomodeWakefield(modes)

Bases: WakefieldBase

Wakefield represented as a sum of damped sinusoidal pseudomodes.

Models the longitudinal wakefield as:

\[W(z) = \sum_i A_i \cdot e^{d_i \cdot z} \cdot \sin(k_i \cdot z + \phi_i)\]

This form is used to approximate short-range wakefields such as the resistive wall wake.

Parameters:

Name	Type	Description	Default
modes	list of Pseudomode or list of dict	List of pseudomodes. Each can be a Pseudomode object or a dict with keys 'A', 'd', 'k', 'phi'.	required

Notes

• Wakefields are defined for z ≤ 0 (trailing the source particle).
• For z > 0, returns 0 (causality).

Examples:

Single mode::

pm = PseudomodeWakefield([Pseudomode(A=1e15, d=1e4, k=1e5, phi=np.pi/2)])
pm.wake(-10e-6)  # Wake at 10 µm behind source

Multiple modes::

modes = [
    Pseudomode(A=1e15, d=1e4, k=1e5, phi=np.pi/2),
    Pseudomode(A=5e14, d=2e4, k=2e5, phi=np.pi/4),
]
pm = PseudomodeWakefield(modes)

Source code in pmd_beamphysics/wakefields/pseudomode.py

def __init__(self, modes: list) -> None:
    self._modes = []
    for mode in modes:
        if isinstance(mode, Pseudomode):
            self._modes.append(mode)
        elif isinstance(mode, dict):
            self._modes.append(Pseudomode(**mode))
        else:
            raise TypeError(f"Mode must be Pseudomode or dict, got {type(mode)}")

Attributes

pmd_beamphysics.wakefields.PseudomodeWakefield.modes property

modes

List of Pseudomode objects.

pmd_beamphysics.wakefields.PseudomodeWakefield.n_modes property

n_modes

Number of pseudomodes.

Functions

pmd_beamphysics.wakefields.PseudomodeWakefield.impedance

impedance(k)

Compute the impedance Z(k) analytically.

Sums the contribution from each pseudomode using the analytic Fourier transform.

Parameters:

Name	Type	Description	Default
k	float or ndarray	Wavenumber [1/m].	required

Returns:

Name	Type	Description
Z	complex or ndarray	Impedance [Ohm/m].

Source code in pmd_beamphysics/wakefields/pseudomode.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) analytically.

    Sums the contribution from each pseudomode using the analytic
    Fourier transform.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    Z = np.zeros(len(k), dtype=complex)
    for mode in self._modes:
        Z += mode.impedance(k)

    if scalar_input:
        return complex(Z[0])
    return Z

pmd_beamphysics.wakefields.PseudomodeWakefield.particle_kicks

particle_kicks(z, weight, include_self_kick=True)

Compute short-range wakefield energy kicks per unit length.

Uses an O(N) single-pass algorithm for each mode by exploiting the exponential form of the pseudomode wakefield.

Parameters:

Name	Type	Description	Default
z	ndarray	Particle positions [m].	required
weight	ndarray	Particle charges [C].	required
include_self_kick	bool	If True, applies the self-kick. Default is True.	True

Returns:

Name	Type	Description
kicks	ndarray	Wake-induced energy kick per unit length [eV/m].

Source code in pmd_beamphysics/wakefields/pseudomode.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
) -> np.ndarray:
    """
    Compute short-range wakefield energy kicks per unit length.

    Uses an O(N) single-pass algorithm for each mode by exploiting
    the exponential form of the pseudomode wakefield.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        If True, applies the self-kick. Default is True.

    Returns
    -------
    kicks : np.ndarray
        Wake-induced energy kick per unit length [eV/m].
    """
    z = np.asarray(z)
    weight = np.asarray(weight)

    if z.shape != weight.shape:
        raise ValueError(
            f"Mismatched shapes: z.shape={z.shape}, weight.shape={weight.shape}"
        )
    if z.ndim != 1:
        raise ValueError("z and weight must be 1D arrays")

    # Sort particles from tail to head
    ix = z.argsort()
    z_sorted = z[ix].copy()
    z_sorted -= z_sorted.max()  # Offset for numerical stability
    weight_sorted = weight[ix]

    n = len(z)
    delta_E = np.zeros(n)

    # Process each mode with O(N) algorithm
    for mode in self._modes:
        s = mode.d + 1j * mode.k
        c = mode.A * np.exp(1j * mode.phi)

        # O(N) accumulator for this mode
        b = 0.0 + 0.0j

        # Iterate from tail to head
        for i in range(n - 1, -1, -1):
            zi = z_sorted[i]
            qi = weight_sorted[i]

            # Kick from trailing particles
            delta_E[i] -= np.imag(c * np.exp(s * zi) * b)

            # Add this particle to accumulator
            b += qi * np.exp(-s * zi)

        if include_self_kick:
            delta_E -= 0.5 * mode.A * weight_sorted * np.sin(mode.phi)

    # Restore original order
    kicks = np.empty_like(delta_E)
    kicks[ix] = delta_E

    return kicks

pmd_beamphysics.wakefields.PseudomodeWakefield.plot

plot(zmax=0.001, zmin=0, n=200, ax=None)

Plot the pseudomode wakefield.

Source code in pmd_beamphysics/wakefields/pseudomode.py

def plot(self, zmax: float = 0.001, zmin: float = 0, n: int = 200, ax=None):
    """Plot the pseudomode wakefield."""
    zlist = np.linspace(zmin, zmax, n)
    Wz = self.wake(-zlist)

    if ax is None:
        _, ax = plt.subplots()
    ax.plot(zlist * 1e6, Wz * 1e-12)
    ax.set_xlabel(r"Distance behind source $|z|$ (µm)")
    ax.set_ylabel(r"$W(z)$ (V/pC/m)")
    title = f"Pseudomode Wakefield ({self.n_modes} mode"
    if self.n_modes > 1:
        title += "s"
    title += ")"
    ax.set_title(title)

pmd_beamphysics.wakefields.PseudomodeWakefield.to_bmad

to_bmad(type='longitudinal', transverse_dependence='none')

Format pseudomode parameters as Bmad-compatible strings.

Parameters:

Name	Type	Description	Default
type	str	Wake type. Default is "longitudinal".	'longitudinal'
transverse_dependence	str	Transverse dependence. Default is "none".	'none'

Returns:

Type	Description
str	Bmad-formatted string (one line per mode).

Source code in pmd_beamphysics/wakefields/pseudomode.py

def to_bmad(
    self,
    type: str = "longitudinal",
    transverse_dependence: str = "none",
) -> str:
    """
    Format pseudomode parameters as Bmad-compatible strings.

    Parameters
    ----------
    type : str, optional
        Wake type. Default is "longitudinal".
    transverse_dependence : str, optional
        Transverse dependence. Default is "none".

    Returns
    -------
    str
        Bmad-formatted string (one line per mode).
    """
    lines = [mode.to_bmad(type, transverse_dependence) for mode in self._modes]
    return "\n".join(lines)

pmd_beamphysics.wakefields.PseudomodeWakefield.wake

wake(z)

Evaluate the pseudomode wakefield at position z.

Parameters:

Name	Type	Description	Default
z	float or ndarray	Longitudinal position [m]. Negative z is behind the source.	required

Returns:

Name	Type	Description
W	float or ndarray	Wakefield value [V/C/m]. Returns 0 for z > 0.

Source code in pmd_beamphysics/wakefields/pseudomode.py

def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the pseudomode wakefield at position z.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0.
    """
    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    out = np.zeros_like(z, dtype=float)
    mask = z <= 0
    z_masked = z[mask]

    # Sum contributions from all modes
    for mode in self._modes:
        out[mask] += mode(z_masked)

    if scalar_input:
        return float(out[0])
    return out

pmd_beamphysics.wakefields.ImpedanceWakefield

ImpedanceWakefield(impedance_func, wakefield_func=None, k_max=10000000.0)

Bases: WakefieldBase

Wakefield defined through its longitudinal impedance Z(k).

Uses FFT-based convolution in frequency domain for efficient wake potential calculations.

Parameters:

Name	Type	Description	Default
impedance_func	callable	Function that returns complex impedance Z(k) [Ohm/m] for wave number k [1/m]. Should handle arrays.	required
wakefield_func	callable	Function that returns wakefield W(z) [V/C/m] for position z [m]. If not provided, uses numerical cosine transform of impedance.	None
k_max	float	Maximum wavenumber for integration [1/m]. Default is 1e7.	10000000.0

Notes

The wakefield W(z) is defined for z ≤ 0 (trailing particles), with W(z) = 0 for z > 0 (causality). The impedance Z(k) is related to the wakefield by:

\[Z(k) = \frac{1}{c} \int_{-\infty}^{0} W(z) e^{-ikz} dz\]

And the inverse:

\[W(z) = \frac{2c}{\pi} \int_0^{\infty} \text{Re}[Z(k)] \cos(kz) dk \quad (z \le 0)\]

Examples:

::

def my_impedance(k):
    return 100 / (1 + 1j * k * 1e-6)  # Simple resonator
wake = ImpedanceWakefield(my_impedance)

Source code in pmd_beamphysics/wakefields/impedance.py

def __init__(
    self,
    impedance_func: Callable[[np.ndarray], np.ndarray],
    wakefield_func: Callable[[np.ndarray], np.ndarray] = None,
    k_max: float = 1e7,
) -> None:
    self._impedance_func = impedance_func
    self._wakefield_func = wakefield_func
    self._k_max = k_max

Attributes

pmd_beamphysics.wakefields.ImpedanceWakefield.k_max property writable

k_max

Maximum wavenumber for integration.

Functions

pmd_beamphysics.wakefields.ImpedanceWakefield.impedance

impedance(k)

Evaluate the impedance at wavenumber k.

Parameters:

Name	Type	Description	Default
k	float or ndarray	Wavenumber [1/m].	required

Returns:

Name	Type	Description
Z	complex or ndarray	Impedance [Ohm/m].

Source code in pmd_beamphysics/wakefields/impedance.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Evaluate the impedance at wavenumber k.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    return self._impedance_func(k)

pmd_beamphysics.wakefields.ImpedanceWakefield.particle_kicks

particle_kicks(z, weight, include_self_kick=True, n_bins=None)

Compute wakefield-induced energy kicks per unit length.

Uses FFT-based density convolution with interpolation back to particle positions. This is O(N + M log M) where N is the number of particles and M is the number of grid points.

Parameters:

Name	Type	Description	Default
z	ndarray	Particle positions [m].	required
weight	ndarray	Particle charges [C].	required
include_self_kick	bool	Whether to include the self-kick term. Default is True.	True
n_bins	int	Number of bins for the density grid. Default is max(100, N//10).	None

Returns:

Name	Type	Description
kicks	ndarray	Energy kick per unit length at each particle [eV/m].

Source code in pmd_beamphysics/wakefields/impedance.py

def particle_kicks(
    self,
    z: np.ndarray,
    weight: np.ndarray,
    include_self_kick: bool = True,
    n_bins: int = None,
) -> np.ndarray:
    """
    Compute wakefield-induced energy kicks per unit length.

    Uses FFT-based density convolution with interpolation back to
    particle positions. This is O(N + M log M) where N is the number
    of particles and M is the number of grid points.

    Parameters
    ----------
    z : np.ndarray
        Particle positions [m].
    weight : np.ndarray
        Particle charges [C].
    include_self_kick : bool, optional
        Whether to include the self-kick term. Default is True.
    n_bins : int, optional
        Number of bins for the density grid. Default is max(100, N//10).

    Returns
    -------
    kicks : np.ndarray
        Energy kick per unit length at each particle [eV/m].
    """
    from scipy.interpolate import interp1d

    z = np.asarray(z)
    weight = np.asarray(weight)
    n = len(z)

    if n_bins is None:
        n_bins = max(100, n // 10)

    # Bin particles into density distribution
    z_min, z_max = z.min(), z.max()
    z_range = z_max - z_min
    if z_range == 0:
        z_range = 1e-6  # Avoid division by zero for single particle

    # Add padding to avoid edge effects
    padding = 0.1 * z_range
    z_min_padded = z_min - padding
    z_max_padded = z_max + padding

    dz = (z_max_padded - z_min_padded) / n_bins
    z_grid = np.linspace(z_min_padded + dz / 2, z_max_padded - dz / 2, n_bins)

    # Create density histogram (charge per bin)
    density, bin_edges = np.histogram(
        z, bins=n_bins, range=(z_min_padded, z_max_padded), weights=weight
    )
    # Convert to charge density [C/m]
    density = density / dz

    # Compute wake potential via FFT convolution
    wake_potential = self.convolve_density(
        density, dz, include_self_kick=include_self_kick
    )

    # Interpolate wake potential to particle positions
    interp = interp1d(
        z_grid, wake_potential, kind="linear", bounds_error=False, fill_value=0.0
    )
    kicks = -interp(z)  # Negative sign: energy loss for trailing particles

    return kicks

pmd_beamphysics.wakefields.ImpedanceWakefield.wake

wake(z, n_fft=4096)

Evaluate the wakefield at position z.

If a wakefield function was provided, uses it directly. Otherwise, computes via cosine transform of Re[Z(k)]: - Scalars: numerical quadrature (accurate) - Arrays: FFT with interpolation (fast)

Parameters:

Name	Type	Description	Default
z	float or ndarray	Longitudinal position [m].	required
n_fft	int	Number of FFT points for array evaluation. Default is 4096.	4096

Returns:

Name	Type	Description
W	float or ndarray	Wakefield value [V/C/m]. Returns 0 for z > 0.

Source code in pmd_beamphysics/wakefields/impedance.py

def wake(self, z: np.ndarray | float, n_fft: int = 4096) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z.

    If a wakefield function was provided, uses it directly.
    Otherwise, computes via cosine transform of Re[Z(k)]:
    - Scalars: numerical quadrature (accurate)
    - Arrays: FFT with interpolation (fast)

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m].
    n_fft : int, optional
        Number of FFT points for array evaluation. Default is 4096.

    Returns
    -------
    W : float or np.ndarray
        Wakefield value [V/C/m]. Returns 0 for z > 0.
    """
    if self._wakefield_func is not None:
        return self._wakefield_func(z)

    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    if scalar_input:
        # Single point: use quadrature for accuracy
        from scipy.integrate import quad

        z_val = float(z[0])
        if z_val > 0:
            return 0.0

        def integrand(k):
            return np.real(self._impedance_func(k)) * np.cos(k * z_val)

        result, _ = quad(integrand, 0, self._k_max, limit=200)
        return (2 * c_light / np.pi) * result
    else:
        # Array: use FFT for speed
        from scipy.fft import irfft
        from scipy.interpolate import interp1d

        k_max = self._k_max
        dk = k_max / (n_fft - 1)
        k_grid = np.linspace(0, k_max, n_fft)

        Zk_grid = self._impedance_func(k_grid)
        ReZ = np.real(Zk_grid)

        n_full = 2 * (n_fft - 1)
        dz = 2 * np.pi / (n_full * dk)
        z_grid = np.arange(n_full) * dz

        W_grid = irfft(ReZ, n=n_full) * n_full * dk * (c_light / np.pi)

        interp = interp1d(
            z_grid, W_grid, kind="cubic", bounds_error=False, fill_value=0.0
        )

        result = interp(-z)
        result = np.where(z > 0, 0.0, result)

        return result

pmd_beamphysics.wakefields.TabularWakefield

TabularWakefield(z, W, fill_value=0.0, kind='cubic')

Bases: WakefieldBase

Wakefield defined by user-supplied tabular data with interpolation.

Uses cubic spline interpolation to evaluate the wakefield at arbitrary positions between the supplied data points.

Parameters:

Name	Type	Description	Default
z	ndarray	Longitudinal positions [m]. Should be negative (behind source) and sorted in ascending order.	required
W	ndarray	Wakefield values [V/C/m] at each z position.	required
fill_value	float	Value to return outside the interpolation range. Default is 0.	0.0
kind	str	Interpolation method. Default is 'cubic'.	'cubic'

Examples:

::

z_data = -np.linspace(1e-6, 1e-3, 100)
W_data = 1e15 * np.exp(z_data / 100e-6) * np.sin(1e5 * z_data)
wake = TabularWakefield(z_data, W_data)
wake.wake(-50e-6)  # Interpolated wake at 50 µm behind source

Source code in pmd_beamphysics/wakefields/tabular.py

def __init__(
    self,
    z: np.ndarray,
    W: np.ndarray,
    fill_value: float = 0.0,
    kind: str = "cubic",
) -> None:
    z = np.asarray(z)
    W = np.asarray(W)

    if z.shape != W.shape:
        raise ValueError(f"Shape mismatch: z.shape={z.shape}, W.shape={W.shape}")
    if len(z) < 4:
        raise ValueError("Need at least 4 points for cubic interpolation")

    # Store data
    self._z = z
    self._W = W
    self._fill_value = fill_value

    # Create interpolator
    self._interp = interp1d(
        z,
        W,
        kind=kind,
        bounds_error=False,
        fill_value=fill_value,
    )

Attributes

pmd_beamphysics.wakefields.TabularWakefield.W_data property

W_data

Return the W data points.

pmd_beamphysics.wakefields.TabularWakefield.z_data property

z_data

Return the z data points.

Functions

pmd_beamphysics.wakefields.TabularWakefield.impedance

impedance(k)

Compute the impedance Z(k) via numerical FFT.

Uses FFT to compute the Fourier transform of the tabular wake data.

Parameters:

Name	Type	Description	Default
k	float or ndarray	Wavenumber [1/m].	required

Returns:

Name	Type	Description
Z	complex or ndarray	Impedance [Ohm/m].

Source code in pmd_beamphysics/wakefields/tabular.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) via numerical FFT.

    Uses FFT to compute the Fourier transform of the tabular wake data.

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    # Use the stored data for FFT
    z_data = self._z
    W_data = self._W

    # Sort by z (ascending)
    sort_idx = np.argsort(z_data)
    z_sorted = z_data[sort_idx]
    W_sorted = W_data[sort_idx]

    # Z(k) = (1/c) * integral of W(z) * exp(-ikz) dz
    # For each k, compute numerical integral
    Z = np.zeros(len(k), dtype=complex)
    for i, ki in enumerate(k):
        integrand = W_sorted * np.exp(-1j * ki * z_sorted)
        Z[i] = np.trapezoid(integrand, z_sorted) / c_light

    if scalar_input:
        return complex(Z[0])
    return Z

pmd_beamphysics.wakefields.TabularWakefield.wake

wake(z)

Evaluate the wakefield at position z using interpolation.

Parameters:

Name	Type	Description	Default
z	float or ndarray	Longitudinal position [m].	required

Returns:

Name	Type	Description
W	float or ndarray	Interpolated wakefield value [V/C/m]. Returns fill_value outside the data range, and 0 for z > 0 (causality).

Source code in pmd_beamphysics/wakefields/tabular.py

def wake(self, z: np.ndarray | float) -> np.ndarray | float:
    """
    Evaluate the wakefield at position z using interpolation.

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m].

    Returns
    -------
    W : float or np.ndarray
        Interpolated wakefield value [V/C/m]. Returns fill_value
        outside the data range, and 0 for z > 0 (causality).
    """
    z = np.asarray(z)
    scalar_input = z.ndim == 0
    z = np.atleast_1d(z)

    # Apply causality
    result = np.where(z > 0, 0.0, self._interp(z))

    if scalar_input:
        return float(result[0])
    return result

pmd_beamphysics.wakefields.Pseudomode dataclass

Pseudomode(A, d, k, phi)

Single pseudomode parameters for a damped sinusoidal wakefield component.

Represents one term in the sum: W(z) = A·exp(d·z)·sin(k·z + φ)

Parameters:

Name	Type	Description	Default
A	float	Amplitude coefficient [V/C/m].	required
d	float	Exponential decay rate [1/m].	required
k	float	Oscillation wavenumber [1/m].	required
phi	float	Phase offset [rad].	required

Functions

pmd_beamphysics.wakefields.Pseudomode.impedance

impedance(k)

Compute the impedance Z(k) for this pseudomode analytically.

For a pseudomode wakefield defined for z ≤ 0:

\[W(z) = A \cdot e^{dz} \cdot \sin(k_0 z + \phi)\]

The impedance is defined as:

\[Z(k) = \frac{1}{c} \int_{-\infty}^{0} W(z) \cdot e^{-ikz} \, dz\]

Derivation:

Substituting the wakefield:

\[Z(k) = \frac{A}{c} \int_{-\infty}^{0} e^{dz} \sin(k_0 z + \phi) e^{-ikz} \, dz\]

Using \(\sin(\theta) = \frac{e^{i\theta} - e^{-i\theta}}{2i}\):

\[Z(k) = \frac{A}{2ic} \left[ e^{i\phi} \int_{-\infty}^{0} e^{(d + i(k_0 - k))z} dz - e^{-i\phi} \int_{-\infty}^{0} e^{(d - i(k_0 + k))z} dz \right]\]

Since \(d > 0\), both integrals converge:

\[\int_{-\infty}^{0} e^{az} dz = \frac{1}{a} \quad \text{for } \text{Re}(a) > 0\]

Evaluating:

\[Z(k) = \frac{A}{2ic} \left[ \frac{e^{i\phi}}{d + i(k_0 - k)} - \frac{e^{-i\phi}}{d - i(k_0 + k)} \right]\]

Mathematica verification:

(* Define the wakefield and compute impedance numerically *)
W[z_, A_, d_, k0_, phi_] := A Exp[d z] Sin[k0 z + phi]

(* Analytical result *)
Zanalytic[k_, A_, d_, k0_, phi_] :=
    A/(2 I c) (Exp[I phi]/(d + I (k0 - k)) - Exp[-I phi]/(d - I (k0 + k)))

(* Numerical integration (should match) *)
Znumeric[k_, A_, d_, k0_, phi_] :=
    1/c NIntegrate[W[z, A, d, k0, phi] Exp[-I k z], {z, -Infinity, 0}]

(* Test with sample values *)
c = 299792458;
{A, d, k0, phi} = {1*^15, 1*^4, 1*^5, Pi/4};
ktest = 5*^4;
{Zanalytic[ktest, A, d, k0, phi], Znumeric[ktest, A, d, k0, phi]}
(* Both should give the same complex number *)

Parameters:

Name	Type	Description	Default
k	float or ndarray	Wavenumber [1/m].	required

Returns:

Name	Type	Description
Z	complex or ndarray	Impedance [Ohm/m].

Source code in pmd_beamphysics/wakefields/pseudomode.py

def impedance(self, k: np.ndarray | float) -> np.ndarray | complex:
    """
    Compute the impedance Z(k) for this pseudomode analytically.

    For a pseudomode wakefield defined for z ≤ 0:

    $$W(z) = A \\cdot e^{dz} \\cdot \\sin(k_0 z + \\phi)$$

    The impedance is defined as:

    $$Z(k) = \\frac{1}{c} \\int_{-\\infty}^{0} W(z) \\cdot e^{-ikz} \\, dz$$

    **Derivation:**

    Substituting the wakefield:

    $$Z(k) = \\frac{A}{c} \\int_{-\\infty}^{0} e^{dz} \\sin(k_0 z + \\phi) e^{-ikz} \\, dz$$

    Using $\\sin(\\theta) = \\frac{e^{i\\theta} - e^{-i\\theta}}{2i}$:

    $$Z(k) = \\frac{A}{2ic} \\left[
        e^{i\\phi} \\int_{-\\infty}^{0} e^{(d + i(k_0 - k))z} dz
        - e^{-i\\phi} \\int_{-\\infty}^{0} e^{(d - i(k_0 + k))z} dz
    \\right]$$

    Since $d > 0$, both integrals converge:

    $$\\int_{-\\infty}^{0} e^{az} dz = \\frac{1}{a} \\quad \\text{for } \\text{Re}(a) > 0$$

    Evaluating:

    $$Z(k) = \\frac{A}{2ic} \\left[
        \\frac{e^{i\\phi}}{d + i(k_0 - k)}
        - \\frac{e^{-i\\phi}}{d - i(k_0 + k)}
    \\right]$$

    **Mathematica verification:**

    ```mathematica
    (* Define the wakefield and compute impedance numerically *)
    W[z_, A_, d_, k0_, phi_] := A Exp[d z] Sin[k0 z + phi]

    (* Analytical result *)
    Zanalytic[k_, A_, d_, k0_, phi_] :=
        A/(2 I c) (Exp[I phi]/(d + I (k0 - k)) - Exp[-I phi]/(d - I (k0 + k)))

    (* Numerical integration (should match) *)
    Znumeric[k_, A_, d_, k0_, phi_] :=
        1/c NIntegrate[W[z, A, d, k0, phi] Exp[-I k z], {z, -Infinity, 0}]

    (* Test with sample values *)
    c = 299792458;
    {A, d, k0, phi} = {1*^15, 1*^4, 1*^5, Pi/4};
    ktest = 5*^4;
    {Zanalytic[ktest, A, d, k0, phi], Znumeric[ktest, A, d, k0, phi]}
    (* Both should give the same complex number *)
    ```

    Parameters
    ----------
    k : float or np.ndarray
        Wavenumber [1/m].

    Returns
    -------
    Z : complex or np.ndarray
        Impedance [Ohm/m].
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    d = self.d
    k0 = self.k
    A = self.A
    phi = self.phi

    # Z(k) = A/(2ic) * [e^{iφ}/(d + i(k₀-k)) - e^{-iφ}/(d - i(k₀+k))]
    denom1 = d + 1j * (k0 - k)
    denom2 = d - 1j * (k0 + k)

    Z = (A / (2j * c_light)) * (
        np.exp(1j * phi) / denom1 - np.exp(-1j * phi) / denom2
    )

    if scalar_input:
        return complex(Z[0])
    return Z

pmd_beamphysics.wakefields.Pseudomode.to_bmad

to_bmad(type='longitudinal', transverse_dependence='none')

Format as Bmad-compatible string.

Source code in pmd_beamphysics/wakefields/pseudomode.py

def to_bmad(
    self,
    type: str = "longitudinal",
    transverse_dependence: str = "none",
) -> str:
    """Format as Bmad-compatible string."""
    return (
        f"{type} = {{{self.A}, {self.d}, {self.k}, "
        f"{self.phi / (2 * np.pi)}, {transverse_dependence}}}"
    )

Low-level Functions

pmd_beamphysics.wakefields.longitudinal_impedance_round

longitudinal_impedance_round(k, a, sigma0, ctau)

Compute longitudinal impedance Z(k) for round (circular pipe) geometry.

Uses closed-form expression (SLAC-PUB-10707 Eq. 2).

Parameters:

Name	Type	Description	Default
k	float or ndarray	Longitudinal wave number [1/m]	required
a	float	Pipe radius [m]	required
sigma0	float	DC conductivity [S/m]	required
ctau	float	Relaxation distance c·τ [m]	required

Returns:

Name	Type	Description
Zk	complex or np.ndarray of complex	Longitudinal impedance [Ohm/m]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def longitudinal_impedance_round(
    k: float | np.ndarray,
    a: float,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Compute longitudinal impedance Z(k) for round (circular pipe) geometry.

    Uses closed-form expression (SLAC-PUB-10707 Eq. 2).

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    a : float
        Pipe radius [m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    Zk : complex or np.ndarray of complex
        Longitudinal impedance [Ohm/m]
    """
    k = np.asarray(k)
    scalar_input = k.ndim == 0
    k = np.atleast_1d(k)

    prefactor = Z0 / (2 * np.pi * a)
    zeta = surface_impedance(k, sigma0, ctau)

    with np.errstate(divide="ignore", invalid="ignore"):
        Zk = prefactor / (1.0 / zeta - 1j * k * a / 2)
        Zk = np.where(k == 0, 0.0 + 0.0j, Zk)

    if scalar_input:
        return complex(Zk[0])
    return Zk

pmd_beamphysics.wakefields.longitudinal_impedance_flat

longitudinal_impedance_flat(k, a, sigma0, ctau)

Compute longitudinal impedance Z(k) for flat (parallel plate) geometry.

Uses numerical integration over transverse modes (SLAC-PUB-10707 Eq. 52).

Parameters:

Name	Type	Description	Default
k	float or ndarray	Longitudinal wave number [1/m]	required
a	float	Half-gap height between parallel plates [m]	required
sigma0	float	DC conductivity [S/m]	required
ctau	float	Relaxation distance c·τ [m]	required

Returns:

Name	Type	Description
Zk	complex or np.ndarray of complex	Longitudinal impedance [Ohm/m]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def longitudinal_impedance_flat(
    k: float | np.ndarray,
    a: float,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Compute longitudinal impedance Z(k) for flat (parallel plate) geometry.

    Uses numerical integration over transverse modes (SLAC-PUB-10707 Eq. 52).

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    a : float
        Half-gap height between parallel plates [m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    Zk : complex or np.ndarray of complex
        Longitudinal impedance [Ohm/m]
    """

    @np.vectorize
    def _Zk_scalar(k_val):
        if k_val == 0:
            return 0.0 + 0.0j

        def integrand(x):
            return _impedance_integrand_flat(k_val, x, a, sigma0, ctau)

        return quad_vec(integrand, 0, np.inf)[0]

    return _Zk_scalar(k)

pmd_beamphysics.wakefields.wakefield_from_impedance

wakefield_from_impedance(z, Zk_func, k_max=10000000.0, epsabs=1e-09, epsrel=1e-06)

Compute wakefield W(z) from Re[Z(k)] using a cosine transform (quadrature).

For z ≤ 0 (trailing particles):

\[W(z) = \frac{2c}{\pi} \int_0^{k_{\max}} \text{Re}[Z(k)] \cos(kz) \, dk\]

Returns 0 for z > 0 (causality).

Parameters:

Name	Type	Description	Default
z	float or ndarray	Longitudinal position [m]. Negative z is behind the source.	required
Zk_func	callable	Function returning complex impedance Z(k) [Ohm/m]	required
k_max	float	Upper limit of k integration [1/m]. Default is 1e7.	10000000.0
epsabs	float	Absolute tolerance. Default is 1e-9.	1e-09
epsrel	float	Relative tolerance. Default is 1e-6.	1e-06

Returns:

Name	Type	Description
Wz	float or ndarray	Wakefield [V/C/m]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def wakefield_from_impedance(
    z: float | np.ndarray,
    Zk_func: callable,
    k_max: float = 1e7,
    epsabs: float = 1e-9,
    epsrel: float = 1e-6,
) -> float | np.ndarray:
    """
    Compute wakefield W(z) from Re[Z(k)] using a cosine transform (quadrature).

    For z ≤ 0 (trailing particles):

    $$W(z) = \\frac{2c}{\\pi} \\int_0^{k_{\\max}} \\text{Re}[Z(k)] \\cos(kz) \\, dk$$

    Returns 0 for z > 0 (causality).

    Parameters
    ----------
    z : float or np.ndarray
        Longitudinal position [m]. Negative z is behind the source.
    Zk_func : callable
        Function returning complex impedance Z(k) [Ohm/m]
    k_max : float, optional
        Upper limit of k integration [1/m]. Default is 1e7.
    epsabs : float, optional
        Absolute tolerance. Default is 1e-9.
    epsrel : float, optional
        Relative tolerance. Default is 1e-6.

    Returns
    -------
    Wz : float or np.ndarray
        Wakefield [V/C/m]
    """

    def _wakefield_scalar(z_val):
        if z_val > 0:
            return 0.0

        def integrand(k):
            return np.real(Zk_func(k)) * np.cos(k * z_val)

        result, _ = quad(integrand, 0, k_max, epsabs=epsabs, epsrel=epsrel, limit=200)
        return (2 * c_light / np.pi) * result

    z = np.asarray(z)
    if z.ndim == 0:
        return _wakefield_scalar(float(z))

    return np.array([_wakefield_scalar(zi) for zi in z])

pmd_beamphysics.wakefields.wakefield_from_impedance_fft

wakefield_from_impedance_fft(z, Zk_func, k_max=10000000.0, n_fft=8192)

Compute wakefield W(z) from Z(k) using FFT and interpolation.

Much faster than wakefield_from_impedance for arrays.

Parameters:

Name	Type	Description	Default
z	ndarray	Longitudinal positions [m]. Negative z is behind the source.	required
Zk_func	callable	Function returning complex impedance Z(k) [Ohm/m]	required
k_max	float	Upper limit of k integration [1/m]. Default is 1e7.	10000000.0
n_fft	int	Number of FFT points. Default is 8192.	8192

Returns:

Name	Type	Description
Wz	ndarray	Wakefield [V/C/m]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def wakefield_from_impedance_fft(
    z: np.ndarray,
    Zk_func: callable,
    k_max: float = 1e7,
    n_fft: int = 8192,
) -> np.ndarray:
    """
    Compute wakefield W(z) from Z(k) using FFT and interpolation.

    Much faster than `wakefield_from_impedance` for arrays.

    Parameters
    ----------
    z : np.ndarray
        Longitudinal positions [m]. Negative z is behind the source.
    Zk_func : callable
        Function returning complex impedance Z(k) [Ohm/m]
    k_max : float, optional
        Upper limit of k integration [1/m]. Default is 1e7.
    n_fft : int, optional
        Number of FFT points. Default is 8192.

    Returns
    -------
    Wz : np.ndarray
        Wakefield [V/C/m]
    """
    from scipy.fft import irfft
    from scipy.interpolate import interp1d

    z = np.asarray(z)

    dk = k_max / (n_fft - 1)
    k_grid = np.linspace(0, k_max, n_fft)

    Zk_grid = Zk_func(k_grid)
    ReZ = np.real(Zk_grid)

    n_full = 2 * (n_fft - 1)
    dz = 2 * np.pi / (n_full * dk)
    z_grid = np.arange(n_full) * dz

    W_grid = irfft(ReZ, n=n_full) * n_full * dk * (c_light / np.pi)

    interp = interp1d(z_grid, W_grid, kind="cubic", bounds_error=False, fill_value=0.0)

    result = interp(-z)
    result = np.where(z > 0, 0.0, result)

    return result

pmd_beamphysics.wakefields.ac_conductivity

ac_conductivity(k, sigma0, ctau)

Frequency-dependent AC conductivity with relaxation time (Drude model).

\[\sigma(k) = \frac{\sigma_0}{1 - i k c \tau}\]

Parameters:

Name	Type	Description	Default
k	float or ndarray	Longitudinal wave number [1/m]	required
sigma0	float	DC conductivity [S/m]	required
ctau	float	Relaxation distance c·τ [m]	required

Returns:

Name	Type	Description
sigma	complex or np.ndarray of complex	AC conductivity [S/m]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def ac_conductivity(
    k: float | np.ndarray,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Frequency-dependent AC conductivity with relaxation time (Drude model).

    $$\\sigma(k) = \\frac{\\sigma_0}{1 - i k c \\tau}$$

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    sigma : complex or np.ndarray of complex
        AC conductivity [S/m]
    """
    return sigma0 / (1 - 1j * k * ctau)

pmd_beamphysics.wakefields.surface_impedance

surface_impedance(k, sigma0, ctau)

Surface impedance for a conducting wall with AC conductivity.

\[\zeta(k) = (1 - i) \sqrt{\frac{k c}{2 \sigma(k) Z_0 c}}\]

Parameters:

Name	Type	Description	Default
k	float or ndarray	Longitudinal wave number [1/m]	required
sigma0	float	DC conductivity [S/m]	required
ctau	float	Relaxation distance c·τ [m]	required

Returns:

Name	Type	Description
zeta	complex or np.ndarray of complex	Surface impedance [dimensionless]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def surface_impedance(
    k: float | np.ndarray,
    sigma0: float,
    ctau: float,
) -> complex | np.ndarray:
    """
    Surface impedance for a conducting wall with AC conductivity.

    $$\\zeta(k) = (1 - i) \\sqrt{\\frac{k c}{2 \\sigma(k) Z_0 c}}$$

    Parameters
    ----------
    k : float or np.ndarray
        Longitudinal wave number [1/m]
    sigma0 : float
        DC conductivity [S/m]
    ctau : float
        Relaxation distance c·τ [m]

    Returns
    -------
    zeta : complex or np.ndarray of complex
        Surface impedance [dimensionless]
    """
    sigma = ac_conductivity(k, sigma0, ctau)
    return (1 - 1j) * np.sqrt(k * c_light / (2 * sigma * Z0 * c_light))

pmd_beamphysics.wakefields.characteristic_length

characteristic_length(a, sigma0)

Characteristic length scale s₀ for resistive wall wakefield.

From SLAC-PUB-10707 Eq. (5):

\[s_0 = \left( \frac{2 a^2}{Z_0 \sigma_0} \right)^{1/3}\]

Parameters:

Name	Type	Description	Default
a	float	Half-gap height (flat) or radius (round) [m]	required
sigma0	float	DC conductivity [S/m]	required

Returns:

Name	Type	Description
s0	float	Characteristic length [m]

Source code in pmd_beamphysics/wakefields/resistive_wall/base.py

def characteristic_length(a: float, sigma0: float) -> float:
    """
    Characteristic length scale s₀ for resistive wall wakefield.

    From SLAC-PUB-10707 Eq. (5):

    $$s_0 = \\left( \\frac{2 a^2}{Z_0 \\sigma_0} \\right)^{1/3}$$

    Parameters
    ----------
    a : float
        Half-gap height (flat) or radius (round) [m]
    sigma0 : float
        DC conductivity [S/m]

    Returns
    -------
    s0 : float
        Characteristic length [m]
    """
    return (2 * a**2 / (Z0 * sigma0)) ** (1 / 3)