A Monte Carlo investigation of the thermal variation of Dy/Tb superlattice magnetization in comparison with experiment is carried out. Our model consists of an alternate stacking along the c-axis of hexagonal compact Dy and Tb layers with abrupt interfaces. Each site is occupied by a classical Heisenberg spin with exchange interactions up to the 4th neighbors in Dy layers to reproduce the low temperature helimagnetic phase whereas Tb is ferromagnetic. Planar anisotropy and interaction with an applied field are also taken into account. Our results indicate that the low temperature magnetization per spin strongly decreases as the number of bilayers increases. For an applied field of 1.04 kOe, the maximum around 200 K is quite well reproduced. The visualization of the magnetic structures at different temperatures provides an explanation of the complex thermal variation of the magnetization.Introduction Magnetism of bulk rare earth (RE) exhibits many magnetic phases [1][2][3][4]. In particular, because of frustration, helimagnetic phases can be observed. For example, Dy [5] and Tb that crystallize in the hexagonal compact (hcp) structure exhibit two temperature-driven magnetic transitions: a paramagnetic-helimagnetic transition around 180 K (Dy) and 230 K (Tb) and a helimagnetic-ferromagnetic transition around 90 K (Dy) and 220 K (Tb). In the helimagnetic phase, the moments of the (a, b) planes are parallel (ferromagnetic planes) and the angle between the moments of two adjacent planes, which is called the "turn angle", increases with temperature. Then, RE used in particular structures, such as superlattices, should exhibit very interesting magnetic structures depending on the layer thicknesses. Indeed, a complex thermal variation of the magnetization of Dy/Tb superlattices ( Fig. 1) has been observed [6] and still has to be precisely interpreted.Most of the numerical investigations devoted to realistic magnetic multilayers have been carried out using mean field theories (see Refs. [7-10], for example). These techniques, which totally or partially neglect the spin fluctuations, require many approximations in the model to be solved. Concerning RE/RE systems, thermal variation of the magnetization of Dy/Gd superlattices has been calculated and compared with experimental results [11]. To our knowledge, very few investigations on magnetic multilayers have been performed by numerical simulations [12,13] and none of them are concerned with RE/RE systems. In this paper, we have chosen the Monte Carlo method, which is a powerful technique to simulate complex models. Our aim is to build a model in order to reproduce and explain experimental results on Dy/Tb superlattices. In the following section, we describe the model and the simulation technique. The numerical results compared with the experimental ones are presented in Section 3. Our conclusion and perspectives are given in Section 4.