def find_neighbor_throats(self, pores, mode='union', flatten=True):
        r"""
        Returns a list of throats neighboring the given pore(s)

        Parameters
        ----------
        pores : array_like
            Indices of pores whose neighbors are sought

        flatten : boolean, optional
            If ``True`` (default) a 1D array of unique throat indices is
            returned. If ``False`` the returned array contains arrays of
            neighboring throat indices for each input pore, in the order
            they were sent.

        mode : string
            Specifies logic to filter the resulting list.  Options are:

            **'or'** : (default) All neighbors of the input pores.  This is
            also known as the 'union' in set theory or 'any' in boolean logic.
            Both keywords are accepted and treated as 'or'.

            **'xor'** : Only neighbors of one and only one input pore.  This
            is useful for finding the thraots that are not shared by any of the
            input pores.

            **'xnor'** : Neighbors that are shared by two or more input pores.
            This is equivalent to finding all neighbors with 'or', minus those
            found with 'xor', and is useful for finding neighbors that the
            inputs have in common.

            **'and'** : Only neighbors shared by all input pores.  This is also
            known as 'intersection' in set theory and (somtimes) as 'all' in
            boolean logic.  Both keywords are accepted and treated as 'and'.

        Returns
        -------
        If ``flatten`` is ``True``, returns a 1D array of throat indices
        filtered according to the specified mode.  If ``flatten`` is ``False``,
        returns a list of lists, where each list contains the neighbors of the
        corresponding input pores.

        Notes
        -----
        The ``logic`` options are applied to neighboring bonds only, thus it
        is not possible to obtain bonds that are part of the global set but
        not neighbors. This is because (a) the list of global bonds might be
        very large, and (b) it is not possible to return a list of neighbors
        for each input site if global sites are considered.

        Examples
        --------
        >>> import openpnm as op
        >>> pn = op.network.Cubic(shape=[5, 5, 5])
        >>> Ts = pn.find_neighbor_throats(pores=[0, 1])
        >>> print(Ts)
        [  0   1 100 101 200 201]
        >>> Ts = pn.find_neighbor_throats(pores=[0, 1], flatten=False)
        >>> print(Ts)
        [array([  0, 100, 200]), array([  0,   1, 101, 201])]

        """
        pores = self._parse_indices(pores)
        if sp.size(pores) == 0:
            return sp.array([], ndmin=1, dtype=int)
        if 'lil' not in self._im.keys():
            self.get_incidence_matrix(fmt='lil')
        neighbors = topotools.find_neighbor_bonds(sites=pores, logic=mode,
                                                  im=self._im['lil'],
                                                  flatten=flatten)
        return neighbors