docs: def lindblad_terms(gamma, rho, Ne, verbose=1):
      u"""Return the Lindblad terms for decays gamma in matrix form.
  
      >>> from sympy import pprint
      >>> aux = define_frequencies(4, explicitly_antisymmetric=True)
      >>> omega_level, omega, gamma = aux
      >>> gamma = gamma.subs({gamma[2, 0]:0, gamma[3, 0]:0, gamma[3, 1]:0})
      >>> pprint(gamma, use_unicode=True)
      ⎡ 0   -γ₂₁   0     0  ⎤
      ⎢                     ⎥
      ⎢γ₂₁   0    -γ₃₂   0  ⎥
      ⎢                     ⎥
      ⎢ 0   γ₃₂    0    -γ₄₃⎥
      ⎢                     ⎥
      ⎣ 0    0    γ₄₃    0  ⎦
      >>> rho = define_density_matrix(4)
      >>> pprint(lindblad_terms(gamma, rho, 4), use_unicode=True)
      ⎡                -γ₂₁⋅ρ₁₂             -γ₃₂⋅ρ₁₃             -γ₄₃⋅ρ₁₄      ⎤
      ⎢ γ₂₁⋅ρ₂₂        ─────────            ─────────            ─────────     ⎥
      ⎢                    2                    2                    2         ⎥
      ⎢                                                                        ⎥
      ⎢-γ₂₁⋅ρ₂₁                          γ₂₁⋅ρ₂₃   γ₃₂⋅ρ₂₃    γ₂₁⋅ρ₂₄   γ₄₃⋅ρ₂₄⎥
      ⎢─────────  -γ₂₁⋅ρ₂₂ + γ₃₂⋅ρ₃₃   - ─────── - ───────  - ─────── - ───────⎥
      ⎢    2                                2         2          2         2   ⎥
      ⎢                                                                        ⎥
      ⎢-γ₃₂⋅ρ₃₁     γ₂₁⋅ρ₃₂   γ₃₂⋅ρ₃₂                         γ₃₂⋅ρ₃₄   γ₄₃⋅ρ₃₄⎥
      ⎢─────────  - ─────── - ───────  -γ₃₂⋅ρ₃₃ + γ₄₃⋅ρ₄₄   - ─────── - ───────⎥
      ⎢    2           2         2                               2         2   ⎥
      ⎢                                                                        ⎥
      ⎢-γ₄₃⋅ρ₄₁     γ₂₁⋅ρ₄₂   γ₄₃⋅ρ₄₂    γ₃₂⋅ρ₄₃   γ₄₃⋅ρ₄₃                     ⎥
      ⎢─────────  - ─────── - ───────  - ─────── - ───────       -γ₄₃⋅ρ₄₄      ⎥
      ⎣    2           2         2          2         2                        ⎦
  
      Notice that there are more terms than simply adding a decay
      gamma_ij*rho_ij/2 for each coherence.
  
      """
      # We count the necessary Lindblad operators.
      Nterms = 0
      for i in range(Ne):
          for j in range(i):
              if gamma[i, j] != 0:
                  Nterms += 1
  
      L = zeros(Ne)
      counter = 0
      t0 = time()
      for i in range(Ne):
          for j in range(i):
              if gamma[i, j] != 0:
                  counter += 1
                  sig = ket(j+1, Ne)*bra(i+1, Ne)
                  L += gamma[i, j]*lindblad_operator(sig, rho)
                  tn = time()
                  if tn-t0 > 1:
                      aux = "Calculated up to i={}, j={}, or {}/{} = {:2.2f} %."
                      if verbose > 0:
                          print(aux.format(i, j, counter, Nterms,
                                           float(counter+1)/Nterms*100))
                      t0 = tn
      return L