At every step \(j\)  of the potential retrieval procedure, we perform a full Monte Carlo simulation of the system, based on the interaction potentials \(V_\mathrm {same}^{(j)}(r)\)  and \(V_\mathrm {opp}^{(j)}(r)\) .
Given these pair interactions of the system, we run a Monte Carlo simulation to compute its final state at equilibrium [1]}.
The Monte Carlo simulation is carried in a square region of side \(L = 10\)  (\(\lambda = 1\) ) with periodic boundary conditions. We note that we performed preliminary simulations with different sizes of the simulation box (\(L = 5,\,10,\,20\) ) and verified that this size does not affect the final results.
The simulation box contains a number of particles \(N=\rho / L^2 = 314\) , where \(\rho \)  is the expected average density of singularities in isotropic random waves [2]}. To initialize the simulation, such particles are arranged in a square lattice with a pitch \(a \approx L/\sqrt{N}\) , perturbed with additive random displacements \(\Delta x\)  and \(\Delta y\)  uniformly distributed in \([-a/4,a/4]\) .
Subsequently, the simulation is carried by proposing at each Monte Carlo step a displacement \(\Delta r\)  for the particle located at \(r_i\) .
We compute the difference in energy \(\Delta E\)  corresponding to this change
\(\Delta E = \sum _{k\;\in \;\mathrm {system}} \left[V^{(j)}(r_i + \Delta r - r_k) - V^{(j)}(r_i - r_k) \right] ,\) 
