How to use the qutip.tensor.tensor function in qutip

"""
    Returns the Bell state:

        |B00> = 1 / sqrt(2)*[|0>|0>+|1>|1>]
        |B01> = 1 / sqrt(2)*[|0>|0>-|1>|1>]
        |B10> = 1 / sqrt(2)*[|0>|1>+|1>|0>]
        |B11> = 1 / sqrt(2)*[|0>|1>-|1>|0>]

    Returns
    -------
    Bell_state : qobj
        Bell state

    """
    if state == '00':
        Bell_state = tensor(
            basis(2), basis(2))+tensor(basis(2, 1), basis(2, 1))
    elif state == '01':
        Bell_state = tensor(
            basis(2), basis(2))-tensor(basis(2, 1), basis(2, 1))
    elif state == '10':
        Bell_state = tensor(
            basis(2), basis(2, 1))+tensor(basis(2, 1), basis(2))
    elif state == '11':
        Bell_state = tensor(
            basis(2), basis(2, 1))-tensor(basis(2, 1), basis(2))

    return Bell_state.unit()

def ghz_state(N=3):
    """
    Returns the N-qubit GHZ-state.

    Parameters
    ----------
    N : int (default=3)
        Number of qubits in state

    Returns
    -------
    G : qobj
        N-qubit GHZ-state

    """
    state = (tensor([basis(2) for k in range(N)]) +
             tensor([basis(2, 1) for k in range(N)]))
    return state/np.sqrt(2)

def _subsystem_apply_reference(state, channel, mask):
    if isket(state):
        state = ket2dm(state)

    if isoper(channel):
        full_oper = tensor([channel if mask[j]
                            else qeye(state.dims[0][j])
                            for j in range(len(state.dims[0]))])
        return full_oper * state * full_oper.dag()
    else:
        # Go to Choi, then Kraus
        # chan_mat = array(channel.data.todense())
        choi_matrix = super_to_choi(channel)
        vals, vecs = eig(choi_matrix.full())
        vecs = list(map(array, zip(*vecs)))
        kraus_list = [sqrt(vals[j]) * vec2mat(vecs[j])
                      for j in range(len(vals))]
        # Kraus operators to be padded with identities:
        k_qubit_kraus_list = product(kraus_list, repeat=sum(mask))
        rho_out = Qobj(inpt=zeros(state.shape), dims=state.dims)
        for operator_iter in k_qubit_kraus_list:
            operator_iter = iter(operator_iter)

if len(tlist) == 2:
                    return output.states[-1]
                else:
                    return output.states

            elif unitary_mode =='batch':
                u = np.zeros(len(tlist), dtype=object)
                _rows = np.array([(N+1)*m for m in range(N)])
                _cols = np.zeros_like(_rows)
                _data = np.ones_like(_rows, dtype=complex)
                psi0 = Qobj(sp.coo_matrix((_data, (_rows, _cols))).tocsr())
                if td_type[1] > 0 or td_type[2] > 0:
                    H2 = []
                    for k in range(len(H)):
                        if isinstance(H[k], list):
                            H2.append([tensor(qeye(N), H[k][0]), H[k][1]])
                        else:
                            H2.append(tensor(qeye(N), H[k]))
                else:
                    H2 = tensor(qeye(N), H)
                options.normalize_output = False
                output = sesolve(H2, psi0, tlist, [],
                                 args=args, options=options,
                                 _safe_mode=False)
                for k, t in enumerate(tlist):
                    u[k] = sp_reshape(output.states[k].data, (N, N))
                    unit_row_norm(u[k].data, u[k].indptr, u[k].shape[0])
                    u[k] = u[k].T.tocsr()

            else:
                raise Exception('Invalid unitary mode.')

def _opto_liouvillian(N):
    from qutip.tensor import tensor
    from qutip.operators import destroy, qeye
    from qutip.superoperator import liouvillian
    Nc = 5                      # Number of cavity states
    Nm = N                     # Number of mech states
    kappa = 0.3                 # Cavity damping rate
    E = 0.1                     # Driving Amplitude
    g0 = 2.4*kappa              # Coupling strength
    Qm = 1e4                    # Mech quality factor
    gamma = 1/Qm                # Mech damping rate
    n_th = 1                    # Mech bath temperature
    delta = -0.43               # Detuning
    a = tensor(destroy(Nc), qeye(Nm))
    b = tensor(qeye(Nc), destroy(Nm))
    num_b = b.dag()*b
    num_a = a.dag()*a
    H = -delta*(num_a)+num_b+g0*(b.dag()+b)*num_a+E*(a.dag()+a)
    cc = np.sqrt(kappa)*a
    cm = np.sqrt(gamma*(1.0 + n_th))*b
    cp = np.sqrt(gamma*n_th)*b.dag()
    c_ops = [cc,cm,cp]

    return liouvillian(H, c_ops)

sz = sigmaz()

    sx_list = []
    sy_list = []
    sz_list = []

    for n in range(N):
        op_list = []
        for m in range(N):
            op_list.append(si)

        op_list[n] = sx
        sx_list.append(tensor(op_list))

        op_list[n] = sy
        sy_list.append(tensor(op_list))

        op_list[n] = sz
        sz_list.append(tensor(op_list))

    # construct the hamiltonian
    H = 0

    # energy splitting terms
    for n in range(N):
        H += - 0.5 * h[n] * sz_list[n]

    # interaction terms
    for n in range(N-1):
        H += - 0.5 * Jx[n] * sx_list[n] * sx_list[n+1]
        H += - 0.5 * Jy[n] * sy_list[n] * sy_list[n+1]
        H += - 0.5 * Jz[n] * sz_list[n] * sz_list[n+1]

elif unitary_mode =='batch':
                u = np.zeros(len(tlist), dtype=object)
                _rows = np.array([(N+1)*m for m in range(N)])
                _cols = np.zeros_like(_rows)
                _data = np.ones_like(_rows, dtype=complex)
                psi0 = Qobj(sp.coo_matrix((_data, (_rows, _cols))).tocsr())
                if td_type[1] > 0 or td_type[2] > 0:
                    H2 = []
                    for k in range(len(H)):
                        if isinstance(H[k], list):
                            H2.append([tensor(qeye(N), H[k][0]), H[k][1]])
                        else:
                            H2.append(tensor(qeye(N), H[k]))
                else:
                    H2 = tensor(qeye(N), H)
                options.normalize_output = False
                output = sesolve(H2, psi0, tlist, [],
                                 args=args, options=options,
                                 _safe_mode=False)
                for k, t in enumerate(tlist):
                    u[k] = sp_reshape(output.states[k].data, (N, N))
                    unit_row_norm(u[k].data, u[k].indptr, u[k].shape[0])
                    u[k] = u[k].T.tocsr()

            else:
                raise Exception('Invalid unitary mode.')


    elif len(c_op_list) == 0 and H0.issuper:
        # calculate the propagator for the vector representation of the
        # density matrix (a superoperator propagator)

The index of the target qubit.

    Returns
    -------
    gate : qobj
        Quantum object representation of N-qubit gate.

    """

    if N < 1:
        raise ValueError("integer N must be larger or equal to 1")

    if target >= N:
        raise ValueError("target must be integer < integer N")

    return tensor([identity(2)] * (target) + [U] +
                  [identity(2)] * (N - target - 1))

# Generate the correct order for qubits permutation,
    # eg. if N = 5, targets = [3,0], the order is [1,2,3,0,4].
    # If the operator is cnot,
    # this order means that the 3rd qubit controls the 0th qubit.
    new_order = [0] * N
    for i, t in enumerate(targets):
        new_order[t] = i
    # allocate the rest qutbits (not targets) to the empty
    # position in new_order
    rest_pos = [q for q in list(range(N)) if q not in targets]
    rest_qubits = list(range(len(targets), N))
    for i, ind in enumerate(rest_pos):
        new_order[ind] = rest_qubits[i]
    id_list = [identity(dims[i]) for i in rest_pos]
    return tensor([oper] + id_list).permute(new_order)