Confusing number density when initializing particles from an array in cylindrical geometry

[Frostleaf1987]

I hope to create a gaussian number density in cylindrical geometry with a defined peak density from an python array. Therefore, I follow the method in the official tutorials(advanced_beam_driven_wake.py) to set up the array. (The whole code is in the attached file)

input.txt

proton bunch parameters

sigma_x = 20mm # initial longitudinal rms size
sigma_r = 3mm # initial transverse/radial rms size (cylindrical symmetry)
bunch_energy_spread = 0.000 # initial rms energy spread / average energy (not in percent)
bunch_normalized_emittance = 0.00000mm_mrad # initial rms emittance, same emittance for both transverse planes
gamma_bunch = mp_to_megammabeta # gammabeta for proton; initial relativistic Lorentz factor of the bunch
nb= 1e18/ncrit # max density of bunch
npart = 200000 # number of computational macro-particles to model the electron bunch
normalized_species_charge = 1 # For positive muons
weight = 2np.pinp.sqrt(2np.pi)nbsigma_r**2sigma_x/npart

initialize the bunch using numpy arrays

the bunch will have npart particles, so an array of npart elements is used to define the x coordinate of each particle and so on ...

array_position = np.zeros((4,npart)) # positions x,y,z weight
array_momentum = np.zeros((3,npart)) # momenta x,y,z

The proton bunch is supposed at waist. To make it convergent/divergent, transport matrices can be used

array_position[0,:] = np.random.normal(loc=center_bunch_x, scale=sigma_x, size=npart) # generate random number from gaussian distribution for x position
array_position[1,:] = np.random.normal(loc=center_bunch_y, scale=sigma_r, size=npart) # generate random number from gaussian distribution for y position
array_position[2,:] = np.random.normal(loc=center_bunch_y, scale=sigma_r, size=npart)
array_momentum[0,:] = np.random.normal(loc=gamma_bunch, scale=bunch_energy_spreadgamma_bunch, size=npart) # generate random number from gaussian distribution for px position # assumption: bunch defined at waist (zero rms divergence)
array_momentum[1,:] = np.random.normal(loc=0., scale=bunch_normalized_emittance/sigma_rgamma_bunch, size=npart) # generate random number from gaussian distribution for py position
array_momentum[2,:] = np.random.normal(loc=0., scale=bunch_normalized_emittance/sigma_r*gamma_bunch, size=npart) # generate random number from gaussian distribution for pz position

array_position[3,:] = np.multiply(np.ones(npart),weight)

However, the number density obtained from the ParticleBinning diagnostic does not match my expectations.

The related ParticleBinning diagnostic is shown below

DiagParticleBinning(
deposited_quantity = "weight",
every = 1,
time_average = 1,
species = ["bunch"],
axes = [
["moving_x", 0. , Lx , 200 ],
["y" , -Lr , Lr , 200],
]
)

Additionally, when I change the sigmar, I find that the peak density is changed in inverse proportion. Also when I let the weight be direct proportion to sigmar and sigmax, the peak density almost remains unchanged.
I wonder why this is happening and in which way I could initialize the particles as wanted.
update: I have found the weight should be set in direct proportion to sigmar _sigmax nb0 1/omega0(omega0 is the reference omega), but I cannot exactly understand or explain the reason.

[mccoys]

I might not have understood your setup but you apparently set the weight from sizes in mm, whereas everything should be in code units (see the units documentation)

Let me also convert this to a discussion as it is not a bug or a request

[Frostleaf1987]

Sorry, I forgot to mention that all the units like mm have already been defined before in code units . My main point is that since in AMcylindrical geometry the particles is defined in 3D space, according to the definition of weight

, then macro-particle weight=(Number of all the species in code units) / (Number of all the macro-particles)=(2pi)^1.5* sigmax* sigmar^2* npeak/npart. And this formula is supported by and consistent with smilei's official tutorial advanced_beam_driven_wake.py(I have check it over in theory and simulation). However this weight gives a unexpected density, and it seams that the weight should be set in direct proportion to sigmax sigmar npeak omega0^-1 (rather than sigmax* sigmar^2* npeak)to get the expected density.

[Z10Frank]

Hello, if some context info can help, I wrote those lines for the weight assuming a bunch with a Gaussian distribution in all directions, e.g.
n(x,y,z) = n_b exp[-x^2/(2sigma_x**2)] * exp[-y^2/(2sigma_y**2)] * exp[-z^2/(2sigma_z**2)] (1),
where n_b is the peak charge density of the bunch, and the sigma_* are the rms sizes of the bunch in the respective directions.

For cylindrical symmetry, I assumed sigma_y=sigma_z=sigma_r.
In that case, the total charge of the bunch is n_b*2*np.pi*np.sqrt(2*np.pi)*(sigma_r**2)*sigma_x, that you see in the line setting the weight of each macro-particle, which is equal to the total bunch charge divided by the number of bunch macro-particles.
All must be done coherently with the chosen normalisation.

This namelist uses the same principle, defining an electron bunch with 60 pC:
https://github.com/SmileiPIC/TP-M2-GI/blob/main/InputNamelist.py

Through the postprocessing script here:
https://github.com/SmileiPIC/TP-M2-GI/blob/main/Postprocessing_Scripts/Compute_bunch_parameters.py
you can verify the total charge of the bunch summing the weights stored in the TrackParticles and then converting to pC.

The Rho_electronbunch in the Probes and Fields should show the density in (1).