TL;DR: A two-qubit unitary operator $U$ is a product operator, i.e. $U=U_1\otimes U_2$ for some single-qubit unitaries $U_1$ and $U_2$ if and only if $U$ has Schmidt rank one
$$
U = U_1\otimes U_2 \iff \mathrm{Sch}(U) = 1.\tag1
$$
This generalizes to $n$-qubit unitaries via straightforward recursive applications of $(1)$ over $n-1$ partitionings of the subsystems.
Operator Schmidt decomposition
The space of linear operators on a vector space $\mathcal{H}$ is a vector space $L(\mathcal{H})$ equipped with an inner product $\langle A,B\rangle_{HS} = \mathrm{tr}(A^\dagger B)$ known as the Hilbert-Schmidt inner product. Therefore, Schmidt decomposition applies to linear operators. Explicitly, any operator $U$ can be written as
$$
U = \sum_{k=1}^r \lambda_k A_k\otimes B_k\tag2
$$
where $A_k$ are orthonormal operators on the first subsystem, $B_k$ are orthonormal operators on the second subsystem and $\lambda_k$ are positive real numbers. The integer $r=:\mathrm{Sch}(U)$ is called the Schmidt rank of $U$.
See e.g. $6.4.2$ in Nielsen's PhD thesis or this paper for more details on operator variant of Schmidt decomposition.
Unitary product is a product of unitaries
Now, suppose that $U=A\otimes B$ for some operators $A$ and $B$. We'll show that $U=U_1\otimes U_2$ for some unitary operators $U_1$ and $U_2$. We have
$$
I=U^\dagger U = A^\dagger A\otimes B^\dagger B,\tag2
$$
so $\mathrm{tr}(A^\dagger A)>0$ and $\mathrm{tr}(B^\dagger B)>0$. Moreover, $A^\dagger A=\alpha I$ for some positive real number $\alpha$ and $B^\dagger B=\beta I$ for some positive real number $\beta$. Consequently, $U=U_1\otimes U_2$ where $U_1:=\frac{A}{\sqrt{\alpha}}$ and $U_2:=\frac{B}{\sqrt{\beta}}$ are unitary.
Generalization to $n$-qubit unitaries
An $n$-qubit unitary operator $U$ is the tensor product of an $(n-1)$-qubit unitary $U_{1\dots(n-1)}$ and a single-qubit unitary $U_n$ if and only if its Schmidt rank with respect to the partitioning $\{1\dots(n-1)\}\cup\{n\}$ is one. Applying this check recursively to $U_{1\dots(n-1)}$, we can determine whether $U$ is the tensor product of $n$ single-qubit unitaries by computing $n-1$ Schmidt decompositions.
cirq.sub_state_vector
for this; it's ultimately the same logic. $\endgroup$