As of #54 the simplest scalar calls to jitted special functions should work.

However there's no support yet for array-valued inputs:

import numpy as np
from numba import njit
from scipy import special

x = np.linspace(-10, 10, 1000)

@njit
def jitted_j0(x):
    res = special.j0(x[0])  # works after PR #54
    # res = special.j0(x)  # breaks
    return res

print(jitted_j0(x))

This is not obviously a shortcoming, since looping in jitted functions should be alright. So this is just a mild suggestion to consider adding support for array-valued signatures. (This should probably be preceded with some benchmarks to see whether this would help anything performance-wise.)