Yes, you can do this step by step. It might help to rewrite this as operators acting on $L(\mathbb{C}^n\otimes \mathbb{C}^n)$. Then the tensor product decomposition that you want is $U\otimes U^\ast\otimes U\otimes U^\ast$. Also note that you can restrict to $SU(n)$.

Now, you can use Littlewood-Richardson rule to find out the decomposition. The fundamental representation has the highest weight (with Dynkin labels) $(1,0^{n-2})$, also denoted by the Young diagram with one box. Its conjugate has weight $(0^{n-2},1)$ (this is the diagram with one column containing $N-1$ boxes). First taking the tensor product of these two, the LR rules give the decomposition as the highest weights $0^{n-1}$ (i.e., a Young diagram with a single column of $N$ boxes) and $(1,0^{n-3},1)$. The former corresponds to the trivial representation (multiples of identity) and the latter corresponds to the adjoint representation (traceless operators) as OP points out. Now, for the tensor product in the question, we want to take the tensor product of $(0^{n-1}) \oplus (1,0^{n-3},1)$ with itself. We get one copy of the trivial, two copies of the adjoint and a tensor product of the adjoint with itself (which we have to decompose).

For the tensor product of the adjoint with itself, one can use the LR rules again. For an arbitrary $N$, the highest weights are

$(0^{n-1})$ trivial representation.

$(0,1,0^{n-3}),(0,0,1,0^{n-4}),\dots,(0^{n-2},1)$

$(3,0^{n-2})$

$(2,1,0^{n-3}),(2,0,1,0^{n-4}),\dots (2,0^{n-3},1)$

$(0,2,0^{n-3})$

$(0,1,1,0^{n-4}),(0,1,0,1,0^{n-5}),\dots ,(0,1,0^{n-4},1)$

$(1,1,0^{n-3}),(1,0,1,0^{n-4}),\dots,(1,0^{n-3},1)$. Two copies of each of these.

$(1,0^{n-2})$

$(2,0^{n-2})$

The general problem (of which your two examples are special cases) is that of decomposing mixed tensor products of the form $U^{\otimes m}\otimes U^{\ast\otimes m}$. This is done in Benkart et al..

**Update:**

Suppose we want to know how these irreducible representations can be realized as subspaces of the original space. This can be done for some of the representations easily. For the others one has to find the Clebsch-Gordan coefficients.

First, consider the simpler problem $U\otimes U^\ast$. We saw that this decomposes into the trivial and the adjoint. In the tensor product of an irreducible representation and its conjugate, the trivial representation is the subspace spanned by the vector $\sum_i |i\rangle |i\rangle$, where $|i\rangle$ is an orthonormal basis for $\mathbb{C}^n$. Under the isomorphism $U\otimes U^\ast\longrightarrow U (\cdot ) U^\dagger$, this corresponds to the identity operator. The adjoint representation, then, consists of everything orthogonal to this i.e., vectors of the form $\sum_{i,j} a_{i,j} |i,j\rangle$ such that $\sum_i a_{i,i} =0$. These, of course, correspond to traceless matrices.

For the second problem, fix an orthonormal basis $|\tilde i\rangle$ for the space of traceless matrices. One possible basis is

$|i,j\rangle$ for $i\neq j$ ($d^2-d$ of these) and

$|i,i\rangle - |i+1,i+1\rangle$ ($d-1$ of these normalized appropriately).

We obtained one copy of the adjoint by tensoring the adjoint and the trivial. This means that subspace is spanned by $\sum_{i,j} a_{i,j}|i,j\rangle \sum_k |k,k\rangle$ with $\sum_i a_{i,i}=0$. As a matrix, this is $\sum_{i,j,k}a_{i,j} |i,j\rangle\langle k,k|$. The other copy of the adjoint is seen to be $\sum_{i,jk} |k,k\rangle\langle i,j| a_{i,j}^\ast$, with $\sum_i a_{i,i}=0$. These two copies are not orthonormal, but their symmetric and antisymmetric spaces will be.

Now in the decomposition of the adjoint with itself, again the trivial is

$\sum_i |\tilde i,\tilde i\rangle$.

This corresponds to matrices of the form

$\sum_{i\neq j} |i,j\rangle\langle i,j| + \sum_i (|i,i\rangle - |i+1,i+1\rangle)(\langle i,i| - \langle i+1,i+1 |)$.

For the others, if we determine the Clebsch-Gordan coefficients, then we can find the matrices. In other words, we need to know how to write vectors $|v_\lambda\rangle$ in a particular irreducible space $\lambda$ as $|v_\lambda\rangle = \sum_{\tilde i, \tilde j}c_{\tilde i,\tilde j}^\lambda |\tilde i,\tilde j\rangle$. The coefficients $c_{\tilde i,\tilde j}^\lambda$ can be determined systematically for the unitary group. See for instance this paper.