The unloading behaviour of adhesion-free elastic-perfectly plastic spheres following contact presents complex nonlinear features. Analytical models capable of accurately predicting this response have not yet been developed for an extensive range of material properties and initial deformation states, and consequently the use of semi-empirical models requiring calibration is widespread in the practical application of contact laws. In this work, we provide insight into contact behaviour during unloading by conducting a finite element study to characterise this response for a comprehensive range of material properties (1 <= E/sigma(y) <= 1000, 0.0 <= v <= 0.45), for particles that have undergone large deformation prior to unloading (0.01 <= d/R <= 0.5), leading to the following findings: Firstly, an empirical relation capable of accurately determining secant unloading stiffness from material properties and degree of initial deformation was formulated, which was expressed in nondimensional form for maximum generality. An analytical model was also developed to help explain some of the contributing mechanisms identified from the finite element analysis. Secondly, the nonlinearity of the force-displacement curve in unloading was quantified and charted, and physical arguments were advanced to explain the trends revealed. Considering both stiffness and nonlinearity results, it was concluded that a single synthetic measure of initial particle deformation relative to deformation at first yield, which is currently used, is insufficient to characterise unloading response at large displacements. The unloading relations developed can be employed with static and dynamic multi-particle simulation approaches such as the Discrete Element Method (DEM) for more accurate simulation of compaction and flow of dense powder beds, and problems reliant on accurate determination of contact areas after unloading between particles following large deformation. (C) 2020 Elsevier B.V. All rights reserved.