docs: def define_density_matrix(Ne, explicitly_hermitian=False, normalized=False,
                            variables=None):
      r"""Return a symbolic density matrix.
  
      The arguments are
  
      Ne (integer):
          The number of atomic states.
      explicitly_hermitian (boolean):
          Whether to make $\rho_{ij}=\bar{\rho}_{ij}$ for $i<j$
      normalized (boolean):
          Whether to make $\rho_{11}=1-\sum_{i>1} \rho_{ii}$
  
      A very simple example:
      >>> define_density_matrix(2)
      Matrix([
      [rho11, rho12],
      [rho21, rho22]])
  
      The density matrix can be made explicitly hermitian
      >>> define_density_matrix(2, explicitly_hermitian=True)
      Matrix([
      [rho11, conjugate(rho21)],
      [rho21,            rho22]])
  
      or normalized
      >>> define_density_matrix(2, normalized=True)
      Matrix([
      [-rho22 + 1, rho12],
      [     rho21, rho22]])
  
      or it can be made an explicit function of given variables
      >>> from sympy import symbols
      >>> t, z = symbols("t, z", positive=True)
      >>> define_density_matrix(2, variables=[t, z])
      Matrix([
      [rho11(t, z), rho12(t, z)],
      [rho21(t, z), rho22(t, z)]])
  
      """
      if Ne > 9:
          comma = ","
          name = r"\rho"
          open_brace = "_{"
          close_brace = "}"
      else:
          comma = ""
          name = "rho"
          open_brace = ""
          close_brace = ""
  
      rho = []
      for i in range(Ne):
          row_rho = []
          for j in range(Ne):
              if i == j:
                  row_rho += [define_symbol(name, open_brace, comma, i, j,
                                            close_brace, variables,
                                            positive=True)]
              elif i > j:
                  row_rho += [define_symbol(name, open_brace, comma, i, j,
                                            close_brace, variables)]
              else:
                  if explicitly_hermitian:
                      row_rho += [conjugate(define_symbol(name, open_brace,
                                                          comma, j, i,
                                                          close_brace,
                                                          variables))]
                  else:
                      row_rho += [define_symbol(name, open_brace, comma, i, j,
                                                close_brace, variables)]
          rho += [row_rho]
  
      if normalized:
          rho11 = 1-sum([rho[i][i] for i in range(1, Ne)])
          rho[0][0] = rho11
  
      rho = Matrix(rho)
      return rho