We have calculated the real part χ ′ of the nonlinear dielectric susceptibility of amorphous insulators in the kHz range, by using the two-level system model and a nonperturbative numerical quantum approach. At low temperature T , it is first shown that the standard two-level model should lead to a decrease of χ ′ when the measuring field E is raised, since raising E increases the population of the upper level and induces Rabi oscillations cancelling the ones induced from the ground level. This predicted E-induced decrease of χ ′ is at odds with experiments. However, a good agreement with low-frequency experimental nonlinear data is achieved if, in our fully quantum simulations, interactions between defects are taken into account by a new relaxation rate whose efficiency increases as √ E, as was proposed recently by Burin et al. (Phys. Rev. Lett. 86, 5616 (2001)). In this approach, the behavior of χ ′ at low T is mainly explained by the efficiency of this new relaxation channel. This new relaxation rate could be further tested since it is shown that it should lead: i) to a completely new nonlinear behavior for samples whose thickness is ≃ 10 nm; ii) to a decrease of nonequilibrium effects when E is increased. PACS numbers: 61.43.Fs,77.22.Ch,72.20.Ht Amorphous materials exhibit universal anomalous properties at low temperature. In 1971, Zeller and Pohl [1] discovered below 1 K a quasilinear behavior of the specific heat in a number of glasses contrasting with the Debye law of crystalline materials. Anderson, Halperin, Varma [2] and Phillips [3] proposed an explaination based upon the existence of localized two-level systems (TLS). Their origin may be due to the tunneling of atoms or groups of atoms between two equilibrium positions separated by a narrow energy barrier featuring asymmetric two-well potentials. They are assumed randomly distributed in energy splittings and tunneling barriers as a consequence of the structural disorder of these materials. This model has proven to be successful to understand most salient experimental features.The standard TLS model assumes defects do not interact with one another. However, defects are strongly coupled to their environment and can emit or absorb phonons. It leads to an indirect interaction between nearest neighbors via the phonon field [4]. Certain recent failures to explain nonequilibrium data (at a few kHz) [5] underscore the likely involvement of these interactions below 100 mK. However, these nonequilibrium effects are small corrections of the kHz stationnary response, and, up to recently, examples of stationnary susceptibilities strongly affected by interactions were very rare : in the kHz regime, it was argued that the ultra-low-T (T ≃ 1 mK) plateau of the dielectric constant in the linear regime, strongly different from the expected logarithmic increase, resulted from interactions [6]. Very recently, such a conclusion was drawn from internal friction experiments [7].In this work, we show that including interactions in the TLS model with a recently proposed m...