Beam-beam collision
This example shows how to simulate the collision between two ultra-relativistic particle beams. This is representative of what happens at the interaction point of a linear collider. We consider a right-propagating electron bunch colliding against a left-propagating positron bunch.
We turn on the Quantum Synchrotron QED module for photon emission (also known as beamstrahlung in the collider community) and the Breit-Wheeler QED module for the generation of electron-positron pairs (also known as coherent pair generation in the collider community).
To solve for the electromagnetic field we use the nodal version of the electrostatic relativistic solver.
This solver computes the average velocity of each species, and solves the corresponding relativistic Poisson equation (see the WarpX documentation for warpx.do_electrostatic = relativistic for more detail). This solver accurately reproduced the subtle cancellation that occur for some component of the E + v x B terms which are crucial in simulations of relativistic particles.
This example is based on the following paper Yakimenko et al. [9].
Run
The PICMI input file is not available for this example yet.
For MPI-parallel runs, prefix these lines with mpiexec -n 4 ... or srun -n 4 ..., depending on the system.
#################################
########## MY CONSTANTS #########
#################################
my_constants.mc2 = m_e*clight*clight
my_constants.nano = 1.0e-9
my_constants.GeV = q_e*1.e9
# BEAMS
my_constants.beam_energy = 125.*GeV
my_constants.beam_uz = beam_energy/(mc2)
my_constants.beam_charge = 0.14*nano
my_constants.sigmax = 10*nano
my_constants.sigmay = 10*nano
my_constants.sigmaz = 10*nano
my_constants.beam_uth = 0.1/100.*beam_uz
my_constants.n0 = beam_charge / (q_e * sigmax * sigmay * sigmaz * (2.*pi)**(3./2.))
my_constants.omegab = sqrt(n0 * q_e**2 / (epsilon0*m_e))
my_constants.mux = 0.0
my_constants.muy = 0.0
my_constants.muz = -0.5*Lz+3.2*sigmaz
# BOX
my_constants.Lx = 100.0*clight/omegab
my_constants.Ly = 100.0*clight/omegab
my_constants.Lz = 180.0*clight/omegab
# for a full scale simulation use: nx, ny, nz = 512, 512, 1024
my_constants.nx = 64
my_constants.ny = 64
my_constants.nz = 128
# TIME
my_constants.T = 0.7*Lz/clight
my_constants.dt = sigmaz/clight/10.
# DIAGS
my_constants.every_red = 1.
warpx.used_inputs_file = warpx_used_inputs.txt
#################################
####### GENERAL PARAMETERS ######
#################################
stop_time = T
amr.n_cell = nx ny nz
amr.max_grid_size = 128
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
#################################
######## BOUNDARY CONDITION #####
#################################
boundary.field_lo = PEC PEC PEC
boundary.field_hi = PEC PEC PEC
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.load_balance_intervals=100
algo.particle_pusher = vay
#################################
########### PARTICLES ###########
#################################
particles.species_names = beam1 beam2 pho1 pho2 ele1 pos1 ele2 pos2
particles.photon_species = pho1 pho2
beam1.species_type = electron
beam1.injection_style = NUniformPerCell
beam1.num_particles_per_cell_each_dim = 1 1 1
beam1.profile = parse_density_function
beam1.density_function(x,y,z) = "n0 * exp(-(x-mux)**2/(2*sigmax**2)) * exp(-(y-muy)**2/(2*sigmay**2)) * exp(-(z-muz)**2/(2*sigmaz**2))"
beam1.density_min = n0 / 1e3
beam1.momentum_distribution_type = gaussian
beam1.uz_m = beam_uz
beam1.uy_m = 0.0
beam1.ux_m = 0.0
beam1.ux_th = beam_uth
beam1.uy_th = beam_uth
beam1.uz_th = beam_uth
beam1.initialize_self_fields = 1
beam1.self_fields_required_precision = 5e-10
beam1.self_fields_max_iters = 10000
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 = NUniformPerCell
beam2.num_particles_per_cell_each_dim = 1 1 1
beam2.profile = parse_density_function
beam2.density_function(x,y,z) = "n0 * exp(-(x-mux)**2/(2*sigmax**2)) * exp(-(y-muy)**2/(2*sigmay**2)) * exp(-(z+muz)**2/(2*sigmaz**2))"
beam2.density_min = n0 / 1e3
beam2.momentum_distribution_type = gaussian
beam2.uz_m = -beam_uz
beam2.uy_m = 0.0
beam2.ux_m = 0.0
beam2.ux_th = beam_uth
beam2.uy_th = beam_uth
beam2.uz_th = beam_uth
beam2.initialize_self_fields = 1
beam2.self_fields_required_precision = 5e-10
beam2.self_fields_max_iters = 10000
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 = 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
ele1.species_type = electron
ele1.injection_style = none
ele1.self_fields_required_precision = 1e-11
ele1.self_fields_max_iters = 10000
ele1.do_qed_quantum_sync = 1
ele1.qed_quantum_sync_phot_product_species = pho1
ele1.do_classical_radiation_reaction = 0
pos1.species_type = positron
pos1.injection_style = none
pos1.self_fields_required_precision = 1e-11
pos1.self_fields_max_iters = 10000
pos1.do_qed_quantum_sync = 1
pos1.qed_quantum_sync_phot_product_species = pho1
pos1.do_classical_radiation_reaction = 0
ele2.species_type = electron
ele2.injection_style = none
ele2.self_fields_required_precision = 1e-11
ele2.self_fields_max_iters = 10000
ele2.do_qed_quantum_sync = 1
ele2.qed_quantum_sync_phot_product_species = pho2
ele2.do_classical_radiation_reaction = 0
pos2.species_type = positron
pos2.injection_style = none
pos2.self_fields_required_precision = 1e-11
pos2.self_fields_max_iters = 10000
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
#################################
############# QED ###############
#################################
qed_qs.photon_creation_energy_threshold = 0.
qed_qs.lookup_table_mode = builtin
qed_qs.chi_min = 1.e-3
qed_bw.lookup_table_mode = builtin
qed_bw.chi_min = 1.e-2
# for accurate results use the generated tables with
# the following parameters
# note: must compile with -DWarpX_QED_TABLE_GEN=ON
#qed_qs.lookup_table_mode = generate
#qed_bw.lookup_table_mode = generate
#qed_qs.tab_dndt_chi_min=1e-3
#qed_qs.tab_dndt_chi_max=2e3
#qed_qs.tab_dndt_how_many=512
#qed_qs.tab_em_chi_min=1e-3
#qed_qs.tab_em_chi_max=2e3
#qed_qs.tab_em_chi_how_many=512
#qed_qs.tab_em_frac_how_many=512
#qed_qs.tab_em_frac_min=1e-12
#qed_qs.save_table_in=my_qs_table.txt
#qed_bw.tab_dndt_chi_min=1e-2
#qed_bw.tab_dndt_chi_max=2e3
#qed_bw.tab_dndt_how_many=512
#qed_bw.tab_pair_chi_min=1e-2
#qed_bw.tab_pair_chi_max=2e3
#qed_bw.tab_pair_chi_how_many=512
#qed_bw.tab_pair_frac_how_many=512
#qed_bw.save_table_in=my_bw_table.txt
warpx.do_qed_schwinger = 0.
#################################
######### DIAGNOSTICS ###########
#################################
# FULL
diagnostics.diags_names = diag1
diag1.intervals = 0
diag1.diag_type = Full
diag1.write_species = 1
diag1.fields_to_plot = Ex Ey Ez Bx By Bz rho_beam1 rho_beam2 rho_ele1 rho_pos1 rho_ele2 rho_pos2
diag1.format = openpmd
diag1.dump_last_timestep = 1
diag1.species = pho1 pho2 ele1 pos1 ele2 pos2 beam1 beam2
# REDUCED
warpx.reduced_diags_names = ParticleNumber ColliderRelevant_beam1_beam2
ColliderRelevant_beam1_beam2.type = ColliderRelevant
ColliderRelevant_beam1_beam2.intervals = every_red
ColliderRelevant_beam1_beam2.species = beam1 beam2
ParticleNumber.type = ParticleNumber
ParticleNumber.intervals = every_red
Visualize
The figure below shows the number of photons emitted per beam particle (left) and the number of secondary pairs generated per beam particle (right).
We compare different results:
* (red) simplified WarpX simulation as the example stored in the directory /Examples/Physics_applications/beam-beam_collision;
* (blue) large-scale WarpX simulation (high resolution and ad hoc generated tables ;
* (black) literature results from Yakimenko et al. [9].
The small-scale simulation has been performed with a resolution of nx = 64, ny = 64, nz = 128 grid cells, while the large-scale one has a much higher resolution of nx = 512, ny = 512, nz = 1024. Moreover, the large-scale simulation uses dedicated QED lookup tables instead of the builtin tables. To generate the tables within WarpX, the code must be compiled with the flag -DWarpX_QED_TABLE_GEN=ON. For the large-scale simulation we have used the following options:
qed_qs.lookup_table_mode = generate
qed_bw.lookup_table_mode = generate
qed_qs.tab_dndt_chi_min=1e-3
qed_qs.tab_dndt_chi_max=2e3
qed_qs.tab_dndt_how_many=512
qed_qs.tab_em_chi_min=1e-3
qed_qs.tab_em_chi_max=2e3
qed_qs.tab_em_chi_how_many=512
qed_qs.tab_em_frac_how_many=512
qed_qs.tab_em_frac_min=1e-12
qed_qs.save_table_in=my_qs_table.txt
qed_bw.tab_dndt_chi_min=1e-2
qed_bw.tab_dndt_chi_max=2e3
qed_bw.tab_dndt_how_many=512
qed_bw.tab_pair_chi_min=1e-2
qed_bw.tab_pair_chi_max=2e3
qed_bw.tab_pair_chi_how_many=512
qed_bw.tab_pair_frac_how_many=512
qed_bw.save_table_in=my_bw_table.txt
Fig. 6 Beam-beam collision benchmark against Yakimenko et al. [9].