We train graph neural networks on halo catalogues from Gadget N-body
simulations to perform field-level likelihood-free inference of cosmological
parameters. The catalogues contain $\lesssim $ 5,000 halos with masses $\gtrsim {10}^{10}\text{}{h}^{-1}{M}_{\odot}$ in a periodic volume of $(25\text{}{h}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}{)}^{3}$ ; every
halo in the catalogue is characterized by several properties such as position,
mass, velocity, concentration, and maximum circular velocity. Our models, built
to be permutationally, translationally, and rotationally invariant, do not
impose a minimum scale on which to extract information and are able to infer
the values of ${\mathrm{\Omega}}_{\mathrm{m}}$ and ${\sigma}_{8}$ with a mean relative error of
$\sim 6\mathrm{\%}$ , when using positions plus velocities and positions plus masses,
respectively. More importantly, we find that our models are very robust: they
can infer the value of ${\mathrm{\Omega}}_{\mathrm{m}}$ and ${\sigma}_{8}$ when tested using halo
catalogues from thousands of N-body simulations run with five different N-body
codes: Abacus, CUBEP${}^{3}$ M, Enzo, PKDGrav3, and Ramses. Surprisingly, the model
trained to infer ${\mathrm{\Omega}}_{\mathrm{m}}$ also works when tested on thousands of
state-of-the-art CAMELS hydrodynamic simulations run with four different codes
and subgrid physics implementations. Using halo properties such as
concentration and maximum circular velocity allow our models to extract more
information, at the expense of breaking the robustness of the models. This may
happen because the different N-body codes are not converged on the relevant
scales corresponding to these parameters.