All in One View

Last updated on 2026-04-06 | Edit this page

Overview

Questions

🤌 What is WarpX?
🤔 What is a PIC code?
🧐 What can I use WarpX for?

Objectives

💡 Understand the basics of PIC codes
🧑‍💻 Learn about the features of WarpX
🎯 Figure out if WarpX can be useful for you!

Overview of PIC

WarpX is a general purpose open-source high-performance Particle-In-Cell (PIC) code.
The PIC method is a very popular approach to simulate the dynamics of physical systems governed by relativistic electrodynamics. Plasmas ⭐ and beams ☄️, which are often made of charged particles that may travel almost at the speed of light, fall into this category. Additionally, other physical effects can be integrated into the PIC algorithm, such as different quantum 🎲 processes.

Here is a picture that condenses the core idea: there are particles and fields. The particles are approximated as macroparticles (usually each representative of many real particles), while the fields are approximated on a grid in space. The particles and the fields are updated, self-consistently, in a temporal loop.

Here is a more informative image that explains the core algorithmic steps. As the particles travel in space, they generate currents, which in turn generate an electromagnetic field. The electromagnetic field then push the particles via the Lorentz force. Therefore, the current density $\textbf{J}$ and the force $\textbf{F}_L$ are the quantitieis that connect the particles and the fields, or in PIC lingo, the macroparticles and the grid. Hence, the idea is the following. First, define a well-posed initial condition, and then iterate the following: * Interpolate the fields from the grid to the particles’ positions and compute the Lorentz force that acts on each macroparticle, * Advance the position and momenta of the macroparticles using Newton equations, * Project while cumulating the contribution to the current density of each macroparticles, * Solve Maxwell’s equations.

In some cases, one can choose to solve Poisson’s equation instead of Maxwell’s. In that case, the current $\textbf{J}$ calculation is replaced with the charge density $\rho$ calculation. Once $\rho$ is known, the electrostic potential is computed to then find the electric field.

If you want to know more about PIC, here are a few references:

The two bibles on PIC 📚
- C. K. Birdsall and A. B. Langdon. Plasma Physics Via Computer Simulation.
- R. W. Hockney and J. W. Eastwood. Computer simulation using particles.
An old review written by one of the pioneers: John M. Dawson, Particle simulation of plasmas, Rev. Mod. Phys. 55, 403
Browse our docs for many, many more references about advanced algorithms and methods.

Features and applications of WarpX

WarpX is developed and used by a wide range of researchers working in different fields, from beam physics to nuclear fusion and a lot more. To learn about what WarpX can be used for, check out our examples and the scientific publications that acknowledged WarpX.

Checklist

Some cool features of WarpX:

📖 Open-source!

✈️ Runs on GPUs: NVIDIA, AMD, and Intel

🚀 Runs on multiple GPUs or CPUs, on systems ranging from laptops to supercomputers

🤓 Many many advanced algorithms and methods: mesh-refinement, embedded boundaries, electrostatic/electromagnetic/pseudospectral solvers, etc.

💾 Standards: openPMD for input/output data, PICMI for inputs

🤸 Active development and mainteinance: check out our GitHub repo

🗺️ International, cross-disciplinary community: plasma physics, fusion devices, laser-plasma interactions, beam physics, plasma-based acceleration, astrophysics, others?

We will add here more details soon!

Key Points

🔮 The particle-in-cell method is used to simulate the self-consistent dynamics of relativistic charged particles

🚀 WarpX is a open-source high-performance particle-in-cell code

✨ WarpX is used in a variety of scientif domains

Content from Install

Last updated on 2026-09-19 | Edit this page

Overview

Questions

🔧 How can I install and run WarpX?
🕵️ How can I analyze the simulation results?

Objectives

💻 Install WarpX on your local machine either with Conda or from source
👀 Install the visualizations tools in Python and Paraview

Basic dependencies

Just a heads-up before we dive deeper.

Callout

📣 Everything you need to know to use WarpX is in the documentation, check it out!

Via Conda-Forge

First, you need a Conda installation and we will assume that you indeed have one.
If not, follow the instruction at this link.
You can install Conda on most operative systems: Windows, macOS, and Linux.
We will also assume you have some familiarity with the terminal. Once you have Conda on your system, WarpX is available as a package via Conda-Forge.
The installation is a one-liner 😌!

Callout

BASH

conda install -c conda-forge warpx

Ok, maybe two lines if you want to keep your system clean by creating a new environment.

BASH

conda create -n warpx -c conda-forge warpx 
conda activate warpx

Callout

The first warpx (after -n) is the name of the conda environment; the second warpx is the name of the package to be installed from the conda-forge channel.

Now you should have 4 different WarpX binaries in your PATH called warpx.1d, warpx.2d, warpx.3d, warpx.rz.
Each binary for a different dimensionality.

To check this, run:

BASH

which warpx.1d warpx.2d warpx.3d warpx.rz

If you get 3 different paths that look something like:

BASH

/home/<username>/anaconda3/envs/warpx/bin/warpx.xd

then you got this 🙌! You can also import pywarpx in Python.

Caution

Conda’s WarpX is serial! To get a parallel WarpX version, install it from source.

From source

Caution

Coming soon. For now, refer to the main documentation.

Key Points

🎯 WarpX is easy to install via Conda: conda -c conda-forge warpx

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

Content from Set Up a Simulation from Scratch

Last updated on 2026-04-08 | Edit this page

Overview

Questions

How do you go from “I want to simulate X” to a working WarpX input file?

Objectives

Navigate the WarpX repository to find relevant examples. Adapt an existing input file to your physics problem. Iterate on simulation parameters to observe the expected physics.

Introduction

Most simulations don’t start from a blank page. In practice, a physicist starting a new simulation will look for the closest existing example, adapt it to their problem, run it, look at the results, and iterate. This episode walks through that process step by step, using the two-stream instability as the target physics.

Rather than handing you a ready-made input file, we will start from a different example altogether – a uniform plasma – and modify it until we observe the instability. Along the way, we’ll make mistakes, debug them, and learn from the process.

Callout

This episode is complementary to the two-stream instability tutorial, which gives you a polished input file with exercises on parameter selection. Here we focus on the workflow: how to find, adapt, and iterate.

Step 1: What do we want to simulate?

Before touching any code, clarify the physics. We want to simulate the two-stream instability: two populations of electrons streaming through each other with opposite drift velocities, in a periodic box with a neutralizing background. The expected result is exponential growth of electrostatic perturbations and the formation of vortex structures (“cat eyes”) in phase space.

What do we need?

A periodic box with a uniform plasma (electrons)
Two counter-streaming populations with a relative drift velocity
Diagnostics to look at the phase space $(z, u_z)$

Step 2: Browse the WarpX examples

WarpX ships with many examples in its GitHub repository, organized under Examples/Physics_applications/. The documentation also has a gallery of examples with descriptions, input files, and analysis scripts.

Challenge

Find a starting point

Browse the examples in the WarpX repository or the documentation gallery. Which example is the closest to what we want?

Show me the solution

The uniform_plasma example is a good starting point: it already has a periodic box with a uniform electron population, a Gaussian (thermal) momentum distribution, and basic diagnostics. We “just” need to add a drift velocity and duplicate the species.

Step 3: Start from the uniform plasma example

Here is the 2D uniform plasma input file, rendered directly from the WarpX repository:

OUTPUT

#################################
####### GENERAL PARAMETERS ######
#################################
max_step = 10
amr.n_cell =  128 128
amr.max_grid_size = 64
amr.blocking_factor = 32
amr.max_level = 0
geometry.dims = 2
geometry.prob_lo = -20.e-6   -20.e-6    # physical domain
geometry.prob_hi =  20.e-6    20.e-6

#################################
####### Boundary condition ######
#################################
boundary.field_lo = periodic periodic
boundary.field_hi = periodic periodic

#################################
############ NUMERICS ###########
#################################
warpx.serialize_initial_conditions = 1
warpx.verbose = 1
warpx.cfl = 1.0
warpx.use_filter = 0

# Order of particle shape factors
algo.particle_shape = 1

#################################
############ PLASMA #############
#################################
particles.species_names = electrons

electrons.charge = -q_e
electrons.mass = m_e
electrons.injection_style = "NUniformPerCell"
electrons.num_particles_per_cell_each_dim = 2 2
electrons.profile = constant
electrons.density = 1.e25  # number of electrons per m^3
electrons.momentum_distribution_type = "gaussian"
electrons.ux_th  = 0.01 # uth the std of the (unitless) momentum
electrons.uy_th  = 0.01 # uth the std of the (unitless) momentum
electrons.uz_th  = 0.01 # uth the std of the (unitless) momentum

# Diagnostics
diagnostics.diags_names = diag1
diag1.intervals = 10
diag1.diag_type = Full

This simulates a single electron species in a box with thermal fluctuations. No drift, no instability.

Step 4: Add a drift velocity

To create two counter-streaming beams, we need to:

Duplicate the species: create ele1 and ele2 with opposite drift velocities
Add a bulk drift to each species along one direction

Here are the key modifications. Replace the single electrons species block with two species:

particles.species_names = ele1 ele2

ele1.species_type = electron
ele1.injection_style = "NUniformPerCell"
ele1.num_particles_per_cell_each_dim = 2 2
ele1.profile = constant
ele1.density = 1.e25
ele1.momentum_distribution_type = "gaussian"
ele1.ux_th = 0.01
ele1.uy_th = 0.01
ele1.uz_th = 0.01
ele1.uz_m = 0.1

ele2.species_type = electron
ele2.injection_style = "NUniformPerCell"
ele2.num_particles_per_cell_each_dim = 2 2
ele2.profile = constant
ele2.density = 1.e25
ele2.momentum_distribution_type = "gaussian"
ele2.ux_th = 0.01
ele2.uy_th = 0.01
ele2.uz_th = 0.01
ele2.uz_m = -0.1

The parameter uz_m sets the mean of the Gaussian distribution for the $z$-momentum (in units of $m c$), so +0.1 and -0.1 give two beams drifting in opposite directions at $\beta_0 \approx 0.1$.

Challenge

What parameters do we need to check?

Before running, take a look at the rest of the input file. Are there any parameters that might need adjusting for the two-stream case?

Show me the solution

A few things to keep in mind:

max_step = 10: way too few timesteps to see anything happen
Resolution: 128 cells may or may not be enough depending on the plasma parameters
Diagnostics: the uniform plasma example uses the default plotfile format, but we’ll want to switch to openPMD so we can use openPMD-viewer to grab particle data easily

We’ll deal with these one at a time.

Define `nt` as a constant

A practical tip: define the number of timesteps as a my_constants variable at the top of the input file, so you can change it in one place and have everything else follow:

my_constants.nt = 10

max_step = nt

This is especially handy because you’ll also want the diagnostic interval to scale with the total number of steps (see below).

Switch diagnostics to openPMD

The uniform plasma example writes diagnostics in WarpX’s native plotfile format. This works, but it’s much more convenient to use the openPMD standard instead, because then you can use the openPMD-viewer Python library to grab particle and field data with a couple of lines of code.

Replace the diagnostics block with:

diagnostics.diags_names = diag1
diag1.intervals = floor(nt/2)
diag1.diag_type = Full
diag1.fields_to_plot = none
diag1.format = openpmd

A few things to note:

diag1.intervals = floor(nt/2) writes 2 snapshots regardless of how many timesteps you run – so when you increase nt later, the diagnostic output adjusts automatically.
diag1.fields_to_plot = none skips the field data, since we only care about the particle phase space for now.
diag1.format = openpmd enables the use of openPMD-viewer for analysis.

Step 5: Run – and nothing happens

Let’s run the simulation as-is with max_step = 10, just to make sure the code doesn’t crash.

BASH

warpx.2d my_inputs.txt

The simulation runs in a few seconds. Now open a Jupyter notebook and plot the phase space $(z, u_z)$ using openPMD-viewer:

PYTHON

from openpmd_viewer import OpenPMDTimeSeries
import matplotlib.pyplot as plt

ts = OpenPMDTimeSeries("./diags/diag1/")
z, uz = ts.get_particle(["z", "uz"], species="ele1", iteration=ts.iterations[-1])
plt.scatter(z, uz, s=0.1, alpha=0.3)
z, uz = ts.get_particle(["z", "uz"], species="ele2", iteration=ts.iterations[-1])
plt.scatter(z, uz, s=0.1, alpha=0.3)
plt.xlabel("z [m]")
plt.ylabel("uz [m_e c]")
plt.show()

The phase space looks like… two flat lines. Nothing is happening. That’s probably because we only ran 10 timesteps, and the instability needs time to grow from noise.

Callout

This is normal. Not seeing what you expect on the first try is part of the process. The important thing is to understand why and to iterate.

Step 6: Iterate

This is where the real work (and the learning) begins. Each iteration teaches you something about the physics and the numerics.

Increase the number of timesteps

If you defined my_constants.nt as suggested earlier, just increase it: try my_constants.nt = 100 or 200 or more. Both max_step and the diagnostic interval will update automatically.

Reduce the temperature

The instability is driven by the relative drift between the two beams. If the temperature is too high, the thermal spread dominates over the drift and the instability won’t develop. Try setting the thermal spread to zero:

ele1.ux_th = 0.0
ele1.uy_th = 0.0
ele1.uz_th = 0.0

This is a “cold beam” limit. In practice, some finite temperature is physical, but starting cold makes it easier to see the instability develop cleanly.

Switch to 1D

The two-stream instability is inherently a 1D phenomenon (along the drift direction). Running in 2D or 3D adds transverse dimensions that increase the computational cost without changing the essential physics. Change to 1D:

geometry.dims = 1
amr.n_cell = 256
geometry.prob_lo = -20.e-6
geometry.prob_hi = 20.e-6

Increase resolution and particles per cell

With only 2 particles per cell, the noise level is high. Try increasing to 100 or more. Also make sure the grid is fine enough to resolve the plasma skin depth $c/\omega_{pe}$.

Challenge

Put it all together

Combine all the modifications above into a single input file. Run the simulation and plot the phase space at multiple time snapshots. Do you see the instability developing?

Show me the solution

You should see the two initially-flat streams begin to develop a sinusoidal modulation, followed by the formation of vortex (“cat eye”) structures as particles get trapped in the electrostatic potential wells. The number of timesteps needed depends on the drift velocity and density; you may need to adjust max_step until you see saturation.

Step 7: Visualize the result

Once the instability is developing, you can produce a sequence of phase-space plots. Loop over the diagnostic snapshots:

PYTHON

from openpmd_viewer import OpenPMDTimeSeries
import matplotlib.pyplot as plt

ts = OpenPMDTimeSeries("./diags/diag1/")
for i in ts.iterations:
    z, uz = ts.get_particle(["z", "uz"], iteration=i)
    plt.figure()
    plt.scatter(z, uz, s=0.1, alpha=0.3)
    plt.xlabel("z [m]")
    plt.ylabel("uz [m_e c]")
    plt.title(f"iteration {i}")
    plt.savefig(f"phase_space_{i:06d}.png", dpi=150)
    plt.close()

From these frames you can make a video using ffmpeg:

BASH

ffmpeg -framerate 10 -pattern_type glob -i 'phase_space_*.png' \
       -c:v libx264 -pix_fmt yuv420p phase_space.mp4

What’s next?

Now that you’ve seen the full workflow – from browsing examples to iterating on parameters to producing visualizations – you can try the polished two-stream instability tutorial, which provides a carefully prepared input file with exercises on choosing physical and numerical parameters, reduced diagnostics for tracking the growth rate, and analysis notebooks.

The key takeaway is the process:

Clarify the physics you want to simulate
Find the closest existing example
Modify it step by step
Run, visualize, and iterate
Don’t be discouraged when things don’t work on the first try
Ask yourself: is what I’m seeing physics, numerics, or both? This is one of the hardest questions in computational physics…

A note on timescales

The characteristic timescale of the two-stream instability is the inverse plasma frequency $\omega_{pe}^{-1}$, where $\omega_{pe} = \sqrt{n_0 e^2 / (m_e \epsilon_0)}$. The instability growth rate is of order $\omega_{pe}$, so the simulation needs to cover many plasma periods for the instability to develop and saturate. This is why increasing nt is the first thing to try when nothing seems to be happening.

Key Points

The WarpX examples and the documentation gallery are the best starting points for a new simulation.

Setting up a simulation is an iterative process: run, visualize, understand what’s wrong, fix it, repeat.

When the simulation doesn’t show the expected physics, check the basics first: enough timesteps? resolved scales? adequate particles per cell? temperature not too high?

The documentation parameter reference is your constant companion for understanding and modifying input files.

Content from A Two-stream Instability

Last updated on 2026-04-08 | Edit this page

Overview

Questions

🤼 How do two counter-streaming electron beams become unstable?

Objectives

Set up, run, and visualize a 1D two-stream instability simulation with WarpX. Understand how the choice of physical parameters affects the instability growth.

Introduction

The two-stream instability is one of the most fundamental kinetic instabilities in plasma physics. When two populations of charged particles (here, electrons) stream through each other with a relative drift velocity, small perturbations in the charge density can grow exponentially. The free energy stored in the relative drift is converted into electrostatic wave energy until the beams thermalize and the instability saturates.

This example uses a 1D periodic box with two counter-streaming electron populations on a neutralizing ion background.

Setup

Make sure to download the input file.

Whenever you need to prepare an input file, this is where you want to go.

Let’s take a look at the input file

OUTPUT

####################
### MY CONSTANTS ###
####################
# PLASMAS
my_constants.n0 = ... # [m^-3]
my_constants.T0 = ...*q_e # [J]
my_constants.v_te = sqrt(T0 / m_e) # [m/s]
my_constants.beta0 = ... # [-]
my_constants.omega_pe = sqrt(n0*q_e**2/(m_e*epsilon0)) # [1/s]
# BOX
my_constants.Lx = ...*clight/omega_pe
my_constants.nx = ...
my_constants.dx = Lx/nx
# TIME
my_constants.cfl = ...
my_constants.T = .../omega_pe
my_constants.dt = cfl * dx / clight
my_constants.nt = floor(T/dt)

##########################
### GENERAL PARAMETERS ###
##########################
stop_time = T
amr.n_cell = nx nx nx
amr.max_level = 0
geometry.dims = 1
geometry.prob_lo = -0.5*Lx -0.5*Lx  -0.5*Lx
geometry.prob_hi =  0.5*Lx  0.5*Lx   0.5*Lx

##########################
### BOUNDARY CONDITION ###
##########################
boundary.field_lo = periodic periodic periodic
boundary.field_hi = periodic periodic periodic
boundary.particle_lo = periodic periodic periodic
boundary.particle_hi = periodic periodic periodic

################
### NUMERICS ###
################
warpx.cfl = cfl
algo.maxwell_solver = yee
algo.particle_shape = 3
algo.particle_pusher = boris
warpx.use_filter = 1

#################
### PARTICLES ###
#################
particles.species_names = ele1 ele2

ele1.species_type = electron
ele1.injection_style = NRandomPerCell
ele1.num_particles_per_cell = 200
ele1.profile = constant
ele1.density = n0
ele1.momentum_distribution_type = maxwell_boltzmann
ele1.theta_distribution_type = constant
ele1.theta = T0 / (m_e * clight**2)
ele1.beta_distribution_type = parser
ele1.beta_function(x,y,z) = beta0
ele1.bulk_vel_dir = +z

ele2.species_type = electron
ele2.injection_style = NRandomPerCell
ele2.num_particles_per_cell = 200
ele2.profile = constant
ele2.density = n0
ele2.momentum_distribution_type = maxwell_boltzmann
ele2.theta_distribution_type = constant
ele2.theta = T0 / (m_e * clight**2)
ele2.beta_distribution_type = constant
ele2.beta = beta0
ele2.bulk_vel_dir = -z

###################
### DIAGNOSTICS ###
###################
# FULL
diagnostics.diags_names =  particles
particles.intervals = floor(nt/200)
particles.diag_type = Full
particles.species = ele1 ele2
particles.fields_to_plot = none
particles.format = openpmd
particles.openpmd_backend = bp
particles.dump_last_timestep = 1
particles.ele1.variables = w z uz
particles.ele2.variables = w z uz

# REDUCED
warpx.reduced_diags_names = FieldEnergy FieldMaximum
FieldEnergy.type = FieldEnergy
FieldEnergy.intervals = 1
FieldMaximum.type = FieldMaximum
FieldMaximum.intervals = 1

Choosing the parameters

You will notice that some parameters in the input file are set to ... (ellipsis). These are left for you to choose 🗳️! Picking sensible physical and numerical parameters is an important part of setting up a simulation. Here is some guidance.

Challenge

Choose the physical parameters

The physical parameters you need to set are:

Parameter	Symbol	Description
`n0`	$n_0$	Electron density of each beam $[\mathrm{m^{-3}}]$
`T0`	$T_0$	Thermal temperature of each beam $[\mathrm{eV}]$
`beta0`	$\beta_0 = v_d/c$	Drift velocity of each beam normalized to $c$

Hints

A typical plasma prototype is an ionized gas. What is the number density of an ionized gas at standard conditions? Think of a fluorescent lamp, a tokamak, or the solar wind.
The temperature should be low enough that the thermal spread does not wash out the instability: the drift velocity should be larger than the thermal velocity, i.e. $v_d \gg v_{\mathrm{th}}$.
The drift velocity $\beta_0$ controls how fast the instability grows. Non-relativistic drifts ($\beta_0 \ll 1$) are a good starting point.

Challenge

Choose the numerical parameters

The numerical parameters you need to set are:

Parameter	Description
`Lx`	Box length, expressed in units of $c/\omega_{pe}$
`nx`	Number of grid cells
`cfl`	CFL number (must be $< 1$)
`T`	Total simulation time, expressed in units of $\omega_{pe}^{-1}$

Hints

The most unstable wavelength is $\lambda \sim 2\pi\,\beta_0\,(c/\omega_{pe})$. The box must be large enough to fit several of these wavelengths.
The grid must resolve the shortest relevant scale, i.e. $\Delta x \lesssim c/\omega_{pe}$. So nx should be at least comparable to Lx (in units of $c/\omega_{pe}$).
The CFL should be $\lesssim 1$ for stability of the Yee solver.
The simulation must run long enough for the instability to develop and saturate: $\sim 50$–$200\,\omega_{pe}^{-1}$ is a typical range.

Notable details

Two electron species (ele1, ele2) are initialized with equal density n0 and temperature T0, but with opposite drift velocities $\pm \beta_0 c$ along the $z$-axis.
There is no ion background in this input. Because we do not ask WarpX to solve Poisson’s equation at initialization, the fields start at zero and the uniform charge density of the electrons never sources a DC field, effectively acting as if a neutralizing background were present.
The diagnostics write out particle data (z, uz) for both species, which lets you construct the phase space $(z, u_z)$ at each snapshot.
Reduced diagnostics (FieldEnergy, FieldMaximum) record the total field energy and the maximum field value at every timestep, so you can track the exponential growth of the instability without post-processing the full data.

Run

Create a new folder and copy the input file there. Then fill it with your chosen parameters. This is a 1D simulation, so it should run very quickly even in serial.

Let’s run

BASH

warpx.1d input_1d_two_stream_instability.txt

Just like that! 💃

You should see a standard output flashing out a lot of info.
At the end, you should find in your folder:

a subfolder diags/particles/ with the full particle diagnostics
files FieldEnergy.txt and FieldMaximum.txt with the reduced diagnostics, inside diags/reducedfiles/
a file called warpx_used_inputs: a summary of the inputs that were used to run the simulation

If that’s the case, yey! 💯

If the run went wrong, you may find a Backtrace.0.0 file which can be useful for debugging purposes.

Visualize

With Python 🐍

Now that we have the results, we can analyze them using Python.
We will use the openPMD-viewer library to grab the data that the simulation produced in openPMD format. Here you can find a few tutorials on how to use the viewer. If you feel nerdy and/or you need to deal with the data in parallel workflows, you can use the openPMD-api.

Here are some things you can visualize:

Phase space $(z, u_z)$: plot the particle positions and momenta at several time snapshots to see how the beams evolve from smooth streams into characteristic “cat-eye” vortices.
Field energy vs time: load FieldEnergy.txt and plot the total electrostatic energy as a function of time on a semi-log scale. You should see an exponential growth phase followed by saturation.
Growth rate: from the slope of $\ln(E_{\mathrm{field}})$ during the linear phase, extract the growth rate and compare it to the theoretical prediction $\gamma \sim \omega_{pe}$.

With a Jupyter notebook 📓

The notebook includes a phase-space plot $(z, u_z)$ and the evolution of the electric and magnetic field energy on a semi-log scale.

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

Challenge

Questions for analysis

Describe what you observe in the phase space $(z, u_z)$ at different times and in the field energy evolution. Can you identify distinct stages?
How does the growth rate of the field energy depend on the drift velocity $\beta_0$? Try running with two different values and compare.
What happens to the phase space after saturation? Can you identify particle trapping in the electrostatic potential wells?
If you increase the temperature (decrease the ratio $v_d / v_{\mathrm{th}}$), does the instability still develop? At what point is it suppressed?
How does the number of particles per cell affect the noise level and the onset of the instability?

Show me the solution

These are open-ended questions meant to guide your exploration. The key insight is that the instability requires $v_d > v_{\mathrm{th}}$; when thermal effects dominate, Landau damping stabilizes the modes.

Gallery

Phase space plot showing the two counter-streaming electron beams forming vortex structures at saturation

Key Points

💡 The two-stream instability converts the kinetic energy of counter-streaming beams into electrostatic wave energy.

🔬 The instability growth rate depends on the ratio of drift velocity to thermal velocity: $v_d / v_{\mathrm{th}} \gg 1$ is needed for the instability to develop.

📊 Reduced diagnostics (FieldEnergy, FieldMaximum) let you track the instability growth in real time without post-processing full datasets.

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.

Content from A Weibel Instability

Last updated on 2026-09-19 | Edit this page

Overview

Questions

How does a momentum anisotropy in a plasma generate magnetic fields from scratch?

Objectives

Set up, run, and visualize a 2D Weibel instability simulation with WarpX. Observe the spontaneous generation of magnetic fields from counter-streaming electron beams.

Introduction

The Weibel instability (also known as the current filamentation instability) is a fundamental electromagnetic instability that arises when there is an anisotropy in the momentum distribution of a plasma. Unlike the two-stream instability, which generates electrostatic waves, the Weibel instability generates magnetic fields.

The physical mechanism is as follows: when two electron populations stream past each other, any small magnetic perturbation will deflect the electrons, causing them to bunch into current filaments. These current filaments in turn amplify the magnetic field that caused the bunching, leading to exponential growth until the fields are strong enough to isotropize the distribution.

This instability plays a key role in astrophysical plasmas (e.g. collisionless shocks in gamma-ray bursts) and in laser-plasma interactions.

In this episode, we simulate it in 2D with two counter-streaming electron populations and periodic boundaries.

Setup

Make sure to download the input file.

Whenever you need to prepare an input file, this is where you want to go.

Let’s take a look at the input file

OUTPUT

####################
### MY CONSTANTS ###
####################
# PLASMAS
my_constants.n0 = ... # [m^-3]
my_constants.T0 = ...*q_e # [J]
my_constants.beta0 = ... # [-]
my_constants.omega_pe = sqrt(n0*q_e**2/(m_e*epsilon0)) # [1/s]
my_constants.skin_depth = clight / omega_pe # [m]
# BOX
my_constants.Lx = ...*skin_depth
my_constants.nx = ...
my_constants.dx = Lx/nx
# TIME
my_constants.cfl = ...
my_constants.T = .../omega_pe
my_constants.dt = cfl*dx/(sqrt(2)*clight)
my_constants.nt = floor(T/dt)

##########################
### GENERAL PARAMETERS ###
##########################
stop_time = T
amr.n_cell = nx nx nx
amr.max_level = 0
geometry.dims = 2
geometry.prob_lo = -0.5*Lx -0.5*Lx  -0.5*Lx
geometry.prob_hi =  0.5*Lx  0.5*Lx   0.5*Lx

##########################
### BOUNDARY CONDITION ###
##########################
boundary.field_lo = periodic periodic periodic
boundary.field_hi = periodic periodic periodic
boundary.particle_lo = periodic periodic periodic
boundary.particle_hi = periodic periodic periodic

################
### NUMERICS ###
################
warpx.cfl = cfl
algo.maxwell_solver = yee
algo.particle_shape = 3
algo.particle_pusher = boris
warpx.use_filter = 1

#################
### PARTICLES ###
#################
particles.species_names = ele1 ele2

ele1.species_type = electron
ele1.injection_style = NRandomPerCell
ele1.num_particles_per_cell = 80
ele1.profile = constant
ele1.density = n0
ele1.momentum_distribution_type = maxwell_boltzmann
ele1.theta_distribution_type = constant
ele1.theta = T0 / (m_e * clight**2)
ele1.beta_distribution_type = parser
ele1.beta_function(x,y,z) = beta0
ele1.bulk_vel_dir = +y

ele2.species_type = electron
ele2.injection_style = NRandomPerCell
ele2.num_particles_per_cell = 80
ele2.profile = constant
ele2.density = n0
ele2.momentum_distribution_type = maxwell_boltzmann
ele2.theta_distribution_type = constant
ele2.theta = T0 / (m_e * clight**2)
ele2.beta_distribution_type = constant
ele2.beta = beta0
ele2.bulk_vel_dir = -y

###################
### DIAGNOSTICS ###
###################
# FULL
diagnostics.diags_names =  fields
fields.intervals = floor(nt/200)
fields.diag_type = Full
fields.write_species = 0
fields.fields_to_plot = Bx Bz
fields.format = openpmd
fields.openpmd_backend = bp
fields.dump_last_timestep = 1
# REDUCED
warpx.reduced_diags_names = FieldEnergy FieldMaximum
FieldEnergy.type = FieldEnergy
FieldEnergy.intervals = 1
FieldMaximum.type = FieldMaximum
FieldMaximum.intervals = 1

Choosing the parameters

You will notice that some parameters in the input file are set to ... (ellipsis). These are left for you to choose! Designing a good simulation requires making physical and numerical choices –
this is part of the exercise.

Challenge

Choose the physical parameters

The physical parameters you need to set are:

Parameter	Symbol	Description
`n0`	$n_0$	Electron density of each beam $[\mathrm{m^{-3}}]$
`T0`	$T_0$	Thermal temperature of each beam $[\mathrm{eV}]$
`beta0`	$\beta_0 = v_d/c$	Drift velocity of each beam normalized to $c$

Hints

The Weibel instability requires an anisotropy in the momentum space. Here, the anisotropy comes from the relative drift between the two beams. A mildly relativistic drift (e.g. $\beta_0 \sim 0.1$–$0.5$) will produce a clear instability.
The temperature should be small enough that the thermal spread does not wash out the anisotropy. As a rule of thumb, $v_d \gg v_{\mathrm{th}}$ where $v_{\mathrm{th}} = \sqrt{T_0/m_e}$.
A typical plasma prototype is an ionized gas. What is the number density of an ionized gas at standard conditions? Think of a fluorescent lamp, a tokamak, or the solar wind. What matters is that lengths and times are measured in units of $c/\omega_{pe}$ and $\omega_{pe}^{-1}$.

Challenge

Choose the numerical parameters

The numerical parameters you need to set are:

Parameter	Description
`Lx`	Box length (in each direction), expressed in units of $c/\omega_{pe}$
`nx`	Number of grid cells per direction
`cfl`	CFL number (must be $< 1/\sqrt{2}$ in 2D for the Yee solver)
`T`	Total simulation time, expressed in units of $\omega_{pe}^{-1}$

Hints

The Weibel filaments have a characteristic size of a few $c/\omega_{pe}$. The box should be large enough to contain several filaments.
You need to resolve $c/\omega_{pe}$ on the grid: $\Delta x \lesssim c/\omega_{pe}$.
In 2D, the Yee solver CFL condition is $\mathrm{cfl} < 1/\sqrt{2} \approx 0.7$. Note how the input file already accounts for this in the dt formula.
The Weibel instability grows slower than the two-stream instability. You will need a longer simulation: $\sim 100$–$500\,\omega_{pe}^{-1}$.

Notable details

The geometry is 2D (geometry.dims = 2) with periodic boundary conditions in all directions.
Two electron species (ele1, ele2) are initialized with equal density and temperature but with opposite drift velocities along the $y$-axis.
The diagnostics record the magnetic field components Bx and Bz at regular intervals,
allowing you to visualize the growth of filamentary magnetic structures.
Reduced diagnostics (FieldEnergy, FieldMaximum) track the total field energy and the maximum field amplitude at every timestep. Since the system starts with essentially zero magnetic field (only numerical noise), you can directly observe the exponential growth of the Weibel instability in these quantities.

Run

Create a new folder and copy the input file there, after filling in your chosen parameters.

Let’s run in serial

BASH

warpx.2d input_2d_weibel_instability.txt

Easy! 🕺

If you want to speed things up, you can run in parallel:

Let’s run in parallel

BASH

export OMP_NUM_THREADS=2
mpirun -np 4 <path/to/your/build/bin/warpx.2d> input_2d_weibel_instability.txt

Testimonial

Depending on your choice of resolution and box size, this simulation may take from a few minutes to much longer. If it takes too long, try reducing nx or T.

You should see a standard output flashing out a lot of info.
At the end, you should find in your folder:

a subfolder diags/ with the full field and particle diagnostics
files FieldEnergy.txt and FieldMaximum.txt with the reduced diagnostics, inside diags/reducedfiles/
a file called warpx_used_inputs: a summary of the inputs that were used to run the simulation

If that’s the case, yey! 💯

If the run went wrong, you may find a Backtrace.0.0 file which can be useful for debugging purposes. Let me know if the code fails in any way!

Visualize

With Python 🐍

Here are some things you can visualize:

Magnetic field maps: use get_field('B', 'z') (or 'x') to plot 2D maps of the magnetic field at different times. You should see the formation and growth of filamentary structures.
Field energy vs time: load FieldEnergy.txt and plot the total field energy on a semi-log scale as a function of time. You should see a clear exponential growth phase followed by saturation.
Filament spacing: from the 2D field maps, estimate the characteristic spacing of the current filaments and compare it to the skin depth $c/\omega_{pe}$.

With a Jupyter notebook 📓

The notebook includes a magnetic field map ($B_x$) and the evolution of the electric and magnetic field energy on a semi-log scale.

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

Challenge

Questions for analysis

Describe what you observe in the magnetic field maps and in the field energy evolution. Can you identify distinct stages?
How does the magnetic field energy grow in time? Can you identify a linear (exponential) growth phase on a semi-log plot? Estimate the growth rate $\gamma$ from the slope.
What is the characteristic size of the filaments at saturation? How does it compare to $c/\omega_{pe}$?
How does the growth rate depend on $\beta_0$? Try at least two different drift velocities and compare the field energy time histories.
What happens if you increase the temperature $T_0$ while keeping $\beta_0$ fixed? At what point does the instability shut off?
Compare this simulation to the two-stream instability episode. What are the key differences in the fields that are generated (electrostatic vs electromagnetic)?

Show me the solution

These are open-ended questions meant to guide your exploration. The Weibel growth rate in the cold limit scales as $\gamma \sim \omega_{pe}\,\beta_0$, and the characteristic filament size is of order the skin depth. Finite temperature tends to stabilize short-wavelength modes.

Gallery

Magnetic field structures induced by the Weibel instability.

Evolution of the energy stored in the electric and magnetic fields.

Key Points

💡 The Weibel instability generates magnetic fields from a momentum anisotropy – unlike the two-stream instability, which generates electrostatic fields.

🔬 The instability produces current filaments with a characteristic size of order the plasma skin depth $c/\omega_{pe}$.

📊 Reduced diagnostics (FieldEnergy, FieldMaximum) are an efficient way to monitor the instability growth without storing large field dumps.

⚡ Choosing parameters wisely ($v_d \gg v_{\mathrm{th}}$, sufficient resolution and box size) is essential for observing the instability clearly.

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.

Content from A Magnetic Mirror

Last updated on 2026-04-06 | Edit this page

Overview

Questions

🪞 How to simulate the dynamics of charged particles in an external field?

Objectives

🏃 Run and 👀 visualize some protons in a magnetic mirror!

Setup

In this example we will simulate a bunch of protons inside a magnetic mirror machine. The protons are initialized with random positions and velocities. The magnetic field is loaded from a .h5 file. Make sure to download the input file.

Whenever you need to prepare an input file, this is where you want to go. By the way, analytics tell us that this is the most popular page of the documentation 👠!

Let’s take a look at the input file

OUTPUT

##########################
# USER-DEFINED CONSTANTS #
##########################
my_constants.Lx = 2 # [m]
my_constants.Ly = 2 # [m]
my_constants.Lz = 5 # [m]
my_constants.dt = 4.4e-7 # [s]
my_constants.Np = 1000

############
# NUMERICS #
############
geometry.dims = 3
geometry.prob_hi =  0.5*Lx  0.5*Ly Lz
geometry.prob_lo = -0.5*Lx -0.5*Ly 0
amr.n_cell = 40 40 40
max_step = 500
warpx.const_dt = dt

##############
# ALGORITHMS #
##############
algo.particle_shape = 1
amr.max_level = 0
warpx.do_electrostatic = labframe
warpx.grid_type = collocated
warpx.serialize_initial_conditions = 0
warpx.use_filter = 0

##############
# BOUNDARIES #
##############
boundary.field_hi = pec pec pec
boundary.field_lo = pec pec pec
boundary.particle_hi = absorbing absorbing absorbing
boundary.particle_lo = absorbing absorbing absorbing

#############
# PARTICLES #
#############
particles.species_names = protons
protons.charge = q_e
protons.mass = m_p
protons.do_not_deposit = 1 # test particles
protons.initialize_self_fields = 0
protons.injection_style = gaussian_beam
protons.x_rms = 0.1*Lx
protons.y_rms = 0.1*Ly
protons.z_rms = 0.1*Lz
protons.x_m = 0.
protons.y_m = 0.
protons.z_m = 0.5*Lz
protons.npart = Np
protons.q_tot = q_e*Np
protons.momentum_distribution_type = uniform
protons.ux_min = -9.5e-05
protons.uy_min = -9.5e-05
protons.uz_min = -0.000134
protons.ux_max = 9.5e-05
protons.uy_max = 9.5e-05
protons.uz_max = 0.000134

##########
# FIELDS #
##########
# field here is applied on directly the particles!
particles.B_ext_particle_init_style = read_from_file
particles.read_fields_from_path = example-femm-3d.h5

###############
# DIAGNOSTICS #
###############
diagnostics.diags_names = diag1
diag1.diag_type = Full
diag1.fields_to_plot = Bx By Bz
diag1.format = openpmd
diag1.intervals = 1
diag1.proton.variables = ux uy uz w x y z
diag1.species = protons
diag1.write_species = 1

A few notable details:

The protons are test particles because of the parameter protons.do_not_deposit=1. This means that the protons do not deposit their current density, therefore they do not contribute to the fields.
The magnetic field is applied directly to the particles with the particles.B_ext_particle_init_style flag, so in principle the grid is not used at all. For technical reasons, we must define a grid nonetheless.

Now that we have an idea of what the input files looks like, let’s set up our environment. Activate the warpx environment if you need to. Create a new directory with your own copy of the input file. Also, don’t forget to download the field file and place it in the directory where you will run the input.

Run

Challenge

Let’s run the code

How would you do it? 🤷

Show me the solution

BASH

warpx.3d inputs_3d_magnetic_mirror.txt

As simple as that! 😉

You should see a standard output flashing out a lot of info.
At the end, you should find in your folder:

a subfolder called diags: here is where the code stored the diagnostics
a file called warpx_used_inputs: this is a summary of the inputs that were used to run the simulation

If that’s the case, yey! 💯

If the run went wrong, you may find a Backtrace.0.0 file which can be useful for debugging purposes. Let me know if the code fails in any way!

Here we have loaded the field of hte magnetic bottle from a file. You can also you can define an external field analytically.

Visualize

With Python 🐍

Now that we have the results, we can analyze them using Python.
We will use the openPMD-viewer library to grab the data that the simulation produced in openPMD format. Here you can find a few great tutorials on how to use the viewer. If you feel nerdy and/or you need to deal with the data in parallel workflows, you can use the openPMD-api.

As an example for the magnetic bottle simulation, we have developed simple Jupyter notebook where we retrieve the magnetic field and the particle attributes at the end of the simulation. With a little bit more work, we also plot the trajectories of the particles.

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

With Paraview

Now it’s time to produce some pretty cool images and videos! 😎 If you don’t have it, you can download Paraview here. In the diags/diag1 directory you should find a file named paraview.pmd: Paraview can read .pmd files. Just open Paraview and from there open the .pmd file. You should see Meshes and Particles in your pipeline browser (usually on the left). We can zhuzh up the pipeline so that we can visualize the trajectories of the protons in time

This is the pipeline that I have used to produce the visualizations below.

simulation of proton trajectories inside a magnetic mirror

If you make any other 3D visualization with this data, let me know! We can add it here 😉!

And that’s all for now! 👋

Key Points

💡 The external B field is loaded from an openPMD file, while the protons are defined as test particles.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python or you can open .pmd files directly in Paraview.

Content from A FODO Cell

Last updated on 2026-04-06 | Edit this page

Overview

Questions

How to simulate a particle beam travelling through a FODO cell 🧲🛝🧲?

Objectives

Learn how to setup, run, and visualize a particle beam simulation through a FODO (Focusing-Defocusing) cell with WarpX 🎢

Setup

In this example we will simulate a particle beam travelling through a FODO cell. A FODO cell is a periodic focusing structure used in particle accelerators, consisting of alternating Focusing (F) and Defocusing (D) quadrupole magnets, with O (drift) sections in between. The FODO cell is one of the most fundamental building blocks in accelerator physics, providing transverse focusing to keep particle beams confined as they travel through the accelerator. The alternating focusing and defocusing quadrupoles create a net focusing effect in both transverse planes, allowing the beam to be transported over long distances while maintaining its size.

WarpX for Beam Dynamics

WarpX is a Particle-in-Cell (PIC) code that can track particles embedded in external fields. Two important features are:

Space charge effects: You can turn space charge effects on or off using the appropriate input parameters. When enabled, WarpX solves the electromagnetic (or electrostatic) fields self-consistently from the particle distributions, allowing you to study how the beam’s own fields affect its dynamics.
Arbitrary external fields: WarpX allows you to define arbitrary external electromagnetic fields, either analytically or by loading them from files. This makes it possible to model complex accelerator elements like quadrupoles, dipoles, and other magnets.

Caution

WarpX is not necessarily the best tool for all beam dynamics simulations. If space charge effects are negligible and you’re primarily interested in linear beam optics, there are more specialized and computationally efficient tools available, such as ImpactX.

WarpX vs ImpactX

ImpactX is another code from the same ecosystem, the Beam, Plasma & Accelerator Simulation Toolkit (BLAST), designed for beam dynamics. See the ImpactX documentation and repository. Key differences:

WarpX — t-based (time-based): solves Maxwell/Poisson on a grid, advances particles and fields in time. Best for:
- Strong space charge, wakefields, full EM solutions
- Plasma–beam and other collective effects
- When time evolution of fields matters
ImpactX — z-based (position-based): tracks particles along the beamline through lattice elements. Best for:
- Weak or approximated space charge
- Optics design, linear dynamics, envelope tracking
- Fast runs when detailed field evolution is not needed

For this tutorial we use WarpX to set up a FODO simulation with optional space charge.

Callout

Equivalent ImpactX Example

If you’re interested in seeing how a similar FODO cell simulation can be set up using ImpactX, check out the ImpactX FODO cell example with 2D space charge using envelope tracking. This example demonstrates the envelope tracking approach.

Whenever you need to prepare an input file, this is where you want to go.

Let’s take a look at the input file

OUTPUT

####################
### MY CONSTANTS ###
####################
# BEAM
my_constants.mass          = m_p                              # kg
my_constants.current       = 0.5                              # A
my_constants.energy        = 6.7e6 * q_e                      # J
my_constants.velocity      = sqrt(2 * energy / mass)          # m/s
my_constants.gamma         = energy / (mass * clight**2) + 1. # -
my_constants.beta          = velocity / clight                # -
my_constants.beta_x        = 0.737881                         # m
my_constants.beta_y        = 0.737881                         # m
my_constants.beta_t        = 0.25                             # m
my_constants.alpha_x       = +2.4685083                       # -
my_constants.alpha_y       = -2.4685083                       # -
my_constants.alpha_t       = 0.                               # -
my_constants.gamma_x       = (1+alpha_x**2)/beta_x            # m^-1
my_constants.gamma_y       = (1+alpha_y**2)/beta_y            # m^-1
my_constants.gamma_t       = (1+alpha_t**2)/beta_t            # m^-1
my_constants.emitt_x       = 1.0e-6                           # m
my_constants.emitt_y       = 1.0e-6                           # m
my_constants.emitt_t       = 1.0e-6                           # m
my_constants.sigma_x       = sqrt(emitt_x * beta_x)           # m
my_constants.sigma_y       = sqrt(emitt_y * beta_y)           # m
my_constants.sigma_z       = sqrt(emitt_t * beta_t)           # m
my_constants.sigma_ux      = beta * sqrt(emitt_x * gamma_x)   # -
my_constants.sigma_uy      = beta * sqrt(emitt_y * gamma_y)   # -
my_constants.sigma_uz      = beta * sqrt(emitt_t * gamma_t)   # -
my_constants.sigma_t       = sigma_z / velocity               # s
my_constants.charge        = current *sigma_t                 # C
# BOX
my_constants.Lx   = 20*sigma_x  # m
my_constants.Ly   = 20*sigma_y  # m
my_constants.zmin = -20*sigma_z # m
my_constants.zmax = 0
# LATTICE
my_constants.L_drift1      =  7.44e-2 # m
my_constants.L_drift2      = 14.88e-2 # m
my_constants.L_drift3      =  7.44e-2 # m
my_constants.quad_gradient = 38.64    # T/m
my_constants.L_quad        = 6.10e-2  # m
# TIME
my_constants.T = (L_drift1+L_drift2+L_drift3+2*L_quad)/velocity # s
my_constants.dt = sigma_t/4                                     # s
my_constants.nt = floor(T/dt)                                   # -

##########################
### GENERAL PARAMETERS ###
##########################
stop_time                = T
geometry.dims            = 3
geometry.prob_lo         = -0.5*Lx -0.5*Ly zmin
geometry.prob_hi         = +0.5*Lx +0.5*Ly zmax
amr.n_cell               = 64 64 128
amr.max_level            = 0
warpx.limit_verbose_step = 1
warpx.const_dt           = dt
warpx.do_electrostatic   = labframe
warpx.poisson_solver     = fft
warpx.do_moving_window   = 1
warpx.moving_window_dir  = z
warpx.moving_window_v    = beta  # in units of c
algo.particle_shape      = 3
boundary.field_lo        = open open open
boundary.field_hi        = open open open

#################
### PARTICLES ###
#################
particles.species_names         = beam
beam.species_type               = proton
beam.injection_style            = gaussian_beam
beam.x_rms                      = sigma_x
beam.y_rms                      = sigma_y
beam.z_rms                      = sigma_z
beam.x_m                        = 0
beam.y_m                        = 0
beam.z_m                        = 0.5*zmin
beam.npart                      = 1e3
beam.q_tot                      = charge
beam.momentum_distribution_type = gaussian
beam.uz_m                       = beta*gamma
beam.uy_m                       = 0.0
beam.ux_m                       = 0.0
beam.ux_th                      = sigma_ux
beam.uy_th                      = sigma_uy
beam.uz_th                      = sigma_uz

###############
### LATTICE ###
###############
lattice.elements = drift1 quad1 drift2 quad2 drift3
drift1.type = drift
drift1.ds   = L_drift1
quad1.type  = quad
quad1.ds    = L_quad
quad1.dBdx  = -quad_gradient
drift2.type = drift
drift2.ds   = L_drift2
quad2.type  = quad
quad2.ds    = L_quad
quad2.dBdx  = quad_gradient
drift3.type = drift
drift3.ds   = L_drift3

###################
### DIAGNOSTICS ###
###################
warpx.reduced_diags_names = beam_stats
beam_stats.type = BeamRelevant
beam_stats.species = beam
beam_stats.intervals = 1

diagnostics.diags_names = particles_in
particles_in.intervals = floor(nt/20)
particles_in.diag_type = Full
particles_in.species = beam
particles_in.fields_to_plot = none
particles_in.format = openpmd
particles_in.openpmd_backend = bp
particles_in.dump_last_timestep = 1

Some notable details:

Space charge effects: Space charge is computed self-consistently by solving the Poisson equation at each timestep using the updated charge density from the particle distribution. This allows the beam’s own fields to affect its dynamics.

Turning off space charge: To turn off space charge effects, set algo.maxwell_solver = none and comment out the following inputs:

BASH

# warpx.const_dt           = dt
# warpx.do_electrostatic   = labframe
# warpx.poisson_solver     = fft
# boundary.field_lo        = open open open 
# boundary.field_hi        = open open open

External fields: To initialize an arbitrary external field (such as the quadrupole magnets in a FODO cell), check out the external fields documentation for details on how to define fields analytically or load them from files.

Run

First things first. Create a new folder where you copy the input file. The simulation should be small enough that you can run it in serial with the Conda installation of WarpX.

Let’s run in serial

BASH

warpx.3d inputs_3d_fodo_cell.txt

Just like that! 💃 Note that with Conda’s WarpX you can run this anywhere in your filesystem (provided that you copied there the input of course) because WarpX’s executables are in your $PATH.

If you want to make the simulation faster and/or bigger, then you should run with the parallel version of WarpX. The optimal setup to run the simulation depends on your hardware. This is an example that should work on many common laptops, even though it might not be ideal.

Let’s run in parallel

This is just one way of doing it!

BASH

export OMP_NUM_THREADS=2
mpirun -np 4 <path/to/your/build/bin/warpx.3d> inputs_3d_fodo_cell.txt

Testimonial

With the default parameters, it may take a while on a laptop. If you have access to a GPU, you can experience the speed-up yourself. My experiments: 35 seconds on an NVIDIA A100 GPU, more than 1 hour on 4 cores of my 12th Gen Intel(R) Core(TM) i9-12900H.

Visualize

With Python 🐍

Now that we have the results, we can analyze them using Python.
We will use the openPMD-viewer library to grab the data that the simulation produced in openPMD format. Here you can find a few tutorials on how to use the viewer. If you feel nerdy and/or you need to deal with the data in parallel workflows, you can use the openPMD-api.

As an example for the FODO cell simulation, we have developed a simple Jupyter notebook where we retrieve the beam variables and analyze the beam dynamics through the cell.

Here is a video of the beam travelling through the FODO cell:

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

Key Points

🎯 Particle tracking follows individual particles, while envelope tracking follows the beam envelope - particle tracking is more detailed but computationally expensive.

🔬 A FODO cell is a periodic focusing structure with alternating focusing and defocusing quadrupole magnets.

⚡ WarpX allows you to turn space charge effects on/off and define arbitrary external fields.

🚀 WarpX is a full PIC code best suited for strong space charge effects and complex electromagnetic interactions. For cases with negligible space charge, ImpactX may be a more efficient choice.

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.

Content from A Laser-Driven Ion Accelerator

Last updated on 2026-04-08 | Edit this page

Overview

Questions

How does an intense laser pulse accelerate ions from a thin solid target?

Objectives

Set up, run, and visualize a 2D Target Normal Sheath Acceleration (TNSA) simulation with WarpX. Understand the role of key laser and target parameters in determining the maximum ion energy.

Introduction

Target Normal Sheath Acceleration (TNSA) is one of the most studied mechanisms for laser-driven ion acceleration. The physics goes roughly as follows:

An intense, short-pulse laser ($I > 10^{18}$ W/cm$^2$) hits a thin solid-density target.
The laser heats electrons at the front surface to relativistic energies.
These hot electrons traverse the target and escape from the rear surface, setting up a strong electrostatic sheath field ($\sim$ TV/m) at the target-vacuum interface.
This sheath field ionizes and accelerates ions (typically protons from surface contaminants) to multi-MeV energies in the direction normal to the target surface.

In this episode, we simulate this process in 2D using a laser pulse impinging on a composite target: a boron slab (the main target) followed by a thin hydrogen contaminant layer from which protons are accelerated.

Caution

This simulation is significantly more demanding than the instability examples. Depending on the resolution you choose, you may need a GPU or a parallel run to complete it in a reasonable time.

Setup

Make sure to download the input file.

Whenever you need to prepare an input file, this is where you want to go.

Let’s take a look at the input file

OUTPUT

####################
### MY CONSTANTS ###
####################
my_constants.micro = 1e-6
my_constants.femto = 1e-15
my_constants.eV = q_e
# LASER
my_constants.FWHM_I = ...*femto # [s]
my_constants.laser_tpeak = 2*FWHM_I # [s]
my_constants.laser_waist = ...*micro # [m]
my_constants.laser_wavelength = ...*micro # [m]
my_constants.a0 = ... # [-]
# PLASMA
my_constants.n_c = epsilon0 * m_e * (2*pi*clight)**2 /(laser_wavelength*q_e)**2 # [m^-3]
my_constants.ne0 = ...*n_c # [m^-3]
my_constants.Te0 = ...*eV # [J]
my_constants.zmin_targ = 2*FWHM_I*clight # [m]
my_constants.thick_targ = ...*micro # [m]
my_constants.thick_cont = ...*micro # [m]
# BOX
my_constants.Lz = 4*FWHM_I*clight
my_constants.Lx = 10*laser_waist
my_constants.nz = ...
my_constants.nx = ...
my_constants.dx = Lx / nx
my_constants.dz = Lz / nz
# TIME
my_constants.cfl = ...
my_constants.T = 2*Lz/clight
my_constants.dt = cfl*dx/(sqrt(2)*clight)
my_constants.nt = floor(T/dt)

##########################
### GENERAL PARAMETERS ###
##########################
stop_time = T
amr.n_cell = nx nz
amr.max_level = 0
geometry.dims = 2
geometry.prob_lo = -0.5*Lx  0
geometry.prob_hi =  0.5*Lx  Lz

##########################
### BOUNDARY CONDITION ###
##########################
boundary.field_lo = pml pml
boundary.field_hi = pml pml
boundary.particle_lo = absorbing absorbing absorbing
boundary.particle_hi = absorbing absorbing absorbing

###############
### NUMERICS ###
################
algo.particle_shape = 3
algo.maxwell_solver = yee
algo.particle_pusher = boris
algo.current_deposition = esirkepov
warpx.cfl = cfl
warpx.use_filter = 1

#################
### PARTICLES ###
#################
particles.species_names = ion_targ ele_targ ion_cont ele_cont

ion_targ.species_type = boron
ion_targ.injection_style = NRandomPerCell
ion_targ.num_particles_per_cell = 20
ion_targ.momentum_distribution_type = maxwell_boltzmann
ion_targ.theta = Te0 / (10*m_p*clight**2)
ion_targ.zmin = zmin_targ
ion_targ.zmax = zmin_targ + thick_targ
ion_targ.profile = constant
ion_targ.density = ne0/5

ele_targ.species_type = electron
ele_targ.injection_style = NRandomPerCell
ele_targ.num_particles_per_cell = 40
ele_targ.momentum_distribution_type = maxwell_boltzmann
ele_targ.theta = Te0 / (m_e*clight**2)
ele_targ.zmin = zmin_targ
ele_targ.zmax = zmin_targ + thick_targ
ele_targ.profile = constant
ele_targ.density = ne0

ion_cont.species_type = proton
ion_cont.injection_style = NRandomPerCell
ion_cont.num_particles_per_cell = 100
ion_cont.momentum_distribution_type = maxwell_boltzmann
ion_cont.theta = Te0 / (m_p*clight**2)
ion_cont.zmin = zmin_targ + thick_targ
ion_cont.zmax =  zmin_targ + thick_targ + thick_cont
ion_cont.profile = constant
ion_cont.density = 5*n_c

ele_cont.species_type = electron
ele_cont.injection_style = NRandomPerCell
ele_cont.num_particles_per_cell = 50
ele_cont.momentum_distribution_type = maxwell_boltzmann
ele_cont.theta = Te0 / (m_e*clight**2)
ele_cont.zmin = zmin_targ + thick_targ
ele_cont.zmax = zmin_targ + thick_targ + thick_cont
ele_cont.profile = constant
ele_cont.density = 5*n_c

#############
### LASER ###
#############
lasers.names        = laser1
laser1.position     = 0 0 dz
laser1.direction    = 0. 0. 1.
laser1.polarization = 1. 0. 0.
laser1.a0           = a0
laser1.wavelength   = laser_wavelength
laser1.profile      = Gaussian
laser1.profile_waist = laser_waist
laser1.profile_duration = FWHM_I/1.17741
laser1.profile_t_peak = laser_tpeak
laser1.profile_focal_distance = zmin_targ

################
#### DIAGNOSTICS
################
# FULL
diagnostics.diags_names = fields particles
fields.diag_type = Full
fields.fields_to_plot = Ex Ez rho_ele_targ rho_ion_targ rho_ele_cont rho_ion_cont
fields.format = openpmd
fields.intervals = floor(nt/200)
fields.openpmd_backend = bp
fields.write_species = 0
particles.diag_type = Full
particles.format = openpmd
particles.openpmd_backend = bp
particles.species = ele_targ ion_cont
particles.fields_to_plot = none
particles.intervals = floor(nt/200)

# REDUCED
warpx.reduced_diags_names = FieldEnergy FieldMaximum
FieldEnergy.type = FieldEnergy
FieldEnergy.intervals = 1
FieldMaximum.type = FieldMaximum
FieldMaximum.intervals = 1

Choosing the parameters

You will notice that some parameters in the input file are set to ... (ellipsis). These are left for you to choose! This is the most complex input of the tutorial, and the parameter choices will determine whether you get meaningful physics out of the simulation.

Challenge

Choose the laser parameters

Parameter	Symbol	Description
`FWHM_I`	$\tau_{\mathrm{FWHM}}$	Pulse duration (FWHM of intensity) $[\mathrm{s}]$
`laser_waist`	$w_0$	Laser spot size (beam waist) $[\mathrm{m}]$
`laser_wavelength`	$\lambda$	Laser wavelength $[\mathrm{m}]$
`a0`	$a_0$	Normalized vector potential (dimensionless laser amplitude)

Hints

A Ti:Sapphire laser has $\lambda \approx 0.8\,\mu\mathrm{m}$. This is the most common choice.
Typical ultrashort pulses have durations $\tau_{\mathrm{FWHM}} \sim 20$–$100$ fs.
The normalized vector potential $a_0$ controls the laser intensity: $I \approx 1.37 \times 10^{18}\, (a_0 / \lambda[\mu\mathrm{m}])^2$ W/cm$^2$. For TNSA, $a_0 \gtrsim 1$ (i.e. relativistic intensity) is required. Values of $a_0 \sim 1$–$10$ are typical.
The laser waist $w_0$ is typically a few micrometers ($\sim 2$–$10\,\mu\mathrm{m}$).

Challenge

Choose the target parameters

Parameter	Symbol	Description
`ne0`	$n_e / n_c$	Electron density of the main target, in units of the critical density $n_c$
`Te0`	$T_e$	Initial electron temperature $[\mathrm{eV}]$
`thick_targ`		Thickness of the main (boron) target $[\mathrm{m}]$
`thick_cont`		Thickness of the hydrogen contaminant layer $[\mathrm{m}]$

Hints

Solid-density targets have $n_e \sim 100$–$500\, n_c$, where $n_c$ is the critical density for the chosen wavelength. However, running at full solid density is very expensive. You may want to start with a reduced density $\sim 10$–$50\, n_c$ to keep the cost manageable.
The initial temperature is typically a few eV (room temperature ionized material), but the exact value does not matter much since the laser will heat the electrons far beyond this.
The target thickness is usually a few micrometers ($\sim 1$–$10\,\mu\mathrm{m}$). Thinner targets lead to higher maximum proton energies (up to a point).
The contaminant layer is very thin ($\sim 0.1$–$1\,\mu\mathrm{m}$).

Challenge

Choose the numerical parameters

Parameter	Description
`nz`	Number of grid cells in the laser propagation direction
`nx`	Number of grid cells in the transverse direction
`cfl`	CFL number (must be $< 1/\sqrt{2}$ in 2D)

Hints

You need to resolve the laser wavelength: $\Delta z \lesssim \lambda / 20$ is a common rule of thumb, with $\lambda / 30$ or finer being safer.
The transverse resolution can be somewhat coarser than the longitudinal, but should still resolve $c/\omega_{pe}$ for the target density.
A CFL of $\sim 0.7$ (i.e. $1/\sqrt{2}$) is the maximum allowed.
For a first test, use a coarser grid to make sure the simulation runs correctly, and then refine.

Notable details

The target consists of four species: boron ions and electrons for the main target (ion_targ, ele_targ), and protons and electrons for the contaminant layer (ion_cont, ele_cont). The boron is initialized with charge state $Z=5$, so the electron density is $5\times$ the ion density.
The laser is injected from the $z_{\mathrm{min}}$ boundary using the built-in Gaussian laser profile, with the focus placed at the front surface of the target.
PML (Perfectly Matched Layer) boundary conditions are used for the fields to absorb outgoing electromagnetic waves, and absorbing boundary conditions are used for particles leaving the box.
The current_deposition = esirkepov option is used to ensure charge conservation, which is essential for a simulation involving dense plasma and strong fields.
The diagnostics write out both field data (Ex, Ez, charge densities of each species) and particle data (target electrons and contaminant protons).

Run

Create a new folder and copy the input file there, after filling in your chosen parameters.

Caution

This simulation is more expensive than the instability examples. Depending on your resolution, running in serial may take a very long time. Consider running in parallel or on a GPU.

Let’s run in serial

BASH

warpx.2d input_2d_tnsa.txt

Let’s run in parallel

BASH

export OMP_NUM_THREADS=2
mpirun -np 4 <path/to/your/build/bin/warpx.2d> input_2d_tnsa.txt

You should see a standard output flashing out a lot of info.
At the end, you should find in your folder:

a subfolder diags/ with the full field and particle diagnostics
files FieldEnergy.txt and FieldMaximum.txt with the reduced diagnostics, inside diags/reducedfiles/
a file called warpx_used_inputs: a summary of the inputs that were used to run the simulation

If that’s the case, yey! 💯

If the run went wrong, you may find a Backtrace.0.0 file which can be useful for debugging purposes. Let me know if the code fails in any way!

Visualize

With Python 🐍

Here are some things you can visualize:

Electric field maps: plot Ex and Ez at several time snapshots to see the laser pulse entering the box, interacting with the target, and the sheath field forming at the rear surface.
Charge density maps: plot the charge densities (rho_ele_targ, rho_ion_cont, etc.) to see the electron heating, expansion, and proton acceleration.
Proton energy spectrum: from the contaminant proton data (ion_cont), compute the kinetic energy of each proton and build a histogram. The maximum proton energy is the key figure of merit in TNSA experiments.
Field energy vs time: use FieldEnergy.txt to track how the electromagnetic energy in the box evolves as the laser enters, interacts, and exits.

With a Jupyter notebook 📓

The notebook includes a field map with the laser electric field overlaid on the target electron density, and a phase-space plot of the accelerated contaminant protons.

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

Challenge

Questions for analysis

Describe what you observe in the field maps and charge density plots at different times. Can you identify the laser entering, the electron heating, and the proton acceleration?
What is the maximum proton energy you observe? How does it compare to the theoretical TNSA scaling with laser intensity $a_0$?
How does the sheath field $E_z$ at the target rear surface evolve in time? Can you identify the moment when the sheath field is strongest?
Try running the simulation with a thinner (or thicker) target. How does the maximum proton energy change? Why?
What happens if you increase $a_0$ (i.e., increase the laser intensity)? How does the proton energy spectrum change?
Look at the electron density behind the target. Can you see the hot electron population that escapes and creates the sheath?
What role does the contaminant layer play? What would happen if you made it thicker, or changed the contaminant species to something heavier (e.g., carbon)?

Show me the solution

These are open-ended questions meant to guide your exploration. A classical TNSA scaling predicts the maximum proton energy $E_{\max} \propto a_0$ (or $\propto \sqrt{I}$), though the exact scaling depends on target thickness, density, and pulse duration. Thinner targets generally yield higher maximum energies because the hot electrons recirculate more efficiently.

Gallery

2D map showing the laser electric field interacting with the solid-density target and accelerated protons

2D map showing the phase space of the accelerated protons at the end of the simulation

Key Points

💡 TNSA accelerates ions via a strong electrostatic sheath field created by hot electrons escaping the rear surface of a laser-irradiated thin target.

🔬 Key parameters controlling the maximum ion energy are the laser intensity ($a_0$), the target thickness, and the target density.

⚙️ Charge-conserving current deposition (esirkepov) and proper boundary conditions (PML, absorbing) are essential for this type of simulation.

⚡ This is a computationally expensive simulation – be prepared to iterate on resolution and target density to find a balance between physical fidelity and computational cost.

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.

Content from A Beam-Beam Collision

Last updated on 2026-04-06 | Edit this page

Overview

Questions

How to use WarpX for beam-beam simulations of colliders 💥?

Objectives

Learn how to setup, run, and visualize your own beam-beam ☄️☄️ simulation

Setup

In this example we will simulate a bunch of electrons colliding against a bunch of positrons. We have selected the parameters of the C$^3$ linear collider.

Whenever you need to prepare an input file, this is where you want to go.

Let’s take a look at the input file

OUTPUT

#################################
########## MY CONSTANTS #########
#################################
my_constants.mc2   = m_e*clight*clight
my_constants.nano  = 1.0e-9
my_constants.micro = 1.e-6
my_constants.milli = 1.e-3
my_constants.GeV   = q_e*1.e9

# BEAMS
my_constants.beam_energy = 125*GeV
my_constants.beam_gamma  = beam_energy/mc2
my_constants.beam_npart  = 6.24e9
my_constants.nmacropart  = 1e5
my_constants.beam_charge = q_e * beam_npart
#my_constants.sigmax     = 210.0*nano
#my_constants.sigmay     = 3.1*nano
my_constants.sigmaz      = 100.*micro
my_constants.beam_uth    = 0.3/100*beam_gamma
my_constants.mux         = 0.*sigmax
my_constants.muy         = 0.*sigmay
my_constants.muz         = 4*sigmaz
my_constants.emitx       = 900*nano
my_constants.emity       = 20*nano
my_constants.dux         = emitx / sigmax
my_constants.duy         = emity / sigmay
my_constants.betax       = 12*milli
my_constants.betay       = 0.12*milli
my_constants.sigmax      = sqrt( emitx * betax / beam_gamma )
my_constants.sigmay      = sqrt( emity * betay / beam_gamma )

# BOX
my_constants.Lx = 20*sigmax
my_constants.Ly = 20*sigmay
my_constants.Lz = 16*sigmaz
my_constants.nx = 128
my_constants.ny = 128
my_constants.nz = 64
my_constants.dx = Lx/nx
my_constants.dy = Ly/ny
my_constants.dz = Lz/nz

# TIME
my_constants.T = 0.5*Lz/clight
my_constants.dt = T / nz
my_constants.nt = floor(T/dt)

# LUMI DIAG
my_constants.bin_num_1d = 2048
my_constants.bin_center_min_1d = 0.
my_constants.bin_center_max_1d = 2*beam_energy/q_e
my_constants.bin_size_1d = (bin_center_max_1d - bin_center_min_1d) / bin_num_1d
my_constants.bin_edge_min_1d = bin_center_min_1d - 0.5 * bin_size_1d
my_constants.bin_edge_max_1d = bin_center_max_1d + 0.5 * bin_size_1d
my_constants.bin_num_1d_eff = bin_num_1d + 1

#################################
####### GENERAL PARAMETERS ######
#################################
stop_time = T
amr.n_cell = nx ny nz
amr.max_level = 0
geometry.dims = 3
geometry.prob_lo = -0.5*Lx -0.5*Ly -0.5*Lz
geometry.prob_hi =  0.5*Lx  0.5*Ly  0.5*Lz

#################################
######## BOUNDARY CONDITION #####
#################################
boundary.field_lo = open open open
boundary.field_hi = open open open
boundary.particle_lo = Absorbing Absorbing Absorbing
boundary.particle_hi = Absorbing Absorbing Absorbing

#################################
############ NUMERICS ###########
#################################
warpx.do_electrostatic = relativistic
warpx.const_dt = dt
warpx.grid_type = collocated
algo.particle_shape = 3
algo.particle_pusher = vay
warpx.poisson_solver = fft
warpx.use_2d_slices_fft_solver = 1

#################################
########### PARTICLES ###########
#################################
particles.species_names = ele1 pos2 pho1 pho2 ele2 pos1
particles.photon_species = pho1 pho2

ele1.species_type = electron
ele1.injection_style = gaussian_beam
ele1.x_rms = sigmax
ele1.y_rms = sigmay
ele1.z_rms = sigmaz
ele1.x_m = - mux
ele1.y_m = - muy
ele1.z_m = - muz
ele1.npart = nmacropart
ele1.q_tot = -beam_charge
ele1.z_cut = 4
ele1.focal_distance = muz
ele1.momentum_distribution_type = gaussian
ele1.uz_m = beam_gamma
ele1.uy_m = 0.0
ele1.ux_m = 0.0
ele1.ux_th = dux
ele1.uy_th = duy
ele1.uz_th = beam_uth
ele1.initialize_self_fields = 1
ele1.do_qed_quantum_sync = 1
ele1.qed_quantum_sync_phot_product_species = pho1
ele1.do_classical_radiation_reaction = 0

pos2.species_type = positron
pos2.injection_style = gaussian_beam
pos2.x_rms = sigmax
pos2.y_rms = sigmay
pos2.z_rms = sigmaz
pos2.x_m = mux
pos2.y_m = muy
pos2.z_m = muz
pos2.npart = nmacropart
pos2.q_tot = beam_charge
pos2.z_cut = 4
pos2.focal_distance = muz
pos2.momentum_distribution_type = gaussian
pos2.uz_m = -beam_gamma
pos2.uy_m = 0.0
pos2.ux_m = 0.0
pos2.ux_th = dux
pos2.uy_th = duy
pos2.uz_th = beam_uth
pos2.initialize_self_fields = 1
pos2.do_qed_quantum_sync = 1
pos2.qed_quantum_sync_phot_product_species = pho2
pos2.do_classical_radiation_reaction = 0

pho1.species_type = photon
pho1.injection_style = none
pho1.do_qed_breit_wheeler = 1
pho1.qed_breit_wheeler_ele_product_species = ele1
pho1.qed_breit_wheeler_pos_product_species = pos1

pho2.species_type = photon
pho2.injection_style = none
pho2.do_qed_breit_wheeler = 1
pho2.qed_breit_wheeler_ele_product_species = ele2
pho2.qed_breit_wheeler_pos_product_species = pos2

ele2.species_type = electron
ele2.injection_style = none
ele2.do_qed_quantum_sync = 1
ele2.qed_quantum_sync_phot_product_species = pho2
ele2.do_classical_radiation_reaction = 0

pos1.species_type = positron
pos1.injection_style = none
pos1.do_qed_quantum_sync = 1
pos1.qed_quantum_sync_phot_product_species = pho1
pos1.do_classical_radiation_reaction = 0

#################################
############# QED ###############
#################################
qed_qs.photon_creation_energy_threshold = 0.

qed_qs.lookup_table_mode = builtin
qed_qs.chi_min = 1.e-7

qed_bw.lookup_table_mode = builtin
qed_bw.chi_min = 1.e-2

warpx.do_qed_schwinger = 0.

#################################
######### DIAGNOSTICS ###########
#################################
# FULL
diagnostics.diags_names = particles_out particles_in

# particles that exit the simulation domain at any given time
particles_out.dump_last_timestep = 1
particles_out.diag_type = BoundaryScraping
particles_out.format = openpmd
particles_out.openpmd_backend = h5
particles_out.intervals = -1

ele1.save_particles_at_xlo = 1
ele1.save_particles_at_ylo = 1
ele1.save_particles_at_zlo = 1
ele1.save_particles_at_xhi = 1
ele1.save_particles_at_yhi = 1
ele1.save_particles_at_zhi = 1

pos2.save_particles_at_xlo = 1
pos2.save_particles_at_ylo = 1
pos2.save_particles_at_zlo = 1
pos2.save_particles_at_xhi = 1
pos2.save_particles_at_yhi = 1
pos2.save_particles_at_zhi = 1

pho1.save_particles_at_xlo = 1
pho1.save_particles_at_ylo = 1
pho1.save_particles_at_zlo = 1
pho1.save_particles_at_xhi = 1
pho1.save_particles_at_yhi = 1
pho1.save_particles_at_zhi = 1

pho2.save_particles_at_xlo = 1
pho2.save_particles_at_ylo = 1
pho2.save_particles_at_zlo = 1
pho2.save_particles_at_xhi = 1
pho2.save_particles_at_yhi = 1
pho2.save_particles_at_zhi = 1

ele2.save_particles_at_xlo = 1
ele2.save_particles_at_ylo = 1
ele2.save_particles_at_zlo = 1
ele2.save_particles_at_xhi = 1
ele2.save_particles_at_yhi = 1
ele2.save_particles_at_zhi = 1

pos1.save_particles_at_xlo = 1
pos1.save_particles_at_ylo = 1
pos1.save_particles_at_zlo = 1
pos1.save_particles_at_xhi = 1
pos1.save_particles_at_yhi = 1
pos1.save_particles_at_zhi = 1

# particles inside the simulation domain
particles_in.intervals = 1
particles_in.diag_type = Full
particles_in.species = ele1 pos2
particles_in.fields_to_plot = none
particles_in.format = openpmd
particles_in.openpmd_backend = h5
particles_in.dump_last_timestep = 1

# REDUCED
warpx.reduced_diags_names = DiffLumi_ele1_pos2 DiffLumi_ele1_ele2 DiffLumi_pos1_pos2 DiffLumi_pos1_ele2 DiffLumi_pos1_pho2 DiffLumi_ele1_pho2 DiffLumi_pho1_pos2 DiffLumi_pho1_ele2 DiffLumi_pho1_pho2 CollRel_ele1_pos2

CollRel_ele1_pos2.type = ColliderRelevant
CollRel_ele1_pos2.species = ele1 pos2
CollRel_ele1_pos2.intervals = 1

DiffLumi_ele1_pos2.type = DifferentialLuminosity
DiffLumi_ele1_pos2.intervals = nt
DiffLumi_ele1_pos2.species = ele1 pos2
DiffLumi_ele1_pos2.bin_number = bin_num_1d_eff
DiffLumi_ele1_pos2.bin_max = bin_edge_max_1d
DiffLumi_ele1_pos2.bin_min = bin_edge_min_1d

DiffLumi_ele1_ele2.type = DifferentialLuminosity
DiffLumi_ele1_ele2.intervals = nt
DiffLumi_ele1_ele2.species = ele1 ele2
DiffLumi_ele1_ele2.bin_number = bin_num_1d_eff
DiffLumi_ele1_ele2.bin_max = bin_edge_max_1d
DiffLumi_ele1_ele2.bin_min = bin_edge_min_1d

DiffLumi_pos1_pos2.type = DifferentialLuminosity
DiffLumi_pos1_pos2.intervals = nt
DiffLumi_pos1_pos2.species = pos1 pos2
DiffLumi_pos1_pos2.bin_number = bin_num_1d_eff
DiffLumi_pos1_pos2.bin_max = bin_edge_max_1d
DiffLumi_pos1_pos2.bin_min = bin_edge_min_1d

DiffLumi_pos1_ele2.type = DifferentialLuminosity
DiffLumi_pos1_ele2.intervals = nt
DiffLumi_pos1_ele2.species = pos1 ele2
DiffLumi_pos1_ele2.bin_number = bin_num_1d_eff
DiffLumi_pos1_ele2.bin_max = bin_edge_max_1d
DiffLumi_pos1_ele2.bin_min = bin_edge_min_1d

DiffLumi_ele1_pho2.type = DifferentialLuminosity
DiffLumi_ele1_pho2.intervals = nt
DiffLumi_ele1_pho2.species = ele1 pho2
DiffLumi_ele1_pho2.bin_number = bin_num_1d_eff
DiffLumi_ele1_pho2.bin_max = bin_edge_max_1d
DiffLumi_ele1_pho2.bin_min = bin_edge_min_1d

DiffLumi_pos1_pho2.type = DifferentialLuminosity
DiffLumi_pos1_pho2.intervals = nt
DiffLumi_pos1_pho2.species = pos1 pho2
DiffLumi_pos1_pho2.bin_number = bin_num_1d_eff
DiffLumi_pos1_pho2.bin_max = bin_edge_max_1d
DiffLumi_pos1_pho2.bin_min = bin_edge_min_1d

DiffLumi_pho1_pos2.type = DifferentialLuminosity
DiffLumi_pho1_pos2.intervals = nt
DiffLumi_pho1_pos2.species = pho1 pos2
DiffLumi_pho1_pos2.bin_number = bin_num_1d_eff
DiffLumi_pho1_pos2.bin_max = bin_edge_max_1d
DiffLumi_pho1_pos2.bin_min = bin_edge_min_1d

DiffLumi_pho1_ele2.type = DifferentialLuminosity
DiffLumi_pho1_ele2.intervals = nt
DiffLumi_pho1_ele2.species = pho1 ele2
DiffLumi_pho1_ele2.bin_number = bin_num_1d_eff
DiffLumi_pho1_ele2.bin_max = bin_edge_max_1d
DiffLumi_pho1_ele2.bin_min = bin_edge_min_1d

DiffLumi_pho1_pho2.type = DifferentialLuminosity
DiffLumi_pho1_pho2.intervals = nt
DiffLumi_pho1_pho2.species = pho1 pho2
DiffLumi_pho1_pho2.bin_number = bin_num_1d_eff
DiffLumi_pho1_pho2.bin_max = bin_edge_max_1d
DiffLumi_pho1_pho2.bin_min = bin_edge_min_1d

Some notable details:

Poisson solver: because the beams are ultra-relativistic (125 GeV electrons and positrons) and flat, we use an FFT-based electrostatic solver. Specifically, it solves many 2D Poisson equations in the $(x,y)$ plane for each $z$. The full 3D version of this solver is also available with warpx.use_2d_slices_fft_solver = 0.
Resolution: the number of grid cells is reduced to fit in a laptop. For production simulations, make sure you increase the resolution.
Timestep: since the beams travel at the speed of light along $z$ and the simulation frame is the center of mass frame, it makes sense to choose $dt = dz / ( 2 c ) $. However this is not strictly necessary. Sometimes it can be useful to save resources by choosing a larger timestep. Just make sure you resolve ‘’well-enough’’ the shortest timescale that you’re interested in.
QED lookup tables: here we use the built-in ones for convenience. For production runs, make sure to use tables with enough points and set the ranges of the $\chi$ parameter to what you need.
Diagnostics:
- the trajectories of the beam particles. This diagnostic can easily take up too much memory. For simulations with many macroparticles, consider using a field diagnostic.
- all the particles that exit the domain (BoundaryScraping)
- the differential luminosity of every pair of left-ward and right-ward moving species

Run

First things first. Create a new folder where you copy the input file. The simulation in small enough that you should be able to run it in serial with the Conda installation of WarpX.

Let’s run in serial

BASH

warpx.3d inputs_3d_beambeam_C3.txt

Just like that! 💃 Note that with Conda’s WarpX you can run this anywhere in your filesystem (provided that you copied there the input of course) because WarpX’s executable are in your $PATH.

Let’s run in parallel

This is just one way of doing it!

BASH

export OMP_NUM_THREADS=2
mpirun -np 4 <path/to/your/build/bin/warpx.3d> inputs_3d_beambeam_C3.txt

Testimonial

On my laptop’s CPU (12th Gen Intel® Core™ i9-12900H × 20) the serial simulation took ~195 seconds, while the parallel one ~44 seconds on 4 cores!

Visualize

With Python 🐍

As an example for the beam-beam simulation, we have developed simple Jupyter notebook where we retrieve the beams’ particles positions and project them on the $(z,x)$ and $(z,y)$ planes.

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

With Paraview

Caution

Coming soon!

Key Points

💅 There are several details one needs to take care when setting up a beam-beam simulation

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python.

Content from OSSFE 2025 - Using WarpX, a general purpose particle-in-cell code

Last updated on 2026-04-08 | Edit this page

Overview

Questions

🤌 What is WarpX?
🔧 How can I install and run WarpX?
🕵️ How can I analyze the simulation results?

Objectives

💻 Install WarpX on your local machine with Conda
🏃 Run a fusion-relevant example on your local machine: protons in a magnetic mirror!
👀 Visualize the results with Python and Paraview

WarpX, a particle-in-cell code

Welcome to the WarpX tutorial at OSSFE 2025 (March 18, 2025)! 👋

WarpX is a general purpose open-source high-performance Particle-In-Cell (PIC) code.
If you are not familiar with the PIC method, here is a picture that condenses the core idea:

And here is a more informative image that explains the core algorithmic steps.

If you want to know more about PIC, here are a few references:

The two bibles on PIC 📚
- C. K. Birdsall and A. B. Langdon. Plasma Physics Via Computer Simulation.
- R. W. Hockney and J. W. Eastwood. Computer simulation using particles.
An old review written by one of the pioneers: John M. Dawson, Particle simulation of plasmas, Rev. Mod. Phys. 55, 403
Browse our docs for many more references about advanced algorithms and methods.

In this tutorial we will go through the basics of WarpX: installation, running a simple example and visualizing the results. Along the way, we will point to some specific locations in the documentation, for your reference.

Callout

📣 Everything you need to know to use WarpX is in the documentation, check it out!

Checklist

Some cool features of WarpX:

📖 Open-source - we wouldn’t be here otherwise!

✈️ Runs on GPUs: NVIDIA, AMD, and Intel

🚀 Runs on multiple GPUs or CPUs, on systems ranging from laptops to supercomputers

🤓 Many many advanced algorithms and methods: mesh-refinement, embedded boundaries, electrostatic/electromagnetic/pseudospectral solvers, etc.

💾 Standards: openPMD for input/output data, PICMI for inputs

🤸 Active development and mainteinance: check our GitHub repo

🗺️ International, cross-disciplinary community: plasma physics, fusion devices, laser-plasma interactions, beam physics, plasma-based acceleration, astrophysics

Installing WarpX using Conda-Forge

Callout

BASH

conda install -c conda-forge warpx

Ok, maybe two lines if you want to keep your system clean by creating a new environment.

BASH

conda create -n warpx -c conda-forge warpx 
conda activate warpx

Now you should have 4 different WarpX binaries in your PATH called warpx.1d, warpx.2d, warpx.3d, warpx.rz.
Each binary for a different dimensionality.

To check this, run:

BASH

which warpx.1d warpx.2d warpx.3d warpx.rz

If you get 3 different paths that look something like:

BASH

/home/<username>/anaconda3/envs/warpx/bin/warpx.xd

then you got this 🙌! You can also import pywarpx in Python.

A simple example of a magnetic mirror

Whenever you need to prepare an input file, this is where you want to go. By the way, analytics tell us that this is the most popular page of the documentation 👠!

Let’s take a look at the input file

OUTPUT

##########################
# USER-DEFINED CONSTANTS #
##########################
my_constants.Lx = 2 # [m]
my_constants.Ly = 2 # [m]
my_constants.Lz = 5 # [m]
my_constants.dt = 4.4e-7 # [s]
my_constants.Np = 1000

############
# NUMERICS #
############
geometry.dims = 3
geometry.prob_hi =  0.5*Lx  0.5*Ly Lz
geometry.prob_lo = -0.5*Lx -0.5*Ly 0
amr.n_cell = 40 40 40
max_step = 500
warpx.const_dt = dt

##############
# ALGORITHMS #
##############
algo.particle_shape = 1
amr.max_level = 0
warpx.do_electrostatic = labframe
warpx.grid_type = collocated
warpx.serialize_initial_conditions = 0
warpx.use_filter = 0

##############
# BOUNDARIES #
##############
boundary.field_hi = pec pec pec
boundary.field_lo = pec pec pec
boundary.particle_hi = absorbing absorbing absorbing
boundary.particle_lo = absorbing absorbing absorbing

#############
# PARTICLES #
#############
particles.species_names = protons
protons.charge = q_e
protons.mass = m_p
protons.do_not_deposit = 1 # test particles
protons.initialize_self_fields = 0
protons.injection_style = gaussian_beam
protons.x_rms = 0.1*Lx
protons.y_rms = 0.1*Ly
protons.z_rms = 0.1*Lz
protons.x_m = 0.
protons.y_m = 0.
protons.z_m = 0.5*Lz
protons.npart = Np
protons.q_tot = q_e*Np
protons.momentum_distribution_type = uniform
protons.ux_min = -9.5e-05
protons.uy_min = -9.5e-05
protons.uz_min = -0.000134
protons.ux_max = 9.5e-05
protons.uy_max = 9.5e-05
protons.uz_max = 0.000134

##########
# FIELDS #
##########
# field here is applied on directly the particles!
particles.B_ext_particle_init_style = read_from_file
particles.read_fields_from_path = example-femm-3d.h5

###############
# DIAGNOSTICS #
###############
diagnostics.diags_names = diag1
diag1.diag_type = Full
diag1.fields_to_plot = Bx By Bz
diag1.format = openpmd
diag1.intervals = 1
diag1.proton.variables = ux uy uz w x y z
diag1.species = protons
diag1.write_species = 1

Challenge

Let’s run the code

How would you do it? 🤷

Show me the solution

BASH

warpx.3d inputs_3d_magnetic_mirror.txt

You should see a standard output flashing out a lot of info.
At the end, you should find in your folder:

a subfolder called diags: here is where the code stored the diagnostics
a file called warpx_used_inputs: this is a summary of the inputs that were used to run the simulation

If that’s the case, yey! 💯

If the run went wrong, you may find a Backtrace.0.0 file which can be useful for debugging purposes. Let me know if the code fails in any way!

Here we have loaded the field of hte magnetic bottle from a file. You can also you can define an external field analytically.

Data handling and visualizations

With Python 🐍

Let’s take a look at the Jupyter notebook

You can download the notebook and try it yourself. Remember to either run the notebook from the simulation directory or change the corresponding path in the notebook.

With Paraview

This is the pipeline that I have used to produce the visualizations below.

If you make any other 3D visualization with this data, let me know! We can add it here 😉!

And that’s all for now! 👋

Key Points

🚀 WarpX is a open-source high-performance particle-in-cell code

🎯 WarpX is easy to install via Conda: conda -c conda-forge warpx

🔍 The documentation is the first place to look for answers, otherwise check out our issues and discussions and ask there.

📷 To analyze and visualize the simulation results in openPMD format, you can use the openPMD-viewer library for Python or you can open .pmd files directly in Paraview.

Content from UCB 2026 -- Two-Stream Instability Live Tutorial

Last updated on 2026-04-08 | Edit this page

Overview

Questions

Where can I find the material from the live WarpX tutorial (April 8, 2026)?

Objectives

Review the steps we followed during the live tutorial session. Access the input files, notebooks, and visualizations produced in class.

What we did

In this live session, we set up and ran a two-stream instability simulation from scratch using WarpX, following an exploratory approach:

We browsed the WarpX examples on GitHub and picked the uniform plasma example as our starting point.
We modified the input file:
- Duplicated the electron species and added opposite drift velocities (uz_m) to create two counter-streaming beams.
- Defined my_constants.nt so we could change the number of timesteps in one place.
- Introduced a temperature constant T and used it to compute the thermal spread: ux_th = sqrt(T/m_e)/clight (and same for uy_th, uz_th). For an isotropic distribution, the temperature is the same in each direction: $T = T_x = T_y = T_z$, and each component of the thermal spread is set independently.
- Switched diagnostics to openPMD format so we could use openPMD-viewer to grab the data.
We ran the simulation in 2D and made a quick Jupyter notebook to visualize the phase space $(z, u_z)$.
The phase space wasn’t changing – so we increased the number of timesteps.
We kept iterating: set the temperature to zero (cold beams), switched to 1D, increased the resolution and the number of particles per cell.
Eventually, we saw the instability develop: the counter-streaming beams formed vortex structures in phase space.
We made a video of the phase-space evolution.

Reference material

For a detailed walkthrough of the exploratory process we followed in class, see the Set Up a Simulation from Scratch episode.

For a polished version of the two-stream instability simulation, with exercises on parameter selection, reduced diagnostics, and analysis notebooks, see the Two-Stream Instability tutorial.

Files from the session

Input file used in class
Jupyter notebook for phase-space visualization and video generation

Video from class

Here is the phase-space video we produced at the end of the session:

As you can see, the dots are too large and it’s hard to see the fine structure. To improve the visualization, use a smaller marker size in the scatter plot, e.g. s=0.1. You may also notice that the instability takes a while to develop – the early frames are not very interesting.

Challenge

How can you make the instability grow faster?

Try to think of what controls the growth rate. What parameters could you change so that the instability develops sooner?

Show me the solution

A couple of ideas:

Increase the density. The instability growth rate scales with the plasma frequency $\omega_{pe} \propto \sqrt{n_0}$, so a higher density means faster growth.
Use random particle positions instead of uniform ones. With NUniformPerCell, the particles start on a regular grid, which is very quiet – the instability has to grow from roundoff-level noise. Switching to NRandomPerCell introduces more initial noise and can seed the instability earlier. Try changing injection_style from NUniformPerCell to NRandomPerCell and num_particles_per_cell_each_dim to num_particles_per_cell.

Key Points

Setting up a simulation is an iterative process: start from an existing example, modify, run, visualize, and refine.

The WarpX documentation and the examples gallery are the best starting points.

When things don’t look right, check the basics: enough timesteps, adequate resolution, and appropriate temperature.

Content from US-FCC 2026: Beam-Beam and Tracking Tutorial

Last updated on 2026-09-19 | Edit this page

Overview

Questions

How can WarpX, ImpactX, Xsuite, and MAD-X be combined for FCC-ee studies?
How closely do ImpactX, Xsuite, and MAD-X agree for linear optics?
How can a WarpX beam-beam interaction be included in a multi-turn simulation with Xsuite?

Objectives

Set up the software environment used by all three exercises.
Simulate one FCC-ee bunch crossing with WarpX.
Compare linear optics through an FCC-ee lattice with ImpactX, Xsuite, and MAD-X.
Combine a WarpX beam-beam interaction with multi-turn Xsuite tracking.

Overview

Lesson authors: Arianna Formenti (LBNL) and Peter Kicsiny (SLAC).

The three exercises cover different scales of an FCC-ee simulation.

A single bunch crossing. WarpX models the electromagnetic interaction between an electron bunch and a positron bunch, including several QED processes such as beamstrahlung, radiative Bhabha scattering, and incoherent pair generation.
An optics comparison. ImpactX, Xsuite, and MAD-X evaluate the same linear FCC-ee lattice with a common reference particle; ImpactX transports a matched beam covariance through it.
Many turns around the ring. Xsuite tracks the bunches around the FCC-ee lattice and WarpX supplies the beam-beam interaction at the collision point.

All three exercises use one Conda environment. Tutorial 1 is designed to run on a laptop in a few minutes. Its numerical resolution is deliberately modest and should not be used as a production setup.

Installation

We assume you have a Conda installation available in your machine. Download the Conda environment file, open a terminal in the directory containing it, and create the environment:

BASH

conda env create -f environment_setup.yml

Activate it with:

BASH

conda activate usfcc26-warpx-tutorial

The environment contains WarpX, ImpactX, Xsuite, MAD-X through CPyMAD, openPMD-viewer, Jupyter, and the Python packages used by the analysis notebooks.

View the environment file

OUTPUT

name: usfcc26-warpx-tutorial
channels:
  - conda-forge
dependencies:
  - python=3.14
  - warpx=26.09
  - impactx=26.09
  - openmpi
  - numpy
  - pandas
  - matplotlib
  - jupyterlab
  - openpmd-viewer
  - pip
  - pip:
      - xsuite
      - cpymad
      - numba

You can use any notebook interface. For example, start Jupyter Lab from the exercise directory with:

BASH

jupyter lab

When you have finished working, deactivate the environment:

BASH

conda deactivate

To remove it completely at a later date, run:

BASH

conda env remove -n usfcc26-warpx-tutorial

Tutorial 1: One FCC-ee bunch crossing with WarpX

What we will simulate

At the FCC-ee Z pole, a 45.6 GeV electron bunch collides with a 45.6 GeV positron bunch at a full crossing angle of 30 mrad. The bunches are flat, with an rms size of about 8 micrometers horizontally and 35 nanometers vertically at the interaction point. The bunches are 16 mm long.

WarpX computes the beam-beam fields and advances both bunches through the collision. The same simulation also generates beamstrahlung photons, radiative Bhabha photons, and incoherent electron-positron pairs.

Download the files

Create a directory for Tutorial 1 and place all of the following files in it:

The notebook and utility module must remain in the same directory as the input file.

Read the input file

The input file is split into constants, numerical settings, particle species, collision processes, and diagnostics. Values defined under my_constants can be reused throughout the rest of the file.

View the WarpX input file

OUTPUT

####################
### MY CONSTANTS ###
####################
my_constants.mc2   = m_e*clight*clight
my_constants.GeV   = q_e*1.e9
my_constants.pico  = 1.e-12
my_constants.nano  = 1.e-9
my_constants.micro = 1.e-6
my_constants.milli = 1.e-3

# BEAMS
my_constants.energy         = 45.6*GeV
my_constants.energy_eV      = energy / q_e
my_constants.gamma          = energy / mc2
my_constants.npart          = 20.20e10
my_constants.nmacropart     = 1e5
my_constants.charge         = q_e * npart
my_constants.betax          = 90*milli
my_constants.betay          = 0.7*milli
my_constants.sigmax         = 7993.75*nano
my_constants.sigmay         = 35.20*nano
my_constants.emitx          = sigmax*sigmax / betax
my_constants.emity          = sigmay*sigmay / betay
my_constants.emitx_n        = emitx * gamma
my_constants.emity_n        = emity * gamma
my_constants.sigmaz         = 16.7*milli
my_constants.dux            = emitx_n / sigmax
my_constants.duy            = emity_n / sigmay
my_constants.espread        = 1.34e-3 * gamma
my_constants.mux            = 0.0*sigmax
my_constants.muy            = 0.0*sigmay
my_constants.muz            = 0.25*Lz
my_constants.focal_distance = muz
my_constants.crossing_angle = 30e-3
my_constants.rotation_angle = 0.5*crossing_angle
my_constants.sigmax_init    = sigmax*sqrt(1 + (muz/betax)**2)
my_constants.sigmay_init    = sigmay*sqrt(1 + (muz/betay)**2)

# BOX
my_constants.Lx = 8*( sigmax_init*cos(rotation_angle) + sigmaz*sin(rotation_angle) )
my_constants.Ly = 8*sigmay_init
my_constants.Lz = 16*sigmaz*cos(rotation_angle)
my_constants.nx = 64    # 512
my_constants.ny = 64    # 2048
my_constants.nz = 128   # 512
my_constants.dx = Lx/nx
my_constants.dy = Ly/ny
my_constants.dz = Lz/nz

# TIME
my_constants.T = 0.5*Lz/clight
my_constants.dt = T / nz
my_constants.nt = floor(T/dt)

# VIRTUAL PHOTONS
# minimum energy
my_constants.hwmin = 1e-4 * mc2 * mc2 / energy

# COLLISIONS
my_constants.sigma_kn_max = 6.65e-29 # m^2, maximum total Compton cross section (Klein-Nishina)
my_constants.sigma_bw_max = 1.7e-29 # m^2, maximum total Breit-Wheeler  cross section
my_constants.probability_estimate =   npart / nmacropart * sigma_bw_max * 2 * clight * dt  / (dx * dy * dz)
my_constants.probability_target_value = 0.01 # we want the event to have this probability
my_constants.multiplier_bw = probability_target_value / probability_estimate
my_constants.probability_threshold = 0.1
# The sampling probability scales as cross section times event multiplier.
# Reduce the Klein-Nishina multiplier so its approximate maximum probability
# matches the Breit-Wheeler target used above.
my_constants.multiplier_kn = sigma_bw_max / sigma_kn_max * multiplier_bw

# DIAGNOSTICS
my_constants.bin_num_1d = 512
my_constants.bin_center_min_1d = 0.
my_constants.bin_center_max_1d = 2.1*energy_eV
my_constants.bin_size_1d = (bin_center_max_1d - bin_center_min_1d) / bin_num_1d
my_constants.bin_edge_min_1d = bin_center_min_1d - 0.5 * bin_size_1d
my_constants.bin_edge_max_1d = bin_center_max_1d + 0.5 * bin_size_1d
my_constants.bin_num_1d_eff = bin_num_1d + 1
warpx.used_inputs_file = warpx_used_inputs.txt

##########################
### GENERAL PARAMETERS ###
##########################
stop_time = T
amr.n_cell = nx ny nz
amr.max_level = 0
geometry.dims = 3
geometry.prob_lo = -0.5*Lx -0.5*Ly -0.5*Lz
geometry.prob_hi =  0.5*Lx  0.5*Ly  0.5*Lz

###########################
### BOUNDARY CONDITIONS ###
###########################
boundary.field_lo = open open open
boundary.field_hi = open open open
boundary.particle_lo = Absorbing Absorbing Absorbing
boundary.particle_hi = Absorbing Absorbing Absorbing

################
### NUMERICS ###
################
warpx.do_electrostatic = relativistic
warpx.const_dt = dt
warpx.grid_type = collocated
algo.particle_shape = 3
algo.particle_pusher = vay
warpx.poisson_solver = fft
warpx.use_2d_slices_fft_solver = 1
warpx.random_seed = 20260908

#################
### PARTICLES ###
#################
particles.species_names = beam1 beam2 pho1 pho2 ele_bh pos_bh ele_bw pos_bw ele_ll pos_ll vpho1 vpho2 pho1_bha pho2_bha beam1_bha beam2_bha test1 test2

# Beams
beam1.species_type = electron
beam1.injection_style = gaussian_beam
beam1.x_rms = sigmax
beam1.y_rms = sigmay
beam1.z_rms = sigmaz
beam1.x_m = - mux
beam1.y_m = - muy
beam1.z_m = - muz
beam1.npart = nmacropart
beam1.q_tot = -charge
beam1.focal_distance = focal_distance
beam1.do_gaussian_beam_rotation = 1
beam1.do_gaussian_beam_rotation_momenta = 0
beam1.gaussian_beam_rotation_angle = rotation_angle
beam1.gaussian_beam_rotation_axis = 0 1 0
beam1.momentum_distribution_type = gaussian
beam1.uz_m = gamma
beam1.uy_m = 0.0
beam1.ux_m = 0.0
beam1.ux_th = dux
beam1.uy_th = duy
beam1.uz_th = espread
beam1.initialize_self_fields = 1
beam1.do_qed_quantum_sync = 1 # BS
beam1.qed_quantum_sync_phot_product_species = pho1
beam1.do_classical_radiation_reaction = 0
beam1.do_qed_virtual_photons = 1
# Keep the baseline consistent with the analytical Bhabha comparison, which
# neglects the finite beam-size effect. Set this and the beam2 flag to 1 for a
# separate finite-beam-size study.
beam1.qed_virtual_photons_do_beam_size_effect = 0
beam1.qed_virtual_photon_species_name = vpho1

beam2.species_type = positron
beam2.injection_style = gaussian_beam
beam2.x_rms = sigmax
beam2.y_rms = sigmay
beam2.z_rms = sigmaz
beam2.x_m = mux
beam2.y_m = muy
beam2.z_m = muz
beam2.npart = nmacropart
beam2.q_tot = charge
beam2.focal_distance = focal_distance
beam2.do_gaussian_beam_rotation = 1
beam2.do_gaussian_beam_rotation_momenta = 0
beam2.gaussian_beam_rotation_angle = -rotation_angle
beam2.gaussian_beam_rotation_axis = 0 1 0
beam2.momentum_distribution_type = gaussian
beam2.uz_m = -gamma
beam2.uy_m = 0.0
beam2.ux_m = 0.0
beam2.ux_th = dux
beam2.uy_th = duy
beam2.uz_th = espread
beam2.initialize_self_fields = 1
beam2.do_qed_quantum_sync = 1
beam2.qed_quantum_sync_phot_product_species = pho2
beam2.do_classical_radiation_reaction = 0
beam2.do_qed_virtual_photons = 1
beam2.qed_virtual_photons_do_beam_size_effect = 0
beam2.qed_virtual_photon_species_name = vpho2

# Test particles
test1.species_type = electron
test1.injection_style = gaussian_beam
test1.x_rms = sigmax
test1.y_rms = sigmay
test1.z_rms = sigmaz
test1.x_m = - mux
test1.y_m = - muy
test1.z_m = - muz
test1.npart = 100
test1.q_tot = -charge
test1.focal_distance = focal_distance
test1.do_gaussian_beam_rotation = 1
test1.do_gaussian_beam_rotation_momenta = 0
test1.gaussian_beam_rotation_angle = rotation_angle
test1.gaussian_beam_rotation_axis = 0 1 0
test1.momentum_distribution_type = gaussian
test1.uz_m = gamma
test1.uy_m = 0.0
test1.ux_m = 0.0
test1.ux_th = dux
test1.uy_th = duy
test1.uz_th = espread
test1.do_not_deposit = 1

test2.species_type = positron
test2.injection_style = gaussian_beam
test2.x_rms = sigmax
test2.y_rms = sigmay
test2.z_rms = sigmaz
test2.x_m = mux
test2.y_m = muy
test2.z_m = muz
test2.npart = 100
test2.q_tot = charge
test2.focal_distance = focal_distance
test2.do_gaussian_beam_rotation = 1
test2.do_gaussian_beam_rotation_momenta = 0
test2.gaussian_beam_rotation_angle = -rotation_angle
test2.gaussian_beam_rotation_axis = 0 1 0
test2.momentum_distribution_type = gaussian
test2.uz_m = -gamma
test2.uy_m = 0.0
test2.ux_m = 0.0
test2.ux_th = dux
test2.uy_th = duy
test2.uz_th = espread
test2.do_not_deposit = 1

# Beamstrahlung photons
pho1.species_type = photon
pho1.injection_style = none

pho2.species_type = photon
pho2.injection_style = none

# Virtual photons
vpho1.qed_virtual_photons_min_energy = hwmin
vpho1.qed_virtual_photons_multiplier = 1
vpho1.species_type = photon
vpho1.injection_style = none
vpho1.do_not_push = 1

vpho2.qed_virtual_photons_min_energy = hwmin
vpho2.qed_virtual_photons_multiplier = 1
vpho2.species_type = photon
vpho2.injection_style = none
vpho2.do_not_push = 1

# Incoherent pairs
ele_bw.species_type = electron
ele_bw.injection_style = none
ele_bw.do_not_deposit = 1

pos_bw.species_type = positron
pos_bw.injection_style = none
pos_bw.do_not_deposit = 1

ele_ll.species_type = electron
ele_ll.injection_style = none
ele_ll.do_not_deposit = 1

pos_ll.species_type = positron
pos_ll.injection_style = none
pos_ll.do_not_deposit = 1

ele_bh.species_type = electron
ele_bh.injection_style = none
ele_bh.do_not_deposit = 1

pos_bh.species_type = positron
pos_bh.injection_style = none
pos_bh.do_not_deposit = 1

# Bhabha photons
pho1_bha.species_type = photon
pho1_bha.injection_style = none
pho1_bha.do_qed_breit_wheeler = 0

pho2_bha.species_type = photon
pho2_bha.injection_style = none
pho2_bha.do_qed_breit_wheeler = 0

# Backscattered beam particles
beam1_bha.species_type = electron
beam1_bha.injection_style = none
beam1_bha.do_not_deposit = 1

beam2_bha.species_type = positron
beam2_bha.injection_style = none
beam2_bha.do_not_deposit = 1

##############
### SF-QED ###
##############

qed_qs.photon_creation_energy_threshold = 0.

qed_qs.lookup_table_mode = builtin
qed_qs.chi_min = 1.e-8

warpx.do_qed_schwinger = 0.

##################
### COLLISIONS ###
##################
collisions.collision_names = ll bh1 bh2 lbw bhabha1 bhabha2

lbw.species = pho1 pho2
lbw.type = linear_breit_wheeler
lbw.product_species = ele_bw pos_bw
lbw.event_multiplier = multiplier_bw
lbw.probability_threshold = probability_threshold
lbw.probability_target_value = probability_target_value

ll.species = vpho1 vpho2
ll.type = linear_breit_wheeler
ll.product_species = ele_ll pos_ll
ll.event_multiplier = multiplier_bw
ll.probability_threshold = probability_threshold
ll.probability_target_value = probability_target_value

bh1.species = vpho1 pho2
bh1.type = linear_breit_wheeler
bh1.product_species = ele_bh pos_bh
bh1.event_multiplier = multiplier_bw
bh1.probability_threshold = probability_threshold
bh1.probability_target_value = probability_target_value

bh2.species = pho1 vpho2
bh2.type = linear_breit_wheeler
bh2.product_species = ele_bh pos_bh
bh2.event_multiplier = multiplier_bw
bh2.probability_threshold = probability_threshold
bh2.probability_target_value = probability_target_value

bhabha1.species = vpho2 beam1
bhabha1.type = linear_compton
bhabha1.product_species = pho1_bha beam1_bha
bhabha1.event_multiplier = multiplier_kn
bhabha1.probability_threshold = probability_threshold
bhabha1.probability_target_value = probability_target_value

bhabha2.species = vpho1 beam2
bhabha2.type = linear_compton
bhabha2.product_species = pho2_bha beam2_bha
bhabha2.event_multiplier = multiplier_kn
bhabha2.probability_threshold = probability_threshold
bhabha2.probability_target_value = probability_target_value

###################
### DIAGNOSTICS ###
###################
# These will contain per particle coordinates, momenta and weights
diagnostics.diags_names = particles_in particles_out trajectories

trajectories.intervals = 1
trajectories.write_species = 1
trajectories.diag_type = Full
trajectories.species = test1 test2
trajectories.fields_to_plot = none
trajectories.format = openpmd
trajectories.openpmd_backend = bp
trajectories.dump_last_timestep = 1
trajectories.test1.additional_variables = Ex Ey Ez Bx By Bz
trajectories.test2.additional_variables = Ex Ey Ez Bx By Bz

particles_in.intervals = floor(nz/4)
particles_in.write_species = 1
particles_in.diag_type = Full
particles_in.species = beam1 beam2 pho1 pho2 ele_bh pos_bh ele_bw pos_bw ele_ll pos_ll pho1_bha pho2_bha beam1_bha beam2_bha
particles_in.fields_to_plot = none
particles_in.format = openpmd
particles_in.openpmd_backend = bp
particles_in.dump_last_timestep = 1

particles_out.intervals = -1
particles_out.diag_type = BoundaryScraping
particles_out.format = openpmd
particles_out.openpmd_backend = bp
particles_out.dump_last_timestep = 1

beam1.save_particles_at_xlo = 1
beam1.save_particles_at_ylo = 1
beam1.save_particles_at_zlo = 1
beam1.save_particles_at_xhi = 1
beam1.save_particles_at_yhi = 1
beam1.save_particles_at_zhi = 1

beam2.save_particles_at_xlo = 1
beam2.save_particles_at_ylo = 1
beam2.save_particles_at_zlo = 1
beam2.save_particles_at_xhi = 1
beam2.save_particles_at_yhi = 1
beam2.save_particles_at_zhi = 1

pho1.save_particles_at_xlo = 1
pho1.save_particles_at_ylo = 1
pho1.save_particles_at_zlo = 1
pho1.save_particles_at_xhi = 1
pho1.save_particles_at_yhi = 1
pho1.save_particles_at_zhi = 1

pho2.save_particles_at_xlo = 1
pho2.save_particles_at_ylo = 1
pho2.save_particles_at_zlo = 1
pho2.save_particles_at_xhi = 1
pho2.save_particles_at_yhi = 1
pho2.save_particles_at_zhi = 1

ele_ll.save_particles_at_xlo = 1
ele_ll.save_particles_at_ylo = 1
ele_ll.save_particles_at_zlo = 1
ele_ll.save_particles_at_xhi = 1
ele_ll.save_particles_at_yhi = 1
ele_ll.save_particles_at_zhi = 1

pos_ll.save_particles_at_xlo = 1
pos_ll.save_particles_at_ylo = 1
pos_ll.save_particles_at_zlo = 1
pos_ll.save_particles_at_xhi = 1
pos_ll.save_particles_at_yhi = 1
pos_ll.save_particles_at_zhi = 1

ele_bh.save_particles_at_xlo = 1
ele_bh.save_particles_at_ylo = 1
ele_bh.save_particles_at_zlo = 1
ele_bh.save_particles_at_xhi = 1
ele_bh.save_particles_at_yhi = 1
ele_bh.save_particles_at_zhi = 1

pos_bh.save_particles_at_xlo = 1
pos_bh.save_particles_at_ylo = 1
pos_bh.save_particles_at_zlo = 1
pos_bh.save_particles_at_xhi = 1
pos_bh.save_particles_at_yhi = 1
pos_bh.save_particles_at_zhi = 1

ele_bw.save_particles_at_xlo = 1
ele_bw.save_particles_at_ylo = 1
ele_bw.save_particles_at_zlo = 1
ele_bw.save_particles_at_xhi = 1
ele_bw.save_particles_at_yhi = 1
ele_bw.save_particles_at_zhi = 1

pos_bw.save_particles_at_xlo = 1
pos_bw.save_particles_at_ylo = 1
pos_bw.save_particles_at_zlo = 1
pos_bw.save_particles_at_xhi = 1
pos_bw.save_particles_at_yhi = 1
pos_bw.save_particles_at_zhi = 1

pho1_bha.save_particles_at_xlo = 1
pho1_bha.save_particles_at_ylo = 1
pho1_bha.save_particles_at_zlo = 1
pho1_bha.save_particles_at_xhi = 1
pho1_bha.save_particles_at_yhi = 1
pho1_bha.save_particles_at_zhi = 1

pho2_bha.save_particles_at_xlo = 1
pho2_bha.save_particles_at_ylo = 1
pho2_bha.save_particles_at_zlo = 1
pho2_bha.save_particles_at_xhi = 1
pho2_bha.save_particles_at_yhi = 1
pho2_bha.save_particles_at_zhi = 1

beam1_bha.save_particles_at_xlo = 1
beam1_bha.save_particles_at_ylo = 1
beam1_bha.save_particles_at_zlo = 1
beam1_bha.save_particles_at_xhi = 1
beam1_bha.save_particles_at_yhi = 1
beam1_bha.save_particles_at_zhi = 1

beam2_bha.save_particles_at_xlo = 1
beam2_bha.save_particles_at_ylo = 1
beam2_bha.save_particles_at_zlo = 1
beam2_bha.save_particles_at_xhi = 1
beam2_bha.save_particles_at_yhi = 1
beam2_bha.save_particles_at_zhi = 1

# REDUCED
warpx.reduced_diags_names = DiffLumi_beam1_beam2 ColliderRelevant

DiffLumi_beam1_beam2.type = DifferentialLuminosity
DiffLumi_beam1_beam2.intervals = nt
DiffLumi_beam1_beam2.species = beam1 beam2
DiffLumi_beam1_beam2.bin_number = bin_num_1d_eff
DiffLumi_beam1_beam2.bin_max = bin_edge_max_1d
DiffLumi_beam1_beam2.bin_min = bin_edge_min_1d

ColliderRelevant.species = beam1 beam2
ColliderRelevant.type = ColliderRelevant

The laptop setup uses:

100000 macroparticles in each primary bunch
100 passive test particles for each beam
a 64 x 64 x 128 mesh
128 time steps
a two-dimensional-slices FFT Poisson solver
the built-in quantum-synchrotron lookup table
a fixed random seed for reproducibility

The mesh is much coarser than a production FCC-ee simulation. It is sufficient for learning the workflow and recovering the main features of the collision. The built-in QED table also keeps the download small but has low resolution. A quantitative physics study would require a higher-resolution QED table, together with mesh and macroparticle convergence scans.

The Conda Forge WarpX package used by this lesson includes QED support and the built-in tables, but it does not currently install the standalone qed_table_generator executable. To generate a higher-resolution table, build WarpX from source with WarpX_QED_TOOLS=ON, following the WarpX QED table-tool instructions. The resulting qed_table_generator can create a quantum-synchrotron table that WarpX reads with qed_qs.lookup_table_mode = load and qed_qs.load_table_from = /path/to/table.

The primary species and their main products are:

Species	Role
`beam1`, `beam2`	Initial electron and positron bunches
`test1`, `test2`	Passive electron and positron probes used for trajectories
`pho1`, `pho2`	Beamstrahlung photons emitted by the two bunches
`vpho1`, `vpho2`	Virtual photons used by the equivalent-photon model
`ele_ll`, `pos_ll`	Landau-Lifshitz pair products
`ele_bh`, `pos_bh`	Breit-Heitler pair products
`ele_bw`, `pos_bw`	Breit-Wheeler pair products
`pho1_bha`, `pho2_bha`	Photons from radiative Bhabha scattering
`beam1_bha`, `beam2_bha`	Scattered primary particles from Bhabha events

Beamstrahlung is enabled directly on each primary beam through the strong-field QED modules. The finite beam-size correction in the equivalent-photon model is disabled in the baseline input. This makes the radiative Bhabha result consistent with the analytical comparison below, which also neglects that correction; it does not disable beamstrahlung. To study finite-beam-size suppression, set qed_virtual_photons_do_beam_size_effect = 1 for both primary beams and run the modified input in a separate directory so that the baseline diagnostics are retained. The incoherent-pair and radiative Bhabha channels are listed in collisions.collision_names. Event multipliers increase the number of sampled rare events, while particle weights preserve the physical yield. Because the sampling probability scales with cross section times event multiplier, the Klein–Nishina multiplier is reduced by sigma_bw_max / sigma_kn_max relative to the Breit–Wheeler multiplier. This keeps their approximate maximum sampling probabilities near the same target.

Notice that the radiative Bhabha products are stored in new species. This makes the emitted photons and scattered primary particles easy to analyze without mixing them into the original bunch species.

Diagnostics

The simulation writes three openPMD particle diagnostics:

diags/particles_in contains snapshots of particles that are still inside the simulation domain
diags/particles_out records particles when they cross a domain boundary
diags/trajectories records the two test-particle species at every time step

The test species start with the same distributions as the primary beams, but do_not_deposit = 1 prevents their charge and current from contributing to the field solve. They respond to the self-consistent beam-beam fields without changing those fields, making them useful single-particle probes.

The ColliderRelevant and DifferentialLuminosity reduced diagnostics are written to diags/reducedfiles.

The analysis utility module combines the particles that remain in the box with the boundary records. For escaped particles, it reconstructs their positions at the final simulation time from their scrape time and momentum:

\[ \boldsymbol{x}(t_{\mathrm{final}}) =\boldsymbol{x}(t_{\mathrm{scrape}}) +\boldsymbol{v}\left(t_{\mathrm{final}}-t_{\mathrm{scrape}}\right). \]

For photons, $\boldsymbol{v}=c\boldsymbol{p}/|\boldsymbol{p}|$. For massive particles, the utility obtains the velocity from $\boldsymbol{v}=\boldsymbol{p}/(\gamma m_e)$. This reconstruction lets the later spectrum and pair plots include products that have already left the mesh instead of silently discarding them.

Run WarpX

Activate the tutorial environment, change to the directory containing the three downloaded files, and run:

BASH

conda activate usfcc26-warpx-tutorial
warpx.3d tutorial_1_input.txt

On the laptop used to prepare this tutorial, the run takes about three to four minutes and produces roughly 80 MB of diagnostics. Runtime will vary with the processor and storage system.

When the run finishes, your directory should contain:

tutorial_1_input.txt
tutorial_1_plots.ipynb
tutorial_1_utils.py
diags/
warpx_used_inputs.txt

warpx_used_inputs.txt is generated by WarpX. It records the fully evaluated input and is useful when an input value is defined through an expression.

Analyze the collision

Open the analysis notebook with your preferred Jupyter interface and run the cells in order. The notebook expects the diags directory in the current Tutorial 1 directory.

View the analysis notebook

The utility module contains the routines used to read particle data, reconstruct boundary particles, integrate the hourglass luminosity, and calculate the radiative Bhabha cross section and lifetime.

View the analysis utility module

OUTPUT

import os
import re
import warnings

import numpy as np
from openpmd_viewer import OpenPMDTimeSeries
from scipy.constants import c, e, m_e
from scipy.integrate import quad, trapezoid


def extract_macroparticles(species_list, sim_folder=".", diags_name="diags", step=-1):
    """
    Read WarpX coordinates located at 'sim_folder' in 'diags_name'
    corresponding to simulation timestep 'step'.
    """
    x_list = []
    y_list = []
    z_list = []
    w_list = []
    ux_list = []
    uy_list = []
    uz_list = []

    if step==-1:
        folder_list = [
            diags_name+'/particles_in',
            diags_name+'/particles_out/particles_at_xlo',
            diags_name+'/particles_out/particles_at_xhi',
            diags_name+'/particles_out/particles_at_ylo',
            diags_name+'/particles_out/particles_at_yhi',
            diags_name+'/particles_out/particles_at_zlo',
            diags_name+'/particles_out/particles_at_zhi',
            ]
    else:
        folder_list = [ diags_name+'/particles_in', ]
    # Loop through the files that contain particles collected at the edges and in the box
    for folder_name in folder_list:
        read_dir = os.path.join(sim_folder, folder_name)
        if os.path.isdir(read_dir):
            series = OpenPMDTimeSeries(read_dir)
            time = series.t[step]
            iteration = series.iterations[step]
            dt =  series.t[-1] / series.iterations[-1]

            for species in species_list:
                x, y, z, ux, uy, uz, w = series.get_particle( ['x', 'y', 'z', 'ux', 'uy', 'uz', 'w'], iteration=iteration, species=species )

                if ("particles_at" in folder_name):
                    it_scrape, = series.get_particle( ['stepScraped', ], iteration=iteration, species=species )
                    t_scrape = it_scrape * dt

                    momentum_squared = ux**2 + uy**2 + uz**2
                    if "pho" in species:
                        momentum_magnitude = np.sqrt(momentum_squared)
                        vx = c * ux / momentum_magnitude
                        vy = c * uy / momentum_magnitude
                        vz = c * uz / momentum_magnitude
                    else:
                        gamma = np.sqrt(1.0 + momentum_squared / (m_e * c) ** 2)
                        vx = ux / (gamma * m_e)
                        vy = uy / (gamma * m_e)
                        vz = uz / (gamma * m_e)

                    time_since_scrape = time - t_scrape
                    x = x + vx * time_since_scrape
                    y = y + vy * time_since_scrape
                    z = z + vz * time_since_scrape

                # convert from SI [kg m s-1] to [eV/c]
                conversion_factor = c/e

                x_list=np.append(x_list, x)
                y_list=np.append(y_list, y)
                z_list=np.append(z_list, z)
                w_list=np.append(w_list, w)
                ux_list=np.append(ux_list, ux * conversion_factor)
                uy_list=np.append(uy_list, uy * conversion_factor)
                uz_list=np.append(uz_list, uz * conversion_factor)

    # x y z [m], ux, uy, uz [eV/c], w=weight (no dimension)
    return np.asarray(x_list), np.asarray(y_list), np.asarray(z_list), np.asarray(ux_list), np.asarray(uy_list), np.asarray(uz_list), np.asarray(w_list)

def get_Ecom(filename):
    """
    Return 1 numpy array:
    - the center-of-mass energy (in eV)
    """
    with open(filename) as f:
        # First line: header, contains the energies
        line = f.readline()
        Ecom = np.array( list(map( float, re.findall('=(.*?)\\(', line) )) )
    return Ecom

def get_dL_dEcom(filename):
    """
    Return the cumulative differential luminosity and its energy integral.

    Returns:
    - the center-of-mass energy [eV]
    - differential luminosity [m^-2 eV^-1]
    - luminosity integrated over center-of-mass energy [m^-2]
    """
    Ecom = get_Ecom(filename)  # eV
    # ``ndmin=2`` also handles a diagnostic containing only its final row.
    values = np.loadtxt(filename, ndmin=2)[-1, 2:]
    if values.size == Ecom.size + 1:
        # Current WarpX appends the total luminosity after the differential
        # energy bins. Keep it out of the plotted spectrum.
        dL_dEcom = values[:-1]  # m^-2 eV^-1
        Ltot = values[-1]  # m^-2
    elif values.size == Ecom.size:
        # Compatibility with older output that omitted the total column.
        warnings.warn(
            "DifferentialLuminosity has no final total-luminosity column; "
            "integrating the energy bins. Check the WarpX version.",
            RuntimeWarning,
            stacklevel=2,
        )
        dL_dEcom = values
        Ltot = trapezoid(dL_dEcom, Ecom)
    else:
        raise ValueError(
            "Unexpected DifferentialLuminosity layout: "
            f"found {values.size} values for {Ecom.size} energy bins"
        )

    return Ecom, dL_dEcom, Ltot


def luminosity_per_bx_hourglass(
    bunch_intensity,
    sigma_x,
    sigma_y,
    sigma_z,
    phi,
    beta_x_star,
    beta_y_star,
    epsabs=0.0,
    epsrel=1.0e-10,
):
    """Return luminosity per bunch crossing, including the hourglass effect.

    This evaluates the longitudinal overlap of two identical Gaussian bunches.
    The transverse beam sizes evolve around the interaction-point waist as
    ``sigma(s) = sigma_star * sqrt(1 + (s / beta_star)**2)``. The crossing
    angle ``phi`` is the half crossing angle.

    Parameters are in SI units (meters and radians); the returned luminosity
    per bunch crossing is in ``m^-2``.
    """
    positive_parameters = {
        "bunch_intensity": bunch_intensity,
        "sigma_x": sigma_x,
        "sigma_y": sigma_y,
        "sigma_z": sigma_z,
        "beta_x_star": beta_x_star,
        "beta_y_star": beta_y_star,
    }
    for name, value in positive_parameters.items():
        if value <= 0.0:
            raise ValueError(f"{name} must be positive, got {value!r}")

    def integrand(s):
        sigma_x_s = sigma_x * np.sqrt(1.0 + (s / beta_x_star) ** 2)
        sigma_y_s = sigma_y * np.sqrt(1.0 + (s / beta_y_star) ** 2)
        longitudinal_overlap = np.exp(-(s / sigma_z) ** 2)
        crossing_angle_reduction = np.exp(-((phi * s) / sigma_x_s) ** 2)
        return longitudinal_overlap * crossing_angle_reduction / (
            sigma_x_s * sigma_y_s
        )

    overlap_integral, _ = quad(
        integrand,
        -np.inf,
        np.inf,
        epsabs=epsabs,
        epsrel=epsrel,
        limit=200,
    )
    normalization = 4.0 * np.pi * np.sqrt(np.pi) * sigma_z
    return bunch_intensity**2 * overlap_integral / normalization


def integrand_qed(y, beam_energy, electron_mass):
    """Return the radiative-Bhabha cross-section integrand.

    ``y`` is the emitted photon's fractional beam energy. ``beam_energy`` and
    ``electron_mass`` must use the same energy unit; the tutorial uses GeV.
    This expression neglects the beam-size effect.
    """
    if not 0.0 < y < 1.0:
        raise ValueError(f"y must lie strictly between 0 and 1, got {y!r}")
    if beam_energy <= 0.0 or electron_mass <= 0.0:
        raise ValueError("beam_energy and electron_mass must be positive")

    photon_spectrum = (4.0 / 3.0 + y**2 - 4.0 * y / 3.0) / y
    logarithmic_factor = (
        2.0 * np.log(4.0 * beam_energy**2 / electron_mass**2)
        + 2.0 * np.log((1.0 - y) / y)
        - 1.0
    )
    return photon_spectrum * logarithmic_factor


def beam_lifetime(
    cross_section,
    luminosity_per_bx,
    bunch_population,
    n_interaction_points,
    revolution_frequency,
    n_bunches=1,
):
    """Return the beam lifetime in hours.

    ``cross_section * luminosity_per_bx`` must be dimensionless. For example,
    use a cross section in mbarn with luminosity converted to mbarn^-1 per
    bunch crossing, as done in the tutorial notebook.
    """
    positive_parameters = {
        "cross_section": cross_section,
        "luminosity_per_bx": luminosity_per_bx,
        "bunch_population": bunch_population,
        "n_interaction_points": n_interaction_points,
        "revolution_frequency": revolution_frequency,
        "n_bunches": n_bunches,
    }
    for name, value in positive_parameters.items():
        if value <= 0.0:
            raise ValueError(f"{name} must be positive, got {value!r}")

    loss_rate = (
        cross_section
        * luminosity_per_bx
        * n_interaction_points
        * revolution_frequency
        * n_bunches
    )
    return bunch_population / loss_rate / 3600.0

What the notebook checks

The notebook does more than plot the output. It turns the particle and reduced diagnostics into beam-physics quantities and, where a compact analytical model is available, compares the simulation with that model. The seven stages below explain what is calculated and what each comparison can tell us.

1. Derived beam parameters

The notebook first reproduces the quantities needed by the later estimates. The relativistic factor and the rms angular divergences at the interaction point are

\[ \gamma=\frac{E_{\mathrm{beam}}}{m_ec^2}, \qquad \sigma_{x'}^*=\frac{\sigma_x^*}{\beta_x^*}, \qquad \sigma_{y'}^*=\frac{\sigma_y^*}{\beta_y^*}. \]

The 30 mrad value in the input is the full crossing angle, so the formulae use the half angle $\phi=15$ mrad. The Piwinski angle and the corresponding effective horizontal overlap size are

\[ \Phi=\frac{\sigma_z}{\sigma_x^*}\tan\phi, \qquad \sigma_{x,\mathrm{eff}}^*=\sigma_x^*\sqrt{1+\Phi^2}. \]

For the supplied parameters, $\Phi\simeq31.34$: the crossing angle therefore dominates the effective horizontal overlap.

The crossing angle increases the effective horizontal overlap of the two bunches.

The notebook evaluates the linearized beam-beam parameters

\[ \xi_x=\frac{N_b r_e\beta_x^*} {2\pi\gamma\sigma_{x,\mathrm{eff}}^* \left(\sigma_{x,\mathrm{eff}}^*+\sigma_y^*\right)}, \]

\[ \xi_y=\frac{N_b r_e\beta_y^*} {2\pi\gamma\sigma_y^* \left(\sigma_{x,\mathrm{eff}}^*+\sigma_y^*\right)}, \]

where $N_b$ is the bunch population and $r_e$ is the classical electron radius. The resulting values are approximately $\xi_x=0.0015$ and $\xi_y=0.0805$.

2. Collision snapshots

At each stored iteration, the notebook plots both primary beams and their beamstrahlung photons in the $(z,x)$ and $(z,y)$ planes. Coordinates are normalized as $z/\sigma_z$, $x/\sigma_x^*$, and $y/\sigma_y^*$ so that the very different horizontal, vertical, and longitudinal scales can be compared in the same figure. These views check the crossing geometry and direction of motion: beam 1 travels toward positive $z$, beam 2 toward negative $z$, and the photon products follow the parent beams. Plotting uses a stride for speed; the underlying diagnostic data are not down-sampled.

The notebook next creates one ParticleTracker for test1 and another for test2. Both trackers select their particles at the first trajectory output. With preserve_particle_index=True, a given particle ID retains the same array index at every later iteration; if a particle is absent, its entry is filled with NaN instead of shifting all subsequent trajectories. Plotting each column of the coordinate arrays therefore traces one physical particle through the collision. The 3D view shows the full $(z,x,y)$ paths, while the $(z,x)$ and $(z,y)$ projections make the small transverse motion easier to read. These passive trajectories reveal the crossing geometry and accumulated beam-beam deflection without the sampling noise of the full bunches.

3. Horizontal beam-beam kick

To measure the kick, the notebook selects the beam-1 macroparticles present in the final diagnostic and uses their unique particle IDs to retrieve those same particles from the initial diagnostic. This avoids pairing unrelated rows if the openPMD particle order changes. It then converts momenta into ultrarelativistic trajectory angles,

\[ x'=\frac{p_x}{p_z}, \qquad \Delta x'=x'_{\mathrm{end}}-x'_{\mathrm{start}}. \]

Close to the opposing bunch axis, a Gaussian beam acts like a linear focusing lens:

\[ \Delta x'(x)\simeq-\frac{4\pi\xi_x}{\beta_x^*}x. \]

Because $\sigma_{x'}^*=\sigma_x^*/\beta_x^*$, normalizing both axes gives a particularly direct test,

\[ \frac{\Delta x'}{\sigma_{x'}^*} \simeq -4\pi\xi_x\frac{x}{\sigma_x^*}. \]

The notebook fits a straight line over $|x|<10\sigma_x^*$. This wide-looking interval still samples the central part of the force because $\sigma_{x,\mathrm{eff}}^*\gg\sigma_x^*$ for this crossing angle. With the supplied laptop input, the fitted normalized slope is about -0.0171, while the analytical slope $-4\pi\xi_x$ is about -0.0183, a difference of roughly 7 percent. This tests the combined field solve and particle push at the chosen resolution. It is not a convergence test, and it considers only particles that remain in the domain through the final recorded iteration.

4. Photon spectra and the loss threshold

The helper combines photons still inside the domain with photons found in the boundary-scraping diagnostics. Since openPMD returns momentum in this analysis as eV/$c$, the photon energy in GeV is calculated from

\[ E_\gamma[\mathrm{GeV}] =10^{-9}\sqrt{p_x^2+p_y^2+p_z^2}. \]

The beamstrahlung and radiative Bhabha histograms are weighted histograms. In an energy bin $j$, the physical photon yield is

\[ N_{\gamma,j}=\sum_{i\in j}w_i, \]

not the number of macroparticle records. This distinction matters because the rare-event multipliers deliberately create extra sampled events and compensate through the macroparticle weights $w_i$.

For the Bhabha spectrum, the notebook also marks the energy associated with a 1 percent ring momentum acceptance. If $\delta_{\mathrm{acc}}=0.01$, a photon is counted as causing the loss of its emitting primary when

\[ E_\gamma>\delta_{\mathrm{acc}}E_{\mathrm{beam}}. \]

This is an approximate loss model: the actual ring acceptance depends on the lattice and on where the off-momentum particle travels.

5. Luminosity and the hourglass effect

The DifferentialLuminosity diagnostic supplies the cumulative center-of-mass-energy spectrum $d\mathcal{L}/dE_{\mathrm{com}}$, in m$^{-2}$ eV$^{-1}$. In current WarpX output, the last value on each data row is the total luminosity per bunch crossing,

\[ \mathcal{L}_{E}=\mathcal{L}_{\mathrm{total}}. \]

The utility separates that last value from the preceding energy bins before plotting. For compatibility with an older output layout that omitted the total, it emits a warning and evaluates the numerical energy integral

\[ \mathcal{L}_{E} \simeq\sum_j \left.\frac{d\mathcal{L}}{dE_{\mathrm{com}}}\right|_j\Delta E_j. \]

The ColliderRelevant diagnostic supplies $d\mathcal{L}/dt$, which gives an independent estimate through a trapezoidal time integral,

\[ \mathcal{L}_{\mathrm{WarpX}} \simeq\sum_k\frac{1}{2} \left[ \left.\frac{d\mathcal{L}}{dt}\right|_k+ \left.\frac{d\mathcal{L}}{dt}\right|_{k+1} \right](t_{k+1}-t_k). \]

Comparing $\mathcal{L}_{E}$ with $\mathcal{L}_{\mathrm{WarpX}}$ first checks the consistency of the two WarpX diagnostics. Their energy and time binning are different, so they are not expected to be bitwise identical; a material difference should prompt checks of the diagnostic binning and numerical resolution.

The first analytical estimate treats the transverse sizes as constant and accounts for the crossing angle only through the effective horizontal size:

\[ \mathcal{L}_0 =\frac{N_b^2}{4\pi\sigma_{x,\mathrm{eff}}^*\sigma_y^*}. \]

FCC-ee has a strong hourglass effect because the beta functions, especially $\beta_y^*$, are short compared with the bunch length. The utility module therefore also evaluates the longitudinal overlap numerically. At a distance $s$ from the interaction-point waist,

\[ \sigma_u(s)=\sigma_u^* \sqrt{1+\left(\frac{s}{\beta_u^*}\right)^2}, \qquad u\in\{x,y\}, \]

and the luminosity model used by the notebook is

\[ \mathcal{L}_{\mathrm{HG}} =\frac{N_b^2}{4\pi\sqrt{\pi}\sigma_z} \int_{-\infty}^{\infty} \frac{ \exp\!\left[-(s/\sigma_z)^2\right] \exp\!\left[-(\phi s/\sigma_x(s))^2\right] }{\sigma_x(s)\sigma_y(s)}\,ds. \]

The two exponential factors describe longitudinal bunch overlap and the crossing-angle reduction, respectively. The notebook reports $\mathcal{L}_{\mathrm{WarpX}}$, $\mathcal{L}_0$, $\mathcal{L}_{\mathrm{HG}}$, and the ratio $\mathcal{L}_{\mathrm{WarpX}}/\mathcal{L}_{\mathrm{HG}}$. At the deliberately coarse laptop resolution, a discrepancy should be interpreted as motivation for mesh, time-step, and macroparticle convergence studies rather than as a precision prediction.

6. Radiative Bhabha cross section and beam lifetime

For each beam, the weighted number of photons above the momentum-acceptance threshold is

\[ N_{\mathrm{loss}}= \sum_{E_{\gamma,i}>\delta_{\mathrm{acc}}E_{\mathrm{beam}}}w_i. \]

Dividing this yield by the simulated luminosity per crossing gives the WarpX cross-section estimate,

\[ \sigma_{\mathrm{Bhabha}}^{\mathrm{WarpX}} =\frac{N_{\mathrm{loss}}}{\mathcal{L}_{\mathrm{WarpX}}}. \]

The notebook compares it with the following QED estimate without the beam-size effect. Defining $y=E_\gamma/E_{\mathrm{beam}}$ and $\delta=\delta_{\mathrm{acc}}$,

\[ \sigma_{\mathrm{Bhabha}}^{\mathrm{QED}} =\frac{2\alpha^3}{m_e^2} \int_\delta^1 \frac{4/3+y^2-4y/3}{y} \left[ 2\ln\!\left(\frac{4E_{\mathrm{beam}}^2}{m_e^2}\right) +2\ln\!\left(\frac{1-y}{y}\right)-1 \right]dy. \]

Energies and masses are inserted in GeV, and the result in GeV$^{-2}$ is converted with $1\ \mathrm{GeV}^{-2}=0.389\ \mathrm{mbarn}$. Since the finite beam-size correction is disabled in the baseline WarpX input, this is a like-for-like comparison of two calculations that neglect that effect. Residual differences can come from the equivalent-photon and linear-Compton implementation, Monte Carlo statistics, numerical resolution, and the simple loss criterion above. Repeating the WarpX calculation with the correction enabled separately demonstrates its influence.

Finally, the revolution frequency is $f_{\mathrm{rev}}=c/C$, where $C$ is the ring circumference. For one bunch, the lifetime inferred from either cross section is

\[ \tau=\frac{N_b} {\sigma_{\mathrm{Bhabha}}\mathcal{L}_{\mathrm{bx}} n_{\mathrm{IP}}f_{\mathrm{rev}}}, \]

with $n_{\mathrm{IP}}=4$ and $\mathcal{L}_{\mathrm{bx}}=\mathcal{L}_{\mathrm{WarpX}}$. The notebook reports the electron and positron WarpX estimates and compares the positron result with the QED estimate.

7. Incoherent-pair distributions

The final stage combines the weights and coordinates of the electrons and positrons from the Landau-Lifshitz, Breit-Heitler, and Breit-Wheeler channels. The simulation keeps the colliding beams head-on in momentum space, so the pair momenta are first rotated into the physical crossing-angle frame. Products with $p_z>0$ are rotated by $+\phi$ about $y$, and products with $p_z<0$ by $-\phi$. For each rotated momentum it then calculates the polar angle

\[ \theta=\operatorname{atan2} \left(\sqrt{p_x^2+p_y^2},p_z\right), \]

and forms the weighted angular density in bin $j$,

\[ \left.\frac{dN}{d\theta}\right|_j \simeq\frac{1}{\Delta\theta_j}\sum_{i\in j}w_i. \]

The accompanying $(z,x)$ and $(x,y)$ projections show where the products are located at the final reconstructed time. These are useful first views of potential detector backgrounds, but the three production channels are summed and no detector geometry or transport through the accelerator magnets is included.

Callout

The close agreement of one observable does not establish numerical convergence. Increase the mesh resolution and macroparticle count before drawing quantitative conclusions from the QED yields or luminosity.

Tutorial 2: Linear optics comparison

What we will compare

This exercise evaluates one turn of the FCC-ee Z-pole fccee_p_ring sequence for a 45.6 GeV electron reference particle with ImpactX, Xsuite, and MAD-X. Rather than including beam-beam forces, radiation, or space charge, it isolates the single-particle linear optics. The main observables are the horizontal and vertical beta functions, $\beta_x(s)$ and $\beta_y(s)$, evaluated along the same lattice.

All three calculations start from fccee_z.madx. This is important: agreement would be much less meaningful if the codes read independently maintained versions of the lattice. Create a Tutorial 2 directory containing:

The notebook runs MAD-X through CPyMAD, imports the live MAD-X sequence into Xsuite, reads the ImpactX diagnostic, and compares all three results. The first run saves the converted Xsuite line as fccee_p_ring.json; later runs load this file instead of repeating the conversion. Keep the three supplied files together and execute them from that directory.

The matched beam model

At the interaction point, the geometric transverse emittances are constructed from the nominal rms sizes and beta functions,

\[ \varepsilon_x=\frac{(\sigma_x^*)^2}{\beta_x^*}, \qquad \varepsilon_y=\frac{(\sigma_y^*)^2}{\beta_y^*}. \]

For either transverse plane $u$, the matched covariance matrix has the Twiss form

\[ \Sigma_u=\varepsilon_u \begin{pmatrix} \beta_u & -\alpha_u\\ -\alpha_u & \gamma_u \end{pmatrix}, \qquad \gamma_u=\frac{1+\alpha_u^2}{\beta_u}. \]

The dispersive part of the horizontal beam size is included through the MAD-X dispersion. For an rms relative momentum spread $\sigma_\delta$,

\[ \sigma_x^2=\varepsilon_x\beta_x+(D_x\sigma_\delta)^2. \]

The script initializes ImpactX with the matched $\beta$, $\alpha$, $D$, and $D'$ values at the beginning of the sequence. For the longitudinal plane it uses

\[ \varepsilon_t=\sigma_z\sigma_\delta, \qquad \beta_t=\frac{\sigma_z}{\sigma_\delta}. \]

These definitions ensure that a disagreement later in the ring is testing the lattice translation and transport maps, rather than an intentionally mismatched initial beam.

Run the ImpactX calculation

Activate the common tutorial environment, change to the Tutorial 2 directory, and run:

BASH

conda activate usfcc26-warpx-tutorial
python fcc_impactx.py

The tested script uses envelope tracking by default:

PYTHON

DO_PARTICLE_TRACKING = False

ImpactX propagates the covariance matrix instead of sampling it with macroparticles. This removes Monte Carlo noise from the comparison and is much faster than tracking the script’s configured $10^7$ macroparticles. Space charge is disabled, the lattice is loaded from fccee_z.madx with one slice per element, and slice_step_diagnostics = True records the beam parameters along the ring.

After a successful run, the file used by the notebook is:

diags/reduced_beam_characteristics.0.0

It contains $s$, beam sizes, emittances, Twiss parameters, and dispersions at the ImpactX diagnostic locations.

Run the three-code comparison

Open the Tutorial 2 analysis notebook from the same directory and run its cells in order. The notebook performs the following steps:

CPyMAD loads fccee_z.madx, selects the fccee_p_ring sequence, defines a 45.6 GeV electron reference beam, and runs MAD-X TWISS.
If fccee_p_ring.json is absent, Xsuite constructs a thick-element Line from that in-memory MAD-X sequence and saves it. Otherwise, it loads the existing line. It then runs a four-dimensional Twiss calculation.
Pandas reads the ImpactX reduced diagnostic, and the notebook maps the common quantity names—for example, MAD-X betx, ImpactX beta_x, and Xsuite betx—onto the same plot.
The beta functions are compared over the full 90.6 km ring and over the first 2 km, where individual optics features are easier to inspect.

For a beam matrix transported by a linear map $R(s)$,

\[ \Sigma(s)=R(s)\Sigma(0)R(s)^{\mathsf T}, \qquad \beta_u(s)=\frac{\Sigma_{uu}(s)}{\varepsilon_u}. \]

This is the quantity produced by the ImpactX envelope calculation and overlaid with the Twiss functions from the other two codes. Because the codes record values at different longitudinal locations, any quantitative comparison would first require an explicit matching or interpolation convention.

Horizontal beta functions from MAD-X, ImpactX, and Xsuite over the first 2 km of the FCC-ee lattice.

Vertical beta functions from MAD-X, ImpactX, and Xsuite over the first 2 km of the FCC-ee lattice.

The curves from the tested setup overlap closely in both planes. The final ImpactX beta functions also return to within about 0.1 percent of their initial values after one turn. Inspect the interaction-region peaks and the one-turn closure separately: agreement at low-beta locations is a more sensitive check than agreement in slowly varying sections.

Callout

Small differences can come from element slicing, coordinate conventions, or the treatment of fringe fields and higher-order terms. This exercise uses one slice per element and a four-dimensional Xsuite Twiss calculation, so it checks the selected linear model; it does not establish equivalence for nonlinear or synchrotron motion.

Tutorial 3: Multi-turn WarpX and Xsuite simulation

What the coupled model does

This exercise connects two deliberately different descriptions of the collider. Xsuite transports the bunches through a linear map representing one FCC-ee superperiod, while WarpX resolves the collective electron-positron collision at the following interaction point. The same macroparticles are passed back and forth, so the beam-beam kick from one encounter affects every later encounter.

One iteration follows this sequence:

Xsuite advances both bunches through one linear superperiod.
The adapter converts the two Xsuite coordinate systems into one laboratory frame and writes one openPMD file per bunch.
A separate Python process starts WarpX, loads those particles, and advances them through a single head-on collision.
The adapter combines particles still inside the WarpX box with particles recorded at its absorbing boundaries.
Particle IDs are matched, the coordinates are converted back, and the Xsuite particle arrays are updated.
Centroids, rms sizes, covariances, and emittances are recorded before the next iteration.

The separation of responsibilities is important. Xsuite does not apply a second beam-beam kick, and WarpX does not model the 90.6 km arc.

Download the files

Create a Tutorial 3 directory containing:

Keep the six files together. The driver renders a run-specific WarpX input in the job directory; it never edits the downloaded template.

Before running the coupling

Tutorial 3 needs numpy, scipy, pandas, matplotlib, xobjects, xpart, xtrack, openpmd-api, openpmd-viewer, and a Python WarpX build with openPMD, FFT, and QED support. The QED capability is only exercised by the --beamstrahlung case, but the supplied input declares the photon product species and quantum-synchrotron table.

Run the driver from the directory containing the six downloaded files. You can inspect all command-line options with

BASH

python exec_tutorial_3.py --help

The --device gpu option selects an Xsuite ContextCupy; WarpX itself uses the backend of the installed WarpX package. The default direct launcher and one rank are intentional. The current adapter can only audit its fallback particle-ID mapping in a single-rank collision.

Simplified FCC-ee model

The reference momentum is 182.5 GeV/$c$, corresponding to the FCC-ee ttbar working point. For an IP with $alpha=0$, the initialized Gaussian widths are

\[ \sigma_x=\sqrt{\varepsilon_x\beta_x}, \qquad \sigma_{p_x}=\sqrt{\frac{\varepsilon_x}{\beta_x}}, \]

and likewise in $y$. The exercise sets $\varepsilon_x=\varepsilon_y=1.59\,\mathrm{nm}$ and $\beta_x=\beta_y=1\,\mathrm{m}$, producing round $39.9\,\mu\mathrm{m}$ beams. This is a teaching configuration, not the nominal flat FCC-ee beam.

The full ring is divided into four identical superperiods. The transverse phase advance used by each Xsuite map is therefore

\[ Q_{x,\mathrm{sp}} =\operatorname{frac}\!\left(\frac{398.148}{4}\right)=0.537, \]

with the same intentionally symmetric value in $y$. The longitudinal tune is $Q_{s,\mathrm{sp}}=0.091/4=0.02275$. Consequently, the command-line --iterations value counts superperiods: 20 iterations represent five full revolutions, not 20 revolutions.

The Xsuite element explicitly selects longitudinal_mode="linear_fixed_qs". This matters because qs and bets are parameters of that particular LineSegmentMap mode; leaving the mode implicit makes the intended synchrotron rotation dependent on constructor behavior.

Synchrotron radiation in the arc, the crossing angle, crab-waist sextupoles, and the nonlinear lattice are omitted. Beamstrahlung in WarpX is optional so that its longitudinal effect can be isolated.

The initialized values $\sigma_z=2.17$ mm and $\sigma_\delta=1.92\times10^{-3}$ are the beamstrahlung-broadened values retained from the original example. They describe an equilibrium only when the competing ring processes are also modeled. Turning on a fresh beamstrahlung kick at every IP while omitting arc damping, quantum excitation, and synchronous-energy compensation does not preserve that equilibrium; this tutorial should therefore be read as a transient coupling study.

The coordinate handoff

This is the numerically delicate part of the exercise. Xsuite stores

\[ p_x=\frac{P_x}{P_0},\qquad p_y=\frac{P_y}{P_0},\qquad \delta=\frac{P-P_0}{P_0}. \]

Therefore $p_x$ is not exactly the geometric slope. At a field-free handoff plane,

\[ \frac{P_s}{P_0} =\sqrt{(1+\delta)^2-p_x^2-p_y^2}, \qquad x'=\frac{P_x}{P_s} =\frac{p_x}{\sqrt{(1+\delta)^2-p_x^2-p_y^2}}. \]

Using $p_x$ directly as $x'$ introduces an energy-dependent drift error. It can be almost invisible without radiation and then become correlated with beamstrahlung energy loss.

WarpX stores normalized mechanical momenta $\boldsymbol{u}=\boldsymbol{P}/(m_ec)$. With $(\beta\gamma)_0=P_0/(m_ec)$, the forward conversion is

\[ u_x=(\beta\gamma)_0p_x, \qquad u_y=(\beta\gamma)_0p_y, \qquad u_s=(\beta\gamma)_0 \sqrt{(1+\delta)^2-p_x^2-p_y^2}. \]

After WarpX, the relative momentum deviation must be reconstructed from the full momentum magnitude,

\[ \delta =\frac{\sqrt{u_x^2+u_y^2+u_s^2}}{(\beta\gamma)_0}-1, \]

not from $u_s$ alone. The polished adapter uses these relations in both directions and tests an algebraic round trip before every collision. Its openPMD writer also attaches the required SI dimensions and unit scales to the position and momentum records; a bare array plus a momentum scale is not a complete external-particle description.

Xsuite defines the longitudinal coordinate as $\zeta=s-\beta_0ct$. There is consequently a distinction between Xsuite’s coordinates on the IP plane and the simultaneous particle snapshot required by WarpX. At the initial common time $t_-=-f/(\beta_0c)$, the local longitudinal position is

\[ s_- = \frac{u_s}{\gamma\beta_0}(\zeta-f), \]

and the transverse positions are drifted by $x_-=x+s_-u_x/u_s$ and $y_-=y+s_-u_y/u_s$. At the final common time $t_+=f/(\beta_0c)$, the inverse longitudinal transformation is

\[ \zeta=\frac{\beta_0\gamma}{u_s}s_+-f. \]

The earlier script implicitly replaced the particle velocity $cu_s/\gamma$ by $\beta_0c$. That is an excellent numerical approximation for the reference particle at 182.5 GeV, but it is unnecessary and ceases to be an exact round trip after radiation changes the momentum.

The WarpX bunch centers start a distance $f=4\sigma_z$ before the IP and end the same distance after it, so the collision duration is

\[ T=\frac{2f}{\beta_0c}. \]

Beam 2 requires additional sign changes because its local $s$ and $x$ axes point opposite to those of beam 1. The present adapter explicitly rejects a nonzero crossing angle: supporting one correctly would require a consistent transformation of positions, momenta, and time, rather than rotating positions alone.

Run the reference and coupled cases

Run the cases below separately. Giving them the same initial seed makes their initial bunches identical, while the driver assigns a new deterministic WarpX seed to each collision.

Case	Important option	Output directory	Purpose
Xsuite arc only	`--no-warpx`	`outputs_without_warpx`	Verify the linear map and unperturbed tunes
Classical collision	none	`outputs_with_warpx`	Test the coordinate handoff and collective beam-beam kick
Collision with beamstrahlung	`--beamstrahlung`	`outputs_with_warpx_beamstrahlung`	Add stochastic photon emission and longitudinal energy loss

Choose the run scale according to the question being asked:

Run level	Macroparticles	Superperiods	Purpose
Smoke test	1,000	1	Verify installation and one Xsuite–WarpX handoff
Tutorial	10,000	20 coupled; 256 arc-only	Inspect moments and coupling checks; obtain a modest-resolution reference FFT
High-statistics spectrum	1,000,000	1,024	Resolve coherent peaks with lower centroid noise on an appropriate GPU or HPC system

The smoke test is:

BASH

python exec_tutorial_3.py --macroparticles 1000 --iterations 1

It establishes that the workflow runs, but one sample cannot say anything about tunes or coherent modes. For the tutorial-scale analysis, first generate the inexpensive arc-only reference:

BASH

python exec_tutorial_3.py --no-warpx --macroparticles 10000 --iterations 256

Then run the classical WarpX collision. A short laptop test is:

BASH

python exec_tutorial_3.py --macroparticles 10000 --iterations 20

It writes outputs_with_warpx. To test the longitudinal coupling with quantum-synchrotron emission enabled, use:

BASH

python exec_tutorial_3.py --macroparticles 10000 --iterations 20 --beamstrahlung

The 20-superperiod coupled runs are intended to expose handoff or gross beam dynamics problems; their Fourier-bin spacing is $1/20=0.05$, which is too coarse for a meaningful coherent tune measurement. The 256-superperiod arc-only run has a finer spacing of about $0.0039$. A spectrum-quality coupled comparison should use hundreds to 1024 superperiods and enough macroparticles to suppress centroid noise.

That run writes outputs_with_warpx_beamstrahlung. Each collision starts a new WarpX process because a WarpX simulation cannot simply be reinitialized in the same Python process. The driver uses a different deterministic random seed for every collision. Reusing the same WarpX seed at every iteration would repeat the same Monte Carlo stream and could turn stochastic photon emission into an artificial coherent excitation.

The openPMD handoff deliberately does not write opticalDepthQSR. This is a WarpX QED runtime attribute that WarpX initializes stochastically for the imported electron and positron species. The previous openPMD file contained a zero-valued record, but the documented external_file interface does not include this runtime attribute and current WarpX initializes it separately. That record was therefore probably ignored or overwritten, so it is not an established cause of the observed oscillation. If a version did honor the supplied value, zero would mean that the sampled emission threshold had already been reached, rather than “no accumulated probability.”

The default direct launcher is appropriate for the single-rank Conda WarpX package. --launcher srun is available for an appropriately compiled HPC installation, but this teaching adapter currently enforces one WarpX rank so that its fallback particle-ID mapping remains auditable.

The run-specific working files are placed under --job-dir, which defaults to the current directory:

Path	Contents
`outputs_*/moments_b1.csv`, `moments_b2.csv`	One row of bunch moments after each superperiod and collision
`outputs_*/handoff_checks.csv`	Before/after WarpX comparisons for every primary particle
`warpx_read/beam1.bp`, `beam2.bp`	Simultaneous openPMD snapshots written by the adapter
`warpx_runtime_input.txt`	Rendered input used for the most recent collision
`warpx_used_inputs_tutorial_3.txt`	Parameters reported by WarpX for the most recent collision
`warpx_dump/diags`	Primary particles from the most recent WarpX call

By default the previous collision diagnostic is removed before the next call. Use --keep-diags only when the individual openPMD outputs are needed; it can consume substantial storage in a long run.

Optional high-statistics reference data

Generated 1-million-macroparticle, 1024-superperiod results do not belong in the lesson’s Git history. The reference-data instructions define the external archive layout, commands, and provenance information. If an archive is extracted as tutorial_3/reference_data/, the notebook detects its outputs_* directories automatically and plots them alongside any local runs.

Only the compact moment CSV files and handoff checks should be published in the archive. Per-collision openPMD files are large and are unnecessary for the notebook. The older coords_b*.txt results supplied with the previous adapter use a different schema and should not be represented as reference results from the revised coordinate handoff.

Coupling checks

Every coupled run writes handoff_checks.csv. For each bunch and collision it records:

the numbers of primary particles sent to and recovered from WarpX
whether WarpX preserved the supplied IDs or the documented single-rank input-order fallback was needed
the maximum algebraic position and momentum round-trip error
$\langle\zeta\rangle$, $\sigma_z$, $\langle\delta\rangle$, and $\sigma_\delta$ immediately before and after WarpX
the longitudinal rms emittance before and after WarpX
the mean energy change and rms instantaneous change in $\zeta$

The longitudinal rms emittance is

\[ \varepsilon_\zeta =\sqrt{ \sigma_\zeta^2\sigma_\delta^2 -\operatorname{Cov}(\zeta,\delta)^2 }. \]

These quantities separate a handoff error from longitudinal dynamics. With beamstrahlung off, the algebraic round-trip errors should be near floating point precision, and every primary should be recovered from either the full or boundary diagnostic. With beamstrahlung on, $\langle\Delta\delta\rangle<0$ and an increase in energy spread are expected; a large instantaneous $\zeta$ jump during the short WarpX call is more suspicious.

The three runs should have distinguishable signatures:

In the arc-only run, a distribution matched to bets should preserve its longitudinal rms emittance and approximately preserve $\sigma_z$ and $\sigma_\delta$. Its centroid spectrum locates the unperturbed map tune.
In the classical WarpX run, the transverse moments and coherent spectrum can change, but there should be no systematic radiative decrease in $\langle\delta\rangle$. This is the most sensitive end-to-end handoff test.
In the beamstrahlung run, the collision should first change $\langle\delta\rangle$, $\sigma_\delta$, and possibly $\varepsilon_\zeta$. A later change in $\sigma_z$ during Xsuite transport is the expected response of a longitudinally mismatched bunch.

If a run does not follow this pattern, use the following checks before interpreting the physics:

Confirm that n_returned == n_sent for both beams. A difference means that the final and boundary diagnostics did not recover every primary.
Inspect both algebraic round-trip errors. They test the adapter itself and should not depend on whether beamstrahlung is enabled.
Compare zeta_std_before_m and zeta_std_after_m. A large change here occurs inside the WarpX call; a change only in the following row occurs in the Xsuite arc.
Compare runs with different --seed values. A feature that changes phase or disappears can be a finite-macroparticle or Monte Carlo fluctuation.
Repeat with more longitudinal cells and macroparticles before drawing a quantitative conclusion from the size of the effect.

Why beamstrahlung can produce a real bunch-length oscillation

Even a perfect handoff does not imply constant $\sigma_z$. Beamstrahlung changes the momentum distribution locally at the IP but does not simultaneously rematch the bunch length and $\zeta$-$\delta$ covariance. The next linear synchrotron map transports the covariance as

\[ \Sigma_s(n+1)=R_s\Sigma_s(n)R_s^{\mathsf T}. \]

For the usual phase-space rotation, the new bunch-length variance contains contributions from both the old bunch length and energy spread,

\[ \sigma_{\zeta,n+1}^2 =\cos^2\mu_s\,\sigma_{\zeta,n}^2 +\beta_s^2\sin^2\mu_s\,\sigma_{\delta,n}^2 +2\beta_s\sin\mu_s\cos\mu_s \operatorname{Cov}(\zeta,\delta), \]

up to the sign convention of the longitudinal map. Thus a beamstrahlung increase in $\sigma_\delta$ naturally becomes an oscillation of $\sigma_z$. In this simplified model there is no radiation damping, quantum excitation, consistent compensation of the mean radiative loss, or self-consistent beamstrahlung equilibrium to damp or rematch that oscillation.

The most useful diagnostic sequence is:

Check whether $\sigma_z$ jumps across the WarpX call itself.
Check whether beamstrahlung first changes $\sigma_\delta$ and $\varepsilon_\zeta$.
Check whether $\sigma_z$ changes mainly after the next Xsuite arc map.
Repeat with several initial seeds to rule out a correlated Monte Carlo excitation.

Analyze the coherent modes

Open the Tutorial 3 notebook and run its cells after producing one or more output directories. It plots the beam moments, the handoff checks, and windowed centroid spectra.

Evolution of normalized beam moments for the arc-only, classical WarpX, and WarpX-with-beamstrahlung runs.

The following diagram illustrates the physical transverse displacement of the two beams at successive collision samples. In the $\sigma$ mode the beams move together, whereas in the $\pi$ mode they move oppositely. The diagram’s “Turn” labels are schematic: this tutorial records one superperiod and one IP encounter per sample, so four samples correspond to one full-ring turn.

Physical motion of two colliding bunches in the coherent sigma and pi modes at successive collision encounters.

Without beam-beam interaction, the centroid has one peak at the unperturbed superperiod tune. Two identical interacting beams instead have two coherent normal modes:

the $\sigma$ mode, in which the physical beam centroids move in phase and $Q_\sigma\simeq Q_0$
the $\pi$ mode, in which their physical centroids move oppositely and the beam-beam force shifts the coherent tune

Because the local horizontal axes point in opposite directions,

\[ x_\sigma=x_1-x_2, \qquad x_\pi=x_1+x_2. \]

The vertical axes have the same orientation, so $y_\sigma=y_1+y_2$ and $y_\pi=y_1-y_2$. For $N$ recorded superperiods, the FFT bin spacing is $\Delta Q=1/N$. Twenty iterations are useful for checking the workflow but are too short for a clean mode measurement; use hundreds of iterations for the spectrum.

Callout

This adapter assumes head-on beams, field-free handoff planes, ultra-relativistic primary particles that all appear in either the final or boundary diagnostic, and a single WarpX rank. Boundary-scraped particles are propagated ballistically to the common final time. These assumptions are checked or stated explicitly, but they remain part of the model.

Key Points

WarpX resolves one collision in detail and writes both particle snapshots and collider-specific reduced diagnostics.
The laptop input prioritizes runtime and teaching value over numerical convergence.
Matching particle IDs between diagnostics makes the normalized kick slope a direct check of the linear beam-beam model.
Keeping different QED products in separate species makes their spectra and weighted yields easier to inspect.
The luminosity, Bhabha cross-section, lifetime, and pair-distribution checks connect the openPMD and reduced diagnostics to accelerator observables.
Using one MAD-X sequence and one matched covariance model makes the Tutorial 2 beta-function comparison sensitive to differences in lattice translation and linear transport.
The Tutorial 3 adapter converts normalized momentum and geometric slope separately; this prevents beamstrahlung energy loss from contaminating the transverse drift.
A beamstrahlung-induced $\sigma_z$ oscillation is not automatically a coupling error: a local increase in energy spread creates a longitudinal mismatch that the arc map rotates into bunch length.
Per-collision handoff logs distinguish an instantaneous coordinate jump from subsequent synchrotron motion.

Content from LLNL HPC Innovation Center 2026: WarpX/ImpactX Tutorial

Last updated on 2026-09-19 | Edit this page

Overview

Questions

What happens when two groups of electrons move through each other?
How can a laser create a wave that accelerates electrons?
Why do the same magnets affect beams of different energies differently?
How could we pass a beam from one simulation to another?

Objectives

Set up the software environment shared by all three exercises.
Run a small one-dimensional (1D) simulation and plot how the electrons move.
Run a small three-dimensional (3D) plasma simulation and plot electron density, electric fields and the distribution of electron energies.
Use ImpactX to follow two example electron beams through accelerator magnets.
Recognize when particles are lost and when a simulation stops giving valid results.
Describe the extra steps needed to pass a beam from WarpX to ImpactX.

Overview

This tutorial starts with a small warm-up, followed by two independent visual examples: a plasma wake in WarpX and beam transport in ImpactX. An optional section explains how the two codes can be coupled for a more detailed study.

We will use two simulation tools. WarpX models charged particles and the electric and magnetic fields that act on them. Its particle-in-cell (PIC) method represents many real particles with each simulated particle and computes fields on a grid. ImpactX follows a beam through a sequence of accelerator components, such as magnets.

Two streams of electrons (two_stream_instability). Start with a quick 1D run: two groups of electrons move in opposite directions and develop a wave. This also checks that your simulation and plotting tools work. The standalone Two-Stream Instability episode lets you explore the parameters in more detail.
A laser-driven plasma wake (htu/lwfa_warpx). Watch a short laser pulse push electrons aside and leave a wave behind it. Plot where the electrons are, the electric field and their energies.
An electron beam passing through magnets (htu/beamline_impactx). Compare two computer-generated beams with energies of 100 MeV and 20 MeV. MeV means million electronvolts, a unit of particle energy. The magnets come from a model of the Hundred Terawatt Undulator (HTU) experiment. Watch the beam shapes change and see which particles reach the end.

The wakefield and magnet examples run independently; neither needs output from the other. Coupling them is an optional extension.

WarpX reads the settings from a plain-text input file using a short Python script. ImpactX uses a Python script to define the sequence of magnets, called a lattice, following the ImpactX HTU beamline example. You can run both examples from the supplied notebooks without writing code.

Use the WarpX notebook and the ImpactX notebook independently. Each focuses on running one example and visualizing its result. Optional coupling is explained at the end of this lesson.

Access and installation 💻

Registered participants: open the AWS link from Slack

Registered participants will receive a link in Slack that opens the AWS-hosted tutorial session directly in their browser. Open that link to start JupyterLab; you do not need to install software on your own computer. If you registered but cannot find your link, ask the organizers in Slack.

Open your session link from Slack.
In JupyterLab’s file browser, navigate to warpx-tutorials/episodes/files/llnl-hpc-2026/.
Start with two_stream_instability/two_stream_instability_plots.ipynb after running the warm-up below, or open htu/wakefield.ipynb and htu/htu_transport.ipynb for the two main examples.
Select WarpX GPU for the wakefield notebook and WarpX CPU for the ImpactX notebook using Kernel > Change Kernel. In the hosted wakefield notebook, leave warpx_executable = None to use its GPU Python kernel.

The tutorial files and software are provided on the instance. No download is needed there. If you do not have a session link, or want to work later on your own machine, keep reading for Docker and Conda instructions.

Files for your own machine

With Git installed, open a terminal and run these commands to download the tutorial folder without the other lessons’ contents:

BASH

git clone --depth 1 --filter=blob:none --sparse \
  https://github.com/BLAST-WarpX/warpx-tutorials.git
cd warpx-tutorials
git sparse-checkout set episodes/files/llnl-hpc-2026
cd episodes/files/llnl-hpc-2026

This is called sparse checkout. It keeps the tutorial’s directory structure and includes files in its parent directories. You do not need a packaging script. The folder contains the notebooks, input files, plotting helpers and optional converter. Open htu/wakefield.ipynb or htu/htu_transport.ipynb. Simulation output files are created when you run the examples.

These commands download the files; use one of the installation options below to get the software. AWS and Docker users already have the tutorial files.

Alternative A: run the tutorial Docker image yourself

The image contains WarpX, ImpactX and the Python analysis tools, served as a browser-based JupyterLab session. To run it locally:

BASH

docker pull ghcr.io/blast-warpx/warpx-tutorials/tutorial:latest
docker run --rm -p 127.0.0.1:3000:3000 ghcr.io/blast-warpx/warpx-tutorials/tutorial:latest

Open http://localhost:3000/lab, launch a Terminal, and go to:

BASH

cd ~/warpx-tutorials/episodes/files/llnl-hpc-2026

For a local NVIDIA GPU with the NVIDIA Container Toolkit installed, add --gpus all to the docker run command. Select the notebook kernels as above. In a terminal, switch to the GPU environment with:

BASH

source /opt/venv-gpu/bin/activate

The files are included in the image. See the tutorial container README for GPU access and other image details.

Alternative B: install locally with Conda

Use this if you want to run the tutorial on your own machine outside of the Docker image – e.g. a laptop with no Docker available, or a cluster where you’d rather use your own Conda installation.

We assume you have a Conda installation available on your machine. Download the Conda environment file, open a terminal in the directory containing it, and create the environment:

BASH

conda env create -f environment_setup.yml

Activate it with:

BASH

conda activate llnlhpc26-warpx-tutorial

The environment installs WarpX and ImpactX from Conda Forge, together with the Python packages used by the simulations and the analysis notebook. This tutorial has been tested with WarpX 26.09 and ImpactX 26.09. GPU execution requires a CUDA-enabled WarpX build; installing this environment alone does not establish GPU support. The warm-up and ImpactX example also run on CPU.

OUTPUT

name: llnlhpc26-warpx-tutorial
channels:
  - conda-forge
dependencies:
  - python=3.14
  - warpx=26.09
  - impactx=26.09
  - openmpi
  - numpy
  - scipy
  - pandas
  - matplotlib
  - jupyterlab
  - openpmd-viewer

When you have finished working, deactivate the environment:

BASH

conda deactivate

To remove it completely at a later date, run:

BASH

conda env remove -n llnlhpc26-warpx-tutorial

Tutorial 1: two-stream instability warm-up

What we will simulate

Two groups of electrons move through each other in opposite directions at one tenth of the speed of light ($\pm 0.1c$, where $c$ is the speed of light). We follow motion along one direction, $z$. The simulation box is periodic: an electron leaving one end re-enters at the other.

Small variations in electron density create an electric field. That field changes the electron motion, which can make the density variations grow. This feedback is the two-stream instability. As the wave grows, it can trap electrons, making them oscillate within the wave. Eventually the rapid growth levels off, or saturates.

All settings are supplied. This small run lets you practice running a simulation and reading its plots before trying the 3D example. For a closer look at the physics and how to choose the settings, see the standalone Two-Stream Instability episode.

Download the files

If you are running from the tutorial Docker image, this is already present at ~/warpx-tutorials/episodes/files/llnl-hpc-2026/two_stream_instability/ – cd there and skip to Read the input file.

Otherwise, create a two_stream_instability directory and place the following files in it:

Read the input file

The Python script below loads the settings and starts WarpX:

OUTPUT

#!/usr/bin/env python3
"""Run the 1D two-stream instability warm-up from a raw WarpX input file."""

from pywarpx import warpx

warpx.load_inputs_file("two_stream_instability_input.txt")
warpx.step()
warpx.finalize()

View the WarpX input file

OUTPUT

# Two-stream instability warm-up: a 1D electrostatic PIC problem, used here
# as a quick, cheap first run before the heavier LWFA and beamline stages of
# this tutorial. Two cold-ish electron populations drift through each other
# at +/-beta0*c; the counter-streaming free energy grows an electrostatic
# wave until the beams trap each other and saturate.
#
# All parameters below are fixed (unlike the companion "A Two-stream
# Instability" episode, where you choose them yourself) so this runs in a
# few seconds and gets everyone to a working WarpX + Python + Jupyter
# pipeline before the LWFA stage.

####################
### MY CONSTANTS ###
####################
# PLASMA
my_constants.n0 = 1.e17          # electron density of each beam [m^-3]
my_constants.T0 = 5.*q_e         # thermal energy of each beam [J], from 5 eV
my_constants.v_te = sqrt(T0 / m_e)
my_constants.beta0 = 0.1         # drift velocity of each beam, in units of c
my_constants.omega_pe = sqrt(n0*q_e**2/(m_e*epsilon0))

# BOX: sized to fit ~16 wavelengths of the fastest-growing mode,
# lambda ~ 2*pi*beta0*(c/omega_pe)
my_constants.Lx = 10.*clight/omega_pe
my_constants.nx = 256
my_constants.dx = Lx/nx

# TIME: long enough for the instability to grow and saturate
my_constants.cfl = 0.9
my_constants.T = 100./omega_pe
my_constants.dt = cfl * dx / clight
my_constants.nt = floor(T/dt)

##########################
### GENERAL PARAMETERS ###
##########################
stop_time = T
amr.n_cell = nx nx nx
amr.max_level = 0
geometry.dims = 1
geometry.prob_lo = -0.5*Lx -0.5*Lx  -0.5*Lx
geometry.prob_hi =  0.5*Lx  0.5*Lx   0.5*Lx

##########################
### BOUNDARY CONDITION ###
##########################
boundary.field_lo = periodic periodic periodic
boundary.field_hi = periodic periodic periodic
boundary.particle_lo = periodic periodic periodic
boundary.particle_hi = periodic periodic periodic

################
### NUMERICS ###
################
warpx.cfl = cfl
algo.maxwell_solver = yee
algo.particle_shape = 3
algo.particle_pusher = boris
warpx.use_filter = 1

#################
### PARTICLES ###
#################
particles.species_names = ele1 ele2

ele1.species_type = electron
ele1.injection_style = NRandomPerCell
ele1.num_particles_per_cell = 200
ele1.profile = constant
ele1.density = n0
ele1.momentum_distribution_type = maxwell_boltzmann
ele1.theta_distribution_type = constant
ele1.theta = T0 / (m_e * clight**2)
ele1.beta_distribution_type = parser
ele1.beta_function(x,y,z) = beta0
ele1.bulk_vel_dir = +z

ele2.species_type = electron
ele2.injection_style = NRandomPerCell
ele2.num_particles_per_cell = 200
ele2.profile = constant
ele2.density = n0
ele2.momentum_distribution_type = maxwell_boltzmann
ele2.theta_distribution_type = constant
ele2.theta = T0 / (m_e * clight**2)
ele2.beta_distribution_type = constant
ele2.beta = beta0
ele2.bulk_vel_dir = -z

###################
### DIAGNOSTICS ###
###################
# FULL
diagnostics.diags_names = particles
particles.intervals = floor(nt/200)
particles.diag_type = Full
particles.species = ele1 ele2
particles.fields_to_plot = none
particles.format = openpmd
particles.openpmd_backend = bp
particles.dump_last_timestep = 1
particles.ele1.variables = w z uz
particles.ele2.variables = w z uz

# REDUCED
warpx.reduced_diags_names = FieldEnergy FieldMaximum
FieldEnergy.type = FieldEnergy
FieldEnergy.intervals = 1
FieldMaximum.type = FieldMaximum
FieldMaximum.intervals = 1

The labels ele1 and ele2 distinguish the two electron groups. They have the same density and temperature, but move in opposite directions. Their random thermal motion is small compared with their initial streaming speed.

We assume a stationary positive background balances the initially uniform negative charge. The input does not track ions explicitly, and the electric field starts at zero. The saved diagnostics are simply measurements from the simulation: electron positions and momenta, plus field energies and maximum field values. The notebook uses these to show how the instability grows.

Run WarpX

BASH

python run_two_stream_instability.py

(If you are running locally with Conda, activate the environment first: conda activate llnlhpc26-warpx-tutorial.)

This small 1D simulation should finish in seconds on a single CPU process. When it is done, your directory should contain a diags/ subfolder with the saved particle data and field measurements.

Analyze the results

If you are running from the tutorial Docker image, the notebook is already present at ~/warpx-tutorials/episodes/files/llnl-hpc-2026/two_stream_instability/two_stream_instability_plots.ipynb. Otherwise, download the analysis notebook into the two_stream_instability directory, alongside the input file and driver script.

Open it with Jupyter:

BASH

jupyter lab two_stream_instability_plots.ipynb

You can also preview the warm-up notebook in your browser.

What to look for 📊

A phase-space plot shows position on one axis and momentum on the other. Here the horizontal coordinate is $z$, and the vertical coordinate is $u_z = p_z/(m_e c)$: electron momentum along $z$, divided by the electron mass and the speed of light. Positive and negative values indicate opposite directions of motion. At these initial speeds, $u_z$ is approximately $v_z/c$.

Electron position versus longitudinal momentum: the two initially separate streams have rolled into cat-eye-shaped vortices as electrons become trapped in the wave.

This reference plot is reused from the Two-Stream Instability episode; your run may show different fine details. Initially, the two groups form nearly horizontal bands. As the instability grows, the bands bend and roll into the loops often called cat-eye vortices. These are loops in position–momentum space, not circular paths in the physical simulation box.

The notebook plots the last saved particle snapshot and the electric and magnetic field energies over time. The energy plot uses a logarithmic vertical axis: exponential growth appears as a straight rising section. Look for the electric field energy to rise and then level off.

💡 Try reading the plots together: the phase-space plot shows how the electrons move; the field-energy plot shows how much energy the growing wave stores. Can you identify evidence of the instability in both?

Once both plots work, continue to the laser-driven wake below.

Tutorial 2: watch a wake in WarpX

A plasma contains free electrons and ions. A strong laser pulse pushes electrons away from its path, leaving the heavier ions behind. The electrons are pulled back, creating a wake behind the pulse. Its electric field can accelerate electrons, rather like a water wave carrying a surfer. This is laser-wakefield acceleration (LWFA).

Three electron-density slices from the reduced WarpX run show a cavity forming and evolving behind the laser at 52, 103 and 155 femtoseconds.

Snapshots from the completed reduced 3D simulation. Read left to right: the wake forms, crosses the density drop, and continues through the lower-density plasma. All panels use the same density scale. The horizontal coordinate $z-ct$ moves at the speed of light, keeping the wake in view; dark regions contain fewer plasma electrons. A femtosecond is $10^{-15}$ seconds.

Open wakefield.ipynb from the htu folder. Run the simulation, then plot a central slice of electron density, the electric field along the direction of travel and the energy spectrum (how many electrons have each energy). Choose a different saved step to watch the wake develop and pass through the density drop. The notebook also puts three density snapshots side by side so the evolution is visible without changing the step manually.

The small 3D model uses a short laser and a dense, short plasma target. It illustrates wakefield dynamics, not the measured HTU beam. It does not need to produce a 100 MeV bunch for the next example.

Files and run commands

Keep these files together (the sparse-checkout commands above do this for you):

From htu/lwfa_warpx/, in the GPU Python environment:

BASH

python run_lwfa_warpx.py
python plot_wakefield.py diags/diag1 --evolution --output wake-evolution.png
python plot_wakefield.py diags/diag1 --iteration 800 --output wake-snapshot.png

If you built a native executable, you can instead run warpx.3d lwfa_warpx_input.txt (use its full path if necessary). The two plotting commands save a density sequence and one field/spectrum snapshot. The reference run took about 8 minutes 25 seconds on an RTX A2000 8 GB laptop GPU; runtime on the provided AWS instance may differ.

In the input file, try changing a0 (a dimensionless measure of laser strength), or set n_up = n_down to remove the drop in plasma density. Rerun the simulation. What changes in the electron-depleted region, the electric field and the electron energies?

Tutorial 3: independent beam transport in ImpactX

Magnets steer and focus an electron beam. Their effect depends on particle momentum, so a set of magnets that works for one energy may spread out a beam at another energy. Here we keep the magnets fixed and change the beam energy. The PMQ triplet is a group of three permanent-magnet quadrupoles used to focus the beam. The chicane is a sequence of bending magnets that takes the beam along a sideways detour.

Particle-density maps at the source, after the PMQ triplet and in the chicane show the 100 MeV beam remaining much narrower than the 20 MeV beam.

Each row follows one synthetic beam downstream. Read the axis scales: source panels use micrometers; downstream panels use millimeters, with a separate range for each panel. Colors show the amount of electric charge in each small square of the plot, and labels show the charge still in the beam. The large spread of the 20 MeV beam shows that these magnet settings do not suit its energy. The model only removes particles at the magnet openings specified in the input; it does not include every wall of the real beam pipe.

Open htu_transport.ipynb. Two computer-generated beams, at 100 MeV and 20 MeV, pass through the same HTU magnets. Each starts with a Gaussian (bell-shaped) distribution. Neither beam is taken from the WarpX run. The notebook shows cross-sections perpendicular to the direction of travel, the charge remaining and the beam size at successive positions along the beamline.

Files and run commands

ImpactX notebook
ImpactX input: the beam settings, magnets and simulation in one file.

The notebook runs this input at both energies and contains all the analysis: beam-density snapshots, remaining charge and beam size. Change energies_MeV in the notebook to [100, 50] to try another pair.

To run just one simulation from htu/, use:

BASH

python beamline_impactx/input_impactx.py --energy-MeV 100 --output runs/100MeV

Choose a new output directory for each run. Use the notebook to run and plot both beams together; it saves beam_density.png, comparison.png and summary.json in the run folder.

Both beams start with the same size, spread in travel directions and fractional spread in energy. The magnet fields are also kept fixed. Particles are removed when they hit a modeled aperture, the opening through a magnet. The model includes known openings in the PMQs and in the focusing magnets of the undulator (a device with alternating magnetic fields). Other pipe dimensions are unspecified.

These are teaching examples, not measured HTU beams or beams specially prepared for these magnets. We also neglect space charge, the electric repulsion between electrons in the beam.

Charge survival and horizontal and vertical beam sizes versus distance; the 20 MeV curves stop where tracking becomes invalid.

A falling charge curve means particles have hit modeled apertures. If the calculation produces invalid particle coordinates, the plot stops at a dotted line and the final fraction reaching the end is marked unknown. The calculation has broken down at that point; we cannot conclude that all the electrons were physically lost.

Optional: couple WarpX to ImpactX

To connect the examples, we would replace the computer-generated ImpactX beam with electrons from WarpX. This requires a simulation with enough spatial and time resolution to describe the accelerated beam reliably. First inspect the positions and momenta and identify a compact group of electrons moving forward together after leaving the plasma. This group is called a bunch.

The converter uses the coordinate utilities to select electrons above an energy threshold and change coordinate conventions. WarpX saves particles at the same time; ImpactX describes their arrival at a shared position along the beamline, relative to a reference particle. The converter preserves how much real charge each simulated particle represents.

From htu/, for standalone WarpX output:

BASH

python coupling/warpx_to_beamline_impactx.py \
  lwfa_warpx/diags/diag1 warpx_bunch.npz --energy-cut-MeV 20
python beamline_impactx/input_impactx.py --bunch warpx_bunch.npz --output coupled_test

Use the diagnostic path of your actual run. The converter selects the final snapshot by default. Choose the energy cut from the observed bunch. The tiny teaching example may have no suitable electrons above the illustrative 20 MeV cut; the converter then stops. Lowering the cut just to produce a file does not establish a usable bunch.

Converting the file does not adjust the magnets. Their strengths and positions need to suit the incoming beam’s energy, size and spread in travel directions. A low-energy beam sent through magnets set for 100 MeV can hit the openings or move at angles too large for the tracking model to describe reliably. See the coupling instructions.

Parameter	Symbol	Description
`n0`	\(n_0\)	Electron density of each beam \([\mathrm{m^{-3}}]\)
`T0`	\(T_0\)	Thermal temperature of each beam \([\mathrm{eV}]\)
`beta0`	\(\beta_0 = v_d/c\)	Drift velocity of each beam normalized to \(c\)

Parameter	Description
`Lx`	Box length, expressed in units of \(c/\omega_{pe}\)
`nx`	Number of grid cells
`cfl`	CFL number (must be \(< 1\))
`T`	Total simulation time, expressed in units of \(\omega_{pe}^{-1}\)

Parameter	Symbol	Description
`FWHM_I`	\(\tau_{\mathrm{FWHM}}\)	Pulse duration (FWHM of intensity) \([\mathrm{s}]\)
`laser_waist`	\(w_0\)	Laser spot size (beam waist) \([\mathrm{m}]\)
`laser_wavelength`	\(\lambda\)	Laser wavelength \([\mathrm{m}]\)
`a0`	\(a_0\)	Normalized vector potential (dimensionless laser amplitude)

Parameter	Symbol	Description
`ne0`	\(n_e / n_c\)	Electron density of the main target, in units of the critical density \(n_c\)
`Te0`	\(T_e\)	Initial electron temperature \([\mathrm{eV}]\)
`thick_targ`		Thickness of the main (boron) target \([\mathrm{m}]\)
`thick_cont`		Thickness of the hydrogen contaminant layer \([\mathrm{m}]\)

All in One View

Overview

Questions

Objectives

Overview of PIC

Features and applications of WarpX

Overview

Questions

Objectives

Basic dependencies

Via Conda-Forge

BASH

BASH

BASH

BASH

From source

Overview

Questions

Objectives

Introduction

Step 1: What do we want to simulate?

Step 2: Browse the WarpX examples

Find a starting point

Show me the solution

Step 3: Start from the uniform plasma example

OUTPUT

Step 4: Add a drift velocity

What parameters do we need to check?

Show me the solution

Define nt as a constant

Switch diagnostics to openPMD

Step 5: Run – and nothing happens

BASH

PYTHON

Step 6: Iterate

Increase the number of timesteps

Reduce the temperature

Switch to 1D

Increase resolution and particles per cell

Put it all together

Show me the solution

Step 7: Visualize the result

PYTHON

BASH

What’s next?

A note on timescales

Overview

Questions

Objectives

Introduction

Setup

Let’s take a look at the input file

OUTPUT

Choosing the parameters

Choose the physical parameters

Hints

Choose the numerical parameters

Hints

Notable details

Run

Let’s run

BASH

Visualize

With Python 🐍

With a Jupyter notebook 📓

Let’s take a look at the Jupyter notebook

Questions for analysis

Show me the solution

Gallery

Overview

Questions

Objectives

Introduction

Setup

Let’s take a look at the input file

OUTPUT

Choosing the parameters

Choose the physical parameters

Hints

Choose the numerical parameters

Define `nt` as a constant