  Ultra-large scales close to the cosmological horizon will be probed by the
upcoming observational campaigns. They hold the promise to constrain
single-field inflation as well as general relativity, but in order to include
them in the forthcoming analyses, their modelling has to be robust. In
particular, general relativistic effects may be mistaken for primordial
signals, and no consensus has emerged either from analytical modelling nor from
the numerical route, obstructed by the large dynamical range to be simulated.
In this work, we present a numerical technique to overcome the latter
limitation: we compute the general relativistic displacement field with the
N-body relativistic code gevolution and combine it with the accurate Newtonian
simulation Gadget-4. This combination leads to an effective simulation
reproducing the desired behaviour at the level of the matter power spectrum and
bispectrum. We then measure, for the first time in a simulation, the
relativistic scale-dependent bias in Poisson gauge; at redshift $z=0$, we find
$b_1^{\mathrm{GR}}=-5.7 \pm 1.7$. Our results at the field level are only valid
in the Poisson gauge and need to be complemented with a relativistic ray
tracing algorithm to compute the number count observable.
