The effects of subgrid scale (SGS) motions on the dispersion of heavy particles raise a challenge to the large-eddy method of simulation (LES). As a necessary first step, we propose the use of a stochastic differential equation (SDE) to represent the SGS contributions to the relative dispersions of heavy particles in LES of isotropic turbulence. The main difficulty is in closing the SGS-SDE model whilst accounting for the effects of particle inertia, filter width and gravity. The physics of the interaction between heavy particles and SGS turbulence is explored using the filtered direct numerical simulation method. It is found in the present work that (i) the ratio of the SGS Lagrangian and Eulerian timescales is different from that of the full-scale Lagrangian and Eulerian timescales. The ratios are also dependent on filter widths. (ii) In the absence of gravity, the SGS timescale seen by heavy particles non-monotonically changes with particle Stokes number and has a maximum at particle Stokes number (St = τ p /δT E ) near 0.5. (iii) In the presence of gravity, a similarity law exists between the SGS Lagrangian correlation function seen by a heavy particle within a timedelay τ and the SGS spatial correlation function with the displacement w τ , where w is the average settling velocity of a heavy particle. The joint effects of particle inertia and gravity are accounted for using the elliptic model for pair fluctuations to the dispersion of particles [9,10]. Chibbaro and Minier [10] addressed the importance of introducing specific features related to the near-wall coherent structure in pipe flows based on the Langevin probability density function (PDF) method. In the SDE, the Lagrangian integral timescale of fluid velocity seen by heavy particles (the fluid-particle interaction timescale), T Lp , is one of the key parameters and it determines the dispersion rates of heavy particles [11][12][13][14]. Direct numerical simulation (DNS) shows that T Lp does not vary monotonically with particle Stokes number St K from T L to T E , but has an 'N' shape with a maximum near St K = 1 due to the relative motion of particles near rotational vortical structures [15][16][17], where St K ≡ τ p /τ K , τ p is the particle Stokes response time and τ p = ρ p d 2 p /18µ, d p is the particle diameter, µ the dynamical viscosity, ρ p particle density and τ K the turbulent Kolmogorov timescale.The RANS method finds wide application in engineering calculations [9] and the DNS method has manifested its power in fundamental research [17][18][19][20]. However, DNS is still limited to flows at relatively low or moderate Reynolds numbers. RANS does not predict well the unsteady structures in turbulent flows which are of importance for turbulent dispersion of heavy particles. Motivated by the limitations of RANS and DNS, large-eddy simulation (LES) is becoming a potential method for future engineering calculations. In LES, large-scale velocitiesū(x, t) are resolved explicitly using the filtered Navier-Stokes equations. The smallscale velocities...