The reversible submonolayer heteroepitaxial growth is studied by means of kinetic Monte Carlo simulations of a simple model with short-and long-range interactions. The Green function of a film on a deformable substrate allows a fast computation of elastic energy barriers and the investigation of the islands statistical properties as functions of strain and the diffusion to deposition flux ratio. We find an unexpected decrease of the island density due to elastic interactions which can be understood as a decrease of the effective adatom binding energy coupled to kinetic effects. The island size distribution satisfies the usual scaling and does not significantly depend on elasticity.Aggregation phenomena are ubiquitous in physics, chemistry, or biology. A paradigm arising in surface science deals with island formation during molecular beam epitaxy, which involves adatom deposition and diffusion, island nucleation, growth, and coalescence. In heteroepitaxial systems, elasticity is an extra crucial effect which can lead to a transition between two-͑2D͒ and three-͑3D͒ dimensional islands and quantum dot formation. 1 The progress in experimental techniques for surface scanning, together with application for electronic devices fabrication has shed some new light on this field. 2-4 The understanding of the density, mean size, and size dispersion of submonolayer strained islands is then of significant importance as they are templates for the subsequent surface evolution. Moreover, beyond aggregation mechanisms, the submonolayer regime may be used to measure adatom diffusion coefficients otherwise hard to measure. 5,6 To tackle this far-from-equilibrium many-body problem, different approaches are available without elasticity in the irreversible case where adatoms never detach from islands. 7,8 Scaling analysis 9,10 first predict that the number N s of islands with s atoms is merelygiven by the scaling function f and the only characteristic island size ͗s͘. They depend solely on the ratio of the adatom diffusion constant D to the deposition flux F, and ͗s͘ = ⌰ z ͑D / F͒ , where ⌰ is the deposited coverage ⌰ = Ft, and z and -two scaling exponents. More detailed information are then found from rate equations which treat the system in a mean-field way 11-14 and compare favorably with kinetic Monte Carlo ͑KMC͒ simulations. However, even if the irreversible hypothesis is well suited for metallic materials at low temperatures, it fails for semiconductors or at high temperatures when E N , the atom binding energy to islands, is no longer large compared to the thermal energy k B T. A simple solid-on-solid ͑SOS͒ model accounting for reversible aggregation 15 reveals that the system depends both on D / F and E N , and exhibits a different regime where the island density saturates before coalescence occurs. In another de-scription, islands with more than the critical nucleus i atoms do not loose matter and are characterized by both i and D / F ͑see Ref. 13 and 16͒. For strained systems with a misfit between the film and substrate, the...