4 20 40 60 80 100 120 140 T  (K) (a.u.) 0 0.2 0.4 0.6 0 0.5 1 1.5 2 c v  / k B 20 40 60 80 100 120 140 T  (K) 0.1 0.2 0.3 0.4 0.5 0.6 u 4 132 134 136 T (K) (a.u.) (a) (c) (b) (d) u 4  (M) 0.62 0.64 0.66 130 131 132 133 134 135 136 137 138 T  (K) 0.5 0.55 0.6 0.65 u 4  (M HM ) 0 0.2 0.4 0.6 0.8 q$_{z}$ S(q z  )  (a.u.) (a) (b) FIG. 4: (color online) Thermodynamic quantities obtained for thickness n = 8 in the temperature range 0-150 K. Colors and symbols as in Fig. 2. (a) : specific heat; (b) : chirality order parameter. (c) : susceptibility χ$_{κ}$ . (d) : Binder cumulant for κ . FIG. 3: (color online) Binder cumulants at thickness n = 16, colors as in Fig. 2. (a) : Binder cumulant for the order parameter defined in Eq. (3). (b) : Binder cumulant extracted from the integral of the structure factor (see Sec. II). Inset: structure factor for L = 64 between T = 131 K (upper curve) and T = 140 K (lower), with 1 K temperature step. 133 . 1(3) K (see Fig. 3b). The crossing points of the Binder’s cumulants of the helical order parameter imme- diately appear to be located, within the error bars, at the same temperature of those for the average magnetization previously discussed. In addition, it is worthwhile to ob- serve that the peak evolution of S (0 , 0 , q$_{z}$ ), in particular close to T$_{N}$ (16) (inset of Fig. 3b), displays the typical behaviour expected for an helical structure. We can thus conclude that for n = 16, as it is commonly observed in bulk samples, the establishment of the in-plane order coincides with onset of the perpendicular helical arrange- ment at T$_{N}$ (16). However, due to helix distortion in the surface regions, the maximum of S (0 , 0 , q$_{z}$ ) stabilizes at values of q$_{z}$ sensibly smaller (e.g. Q$_{z}$ ( T$_{N}$ (16)) ≈ 16 $^{◦}$, and Q$_{z}$ ( T = 10 K ) ≈ 28 $^{◦}$) with respect to the bulk one ( Q Ho z = 30 . 5 $^{◦}$). The results obtained by MC simulations of the model introduced in Sec. II will be presented starting from n = 16, i.e. the highest investigated film thickness which still displays a bulk-like behaviour. In Fig. 2 the spe- cific heat for samples with n = 16 and lateral dimension L = 24 , 32 , 48 , 64 is shown. The location of the specific heat maximum shows a quite definite evolution toward the bulk transition temperature, T Ho N ≃ 132 K 10 (it is worthwhile to note that for this XY model the mean field theory predicts a critical temperature T Ho N,MF ≃ 198 K). The intensity of the maximum of c$_{v}$ has been analyzed by the MH technique for the same lateral dimensions (see inset of Fig. 2): it clearly appears as it increases with L in a smooth way. III. RESULTS The MC simulations outcomes for n = 16 we just pre- sented appear quite different with respect to those ob- tained at the same thickness for the model with six cou- pling constants along the z direction $^{14,15}$. Indeed, for the J$_{1}$ - J$_{2}$ model here investigated, we observe that all layers order at the same temperature, and we do not find any hint of the block-phase, with inner disordered planes intercalated to antiparallel quasi -FM four-layer blocks, previously observed; sample MC runs we made using the same hcp lattice employed in Refs. 14,15 shows that the presence or absence of the block phase is not related to the lattice geometry, but it is a consequence of the inter- action range only. The Binder cumulant for the average order parameter defined in Eq. (3) was obtained close to the c$_{v}$ peak and is reported in Fig. 3a; its analysis leads to an estimate of the critical temperature of the sample (given by the location of the common crossing point of the different curves re- ported in the figure) of T$_{C}$ (16) = 133 . 2(5) This value can be considered in a rather good agreement with the exper- imental ordering temperature of Holmium T Ho N , the rel- ative difference being about 1%. Even such a mismatch between T Ho N and T$_{C}$ (16) could be completely eliminated by slightly adjusting the in-plane coupling constant J$_{0}$ , but, as discussed in Sec. II, we shall preserve the value reported in Refs. 13, and 12 in order to allow for a correct comparison with the results reported in those papers. We now move to describe and discuss MC simulation data for thinner samples. A graphical synthesis of the results obtained for n = 8 in reported in Fig. 4a-d. The specific heat c$_{v}$ , shown in Figs. 4a, reveals very small finite-size effects, which, however, cannot be unambigu- ously detected for the largest lattice size ( L = 64), as they fall comfortably within the error range. Surpris- ingly, the specific heat maximum is located close to the bulk transition temperature as found for n = 16, and The development of the helical arrangement of magne- tization along the film growth direction was investigated by looking at the integral of the structure factor S ( ⃗ q ) along the z -direction, i.e. by taking ⃗ q = (0 , 0 , q$_{z}$ ), and making again use of the cumulant analysis in order to locate the helical transition temperature at T$_{N}$ (16) =