Virtual photons (#5595)

This PR adds the option of creating a virtual photon population
associated to a given electron or positron species.
The virtual photons are generated from scratch at every timestep.
They are initialized in the same position as the parent particle.
Their momenta is sampled from the theoretical analytical spectrum (see
equation 99.16 in Berestetskii, Pitaevskii, Lifshitz. Quantum
Electrodynamics: Volume 4.)

The virtual photons can be used in linear Breit-Wheeler or linear
Compton collisions. This is used to simulate incoherent beam-beam
effects at the interaction point of colliders.

To do:
- [x] create a virtual photon species with the correct number of
particles and with attributes same as primary particles
- [x] draw the virtual photon energy from the correct probability
distribution
- [x] add infrastructure
- [x] add docs and test

Note: As of this PR, the virtual photons are technically not virtual at
all, as their virtuality is zero. A more accurate name would be
"equivalent quanta". However, follow-up PRs will soon add virtuality
calculations.

🧨

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Luca Fedeli <[email protected]>
Co-authored-by: Edoardo Zoni <[email protected]>
Co-authored-by: Edoardo Zoni <[email protected]>
Co-authored-by: Neïl Zaim <[email protected]>
Co-authored-by: Andrew Myers <[email protected]>
Co-authored-by: Remi Lehe <[email protected]>

Author: 7 people

Docs/source/usage/parameters.rst

@@ -1512,6 +1512,24 @@ Particle initialization
     in each dimension with a matched shape function and filtering used for current deposition.
     This is required when using the electron energy solver with electron-ion temperature relaxation.
 
+* ``<species>.do_qed_virtual_photons`` (`boolean`) optional (default `false`)
+    Create a population of virtual photons associated with ``<species>``.
+    It only works if ``<species>`` is an electron or a positron species.
+    The virtual photon species has to be created as a regular photon species in the input file.
+    Virtual photons are created from scratch at each timestep at the same position as the parent particle.
+    This implies that different primary species must have different virtual photon species.
+    The energy of the virtual photons is sampled from their spectrum (see :cite:t:`LandauVol4` section 99 for more details).
+    The momentum of the virtual photons is parallel to that of the parent particle.
+    This feature also requires the following input parameters:
+      * ``<species>.qed_virtual_photon_species_name`` (`string`) name of the virtual photon species associated with the current lepton species.
+      * ``<virtual_photon_species>.qed_virtual_photon_min_energy`` (`float`, in Joules) minimum energy of the virtual photons
+      * ``<virtual_photon_species>.qed_virtual_photon_multiplier`` (`int`), sampling factor for the virtual photons.
+        A sampling factor of ``f`` means that the number of virtual photons is multiplied by ``f``, while their weights are divided by ``f``.
+    The virtual photons can undergo collisions via the linear Breit-Wheeler or linear Compton processes.
+    This is useful to model incoherent beam-beam effects in colliders (e.g. pair generation, radiative Bhabha scattering).
+    This QED feature is separated from the strong-field QED modules (quantum synchrotron and non-linear Breit-Wheeler).
+    It requires WarpX to be compiled with ``WarpX_QED=ON`` (CMake) or ``QED=TRUE`` (GNU Make).
+
 .. _running-cpp-parameters-fluids:
 
 Cold Relativistic Fluid initialization

Examples/Tests/CMakeLists.txt

@@ -84,3 +84,4 @@ add_subdirectory(single_particle)
 add_subdirectory(space_charge_initialization)
 add_subdirectory(subcycling)
 add_subdirectory(vay_deposition)
+add_subdirectory(virtual_photons)

Examples/Tests/virtual_photons/CMakeLists.txt

@@ -0,0 +1,12 @@
+# Add tests (alphabetical order) ##############################################
+#
+
+add_warpx_test(
+    test_3d_virtual_photons  # name
+    3  # dims
+    1  # nprocs
+    inputs_test_3d_virtual_photons  # inputs
+    analysis_virtual_photons.py  # analysis
+    "analysis_default_regression.py --path diags/diag1"  # checksum
+    OFF  # dependency
+)

Examples/Tests/virtual_photons/analysis_default_regression.py

@@ -0,0 +1 @@
+../../analysis_default_regression.py

Examples/Tests/virtual_photons/analysis_virtual_photons.py

@@ -0,0 +1,111 @@
+#!/usr/bin/env python3
+
+# This test generates the population of virtual photons
+# of one high-energy electron.
+# The total number and spectrum of the virtual photons are
+# compared to the theoretical prediction.
+# Checks that the photons are in the same position of the electron.
+
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy import log
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.constants import alpha, c, eV, m_e, pi
+
+###########################
+### ENERGY AND SPECTRUM ###
+###########################
+
+# useful constants
+GeV = 1e9 * eV
+
+# electron energy
+energy = 125 * GeV
+
+# virtual photons min energy
+hw_min = 1e-12 * m_e * c**2
+
+# min fractional energy of the virtual photon wrt electron energy
+ymin = hw_min / energy
+
+#############
+### WarpX ###
+#############
+
+series = OpenPMDTimeSeries("./diags/diag1/")
+sampling_factor = 1e7
+uz_virtual_photons, w_virtual_photons = series.get_particle(
+    ["uz", "w"], species="virtual_photons", iteration=1
+)
+w_electrons = series.get_particle(["w"], species="beam", iteration=1)
+
+# fractional photon energy (photon energy / electron energy)
+y_warpx = uz_virtual_photons * c / energy
+
+# bins for the fractional photon energy
+y = np.geomspace(ymin, 1, 401)
+
+# number of virtual photons per electron obtained with WarpX
+N_warpx = np.sum(w_virtual_photons) / np.sum(w_electrons)
+
+# spectrum of the virtual photons per electron
+H, b = np.histogram(y_warpx, bins=y, weights=w_virtual_photons)
+db = np.diff(b)
+b = 0.5 * (b[1:] + b[:-1])
+dN_dy_warpx = H / db / np.sum(w_electrons)
+
+##############
+### Theory ###
+##############
+
+y = b
+# spectrum of virtual photons for one electron
+dN_dy_theory = alpha / pi / y * (-2 * log(y))
+# dN_dy[dN_dy < 0] = 0.0
+
+# number of virtual photons for one electron from theory
+N_theory = alpha / pi * log(ymin) ** 2
+
+############
+### Plot ###
+############
+
+fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(5, 4), dpi=200)
+ax.plot(y, dN_dy_theory, color="black", lw=6, label="theory")
+ax.plot(y, dN_dy_warpx, color="dodgerblue", lw=4, label="WarpX")
+ax.legend()
+ax.set_yscale("log")
+ax.set_xscale("log")
+ax.set_xlabel("Fractional photon energy")
+ax.set_ylabel("dN/dy")
+ax.set_title("Virtual photons spectrum")
+fig.savefig("spectrum_virtual_photons.png")
+
+#############
+### Error ###
+#############
+
+number_rel_error = np.abs(N_warpx - N_theory) / N_theory
+spectrum_rel_error = np.abs(dN_dy_warpx - dN_dy_theory) / dN_dy_theory
+
+print("Number of virtual photons per electron:")
+print(f"From simulation : {N_warpx}")
+print(f"From theory     : {N_theory}")
+print(f"Relative error  : {number_rel_error:.4%}")
+
+print("Spectrum of virtual photons per electron:")
+print(f"Max relative error: {spectrum_rel_error.max()}")
+
+assert number_rel_error < 0.02
+assert (spectrum_rel_error < 0.04).all()
+
+################
+### Position ###
+################
+
+x, y, z = series.get_particle(["x", "y", "z"], species="virtual_photons", iteration=1)
+x_e, y_e, z_e = series.get_particle(["x", "y", "z"], species="beam", iteration=1)
+
+assert np.unique(x) == x_e
+assert np.unique(y) == y_e
+assert np.unique(z) == z_e

Examples/Tests/virtual_photons/inputs_test_3d_virtual_photons

@@ -0,0 +1,104 @@
+#################################
+########## MY CONSTANTS #########
+#################################
+my_constants.mc2_eV = m_e*clight*clight/q_e
+
+# BEAMS
+my_constants.beam_energy_eV = 125.e9
+my_constants.beam_gamma = beam_energy_eV/(mc2_eV)
+my_constants.beam_N = 1.2e10
+my_constants.sigmax = 500e-9
+my_constants.sigmay = 10e-9
+my_constants.sigmaz = 300e-3
+my_constants.muz = 4*sigmaz
+my_constants.nmacropart = 1
+
+# VIRTUAL PHOTONS
+# minimum energy
+my_constants.hwmin = 1e-12*m_e*clight*clight
+
+# BOX
+my_constants.Lx = 8*sigmax
+my_constants.Ly = 8*sigmay
+my_constants.Lz = 16*sigmaz
+my_constants.nx = 32
+my_constants.ny = 32
+my_constants.nz = 32
+
+# TIME
+my_constants.dt = sigmaz/clight/10.
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+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
+warpx.numprocs = 1 1 1
+
+#################################
+######## BOUNDARY CONDITION #####
+#################################
+boundary.particle_lo = Absorbing Absorbing Absorbing
+boundary.particle_hi = Absorbing Absorbing Absorbing
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.maxwell_solver = none
+warpx.const_dt = dt
+algo.particle_shape = 3
+algo.particle_pusher = vay
+
+#################################
+########### PARTICLES ###########
+#################################
+particles.species_names = beam virtual_photons
+particles.photon_species = virtual_photons
+
+beam.species_type = electron
+beam.injection_style = gaussian_beam
+beam.x_rms = sigmax
+beam.y_rms = sigmay
+beam.z_rms = sigmaz
+beam.x_m = 0
+beam.y_m = 0
+beam.z_m = -muz
+beam.npart = nmacropart
+beam.q_tot = -beam_N*q_e
+beam.z_cut = 4
+beam.momentum_distribution_type = gaussian
+beam.uz_m = beam_gamma
+beam.uy_m = 0.0
+beam.ux_m = 0.0
+beam.ux_th = 0.0
+beam.uy_th = 0.0
+beam.uz_th = 0
+beam.do_not_deposit = 1
+beam.do_not_push = 1
+beam.do_qed_virtual_photons = 1
+beam.qed_virtual_photon_species_name = virtual_photons
+
+virtual_photons.species_type = photon
+virtual_photons.injection_style = none
+virtual_photons.do_not_push = 1
+virtual_photons.do_not_gather = 1
+virtual_photons.do_not_deposit = 1
+virtual_photons.qed_virtual_photons_min_energy = hwmin
+virtual_photons.qed_virtual_photons_multiplier = 10000000
+
+#################################
+######### DIAGNOSTICS ###########
+#################################
+# FULL
+diagnostics.diags_names = diag1
+
+diag1.intervals = -1
+diag1.diag_type = Full
+diag1.write_species = 1
+diag1.species = beam virtual_photons
+diag1.format = openpmd
+diag1.virtual_photons.variables = uz w x y z

Regression/Checksum/benchmarks_json/test_3d_virtual_photons.json

@@ -0,0 +1,29 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 0.0,
+    "Ey": 0.0,
+    "Ez": 0.0,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0
+  },
+  "beam": {
+    "particle_position_x": 2.748730894777248e-07,
+    "particle_position_y": 1.5827522919751402e-08,
+    "particle_position_z": 1.122721937612553,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 6.680357490847885e-17,
+    "particle_weight": 12000000000.0
+  },
+  "virtual_photons": {
+    "particle_position_x": 10.23535247755454,
+    "particle_position_y": 0.5893638996747153,
+    "particle_position_z": 41806401.59275659,
+    "particle_momentum_z": 0.0,
+    "particle_weight": 44683977600.0
+  }
+}