diff --git a/doc/source/run/parameters.rst b/doc/source/run/parameters.rst
index 48e74492ce..6b60af4911 100644
--- a/doc/source/run/parameters.rst
+++ b/doc/source/run/parameters.rst
@@ -92,6 +92,10 @@ Plasma parameters
 * ``plasma.radius`` (`float`) optional (default `infinity`)
     Radius of the plasma. Set a value to run simulations in a plasma column.
 
+* ``plasma.channel_radius`` (`float`) optional (default `0.`)
+    Channel radius of a parabolic plasma profile. The plasma density is set to
+    :math:`\mathrm{plasma.density} * (1 + r^2/\mathrm{plasma.channel\_radius}^2)`.
+
 * ``plasma.max_qsa_weighting_factor`` (`float`) optional (default `35.`)
     The maximum allowed weighting factor :math:`\gamma /(\psi+1)` before particles are considered
     as violating the quasi-static approximation and are removed from the simulation.
diff --git a/src/particles/PlasmaParticleContainer.H b/src/particles/PlasmaParticleContainer.H
index 93c87937ac..a2bb601f82 100644
--- a/src/particles/PlasmaParticleContainer.H
+++ b/src/particles/PlasmaParticleContainer.H
@@ -58,6 +58,8 @@ public:
      *  the quasi-static approximation and is removed */
     amrex::Real m_max_qsa_weighting_factor {35.};
     amrex::Real m_radius {std::numeric_limits<amrex::Real>::infinity()}; /**< radius of the plasma */
+    /** defines the channel radius of a parabolic plasma profile */
+    amrex::Real m_channel_radius {std::numeric_limits<amrex::Real>::infinity()};
     amrex::IntVect m_ppc {0,0,1}; /**< Number of particles per cell in each direction */
     amrex::RealVect m_u_mean {0,0,0}; /**< Avg momentum in each direction normalized by m*c */
     amrex::RealVect m_u_std {0,0,0}; /**< Thermal momentum in each direction normalized by m*c */
diff --git a/src/particles/PlasmaParticleContainer.cpp b/src/particles/PlasmaParticleContainer.cpp
index b88efa612d..959314ad58 100644
--- a/src/particles/PlasmaParticleContainer.cpp
+++ b/src/particles/PlasmaParticleContainer.cpp
@@ -8,6 +8,9 @@ PlasmaParticleContainer::PlasmaParticleContainer (amrex::AmrCore* amr_core)
     amrex::ParmParse pp("plasma");
     pp.query("density", m_density);
     pp.query("radius", m_radius);
+    pp.query("channel_radius", m_channel_radius);
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_channel_radius != 0,
+                                     "The plasma channel radius must not be 0");
     pp.query("max_qsa_weighting_factor", m_max_qsa_weighting_factor);
     amrex::Vector<amrex::Real> tmp_vector;
     if (pp.queryarr("ppc", tmp_vector)){
@@ -34,7 +37,7 @@ PlasmaParticleContainer::InitData ()
     reserveData();
     resizeData();
 
-    InitParticles(m_ppc,m_u_std, m_u_mean, m_density, m_radius);
+    InitParticles(m_ppc, m_u_std, m_u_mean, m_density, m_radius);
 
     m_num_exchange = TotalNumberOfParticles();
 }
diff --git a/src/particles/PlasmaParticleContainerInit.cpp b/src/particles/PlasmaParticleContainerInit.cpp
index fc725c3c79..8423c3f7f4 100644
--- a/src/particles/PlasmaParticleContainerInit.cpp
+++ b/src/particles/PlasmaParticleContainerInit.cpp
@@ -90,6 +90,8 @@ InitParticles (const amrex::IntVect& a_num_particles_per_cell,
 
         PhysConst phys_const = get_phys_const();
 
+        const amrex::Real channel_radius = m_channel_radius;
+
         amrex::ParallelFor(tile_box,
         [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
         {
@@ -120,9 +122,10 @@ InitParticles (const amrex::IntVect& a_num_particles_per_cell,
                     y >= a_bounds.hi(1) || y < a_bounds.lo(1) ||
                     x*x + y*y > a_radius*a_radius ) continue;
 
+                const amrex::Real rp = std::sqrt(x*x + y*y);
+
                 amrex::Real u[3] = {0.,0.,0.};
-                ParticleUtil::get_gaussian_random_momentum(u, a_u_mean,
-                                                           a_u_std);
+                ParticleUtil::get_gaussian_random_momentum(u, a_u_mean, a_u_std);
 
                 ParticleType& p = pstruct[pidx];
                 p.id()   = pid + pidx;
@@ -131,8 +134,10 @@ InitParticles (const amrex::IntVect& a_num_particles_per_cell,
                 p.pos(1) = y;
                 p.pos(2) = z;
 
-                arrdata[PlasmaIdx::w        ][pidx] = a_density * scale_fac;
-                arrdata[PlasmaIdx::w0       ][pidx] = a_density * scale_fac;;
+                arrdata[PlasmaIdx::w        ][pidx] =
+                        a_density*(1. + rp*rp/(channel_radius*channel_radius)) * scale_fac;
+                arrdata[PlasmaIdx::w0       ][pidx] =
+                        a_density*(1. + rp*rp/(channel_radius*channel_radius)) * scale_fac;
                 arrdata[PlasmaIdx::ux       ][pidx] = u[0] * phys_const.c;
                 arrdata[PlasmaIdx::uy       ][pidx] = u[1] * phys_const.c;
                 arrdata[PlasmaIdx::psi      ][pidx] = 0.;