2 neighbor outgoing links to a different site chosen at ran- dom. After repeating this process for every link, we are left with a network with a density q of random two- dimensional lattice directed links. Therefore, with this procedure every site will have exactly the same amount of undirected random two-dimensional lattice outgoing links and a varying (random) number of incoming links. Model - We consider now the two-dimensional spin-1 Blume-Capel model on these Poissonian random lattices. The Blume-Capel Model is a generalization of the stan- dard Ising model 10 and was originally proposed for spin- 1 to account for first-order phase transition in magnetic systems $^{11,12}$. The Hamiltonian can be written as where the first sum runs over all nearest-neighbor pairs of sites (points in the Voronoi construction) and the spin- 1 variables S$_{i}$ assume values ± 1 , 0. In eq. (1) J is the exchange coupling and Δ is the single ion anisotropy pa- rameter. The second sum is taken over the N spins on a D -dimensional lattice. The case where S = 1 has been extensively studied by several approximate techniques in two- and three-dimensions and its phase diagram is well established $^{11–17}$. The case S > 1 has also been investi- gated according to several procedures $^{18–24}$. Simulations - The simulations have been performed for Δ = 0, which is the simplest case and where we have some results in the literature to compare with. The lattices comprised different sizes with N ranging from N = 250 , 500 , 1000 , 2000 , 4000 , 8000, and 16000, where N is total umber of sites. For simplicity, the linear length of the system is defined here in terms of the size of a reg- ular two-dimensional lattice L = N 1 / $^{2}$. For each system size quenched averages over the connectivity disorder are approximated by averaging over R = 100 ( N = 500 to 4000), R = 50 ( N = 8000) and R = 25 ( N = 16000) independent realizations. For each simulation we have started with a uniform configuration of spins (the results are, however, independent of the initial configuration). We ran 2 . 52 × 10 6 Monte Carlo steps (MCS) per spin with 1 . 2 × 10 5 configurations discarded for thermaliza- tion using the “perfect” random-number generator $^{25}$. In both cases we have employed the heat bath algorithm and for every 12th MCS, the energy per spin, e = E/N , and magnetization per spin, m = ∑ $_{i}$S$_{i}$ /N , were measure and recorded in a time series file. From these thermody- namic quantities the behavior of the transition can be analyzed. From the series of the energy measurements we can compute, by re-weighting over a controllable temperature interval Δ T , the average energy, the specific heat, and the energy fourth-order cumulant which are respectively given by In order to get the transition temperature, as well as to determine the order of the transition of this model, we apply the finite-size scaling (FSS) procedure. Initially, we search for the minima of the energy fourth-order cumu- lant given by Eq. (4). This quantity gives a qualitative, as well as a quantitative, description of the order of the transition $^{26}$. For instance, it is known27 that this param- eter takes a minimum value B$_{min}$ at the effective transi- tion temperature T$_{c}$ ( N ). One can also show 28 that for a second-order transition, lim$_{N}$ →∞ (2 / 3 − B$_{min}$ ) → 0, even at T$_{c}$ , while at a first-order transition the same limit mea- suring the same quantity is small and (2 / 3 − B$_{min}$ ) = 0. A more quantitative analysis can be carried out through the FSS of the energy fluctuation C$_{max}$ , the sus- ceptibility maxima χ$_{max}$ , and the minima of the Binder cumulants B , U$_{2}$ and U$_{4}$ . If the hypothesis of a first-order phase transition is correct, we should then expect, for large systems sizes, an asymptotic FSS behavior of the form $^{29,30}$, where K = J/k$_{B}$T , with a normalized exchange J = 1, and k$_{B}$ is the Boltzmann constant. In the above equa- tions < ... > stands for thermodynamic averages and [ ... ]$_{av}$ for averages over the different realizations. Sim- ilarly, we can derive from the magnetization measure- ments the average magnetization, the susceptibility, and the magnetic cumulants, Further useful quantities involving both the energy and magnetization are their derivatives H = − J ∑ <i,j> S$_{i}$S$_{j}$ + Δ ∑ i S 2 $_{i}$, (1) B ( K ) = [1 − < e 4 > 3 < e 2 > $_{2}$]$_{av}$ , (4) m ( K ) = [ < | m | > ]$_{av}$ , (5) χ ( K ) = KN [ < m 2 > − < | m | > $^{2}$]$_{av}$ , (6) U$_{2}$ ( K ) = [1 − < m 2 > 3 < | m | > $_{2}$]$_{av}$ , (7) U$_{4}$ ( K ) = [1 − < m 4 > 3 < | m | > $_{2}$]$_{av}$ . (8) d [ < | m | > ]$_{av}$ dK = [ < | m | E > − < | m | >< E > ]$_{av}$ , (9) d ln[ < | m | > ]$_{av}$ dK = [ < | m | E > < | m | > − < E > ]$_{av}$ , (10) d ln[ < | m $^{2}$| > ]$_{av}$ dK = [ < | m $^{2}$| E > < | m $^{2}$| > − < E > ]$_{av}$ . (11) u ( K ) = [ < E > ]$_{av}$ /N, (2) C ( K ) = K $^{2}$N [ < e 2 > − < e > $^{2}$]$_{av}$ , (3) C$_{max}$ = a$_{C}$ + b$_{C}$N + ... (12) χ$_{max}$ = a$_{χ}$ + b$_{χ}$N + ... (13)