docs: def bcs_parameters(n_site, n_fermi, u, t) :
      """Generate the parameters for the BCS ground state, i.e., the
      superconducting gap and the rotational angles in the Bogoliubov
      transformation.
  
       Args:
          n_site: the number of sites in the Hubbard model
          n_fermi: the number of fermions
          u: the interaction strength
          t: the tunneling strength
  
      Returns:
          float delta, List[float] bog_theta
      """
  
      # The wave numbers satisfy the periodic boundary condition.
      wave_num = np.linspace(0, 1, n_site, endpoint=False)
      # The hopping energy as a function of wave numbers
      hop_erg = -2 * t * np.cos(2 * np.pi * wave_num)
      # Finding the Fermi energy
      fermi_erg = hop_erg[n_fermi // 2]
      # Set the Fermi energy to zero
      hop_erg = hop_erg - fermi_erg
  
      def _bcs_gap(x):
          """Defines the self-consistent equation for the BCS wavefunction.
  
          Args:
              x: the superconducting gap
          """
  
          s = 0.
          for i in range(n_site):
              s += 1. / np.sqrt(hop_erg[i] ** 2 + x ** 2)
          return 1 + s * u / (2 * n_site)
  
      # Superconducting gap
      delta = scipy.optimize.bisect(_bcs_gap, 0.01, 10000. * abs(u))
      # The amplitude of the double excitation state
      bcs_v = np.sqrt(0.5 * (1 - hop_erg / np.sqrt(hop_erg ** 2 + delta ** 2)))
      # The rotational angle in the Bogoliubov transformation.
      bog_theta = np.arcsin(bcs_v)
  
      return delta, bog_theta