Super-ohmic bosonic bath correlation function

In quantum Brownian motion, bosonic/harmonic oscillator bath and interaction described by Hamiltonian

$$
H_B = sum_{n}hbaromega_n(b_n^dagger b_n) \
H_I = -sigma_x otimes B
$$

and

$$ B = sum_n kappa_nsqrt{frac{hbar}{2m_nomega_n}}(b_n+b_n^dagger)$$

We have two correlation function to describe time evolution of system (See Breuer p.174)

$$
D(tau) equiv iTr_R([B,B(-tau)]rho_B) \
D_1(tau) equiv Tr_R({B,B(-tau)}rho_B)
$$

And using Spectral Density (approximating infinite mode) we get (setting $hbar=1$)

$$
D(tau)=2int_0^infty domega , J(omega)sin(omegatau) \
D_1(tau)=2int_0^infty domega , J(omega)coth left(frac{hbaromega}{2k_BT} right)cos(omegatau)
$$

In many references (including Breuer) ohmic spectral density with Drude-Lorentz cutoff is used
$$
J(omega)=alpha^2omegafrac{lambda^2}{omega^2+lambda^2}$$

And correlation function can be easily analytically calculated.

My question is, if we define more general spectral density

$$J(omega)=alpha^2omega^somega_{ph}^{1-s}frac{lambda^2}{omega^2+lambda^2}$$

With $S = 1$ (Ohmic), $S > 1$ (Super-Ohmic), $S < 1$ (Sub-Ohmic). In Sub-ohmic region, working the integral for $D(tau)$ in Mathematica gives Hypergeometric function which is more complicated than Ohmic case but still a well behaved function. The problem arise when choosing spectral density as a Super-Ohmic

$$
D(tau)=2alpha^2lambda^2omega_ph^{1-s}int_0^infty frac{omega^s sin(omegatau)}{omega^2+lambda^2}$$

For $S > 1$ the integral is divergent. From this result, is the boson bath "limited" to the ohmic and sub-ohmic spectral density(is there any microscopical model to suggest this) ? How the divergence of correlation function can be interpret in relation to super-ohmic ?