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-04-06 | 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

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-04-08 | 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.

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