Draft: Hamilton-Jacobi Equation in Quantum Mechanics

Hamilton-Jacobi form

In the non-relativistic quantum mechanics (NRQM), Schoedinger Equation dictates the dynamical evolution of the system, $i \hbar \partial_t \psi = -\frac{\hbar^2}{2m}\nabla^2\psi + V\psi.$ Schroedinger equation is a linear second-order partial differential equation (PDE). Differential operator $\hat{\mathbf{p}}=\frac{\hbar}{i}\nabla$ is called momentum, and $\hat H = \frac{\hat{\mathbf{p}}^2}{2m} + V$ Hamiltonian. Let's do a change of variable $\psi = \exp \frac{i}{\hbar}S$. We obtain a non-linear second-order PDE, $$\partial_t S + \frac{1}{2m}(\nabla S)^2 + V =  \frac{i\hbar}{2m}\nabla^2_qS$$ Let's call the complex function $\mathbf{p} = \nabla_q S$ momentum and $H = \frac{\mathbf p^2}{2m} + V$ the Hamiltonian. The above equation can be written as $\partial_t S + H(q, \nabla_q S, t) =  \frac{i\hbar}{2m}\nabla^2 S$. This is the Hamilton-Jacobi form of the wave equation. Function $\mathbf p$ is different from but consistent with the differential operator $\hat{\mathbf{p}}$. It is easy to see, $\hat{\mathbf{p}}\psi = \mathbf{p}\cdot \psi$.

If the Hamiltonian $H$ is time-independent, we separate the time by $S = W - Et$. The resulted equation, $i\hbar \nabla \cdot \mathbf p - \mathbf{p}^2 + 2m(E-V)$, corresponds to the time-independent Schroedinger equation.

Second quantization

The quantum HJ equation differs the classical one by an extra term $\frac{i\hbar}{2m}\nabla \cdot \mathbf{p}$. It has been argued by Roncadelli and Schulman (2007) that this term arises from second-quantization of the classical HJ equation, $$\partial_t \hat S + H(\hat q, \partial_{\hat q} \hat S, t) = 0.$$ Here $\hat S = S(\hat q, \hat Q, t)$ is a quantum mechanical operator that depends on the generalized coordinate operator before and after canonical transformation, $\hat{q}$ and $\hat Q$. In order to get the wave equation, define $S = \langle q |\hat S |Q \rangle$. Obviously, $\langle q |\partial_t \hat S |Q \rangle = \partial_t S$, $\langle q |V(\hat q) |Q \rangle = V(q)$. However, one has to be careful with the quadratic term $\langle q |\hat p^2 |Q \rangle$, as $\hat q$ may not commute with $\hat Q$, which means $\hat p^2$ need to be ordered. In general $\hat p = \sum_i a_i(\hat q, t) b_i(\hat Q,t)$. Therefore,
$$\hat p^2 = \sum_i  \hat p a_i(\hat q, t) b_i(\hat Q,t) = \sum_i a_i(\hat q, t) \hat p b_i(\hat Q,t) + \sum_i [\hat p, a_i(\hat q,t)] b_i (\hat Q,t).$$ The first term is well-ordered. For the second term, recall that $[\hat p, A] = -i\partial_q \hbar A$. So $\hat p^2 = :\hat p^2: - i\hbar \partial_q \hat p$. Sandwiched with $\langle q |$ and $|Q\rangle$, we derive the quantum HJ equation.

Classical limit

The quantum Hamilton-Jacobi equation reduces to the classical Hamilton-Jacobi equation as $\hbar \to 0$. This observation suggests an expansion of the quantum action $S$ around the classical action $S_{cl} \equiv S_0$ with respect to $\hbar$: $$S = S_0 + \frac{\hbar}{i} S_1 + \left(\frac{\hbar}{i}\right)^2 S_2 + \cdots; \\ W = W_0 + \frac{\hbar}{i} W_1 + \left(\frac{\hbar}{i} \right)^2 W_2 + \cdots.$$ For a semi-classical system, $\hbar | S_n | \ll |S_{n-1}|$. Therefore, the quantum theory can be solved from improving the classical action order by order. This method is known as the Eikonal approximation (cf. WKB approximation). The classical action obeys the classical HJ equation, $\partial_t S_{cl} + H(q, \partial_q S_{cl}, t) = 0$. To the first order of $\hbar$, $$\partial_t S_1 + H_1(q, \partial_q S_{cl}, \partial_q S_1, t) = 0,$$ here $H_1 = H(q, p_{cl}+p_1,t) - H(q, p_{cl},t) \simeq \partial_p H \partial_q S_1$.

Fig. The solution surface of $S_1$ is the collection of the integral curve $\dot{\mathbf{q}} = \mathbf{p}_{cl}$, $\dot{S}_1 = -\frac{1}{2}\nabla \cdot \mathbf{p}_{cl}$.

In the time-independent problem, the classical momentum is just $\mathbf p_{cl}^2 = 2m (E - V)$. Up to the first order of $\hbar$, $\mathbf p = \mathbf p_{cl} + \frac{\hbar}{i}\mathbf p_1$. Substitute it to the HJ equation and drop the second order terms, $\mathbf p_{cl} \cdot \nabla S_1 = -\frac{1}{2}\nabla \cdot \mathbf p_{cl}$. $S_1$ is a hyper surface composed of the integral curves $$\dot{\mathbf{q}} = \mathbf{p}_{cl}, \quad \dot{S_1} = -\frac{1}{2}\nabla \cdot \mathbf{p}_{cl}$$  The wavefunction is, $$\psi = \exp\left( \frac{i}{\hbar} \int \mathrm{d} \mathbf{q} \cdot \mathbf{p}_{cl} + S_1 \right)$$  For example, in the 1D case, $S_1 = -\frac{1}{2} \log p_{cl} + c$ and the full solution is, $$\psi = \frac{C_+}{\sqrt{p_{cl}}} \exp\left( \frac{i}{\hbar}\int p_{cl} \mathrm dq \right) + \frac{C_-}{\sqrt{p_{cl}}} \exp\left(- \frac{i}{\hbar}\int p_{cl} \mathrm dq \right).$$
In the classically forbidden region $(E<V)$, $p_{cl} = \pm i \sqrt{2m|E-V|}$ becomes imaginary and the wavefunction developed an exponential-decline factor. This is the phenomenon of quantum tunneling.

Separation of variables

In one-dimension, the time-independent HJ equations, $i\hbar p' = p^2 - 2m(E-V)$, is a Riccati equation. If given a particular solution $p_0$, the general solution is $$p = p_0 + \frac{ \exp\left( \frac{2}{i\hbar} \int^x \mathrm dy p_0(y) \right) }{C + \frac{i}{\hbar} \int^x \mathrm dy \exp\left( \frac{2}{i\hbar} \int^y \mathrm dz p_0(z) \right) }.$$

Some higher-dimensional system can be reduced to one-dimensional by the separation of variables. The notable example is the central potential: $V = V(r)$. We can separate variables by $S = W_r + W_\theta + W_\phi - Et$ and get, $$i\hbar W''_\phi - {W'}_\phi^2 = -m_s^2 \hbar^2; \\ i\hbar \frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left(\sin\theta W'_\theta \right) - \frac{m_s^2\hbar^2}{\sin^2\theta} - {W'}_\theta^2 = -l(l+1) \hbar^2; \\ i\hbar \frac{1}{2m r^2} \frac{\partial}{\partial r}\left( r^2 W'_r \right) - \frac{1}{2m r^2}l(l+1)\hbar^2 = V(r) - E + {W'_r}^2$$ with the boundary condition $W \to W_{cl}$ if $\hbar \to 0, W_\phi(\phi + 2\pi) = W_\phi(\phi)$.

The first equation has the solution $W_\phi(\phi) = \frac{\hbar}{i}\log \cos m_s\phi$ with $m_s = 0, \pm 1, \pm 2, \cdots$.

The second equation is a Riccati equation, that can be rewritten to a linear second order differential equation. The resultant equation is just the Bessel equation. Meanwhile, if we do a change of variable $W'_\theta = p_\theta = i\hbar ( u_\theta + \frac{1}{2}\cot\theta)$, it becomes, $$u'_\theta = u_\theta^2 + \frac{1-4m_s^2}{4\sin^2\theta} + (l+ \frac{1}{2})^2$$  The third equation is also a Riccati equation. With a change of variable $W'_r = p_r = \frac{i\hbar}{2m}( u_r + \frac{1}{r})$, the equation becomes $$u'_r = u_r^2 + \frac{4 m^2}{\hbar^2} (E - V(r)) - \frac{2 m l(l+1)}{r^2}$$

Bohr-Sommerfeld quantization

Recall the angular action $J_q$ is defined as $$J_q = \oint_C \mathrm{d}q \; p_q$$ where $p_q = \frac{\partial}{\partial q}W$. The value of $J_q$ only depends on the analytic structure of the solution $p_q$ of the Hamilton-Jacobi equation. The poles are the "good" singular points that gives finite result for $J_q$. The "bad" singularities include the branch points and the essential singular points are called the critical points. $J_q$ gains contributions from both poles and critical points.

Non-linear first order differential equations $\frac{\mathrm{d}w}{\mathrm{d}z} = F(z, w)$ where $z \in \mathbb{C}$, and $F(z, w)$ is locally analytic can have "internal singularities" or the movable singular points. The location of the singular points depends on integration constant. A leading example is the equation $$\frac{\mathrm{d}w}{\mathrm{d}z} = w^2 \quad \Rightarrow \quad w(z) = \frac{1}{c - z}.$$ Riccati equations do not have movable critical singularities. In fact, in 1884 W. Fuchs showed that Riccati equations are the only class of first order differential equations without movable critical singular points [2].

Riccati equation $w' = w^2 + f(z)$ admits a solution $w(z) = -\frac{v'(z)}{v(z)}$ where $v(z)$ is a solution of the second order linear equation $v'' + f(z) v = 0$. Then, the original solution $w(z) = - \frac{v'(z)}{v(z)} = - \frac{g(z) + (z-z_1) g'(z)}{(z-z_1)g(z)} = -\frac{1}{z-z_1} - \frac{g'(z)}{g(z)} = -\sum_i \frac{1}{z-z_i} - \varphi(z)$ ($\varphi(z)$ only has critical singularities ) can only have first order non-movable singularities. Then the angular action is quantized with the famous Bohr-Sommerfeld condition $$J_q = 2 \pi i n (-1) i \hbar = 2 \pi n \hbar + C, \quad n = 0,1,2,...$$ where $n$ is the number of single poles lying on the real axis of the solution $w(z)$ hence $p_q$. $C$ is constant. With any luck, the solution has no critical singularities and then $C = 0$. A nice feature of this analysis is the manifest of the correspondence principle. In the classic mechanics, $p_q = \sqrt{2m (E-V)}$ that has a branch cut. Whereas in quantum theory in the classical limit $n \to \infty$, the poles on the real axis behaves like the branch cut. Along this line, the ground state is the solute with the minimal pole.


- - -
[1]:http://eqworld.ipmnet.ru/en/solutions/ode/ode0123.pdf
[2]:http://www.encyclopediaofmath.org/index.php/Painlev%C3%A9-type_equations#References
[3]:http://www.risc.jku.at/publications/download/risc_2773/prep2.pdf
[4]:Robert A. Leacock and Michael J. Padgett, Hamilton-Jacobi Theory and the Quantum Action Variable, Phys. Rev. Lett. 50, 3–6 (1983) URL: http://prl.aps.org/abstract/PRL/v50/i1/p3_1
[5]:Marco Roncadelli and L. S. Schulma, Quantum Hamilton-Jacobi Theory, PRL, 99, 170406 (2007), arxiv URL: http://arxiv.org/pdf/0712.0307v1.pdf
[6]: https://michaelberryphysics.files.wordpress.com/2013/07/berry023.pdf

Symmetry, conservation laws and compatible operators

Symmetry

A symmetry of a dynamical system is a transformation that leaves the action $S = \int \mathrm{d}^4 x \mathcal{L}$ invariant. A set of symmetry may form a group, the symmetry group. Most important symmetries in physics are from symmetry groups. Notable examples of the symmetry group include: the spacetime symmetry - Poincaree group and conformal group, the gauge groups, the lattice group and the supersymmetry.

A continuous group (Lie group) can be written in terms of the single parameter exponential map, $g = \exp(-i \theta_r X^r), \; r = 1, 2, ..., n$, where $\theta_r$ are scalars. $X^r$ are called the group generators. The group generators form a Lie algebra. $n$ is called the dimension of the Lie algebra. In quantum mechanics, the group transformation is represented by unitary or anti-unitary operators (Wigner theorem). Therefore, the group generators $X^r$ are Hermitian operators ("real operators"), hence are observables. The Lie algebra, especially the commutation relations between groups generators (commutators) $[X^r, X^s] = iC^{rst} X^t, r,s,t = 1,2,3..., n$ is dictated by the Lie group. If in a group, all the commutators vanish, the group is called an Abelian group, otherwise called Non-Abelian. Most important symmetry groups in physics are Non-Abelian.

It's worth mentioning the Hamiltonian $H$ is a generator of the Poincaree group.

Conservation laws

A quantity $A$ is called conserved if it does not change upon time, $\frac{\mathrm{d}}{\mathrm{d}t} A = 0$. In quantum mechanics, the time evolution is dictated by the Hamiltonian of the system. In Heisenberg picture, $$\frac{\mathrm{d}}{\mathrm{d}t} A = i [ H, A ] + \frac{\partial}{\partial t} A.$$ Similarly, in classical mechanics, $$\frac{\mathrm{d}}{\mathrm{d}t} A = -\{H, A\} + \frac{\partial}{\partial t} A$$ where $\{ A, B\}$ is the Poisson bracket. One may note the similarity between the Poisson bracket and the commutator. Indeed, in the canonical quantization, Poisson bracket is replaced by the commutator $\{ \cdot, \cdot \} \to \frac{1}{i}[ \cdot, \cdot ]$.

Note that a conserved quantity does not necessarily commute with the Hamiltonian, as it may have explicit time dependence init. One has to be especially cautious in a time-independent problem, as we normally only consider the operator at $t = 0$.

Let $U = e^{-i \theta A}$ be a symmetry. Then the principle of time evolution should not change for a transformed state: $U \left.|\psi(t)\right> = V(t) U\left.|\psi(0)\right>$, where $V(t) = \exp(-i H t)$ is the time-evolution operator, $V(t) \left.|\psi(0)\right>$. That implies $U V = V U$. So we conclude $[H, A] = 0$. There are two pitfalls: 1). a symmetry may be represented by anti-unitary operators. In that case, the commutator is non-vanishing. 2). the generator may be explicitly time dependent, especially the symmetry transformation involves time.

Noether's theorem

Symmetry and conservation laws are closely related by Noether's theorem. Noether theorem states, for each differential symmetry of a physical system, there exists a conservation law. More specifically, let $q_i$ be the generalized coordinates (for instance, the quantum fields). Under the infinitesimal symmetry transformation, $q_i \to q_i + \epsilon \Delta q_i$. The Lagrangian must takes the form, $$\frac{\mathrm{d}}{\mathrm{d}t} L = \frac{\mathrm{d}}{\mathrm{d}t} G$$ to keep the action $S = \int \mathrm{d} t L(q, \dot{q}, t)$ invariant. Then, there exists a conserved quantity (often called conserved charge), $$Q = G - \frac{\partial L}{\partial \dot{q}_i} \Delta q_i; \\ \dot{Q} = 0$$ The conserved charges are the generators of the symmetry group.

Boost generator $\vec{K}$ is the conserved charge of boost transformation. But it does not commute with $H$. So clearly, it has explicit time dependence. In fact, in the Poincaree algebra, $$K^i = M^{0i} = x^0 P^i - x^i P^0 = t P^i - x^i H.$$ This is even clearer in non-relativistic quantum mechanics, as the boost at $t = 0$ is simply a translation in the momentum space $$\vec{p} \to \vec{p} + m \vec{v} \\ \vec{x} \to \vec{x} + t \vec{v} \\ t \to t$$ Therefore we must have $\vec{K} = t\vec{P} - \vec{X} M \simeq t \vec{P} - \vec{X} H$. The conserved charge (i.e. $\vec{K}$) is simply the mass center.

Compatible operators

According to linear algebra, two linear operators can be simultaneously diagonalized iff they commute. $A$ and $B$ are simultaneously diagonalizable, means there exists a matrix $S$, both $S^{-1} A S$ and $S^{-1} B S$ are diagonal. Note that matrix $S$ may not be unique. Similarly, a set of linear operators can be simultaneously diagonalized iff each pair commutes. In quantum mechanics, the the of compatible operators are very useful to identify quantum states. Indeed, if one can find a complete set of compatible operators $\{A_1, A_2, \cdots, A_n\}$, the mutal eigenstates can be identified by the eigenvalues, $\left.|\lambda_1, \lambda_2, \cdots \lambda_n \right>$.

In practice, one of the operator is often chosen to be the Hamiltonian. Therefore, Hamiltonian is block diagonal in the mutual eigenstates of a set of its compatible operators. If the compatible operators is complete, the eigenstate of the Hamiltonian is uniquely determined by their eigenvalues. For example, in the non-relativistic Hydrogen atom problem, a set of the compatible operator is $\{ H, J^2, J_z\}$. The eigenstates can be identified by $\left.|n, l, m \right>$, where the eigenvalues are $\{E_n, l(l+1), m\}$. The eigenvalues of a set of compatible operators for the Hamiltonian are often called "good quantum numbers". Note that there may be more than one set of compatible operators for the Hamiltonian.

Clearly, conserved observables are the natural candidates as compatible operators. The leading example is the momentum operator $\vec{P}$. A counterexample is the boost operator $\vec{K}$, as it does not commute with the Hamiltonian $[\vec{K}, H] = -i \vec{P}$. Another counterexample is $J_x$ or $J_y$ if we have already used $J_z$. For a symmetry group, the Casimir operators are also natural candidates. The leading example is the total angular momentum operator $\vec{J}^2$ from the rotation group. According to what we have analyzed above, the generators of a symmetry is a promising candidate for the compatible operator set for the Hamiltonian. In practice, a degeneracy is the sign of extra compatible operators. When a unexpected degeneracy is detected, it is often said an "unknown" symmetry protects the degeneracy.

On the conserved observable in a truncated subspace

In quantum mechanics, a time-independent observable $A$ is a conserved quantity (constant of motion) if, according to Heisenberg equation, it commutes with the Hamiltonian, $$\frac{\mathrm{d}}{\mathrm{d}t} A = i [H, A] \quad \Rightarrow \quad [H, A] = 0.$$ An interesting question is that if $A$ is still a conserved quantity in a truncated subspace.

Let's $P$ be the projection operator for the truncated subspace, $P^2 = P$. In the truncated subspace, the operators $H$ and $A$ are replaced with $P H P$ and $P A P$ respectively. Then, $$[ P H P, P A P ] = P H P A P - P A P H P = P [H, A] P + P [ [P, H], [P, A] ] P = P [ [P, H], [P, A] ] P$$ Therefore, in the truncated space, the truncated observable is not necessarily a conserved quantity. There are two important exceptions:

1. The truncation is consistent with the eigensubspaces of the Hamiltonian, i.e. $[P, H] = 0$, $A$ is still conserved. This means solving the problem first. 
2. The truncation is consistent with the observable $A$, namely $[A, P] = 0$. 

Example: the shell model. The basis space is generated with 3D harmonic oscillator Hamiltonian $G$. A cutoff is performed on its eigensubspaces (all the degenerate eigensubspaces have to be included). Then the angular momentum $[G, J] = 0$, so $[P, J] = 0$. Therefore, in the truncated subspace, the angular momentum is still a good quantum number (a conserved quantity). 

Many-Body Method for Angular Momentum Addition

Let's start from a simple problem: the total spin of two spin-1/2 particles.

The two particles have total spin $j_1$ and $j_2$, as well as spin projection $m_1$ and $m_2$. Similarly, the composite system will also have a total spin $j$ and a spin projection $m$. The total spin is the sum of the two $\mathbf{J} = \mathbf{J}_1 + \mathbf{J}_2$, according to angular momentum conservation. The (scalar) total spin $J = \sqrt{\mathbf{J}^2}$. In some interactions, such as the $\mathbf{L}\cdot\mathbf{S}$ coupling, total $J$ is conserved, instead of $L$ or $S$. It can be seen from $$\mathbf{L}\cdot\mathbf{S} = \frac{1}{2}\left[ \mathbf{J}^2 - \mathbf{L}^2 - \mathbf{S}^2 \right] \\  = \frac{1}{2}\left[ j(j+1) - l(l+1) - s(s+1) \right]$$ A notable example of the $\mathbf{L}\cdot\mathbf{S}$ coupling is the spin-orbit coupling in the hydrogen-like atom caused by the relativistic correction: $$V(r) = \frac{Ze^2}{8\pi\varepsilon_0 m^2 r^3}\mathbf{L}\cdot\mathbf{S}.$$ However, a product state like $\left.|1/2,+1/2\right>\left.|1/2,-1/2\right>$ (often also denoted as $\left.|\uparrow\downarrow\right>$ for simplicity) is usually not a mutual eigenstate of the total spin operator $J$ and spin projection operator $J_z$ due to the non-commuting nature of the angular momentum operator. According to quantum mechanics, it is in fact the superposition of several eigenstates with various total spin $j$ ( and in principle spin projection $m$ ). The exact expression relating the product states and the eigenstates of the total spin operator is given by Clebsch-Gordan (CG) decomposition: $$\left.|j_1,m_1\right>\left.|j_2,m_2\right> = \sum_{j,m} \left( j_1,m_1,j_2,m_2|j,m\right) \; \left.|j,m\right>$$ or $$\left.|j,m\right> = \sum_{j_1,m_1,j_2,m_2} \left(j,m|j_1,m_2,j_2,m_2\right) \;  \left.|j_1,m_1\right>\left.|j_2,m_2\right>$$
In the case of the spin-1/2 particles, there are four eigenstates categorized into a spin singlet: $$\left.|0,0\right> = \frac{1}{\sqrt{2}}\left( \left.| \uparrow\downarrow\right> - \left.|\downarrow\uparrow\right>  \right)$$ and three spin triplets: $$\left.|1,+1\right>      = \left.| \uparrow\uparrow\right>    \\ \left.|1,0\right> = \frac{1}{\sqrt{2}}\left( \left.| \uparrow\downarrow\right> + \left.|\downarrow\uparrow\right>  \right)    \\ \left.|1,-1\right> = \left.| \downarrow\downarrow  \right>$$
CG coefficients give the analytic solutions to the eigenvalue problem of the two-body total angular momentum operator $J$. In a more general many-body system, it may be more convenient to solve the eigenvalue problem directly (including using numeric methods).

The total spin squared operator
$$\mathbf{J}^2 = \left( \sum_a \mathbf{J}_a \right)^2 = \sum_i \mathbf{J}_a^2 + 2 \sum_{a<b} \mathbf{J}_a\cdot \mathbf{J}_b.$$ where $a, b$ enumerate the constituents. Noting that $2\mathbf{J}_a \cdot \mathbf{J}_b =2 J_a^x J_b^x + 2 J_a^y J_b^y +2 J_a^z J_b^z = J_a^+ J_b^- + J_a^- J_b^+ + 2 J_a^z J_b^z$, where $J^\pm = J_x \pm i J_y$ is the raising/lowering operator. Then the total $J$ operator becomes $$\mathbf{J}^2 = \sum_a \mathbf{J}_a^2 + 2 \sum_{a<b} J^z_a \cdot J_b^z + J_a^+ J_b^- + J_a^- J_b^+.$$ The raising/lowering operator acts on spin states as $$J^\pm \left.| j, m\right> = \sqrt{j(j+1) - m(m \pm 1)}\left.|j, m\pm 1\right>.$$ Suppose $\alpha=(\{j_1,m_1\}, \{j_2,m_2\}, \cdots  )$ and $\beta=(\{j'_1,m'_1\}, \{j'_2,m'_2\}, \cdots)$ are many-body basis states. Then the matrix element is $$\left< \beta \left| \mathbf{J^2}\right| \alpha \right> = \delta_{\alpha, \beta} \sum_a  j_a (j_a+1) + \delta_{\alpha, \beta} 2  \sum_{a<b} m_a m_b  \\ + \sum_{a<b} \delta_{\hat{\alpha}_{ab}, \hat{\beta}_{ab}} \delta_{j_a,j'_a} \delta_{j_b,j'_b} \delta_{m_a, m'_a+1} \delta_{m_b, m'_b-1}\sqrt{\left[j_a(j_a+1)-m_a(m_a+1) \right] \left[j_b(j_b+1) - m_b(m_b-1) \right]}$$ where $\hat{\alpha}_{a,b}$ is the list of quantum numbers except the $a$-th and $b$-th.

We have constructed the matrix elements of the total spin squared operator for the many-body basis. The total $j$ states $\left.|j, m\right>$ are merely the eigenstates of the total spin squared operator. The CG coefficients are an analytical solution of the eigenvalue problem. In general, we obtain the total spin $j$ along with the coefficients by diagonalize the total spin operator. With the many-body matrix in hand, the diagonalization can be done numerically.

This method is particularly suitable for many-body system (particle number larger than 2). It is also superior in system with identical particles. In that case, the basis space shrinks to symmetric/anti-symmetric basis.

The many-body method can also be generalized to $SU(N_c)$ problem, in which the corresponding quantity is often termed as "color". The total color squared operator is the Casimir operator $$\lambda^2 = \sum_c^{N_c^2 - 1} \sum_i \lambda^c_i \sum_j \lambda^c_j = \sum_i \sum_c^{N_c^2-1} \lambda_i^c \lambda_i^c + 2 \sum_{i<j} \sum_c^{N^2_c-1} \lambda^c_i \lambda_j^c$$
where $\lambda^c, c=1,2,\cdots N_c^2-1$ is the generalized Gell-Mann matrix. The matrix elements can be constructed.

It should be mentioned that another systematic method is the Young tableau.

It is also not a surprise that we are looking for irreducible representations from fundamental and/or induced representations.

update:
editing