Increasing the density of a plasma beam without its rapid expansion

[Aquios7]

I have a simulation of a plasma beam going across a 3D grid that works at lower beam densities, but when I increase the flux to 8.75e23 (1e17 beam density in my case) the beam expands rapidly beyond some point. Is there a parameter(s) that needs to be adjusted that could be causing the beam to spread radially?

I've added an image of the beam that I'm getting when I try to go to higher density.

[ax3l]

Thanks for the question. I cannot really follow, is this maybe just physical expansion from space charge going up in your higher density beam? Can you share more details about what you expect and set?

[Aquios7]

Thanks for the reply.

For the plasma beam, I assume it is charge and current neutralized. With ions and electrons matching charge and moving at the same speed, I expect to have a beam that stays together as I increase the density over the 2 meter distance.

With the kinetic velocity for both particles at 8.75e6 m/s (400keV ions) and the thermal velocity at 1.38e5 (100eV ions) I wouldn't expect this much spread in the beam at 1 meter distance, moreover the spread gets worse the higher I go in beam density.

Here is how I set up the ion and electron beams:

# electron beam
electron_beam = picmi.UniformFluxDistribution(lower_bound=[xmin, -bwidth/2, -bwidth/2],
                                              upper_bound=[xmin, bwidth/2, bwidth/2],
                                              directed_velocity=[V1, V2, 0.],
                                              rms_velocity=[rms_vel_p]*3,
                                              flux=flux,
                                              flux_normal_axis="x",
                                              surface_flux_position=0.,
                                              flux_direction=1,
                                              flux_tmax=1e-7
                                              )

# proton beam
proton_beam = picmi.UniformFluxDistribution(lower_bound=[xmin, -bwidth/2, -bwidth/2],
                                            upper_bound=[xmin, bwidth/2, bwidth/2],
                                            directed_velocity=[V1, V2, 0.],
                                            rms_velocity=[rms_vel_p]*3,
                                            flux=flux,
                                            flux_normal_axis="x",
                                            surface_flux_position=0.,
                                            flux_direction=1,
                                            flux_tmax=1e-7
                                            )

and here is how I set up the simulation:

# layout of particles (sets of layouts)
pseudoLayout = picmi.PseudoRandomLayout(grid=grid,
                                        n_macroparticles_per_cell=number_sim_macroparticles)

#################################
####### SIMULATION SETUP ########
#################################

sim = picmi.Simulation(
        solver=solver,
        max_steps=max_steps,
        verbose=verbose,
        warpx_serialize_initial_conditions=False,
        warpx_grid_type='staggered',
        warpx_do_dynamic_scheduling=False,
        warpx_use_filter=use_filter,
        time_step_size=dt,
        particle_shape=particle_shape,
        # warpx_amr_restart=restart
        )

sim.add_applied_field(initial_field)

# adding species to the simulation
sim.add_species(
    electrons,
    initialize_self_field=True,
    layout=pseudoLayout)
sim.add_species(
    protons,
    initialize_self_field=True,
    layout=pseudoLayout)

Is there a density to velocity ratio that I don't know about for keeping a neutral beam together?