def nD0_active(N, vzero, dr, L):
    """
    Returns product of particle density n = N/L^2 and active diffusion constant
    D_0 = vzero^2/2*dr.

    Parameters
    ----------
    N : int or float
        Number of particles.
    vzero : float
        Self-propelling velocity.
    dr : float
        Rotation diffusion constant.
    L : float
        Characteristic system length.

    Returns
    -------
    product : float
        n D_0
    """

    return (N*(vzero**2))/(2*dr*(L**2))