docs: 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))