beam expansion over a certain density

[Aquios7]

Hello,
I'm having trouble getting a beam density over 1e18 without the entire beam exploding out in every direction.
This is the beam that I have:

electron beam

electron_beam = picmi.UniformDistribution(density = 1e15,
lower_bound=[-0.1, -0.1, -1],
upper_bound=[+0.1, 0.1), -0.3],
directed_velocity=[0., 0.,9.79e6],
fill_in=True)

proton beam

proton_beam = picmi.UniformDistribution(density = 1e15,
lower_bound=[-0.1, -0.1, -1],
upper_bound=[+0.1, 0.1), -0.3],
directed_velocity=[0., 0., 9.79e6],
fill_in=True)

It's a neutral beam, and should be able to move in the Z direction without expanding rapidly after a few hundred steps.
I hope I've put enough here to figure out a solution. If not I can post more or all of the code for replication

Thank You

EDIT:
It is an EM simulation, and I have tried with and without the initialize_self_field that I believe refers to the initial Poisson solve.

Here are the species used for the electron and proton beams:

electron species setup

electrons = picmi.Species(particle_type='electron',
name='electrons',
initial_distribution=electron_beam,
initialize_self_field=True,
method='Boris')

proton species setup

protons = picmi.Species(particle_type='proton',
name='protons',
initial_distribution=proton_beam,
initialize_self_field=True,
method='Boris')

The simulation I am trying to run is a neutral beam moving into a transverse magnetic field. I wanted to run a case of medium and high beta, which would need to be at ~0.2 and ~2 thermal beta, which requires high density.

Thank you again

Edit 2:
I've added iamges from VisIt on the low, medium, and high density examples, which show the timestep they expand, up to 20k steps.

[ax3l]

Hi @Aquios7, can you post more details on your input set?

Is this an EM or ES simulation?
Do you do an initial Poisson solve to initialize the fields (if EM), etc?

[Aquios7]

Thank you for the reply.
I'll be updating my question as well as answering here. It is an EM simulation, and I have tried with and without the initialize_self_field that I believe refers to the initial Poisson solve.