How to implement filter poles in noise shaping for ΔΣ modulator?
I believe I have actually found an answer myself after even more research ;)
The crucial piece of information is the following: The filter function H(z) shown in many models is not the actual Noise Filter Function. Usually the noise filter should be some kind of high pass and I used this directly as H(z) erroneously.
Instead, the noise is shaped by the noise transfer function NTF(z) = 1-H(z).
Therefore, after designing the highpass noise transfer filter NTF(z), one needs to calculate the filter model via H(z) = 1-NTF(z).
It turns out that, if the NTF is a &quot;poleless&quot; filter (FIR filter) of form NTF(z) = 1 + Sum(b_n z^-n), then the filter function will be very similar: H(z) = 0 - Sum(b_n z^-n) where n starts from 1. The first coefficient will be 0.
I actually implemented this noise shaping directly in the past, by ignoring the leading coefficient 1 and inverting the sign of the filter result at the summing node with the signal. I came to this conclusion, because I picked it up somewhere but didn't understand its origins.
When the NTF contains poles, however, the H(z) will deviate more significantly. The denominator of NTF(z) and H(z) will be in fact identical. But the numerator will be X(z) = Sum(a_n z^-n) - Sum(b_n z^-n), where a_n are the pole coefficients, and b_n are the zero coefficients. As typically, the leading coefficients of the NTF(z) will be both 1 ( it is a minimum phase filter after all), the leading coefficient of H(z) will become 0. This produces the missing first coefficient naturally.
With this properly defined H(z), following the flow chart as given e.g. in Lipshitz 1991 paper, gives correct results. I used +/- 0.5 of the quantizer range as uniform white dither.
Update:
after correcting some minor mistakes, it all works as it should. In Python I use the Inverse Chebyshev Filter model to make the Noise Transfer Function. Then I correct the coefficients, to turn them into the H(z) filter model:
fDAC=20e6
b, a = signal.cheby2(5, 160, 20000, 'hp', fs=fDAC)
b /= b[0] # normalize according to Gerzon/Craven
b = a - b # turn NTF(z) into H(z) for noise shaping

Then the actual noise shaping (left some redundancy in the code, to keep it closer to Lipshitz 1991):
l = 200000
ts = np.arange(l)/fDAC
wav = 0.55 * np.sin(2*np.pi*3200*ts)
errbuf = np.zeros_like(b)
samps = np.zeros_like(wav)
for iSAM in range(l):
    desired = wav[iSAM] - np.dot(b, errbuf)
    samps[iSAM] = np.sign( desired )
    error = samps[iSAM] - desired 
    errbuf[0] = - error
    errbuf[0] = - np.dot(a, errbuf)
    errbuf = np.roll(errbuf, 1)

Note that no additive dither is being used. Even though the quantizer noise is not white itself, the noise shaping is efficient enough to remove the quantizer distortion. Undithered Delta Sigma modulators do have drawbacks tough, e.g. on stationary input, as can be read in numerous sources. If you want dither, it should be added right before quantization, i.e. inside the sign() function

Figure: Example of a 5th order delta sigma modulator. Vertical scale unit is full scale. Horizontal is Frequency in Hz. The orange curve is a designed 5th order inverse Chebyshev Highpass filter as Noise transfer function. Blue is the result with a ~-5dBFS sine wave.
PS: Since I originally set a bounty on this question, which was keeping me up at night: Please feel free to correct me and claim the bounty with a better answer of your own.