def electric_succeptibility(l, Ep, epsilonp, rm, n, rho, unfolding, part=0):
    r"""Return the electric succeptibility for a given field.

    INPUT:

    -  ``l`` - The index labeling the probe field.
    -  ``Ep`` - A list of the amplitudes of all pump fields.
    -  ``epsilonp`` - The polarization vector of the probe field.
    -  ``rm`` -     The below-diagonal components of the position operator \
    in the cartesian basis:
    -  ``n`` - The number density of atoms.
    -  ``rho`` - A density matrix in unfolded format, or a list of such \
    density matrices.
    -  ``unfolding`` - A mapping from matrix element indices to unfolded \
    indices.


    >>> import numpy as np
    >>> from sympy import symbols
    >>> from scipy.constants import physical_constants
    >>> from fast import vapour_number_density

    >>> e_num = physical_constants["elementary charge"][0]
    >>> hbar_num = physical_constants["Planck constant over 2 pi"][0]

    >>> Ne = 2
    >>> Nl = 1
    >>> Ep = [-1.0]
    >>> epsilonp = np.array([[0, 0, 1.0]])
    >>> delta = symbols("delta")

    >>> detuning_knob = [delta]
    >>> gamma = np.array([[0.0, -1.0], [1.0, 0.0]])
    >>> omega_level = np.array([0.0, 100.0])
    >>> rm = [np.array([[0.0, 0.0], [1.0, 0.0]])*hbar_num/e_num
    ...       for p in range(3)]
    >>> xi = np.array([[[0, 1], [1, 0]]])
    >>> theta = phase_transformation(Ne, Nl, rm, xi)
    >>> sweep_steady_state = fast_sweep_steady_state(Ep, epsilonp, gamma,
    ...                                              omega_level, rm, xi,
    ...                                              theta)
    >>> deltas, rho = sweep_steady_state([[-20, 20, 11]])
    >>> n = vapour_number_density(273.15+20, "Rb")
    >>> unfolding = Unfolding(Ne, True, True, True)
    >>> chire = electric_succeptibility(0, Ep, epsilonp, rm, n,
    ...                                 rho, unfolding)
    >>> print(chire)
    [ 4.4824e-09-1.1206e-10j  5.5971e-09-1.7491e-10j  7.4459e-09-3.1024e-10j
      1.1097e-08-6.9356e-10j  2.1449e-08-2.6811e-09j  0.0000e+00-5.9877e-08j
     -2.1449e-08-2.6811e-09j -1.1097e-08-6.9356e-10j -7.4459e-09-3.1024e-10j
     -5.5971e-09-1.7491e-10j -4.4824e-09-1.1206e-10j]

    """
    epsilonm = epsilonp.conjugate()
    rp = np.array([rm[i].transpose().conjugate() for i in range(3)])
    if part == 1:
        op = cartesian_dot_product(rp, epsilonm[0])
        op += cartesian_dot_product(rm, epsilonp[0])
        op = -e_num*n/epsilon_0_num/np.abs(Ep[0])*op
    elif part == -1:
        op = cartesian_dot_product(rm, epsilonp[0])
        op += - cartesian_dot_product(rp, epsilonm[0])
        op = -1j*e_num*n/epsilon_0_num/np.abs(Ep[0])*op
    elif part == 0:
        chire = electric_succeptibility(l, Ep, epsilonp, rm,
                                        n, rho, unfolding, +1)
        chiim = electric_succeptibility(l, Ep, epsilonp, rm,
                                        n, rho, unfolding, -1)
        return chire + 1j*chiim

    return np.real(observable(op, rho, unfolding))