Beam-beam simulations at the interaction point

This section contains example scripts for simulations of beam-beam crossing at the interaction point using WarpX and Guinea-Pig. As an arbitrary choice, we selected the \(C^3\) collider parameters at 125 GeV.

Overview

This simulates the collision between a beam of electrons and a beam of positrons traveling at 125 GeV. As usual, \(x, y, z\) are the horizontal, vertical, and longitudinal coordinates, respectively. The beams are initialized such that their centroids are \(4 \sigma_z\) away from interaction point.

Beamstrahlung and coherent pair production (though not actually important for these beam parameters) are activated. Note that incoherent processes are not enabled in either WarpX and Guinea-Pig. We will add them in future versions of this example.

WarpX

Resources

Some useful WarpX’s links if you get stuck or want to know more:

Build

WarpX is open-source and available on GitHub. You can install WarpX in several different ways (for example from source, using Conda, Spack), and on all the major operative systems (Linux, Windows, MacOS). You can compile it on CPU or on GPU, depending on your machine. Find the instructions here.

Run

Download or copy the input file below.

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

# 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.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 = 256
my_constants.ny = 256
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)

##############################
##### GENERAL PARAMETERS #####
##############################
stop_time = T
amr.n_cell = nx ny nz
amr.max_grid_size = nz
amr.blocking_factor = 2
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

####################
##### 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

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

#####################
##### PARTICLES #####
#####################
particles.species_names = beam1 beam2 pho1 pho2 ele_nlbw1 pos_nlbw1 ele_nlbw2 pos_nlbw2
particles.photon_species = pho1 pho2

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

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

pho1.species_type = photon 
pho1.injection_style = none 
pho1.do_qed_breit_wheeler = 1
pho1.qed_breit_wheeler_ele_product_species = ele_nlbw1
pho1.qed_breit_wheeler_pos_product_species = pos_nlbw1

pho2.species_type = photon 
pho2.injection_style = none 
pho2.do_qed_breit_wheeler = 1
pho2.qed_breit_wheeler_ele_product_species = ele_nlbw2
pho2.qed_breit_wheeler_pos_product_species = pos_nlbw2

ele_nlbw1.species_type = electron 
ele_nlbw1.injection_style = none 
ele_nlbw1.do_qed_quantum_sync = 1
ele_nlbw1.qed_quantum_sync_phot_product_species = pho1
ele_nlbw1.do_classical_radiation_reaction = 0

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

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

pos_nlbw2.species_type = positron
pos_nlbw2.injection_style = none
pos_nlbw2.do_qed_quantum_sync = 1
pos_nlbw2.qed_quantum_sync_phot_product_species = pho2
pos_nlbw2.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 = bound trajs

bound.dump_last_timestep = 1
bound.diag_type = BoundaryScraping
bound.format = openpmd
bound.openpmd_backend = bp
bound.intervals = 1

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

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

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

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

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

pos_nlbw2.save_particles_at_xlo = 1
pos_nlbw2.save_particles_at_ylo = 1
pos_nlbw2.save_particles_at_zlo = 1
pos_nlbw2.save_particles_at_xhi = 1
pos_nlbw2.save_particles_at_yhi = 1
pos_nlbw2.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

trajs.intervals = 32
trajs.diag_type = Full
trajs.species = beam1 beam2 pho1 pho2 ele_nlbw1 pos_nlbw1 ele_nlbw2 pos_nlbw2
trajs.fields_to_plot = none
trajs.format = openpmd
trajs.openpmd_backend = h5
trajs.dump_last_timestep = 1

# REDUCED
warpx.reduced_diags_names = ParticleEnergy RhoMaximum ParticleNumber CollRel_beam1_beam2 DiffLumi_beam1_beam2 DiffLumi_beam1_ele2 DiffLumi_beam1_pos2 DiffLumi_beam1_pho2 DiffLumi_ele1_beam2 DiffLumi_ele1_ele2 DiffLumi_ele1_pos2 DiffLumi_ele1_pho2 DiffLumi_pos1_beam2 DiffLumi_pos1_ele2 DiffLumi_pos1_pos2 DiffLumi_pos1_pho2 DiffLumi_pho1_beam2 DiffLumi_pho1_ele2 DiffLumi_pho1_pos2 DiffLumi_pho1_pho2 

ParticleEnergy.type = ParticleEnergy
ParticleEnergy.intervals = 1

RhoMaximum.type = RhoMaximum
RhoMaximum.intervals = 1

ParticleNumber.type = ParticleNumber
ParticleNumber.intervals = 1

CollRel_beam1_beam2.type = ColliderRelevant
CollRel_beam1_beam2.intervals = 1
CollRel_beam1_beam2.species = beam1 beam2

DiffLumi_beam1_beam2.type = DifferentialLuminosity
DiffLumi_beam1_ele2.type  = DifferentialLuminosity
DiffLumi_beam1_pos2.type  = DifferentialLuminosity 
DiffLumi_beam1_pho2.type  = DifferentialLuminosity 
DiffLumi_ele1_beam2.type  = DifferentialLuminosity 
DiffLumi_ele1_ele2.type   = DifferentialLuminosity 
DiffLumi_ele1_pos2.type   = DifferentialLuminosity 
DiffLumi_ele1_pho2.type   = DifferentialLuminosity 
DiffLumi_pos1_beam2.type  = DifferentialLuminosity 
DiffLumi_pos1_ele2.type   = DifferentialLuminosity 
DiffLumi_pos1_pos2.type   = DifferentialLuminosity
DiffLumi_pos1_pho2.type   = DifferentialLuminosity
DiffLumi_pho1_beam2.type  = DifferentialLuminosity 
DiffLumi_pho1_ele2.type   = DifferentialLuminosity 
DiffLumi_pho1_pos2.type   = DifferentialLuminosity 
DiffLumi_pho1_pho2.type   = DifferentialLuminosity 


DiffLumi_beam1_beam2.bin_max = 2.1*beam_energy/q_e 
DiffLumi_beam1_ele2.bin_max  = 2.1*beam_energy/q_e 
DiffLumi_beam1_pos2.bin_max  = 2.1*beam_energy/q_e  
DiffLumi_beam1_pho2.bin_max  = 2.1*beam_energy/q_e  
DiffLumi_ele1_beam2.bin_max  = 2.1*beam_energy/q_e  
DiffLumi_ele1_ele2.bin_max   = 2.1*beam_energy/q_e  
DiffLumi_ele1_pos2.bin_max   = 2.1*beam_energy/q_e  
DiffLumi_ele1_pho2.bin_max   = 2.1*beam_energy/q_e  
DiffLumi_pos1_beam2.bin_max  = 2.1*beam_energy/q_e  
DiffLumi_pos1_ele2.bin_max   = 2.1*beam_energy/q_e  
DiffLumi_pos1_pos2.bin_max   = 2.1*beam_energy/q_e 
DiffLumi_pos1_pho2.bin_max   = 2.1*beam_energy/q_e 
DiffLumi_pho1_beam2.bin_max  = 2.1*beam_energy/q_e  
DiffLumi_pho1_ele2.bin_max   = 2.1*beam_energy/q_e  
DiffLumi_pho1_pos2.bin_max   = 2.1*beam_energy/q_e  
DiffLumi_pho1_pho2.bin_max   = 2.1*beam_energy/q_e  

DiffLumi_beam1_beam2.bin_min = 0. 
DiffLumi_beam1_ele2.bin_min = 0. 
DiffLumi_beam1_pos2.bin_min = 0.  
DiffLumi_beam1_pho2.bin_min = 0.  
DiffLumi_ele1_beam2.bin_min = 0.  
DiffLumi_ele1_ele2.bin_min = 0.
DiffLumi_ele1_pos2.bin_min = 0.  
DiffLumi_ele1_pho2.bin_min = 0.  
DiffLumi_pos1_beam2.bin_min = 0.  
DiffLumi_pos1_ele2.bin_min = 0.  
DiffLumi_pos1_pos2.bin_min = 0. 
DiffLumi_pos1_pho2.bin_min = 0. 
DiffLumi_pho1_beam2.bin_min = 0.  
DiffLumi_pho1_ele2.bin_min = 0.  
DiffLumi_pho1_pos2.bin_min = 0.  
DiffLumi_pho1_pho2.bin_min = 0.

DiffLumi_beam1_beam2.species = beam1 beam2
DiffLumi_beam1_ele2.species = beam1 ele_nlbw2
DiffLumi_beam1_pos2.species = beam1 pos_nlbw2
DiffLumi_beam1_pho2.species = beam1 pho2
DiffLumi_ele1_beam2.species = ele_nlbw1 beam2
DiffLumi_ele1_ele2.species= ele_nlbw1 ele_nlbw2
DiffLumi_ele1_pos2.species= ele_nlbw1 pos_nlbw2
DiffLumi_ele1_pho2.species= ele_nlbw1 pho2
DiffLumi_pos1_beam2.species = pos_nlbw1 beam2
DiffLumi_pos1_ele2.species= pos_nlbw1 ele_nlbw2  
DiffLumi_pos1_pos2.species= pos_nlbw1 pos_nlbw2
DiffLumi_pos1_pho2.species= pos_nlbw1 pho2
DiffLumi_pho1_beam2.species = pho1 beam2
DiffLumi_pho1_ele2.species= pho1 ele_nlbw2
DiffLumi_pho1_pos2.species= pho1 pos_nlbw2 
DiffLumi_pho1_pho2.species= pho1 pho2 


DiffLumi_beam1_beam2.bin_number = 256
DiffLumi_beam1_ele2.bin_number = 256
DiffLumi_beam1_pos2.bin_number = 256 
DiffLumi_beam1_pho2.bin_number = 256 
DiffLumi_ele1_beam2.bin_number = 256 
DiffLumi_ele1_ele2.bin_number = 256 
DiffLumi_ele1_pos2.bin_number = 256 
DiffLumi_ele1_pho2.bin_number = 256 
DiffLumi_pos1_beam2.bin_number = 256 
DiffLumi_pos1_ele2.bin_number = 256 
DiffLumi_pos1_pos2.bin_number = 256
DiffLumi_pos1_pho2.bin_number = 256
DiffLumi_pho1_beam2.bin_number = 256 
DiffLumi_pho1_ele2.bin_number = 256 
DiffLumi_pho1_pos2.bin_number = 256 
DiffLumi_pho1_pho2.bin_number = 256 

DiffLumi_beam1_beam2.intervals = nt
DiffLumi_beam1_ele2.intervals = nt 
DiffLumi_beam1_pos2.intervals = nt  
DiffLumi_beam1_pho2.intervals = nt  
DiffLumi_ele1_beam2.intervals = nt  
DiffLumi_ele1_ele2.intervals = nt  
DiffLumi_ele1_pos2.intervals = nt  
DiffLumi_ele1_pho2.intervals = nt  
DiffLumi_pos1_beam2.intervals = nt  
DiffLumi_pos1_ele2.intervals = nt  
DiffLumi_pos1_pos2.intervals = nt 
DiffLumi_pos1_pho2.intervals = nt 
DiffLumi_pho1_beam2.intervals = nt  
DiffLumi_pho1_ele2.intervals = nt  
DiffLumi_pho1_pos2.intervals = nt  
DiffLumi_pho1_pho2.intervals = nt

On your local machine, you can run it with:

mpirun -np <n_ranks> ./<warpx_executable> <input_file>

The simulation takes about ~500 seconds on two cores on an ordinary not-too-old laptop.

Guinea-Pig

Build

First, install FFTW3. Then follow the instructions here.

Run

You can download or copy the input file below

See Guinea-Pig inputs
$ACCELERATOR:: C3_250_PS1
{
energy=125.0;
particles=0.624;
charge_sign=-1;
beta_x=12.0;
beta_y=0.12;
emitt_x=0.9;
emitt_y=0.02;
sigma_z=100.0;
scale_step=1.0;
espread=0.003;
which_espread=3;
dist_z=0;
f_rep=1;
n_b=1;
offset_y=0.;
offset_x=0;
waist_y=0;
}


$PARAMETERS:: params
{
n_z=33;
n_t=1;
n_m=100000;
cut_z=4*sigma_z.1;
n_x=256;
n_y=256;
cut_x=20*sigma_x.1;
cut_y=20*sigma_y.1;
pair_q2=1;
beam_size=1;
grids=7;
store_beam=1;
do_pairs=0;
track_pairs=0;
store_pairs=0;
do_photons=1;
store_photons=1;
do_hadrons=0;
do_jets=0;
do_coherent=1;
electron_ratio=1;
photon_ratio=1;
do_eloss=1;
do_espread=1;
rndm_seed=978360;
rndm_load=0;
rndm_save=1;
do_lumi=7;
automatic_grid_sizing=0;
bmt_precession=0;
pair_ecut=0.;
do_dump=1;
}

To run, do:

./<guineapig_executable> C3_250_PS1 params acc.out

Visualization and postprocessing

You can download this Jupyter notebbok that shows how to extract and make a video of the main particle data.

Soon we will add the luminosities.

Jupyter Notebook

Visualizing the beams with Guinea-Pig and WarpX

import os, sys
import numpy as np
import matplotlib.pyplot as plt
from openpmd_viewer import OpenPMDTimeSeries
from scipy.constants import c, micro, nano, pi

Collider parameters

sigmaz = 100*micro
sigmax = 210*nano
sigmay = 3.1*nano
npart = 6.24e9
n0 = npart / (sigmax * sigmay * sigmaz * (2.*pi)**(3./2.))

Simulation parameters

# time
n_iterations = 256
iterations = range(n_iterations)
# box
nx = 512
ny = 512
nz = 512
Lx = 20*sigmax
Ly = 20*sigmay
Lz = 16*sigmaz
gridx = np.linspace(-0.5*Lx, 0.5*Lx, nx+1)
gridy = np.linspace(-0.5*Ly, 0.5*Ly, ny+1)
gridz = np.linspace(-0.5*Lz, 0.5*Lz, nz+1)
dx = gridx[1]-gridx[0]
dy = gridy[1]-gridy[0]
dz = gridz[1]-gridz[0]
# particles
nmacropart = 1e5
w0 = npart / nmacropart # weight

Function to plot a single Guine-Pig step

def one_step_gp(n, gp_dir):
    '''
    inputs
        n: current timestep
        gp_dir: simulation folder
    outputs
        H_zx, H_zy: density of the beams integrated along y and x, resp
    '''
    global w0, dx, dy, dz, gridx, gridy, gridz

    # get beams' data, columns are:
    # Particle Energy [GeV] | x [um] | y [um] | z [um] | x' [urad] | y' [urad]
    data1 = np.loadtxt(os.path.join(gp_dir, 'b1.%d' % n))
    data2 = np.loadtxt(os.path.join(gp_dir, 'b2.%d' % n))

    # get the number of macroparticles in the beams
    N1 = np.shape(data1)[0]
    N2 = np.shape(data2)[0]

    # stack the data together and convert to SI
    x_data = np.hstack((data1[:,1], data2[:,1]))*micro
    y_data = np.hstack((data1[:,2], data2[:,2]))*micro
    z_data = np.hstack((data1[:,3], data2[:,3]))*micro

    weights = np.ones(N1+N2) * w0
    weights[-N2:] = - w0 # assign negative weights to the second beam

    H_zx, bx, bz = np.histogram2d(x_data, z_data, bins=(gridx, gridz), weights=weights)
    H_zy, by, bz = np.histogram2d(y_data, z_data, bins=(gridy, gridz), weights=weights)

    print(np.max(np.abs(H_zx)), np.max(np.abs(H_zy)))

    return H_zx/(dz*dx), H_zy/(dz*dy)

Function to plot a single WarpX step

def one_step_wx(n, wx_series):
    '''
    inputs
        n: current timestep
        wx_dir: simulation folder
    outputs
        H_zx: density of the beams projected on the plane (z,x), integrated along y
        H_zy: density of the beams projected on the plane (z,y), integrated along x
    '''
    global dx, dy, dz, gridx, gridy, gridz

    x1,y1,z1,w1 = wx_series.get_particle(['x','y','z','w'], species='beam1', iteration=n)
    x2,y2,z2,w2 = wx_series.get_particle(['x','y','z','w'], species='beam2', iteration=n)
    w1 = -w1
    X = np.hstack((x1,x2))
    Y = np.hstack((y1,y2))
    Z = np.hstack((z1,z2))
    W = np.hstack((w1,w2))

    H_zx, bx, bz = np.histogram2d(X, Z, bins=(gridx, gridz), weights=W)
    H_zy, by, bz = np.histogram2d(Y, Z, bins=(gridy, gridz), weights=W)

    print(np.max(np.abs(H_zx)), np.max(np.abs(H_zy)))

    return H_zx/(dz*dx), H_zy/(dz*dy)

Generate a video of the colliding beams

!mkdir -p "plots"

gp_dir = "gp"
wx_dir = "wx"

path=os.path.join(wx_dir, 'diags/trajs')
series = OpenPMDTimeSeries(path)

v0 = n0 * np.sqrt(sigmax * sigmay) * 0.05
extent_zx = [gridz[0]/micro, gridz[-1]/micro, gridx[0]/nano, gridx[-1]/nano]
extent_zy = [gridz[0]/micro, gridz[-1]/micro, gridy[0]/nano, gridy[-1]/nano]

plt.rcParams.update({'font.size': 16})

# loop through the timesteps
for n in iterations:

    # prepare canvas
    fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(16,12), dpi=300, sharex='col', sharey='row')

    H_zx, H_zy = one_step_gp(n, gp_dir)
    im=ax[0][0].imshow(H_zx, extent=extent_zx, cmap='seismic', origin='lower', interpolation='nearest', vmin=-v0, vmax=v0)
    im=ax[1][0].imshow(H_zy, extent=extent_zy, cmap='seismic', origin='lower', interpolation='nearest', vmin=-v0, vmax=v0)
    ax[0][0].set_title("Guinea-Pig")


    H_zx, H_zy = one_step_wx(n, series)
    im=ax[0][1].imshow(H_zx, extent=extent_zx, cmap='seismic', origin='lower', interpolation='nearest', vmin=-v0, vmax=v0)
    im=ax[1][1].imshow(H_zy, extent=extent_zy, cmap='seismic', origin='lower', interpolation='nearest', vmin=-v0, vmax=v0)
    ax[0][1].set_title("WarpX")

    fig.subplots_adjust(right=0.85)
    cbar_ax = fig.add_axes([0.88, 0.085, 0.015, 0.87])
    fig.colorbar(im, cax=cbar_ax, label='density [arb. units]')

    for a in ax.reshape(-1):
        a.set_aspect('auto')

    ax[1][0].set_xlabel(r'z [$\mu$m]')
    ax[1][1].set_xlabel(r'z [$\mu$m]')

    ax[0][0].set_ylabel(r'x [nm]')
    ax[1][0].set_ylabel(r'y [nm]')

    plt.savefig(f"plots/img_{n:04d}",dpi=300, bbox_inches='tight')
    plt.close("all")

! ffmpeg -framerate 30 -i 'plots/img_%04d.png' -y out.mp4

from IPython.display import Video, display
display(Video("out.mp4", embed=True))