Matrix representation of total quantum angular momentum so I'm currently learning quantum mechanics in the context of a physical chemistry course. So it's not very mathematical which leads to a lot of little traps that just confuse me.
Let me first introduce the Bra-Ket Notation I'm going to use. The inner product of two complex valued functions $\varphi_m$ and $\varphi_n$ can be written as
$$\int\varphi_m^*\varphi_nd\tau := \langle \varphi_m | \varphi_n\rangle := \langle m | n \rangle \tag{1}$$
whereas $m$ and $n$ are indices for the eigenfunctions and eigenvalues and $\varphi_m^*$ is the complex conjugate.
For an operator $\hat{A}$ we get:
$$\int\varphi_m^*\hat{A}\varphi_nd\tau := \langle \varphi_m | \hat{A}\varphi_n\rangle := \langle m | \hat{A} | n \rangle \tag{2}$$
So that's everything we ever learned about the Dirac-Notation.
Now if we have a complete orthonormal basis consisting of the eigenfunctions $\{ \varphi_n \}$ of a operator $\hat{A}$, we can write $\hat{A}$ as an matrix:
$$ A_{nm} = \int \varphi^*\hat{A}\varphi_md\tau = \langle n|\hat{A}|m\rangle \tag{3} $$
Now, so far so good. Now consider the following problem:
Let $\hat{J}$ be a general angular momentum. Furthermore we denote the spin "electron spin angular momentum" (I guess you could just think of spin) with $\vec{s}$, the "orbital angular momentum" with $\vec{l}$ and the total angular momentum with $\vec{j}$.
We want to calculate the matrices for $\hat{l_x},\hat{l_y},\hat{l_z},\hat{l}^2,\hat{s_x},\hat{s_y},\hat{s_z},\hat{s}^2$.
Since the operators $\hat{l}^2$ and $\hat{l_z}$ commute, we know they share the same basis. (Similar for $\hat{s}^2$ and $\hat{s_z}$). So we can use $|l, m_l\rangle$ respectively $s|m_s\rangle$ as a basis.
We further know that for a general angular Momentum $\hat{J}$ (i.e. especially for $\hat{s_z}, \hat{l_z}$) we have:
$\hat{J}_\pm = \hat{J}_x \pm i\hat{J_y} \tag{3}$
$\hat{J_\pm}|J,M\rangle = \hbar\sqrt{J(J+q) - M(M\pm1)}|J,M\pm 1\rangle \tag{4}$
$\hat{J_z}|J,M\rangle = \hat M|J,M\rangle \tag{5}$
Now, it'd be easy to calculate the matrices for $\hat{l_x},\hat{l_y},\hat{l_z},\hat{l}^2,\hat{s_x},\hat{s_y},\hat{s_z},\hat{s}^2$.
For the spin, we would get a 2 dimensional vector space $V_s$ and for the angular momentum we would get a 3 dimensional vector space $V_l$. We could now ask about the total angular momentum. What space does that one live in? The answer would be: $ V_j = V_s \otimes V_l$ (I hope the notation is correct). It's basis would be
$|l, m_l\rangle s, m_s\rangle \tag{6}$
So for me it's clear, that if we look at the total angular momentum, we basically look at the "union" of both other vector spaces. Now comes the actual question: What they did was the following. First they gave us the following basis:

As you can see, they used
$|l,m_l,s,m_s\rangle = |l, m_l\rangle s, m_s\rangle \tag{6}$
Which already confuses me. Sure it makes sense that a eigenfunction of the total angular momentum basis on the numbers $l, m_l, s, m_s$ but what does the notation $|l,m_l,s,m_s\rangle$ actually notate?
For me $|l, m_l\rangle s, m_s\rangle$ is a "multiplication" of two functions. Which would be an inner product (since we are in a Hilber space), but such an inner product would be donated e.g. as $\langle l, m_l | s, m_s\rangle$ but that obviously is in general not the same since we can't just "flip" a Bra to a Ket like that.
Further more, if we look at the defintion of the Bra-Ket Notation, we suddendly have 4 indices. (Sure, two are fixed and we could shorten it to $|m_l, m_s\rangle$, but there is no mention of such a thing being done.) How should I now apply the Bra-Ket Notation?
Anyway, to get the matrix of my $\hat{l_z}$ operator in the basis $|l, m_l\rangle s, m_s\rangle$ I'd do the following:
We know: $\hat{\ell}_z| l,m_l\rangle = c|l, m_l+1\rangle$ (with $c$ according to (5))
We multiply it with $|s, m_s\rangle$ and get
$\hat{\ell}_z| l,m_l\rangle|s, m_s\rangle = c|l, m_l+1\rangle|s, m_s\rangle$
Then
$\langle s,m_s|\langle l,m_l|\hat{\ell}_z| l',m_l'\rangle|s', m_s'\rangle$
$ = c\langle s,m_s|\langle l,m_l|l', m_l'+1\rangle|s', m_s'\rangle$
$= c\langle s,m_s|\delta_{l,l'}\delta_{m_l,(m_l'+1)}|s', m_s'\rangle$
$= c\delta_{l,l'}\delta_{m_l,(m_l'+1)}\langle s,m_s|s', m_s'\rangle$
$= c\delta_{l,l'}\delta_{m_l,(m_l'+1)}\delta_{s,s'}\delta_{m_s,m_s'}$
$= c\delta_{m_l,(m_l'+1)}\delta_{m_s,m_s'}$
So that works fine. But I don't see how I could do the same thing for the spin operator. I don't know how to handle the Dirac-Notation correctly here.
I hope you see my problem. It was a bit of a big text, so let's summarize my questions:
Q1: What does $|l,m_l,s,m_s\rangle$ denote exactly?
Q2: What does $|l, m_l\rangle s, m_s\rangle$ exactly?
Q3: To get the matrix representation of an operator we'd use (2). I'm confused on how to get to such an expression if we have more than two quantum numbers. How would I get the expression for the matrix?
 A: Like you correctly understood, the total Hilbert space we’re working with is a direct product of the form $V_j=V_l\otimes V_s$. What this means is that to specify one state in the j-space, you need to specify two states, one in each of the constituent spaces.
The states are combined similarly as the spaces by an outer (direct) product. Not inner product. So the state in true sense is actually
$$|l,m_l,s,m_s\rangle = |l\rangle\otimes |m_l\rangle\otimes |s\rangle\otimes |m_s\rangle$$
It’s just like in hydrogen atom where to specify a particular state you need to specify all three things $n,l,m$ And the operators act independently on the relevant states. For example the operator $L_z$ acting on our j-space is actually
$$L_z|l,m_l,s,m_s\rangle = \mathbb{1}\otimes L_z\otimes\mathbb{1}\otimes\mathbb{1}|l,m_l,s,m_s\rangle\\
=\mathbb{1}|l\rangle\otimes L_z|m_l\rangle\otimes \mathbb{1} |s\rangle\otimes \mathbb{1} |m_s\rangle\\
=|l\rangle\otimes m_l|m_l\rangle\otimes |s\rangle\otimes |m_s\rangle\\
\Rightarrow L_z|l,m_l,s,m_s\rangle =m_l|l,m_l,s,m_s\rangle $$
