Independent boson model with an arbitrary finite-dimensional impurity The independent boson model consists of the following Hamiltonian:
$$ H_s = E \sigma^z $$
$$ H_b = \sum_k \omega_k b^{\dagger}_kb_k $$
$$H_{sb} = \sigma^z \sum_k (g_k b_k + g_k^{\ast}b^{\dagger}_k).$$
The model describes a single spin-1/2 impurity with Pauli operators $\sigma^{x,y,z}$ linearly coupled to an infinity of bosonic modes $b_k$. Importantly, the interaction $H_{sb}$ commutes with $H_s$. 
The model is exactly solvable by introducing a state-dependent displacement:
$$ U = \exp \left[ \sigma^z \sum_k (g_k^{\ast}b^{\dagger}_k - g_k b_k)\right],$$
leading to the transformed Hamiltonian
$$U H U^{\dagger} = E^{\prime} \sigma^z + \sum_k \omega_k b^{\dagger}_kb_k + \mathrm{const.}$$
where $E^{\prime}$ is the renormalised impurity energy. Similar tricks allow one to compute time evolution etc. The solutions can be found in detail in Mahan's book Many-Particle Physics.
Note that there exists an equivalence between a spin-1/2 particle and a single fermionic mode, i.e. we can rewrite the above Hamiltonian by replacing $\sigma^z \to c^{\dagger} c$, where $\{c,c^{\dagger}\} = 1$ are fermionic ladder operators. The resulting model is equivalent up to a shift of the equilibrium position of the oscillators. 
However, when $c$ is instead taken to be bosonic, the solution fails. The fermionic/spin solution relies on the fact that $(c^{\dagger}c)^2 = c^{\dagger}c$, which ultimately stems from the fact that the fermionic Hilbert space has 2 states. In contrast, the Hilbert space of a bosonic mode is infinite-dimensional. 
Is the independent boson model always exactly solvable so long as the Hilbert space of the impurity is finite-dimensional?
I mean precisely the following: imagine replacing $\sigma^z$ with $S^z$, the $z$ projection of a spin with total angular momentum $S > 1/2$. Is the model exactly solvable? Constructive answers which describe the form of the solution would be great, or any references to where this problem has already been solved in the literature.
 A: I'm not completely sure why you think that the bosonic mode fails, but it seems to me that the answer is definitely yes. The system is solvable in both the finite-dimensional and the bosonic case; the problem with the bosonic case is that the solution is ugly, because the hamiltonian is ugly.
Take a hamiltonian of the form
$$
H=E_0S+S\sum_k (g_k b_k+g^*b_k^\dagger)+\sum_k\omega_k b_k^\dagger b_k,
$$
where $S$ is so far unspecified. Let $|s⟩$ be an eigenstate of $S$, with $S|s⟩=s|s⟩$, and look for an eigenstate of the form $|\psi⟩=|s⟩|\phi⟩$. This eigenstate ought to exist because of the structure of the hamiltonian, but for now you can simply look at the action of $H$ on states of that form:
\begin{align}
H|\psi⟩
&=
\left[\left(E_0+\sum_k (g_k b_k+g^*b_k^\dagger)\right)s+\sum_k\omega_k b_k^\dagger b_k\right]|s⟩|\phi⟩
\\ &=
\left[E_0s-\sum_k\omega_k \delta_k\delta_k^* +\sum_k\omega_k (b_k+\delta_k)^\dagger(b_k+\delta_k)\right]|s⟩|\phi⟩,
\end{align}
where $\delta_k=sg_k/\omega_k$.
It is now easy to find eigenstates of this form - all you need is to take $|\phi⟩$ as a product of eigenstates of each displaced number operator $(b_k+\delta_k )^\dagger (b_k+\delta_k)$. This is easy to write as
\begin{align}
|\phi⟩
&=
\bigotimes_k D(-\delta_k)|n_k⟩
=
\exp\left(\sum_k(\delta_k^*b_k-\delta_k b_k)\right)|\{n_k\}⟩
\\&=
\exp\left(s\sum_k\frac{g_k^*b_k-g_k b_k}{\omega_k}\right)|\{n_k\}⟩.
\end{align}
This eigenstate then has energy
$$
E=sE_0-s^2\sum_k |g_k|^2/\omega_k +\sum_k\omega_k n_k.
$$
In some ways, you're sort of done. This is enough to provide a basis of eigenstates of $H$, and there is definitely no problem if $S$ is finite-dimensional. On the other hand, I imagine you still want a canonical-transformation formulation of this. To do that, you can write
$$
|\psi⟩=\exp\left(S\sum_k\frac{g_k^*b_k-g_k b_k}{\omega_k}\right)|s⟩|\{n_k\}⟩=U|s⟩|\{n_k\}⟩,
$$
via a unitary transformation
$$
U=\exp\left(S\sum_k\frac{g_k^*b_k-g_k b_k}{\omega_k}\right).
$$
This means that the transformed hamiltonian acts as
\begin{align}
U^\dagger HU|s⟩|\{n_k\}⟩
&=
U^\dagger H|\psi⟩
=
\left[sE_0-s^2\sum_k |g_k|^2/\omega_k +\sum_k\omega_k n_k\right]U^\dagger |\psi⟩
\\ & =
\left[sE_0-s^2\sum_k |g_k|^2/\omega_k +\sum_k\omega_k n_k\right]|s⟩|\{n_k\}⟩
\\ & =
\left[SE_0-S^2\sum_k |g_k|^2/\omega_k +\sum_k\omega_k b_k^\dagger b_k \right]|s⟩|\{n_k\}⟩
\end{align}
on a basis, and therefore has that same action everywhere:
\begin{align}
U^\dagger HU
=
SE_0-S^2\sum_k |g_k|^2/\omega_k +\sum_k\omega_k b_k^\dagger b_k.
\tag1
\end{align}
This works regardless of the dimensionality of $S$, as far as I can tell. I can't completely rule out funky business with the unitarity of $U$ in the infinite-dimensional case but I honestly can't see where the problem would come in from.

That said, the bosonic case $S=a^\dagger a$ is obviously problematic, because the transformed hamiltonian in $(1)$ is unbounded from below. as soon as you have one nonzero $g_k$. However, I think this is a problem with $H$ itself rather than the transformation.
To see this, consider the case where a single $g_k$ is nonzero, so
$$
H=\omega_s a^\dagger a+(gb+g^*b^\dagger)a^\dagger a +\omega_b b^\dagger b.
$$
Consider further the state $|\chi⟩=|n⟩|-n g/\omega_b⟩$, i.e. the number state $|n⟩$ tensor-times a coherent state at $-ng/\omega_b$ for the bath, and calculate its energy expectation value:
\begin{align}
⟨\chi|H|\chi⟩
&=
\omega_s n
+n(g(-ng/\omega_b)+g(-ng/\omega_b))
+\omega_b|ng/\omega_b|^2
=
\omega_s n-n^2g^2/\omega_b.
\end{align}
This is arbitrarily negative for sufficiently large $n$, which proves that there is no ground state with finitely negative energy.
So: the problem with a boson-boson interaction of this form is not that the hamiltonian is not diagonalizable (which it patently is), but that the interaction is unphysical.
