text:
Diego Magaldi

Maria Ulfa

Mai-Phuong Nghiêm

Gjergji Sini

Fabrice Goubard

Thierry Pauporté
email: thierry.pauporte@chimie-paristech.fr

Thanh-Tuân Bui
email: tbui@u-cergy.fr

Keywords: carbazole, hole transporting material, perovskite solar cell, polymer 1. Introduction

Designing organic hole transporting materials (HTMs) for stable perovskite photovoltaic devices remains a great challenge. Herein, we have prepared two small carbazole-based HTMs bearing a polymerizable double bond and one of the corresponding polymers. All compounds can be easily synthetized by a short procedure from largely available commercial products. These compounds have thermal, morphological, optical and electrochemical properties suitable for an application in perovskite solar cells. The side chains at the N-positon on the carbazole have a negligible influence on the opto-electrochemical properties.

However, the thermal properties differ largely between the monomer and its corresponding polymer.

While these materials have the same photovoltaic performance, the polymeric HTM led to more stable devices. This illustrates the higher thermal stability of polymeric HTM versus its corresponding monomer.

In recent years, the development of carbazole-based functional materials for emerging technological applications has received considerable interest. Carbazole offers very diverse chemistry enabling the design of different kind of functional materials having interesting optical and electrochemical properties, good charge transport, high thermal and morphological stability. [START_REF] Calió | Vacuum Deposited Perovskite Solar Cells Employing Dopant-Free Triazatruxene as the Hole Transport Material[END_REF][START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF][START_REF] Bui | Recent Advances in Small Molecular, Non-Polymeric Organic Hole Transporting Materials for Solid-State DSSC[END_REF][START_REF] Bui | Carbazol-N-Yl and Diphenylamino End-Capped Triphenylamine-Based Molecular Glasses: Synthesis, Thermal, and Optical Properties[END_REF][START_REF] Xue | Readily Synthesized Dopant-Free Hole Transport Materials with Phenol Core for Stabilized Mixed Perovskite Solar Cells[END_REF][START_REF] Al Mousawi | Carbazole Derivatives with Thermally Activated Delayed Fluorescence Property as Photoinitiators/Photoredox Catalysts for LED 3D Printing Technology[END_REF][START_REF] Al Mousawi | Azahelicenes as Visible Light Photoinitiators for Cationic and Radical Polymerization: Preparation of Photoluminescent Polymers and Use in High Performance LED Projector 3D Printing Resins[END_REF][START_REF] Bui | Thermally Activated Delayed Fluorescence Emitters for Deep Blue Organic Light Emitting Diodes: A Review of Recent Advances[END_REF][START_REF] Bui | Recent Advances on Organic Blue Thermally Activated Delayed Fluorescence (TADF) Emitters for Organic Light-Emitting Diodes (OLEDs)[END_REF] A huge number of molecular and polymeric carbazole derivatives have been developed for photovoltaics, especially for the emerging perovskite solar cells (PSCs). [START_REF] Sathiyan | Review of Carbazole Based Conjugated Molecules for Highly Efficient Organic Solar Cell Application[END_REF] PSCs have become the most promising renewable energy technology due to their high light-to-electricity conversion efficiency, their solution processability, low cost and possible large scale production. [START_REF] Di Giacomo | Flexible Perovskite Photovoltaic Modules and Solar Cells Based on Atomic Layer Deposited Compact Layers and UV-Irradiated TiO 2 Scaffolds on Plastic Substrates[END_REF] In PSC technology, the hole-transporting materials (HTM) and electron-transporting materials (ETM) implemented for selective contacts and serving for electron and hole extraction are important. [START_REF] Kim | Control of I-V Hysteresis in CH3NH3PbI3 Perovskite Solar Cell[END_REF] HTM and ETM having high charge transport properties with accurate frontier energy levels (EHOMO and ELUMO) compatible with perovskite valence band maximum (VBM) and conduction band minimum (CBM) bands are key for achieving high efficiency. However, the accurate design of efficient HTMs and ETMs for PSC devices, remains a great challenge. An ideal HTM must fulfil some general requirements to work efficiently in PSCs. Its EHOMO must be compatible with the VBM energy of the perovskite in order to provide the suitable driving force for charge transfer. Furthermore, its ELUMO should be significantly higher than that of the CBM of the perovskite absorbers to exhibit electron-blocking properties. Sufficient hole mobility, excellent thermal, morphological and photochemical stability, as well as a high hydrophobic property to protect the perovskite materials from moisture are also important key properties. Organic HTMs can be essentially divided into two subclasses depending on their molecular weights: molecular and polymeric HTMs. [START_REF] Ulfa | Molecular versus Polymeric Hole Transporting Materials for Perovskite Solar Cell Application[END_REF] While molecular HTMs feature well-defined molecular weight and structure with high synthetic reproducibility and easy purification, the polymer ones offer important advantages in terms of processability, thermal and mechanical stability, and higher intrinsic hole mobility. [START_REF] Hsu | Solid-State Dye-Sensitized Solar Cells Based on Spirofluorene (Spiro-OMeTAD) and Arylamines as Hole Transporting Materials[END_REF][START_REF] Holliday | Advances in Charge Carrier Mobilities of Semiconducting Polymers Used in Organic Transistors[END_REF][START_REF] Arias | Materials and Applications for Large Area Electronics: Solution-Based Approaches[END_REF][START_REF] Zaumseil | Electron and Ambipolar Transport in Organic Field-Effect Transistors[END_REF] A huge number of organic HTM have been reported in the literatures. [START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF][START_REF] Bui | Carbazol-N-Yl and Diphenylamino End-Capped Triphenylamine-Based Molecular Glasses: Synthesis, Thermal, and Optical Properties[END_REF][START_REF] Xue | Readily Synthesized Dopant-Free Hole Transport Materials with Phenol Core for Stabilized Mixed Perovskite Solar Cells[END_REF][START_REF] Bui | p-Methoxyphenyl)Amine End-Capped Tri(p-Thiophenylphenyl)Amine Based Molecular Glasses as Hole Transporting Materials for Solid-State Dye-Sensitized Solar Cells[END_REF][START_REF] Bui | Simple 3,6-Bis(Diphenylaminyl)Carbazole Molecular Glasses as Hole Transporting Materials for Hybrid Perovskite Solar Cells[END_REF][START_REF] Cho | Role of LiTFSI in High Tg Triphenylamine-Based Hole Transporting Material in Perovskite Solar Cell[END_REF][START_REF] Daskeviciene | Carbazole-Based Enamine: Low-Cost and Efficient Hole Transporting Material for Perovskite Solar Cells[END_REF][START_REF] Rakstys | Molecular Engineering of Face-on Oriented Dopant-Free Hole Transporting Material for Perovskite Solar Cells with 19% PCE[END_REF][START_REF] Zhang | Two Dimethoxyphenylamine-Substituted Carbazole Derivatives as Hole-Transporting Materials for Efficient Inorganic-Organic Hybrid Perovskite Solar Cells[END_REF] Structurally, most of them contain electronrich tertiary aromatic amine function. [START_REF] Zhao | Organic Hole-Transporting Materials for Efficient Perovskite Solar Cells[END_REF] Among the different HTM families, carbazole derivatives occupy an important place, especially the materials containing the 3,6-bis(di(4-methoxyphenyl)aminyl)carbazole building block. [START_REF] Berton | DMPA-Containing Carbazole-Based Hole Transporting Materials for Perovskite Solar Cells: Recent Advances and Perspectives[END_REF] The methoxy group effectively tune the solubility as well as the optoelectrochemical properties of targeted molecules and it provides a good interaction with perovskite materials' surface, notably with Pb(II) ion. [START_REF] Torres | Surface Effects and Adsorption of Methoxy Anchors on Hybrid Lead Iodide Perovskites: Insights for Spiro-MeOTAD Attachment[END_REF] Initially we employed this building block to develop HTM for solid state dyesensitized solar cells 2 prior to successfully applied them in PSCs. [START_REF] Bui | Simple 3,6-Bis(Diphenylaminyl)Carbazole Molecular Glasses as Hole Transporting Materials for Hybrid Perovskite Solar Cells[END_REF] Recently, we employed this building block to design a dendritic molecular HTM which took advantages of both molecular and polymeric materials. [START_REF] Bui | Design of Dendritic Core Carbazole-Based Hole Transporting Materials for Efficient and Stable Hybrid Perovskite Solar Cells[END_REF] Comparative study of this molecule with its simple model clearly showed that this HTM designing strategy could give PSC device with higher stability. [START_REF] Bui | Simple 3,6-Bis(Diphenylaminyl)Carbazole Molecular Glasses as Hole Transporting Materials for Hybrid Perovskite Solar Cells[END_REF][START_REF] Bui | Design of Dendritic Core Carbazole-Based Hole Transporting Materials for Efficient and Stable Hybrid Perovskite Solar Cells[END_REF][START_REF] Ulfa | Impact of Organic Hole Transporting Material and Doping on the Electrical Response of Perovskite Solar Cells[END_REF] In continuity with this work, in the present study, the same building block was used to design novel polymerizable molecular HTMs (DM1 and DM2, Scheme 1). Scheme 1. Structure of studied HTMs DM1 could be successfully transformed into its corresponding polymer (DM1P, Scheme 1). We thus made a comparative study of the molecular HTMs bearing a polymerizable double bond function and its corresponding polymeric counterparts and evaluated the devices performance and operation stability.

Results and discussion

Materials synthesis

DM1 and DM2 were synthetized using the synthetic route given in Scheme 2 and already employed in our previous works. [START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF][START_REF] Abdallah | Development of New High-Performance Visible Light Photoinitiators Based on Carbazole Scaffold and Their Applications in 3D Printing and Photocomposite Synthesis[END_REF] Briefly, 3,6-dibromocarbazole reacts with alkyl halide or benzyl halide in basic medium leading to N-substituted carbazoles (1-2). 1 and 2 are then subjected to a twofold Pd-catalyzed amination with 4,4'-dimethoxydiphenylamine, leading to the final products (DM1-2) in good yields. The detailed description and the analytical analysis of DM1 and DM2 has been recently reported in a separate article focused on the photopolymerization application. In the next step, we attempted to convert DM1 and DM2 into their corresponding polymers. Different classic alkene polymerization conditions were implemented such as radical (AIBN and BPO) and anionic (tert-butyl lithium, NaOH) ones in different solvents (toluene, cyclohexanone, DMF) but they only led to traces of polymers. Finally, we found that during the synthesis of DM1, if chloroform is added into the reaction medium in the amination step, a polymer (called DM1P) is spontaneously formed. DM1P was also formed when t BuONa was added into a solution of purified DM1 dissolved in a mixture of CHCl3/toluene at room temperature and let to react under stirring for 4 hours. The solvents were then removed to the saturation point and the product was added into MeOH to precipitate the polymer. In both cases, yellow amorphous polymers were obtained in 60-80% of yield. We believe that DM1P is formed by a chain transfer reaction catalyzed by the chloroform and the strong basic conditions. [START_REF] Konishi | Base-Catalyzed Polymerization of p-Styrenesulfonamides. II. N-Methyl and N-Phenyl Derivatives[END_REF][START_REF] Jenkins | Terminology for Reversible-Deactivation Radical Polymerization Previously Called "Controlled" Radical or "Living" Radical Polymerization (IUPAC Recommendations[END_REF][START_REF] Mayo | Chain Transfer in the Polymerization of Styrene: The Reaction of Solvents with Free Radicals 1[END_REF][START_REF] Flory | The Mechanism of Vinyl Polymerizations[END_REF][START_REF] Furuncuoǧlu | Role of Chain Transfer Agents in Free Radical Polymerization Kinetics[END_REF][START_REF] Matyjaszewski | Handbook of Radical Polymerization[END_REF] GPC analysis showed that DM1P had an average molecular weight of ca. 43 kDa (PDI = 2.19). All of these mentioned polymerization methods have been also applied to DM2. However, no polymer formation was observed which could be explained by the lack of the stability of the intermediate radicals due to the absence of a conjugated aromatic ring like the case of DM1. All materials (DM1, DM2, DM1P) were well-soluble in common organic solvent, rending them suitable for solution processing.

Structural, thermal, and morphological properties

The DM1P polymer was studied by density functional theory (DFT) calculations using the Gaussian 09 [START_REF] Frisch | [END_REF] suite of programs. Ground state geometries of DM1P series were fully optimized starting from two different initial arrangements between the median axes of successive pending groups, 120° and 180° arrangements. The wB97WD [START_REF] Chai | Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections[END_REF] functional was used for these calculations in conjunction with the 6-31G(d,p) basis set. Thermogravimetric analysis (TGA) and differential scanning calorimetry (DSC) experiments were conducted to examine the thermal properties of the novel carbazole derivatives. Figure 2 shows the TGA curves of the different compounds while all DSC curves are given in the ESI (Figures S 1, 2, 3). Thermal data are gathered in Table 1. All these compounds were found thermally stable with a thermal decomposition temperatures (Td) higher than 400°C which are high enough for practical applications. The surface roughness was estimated at 12.9 nm for DM1, 6.1 nm for DM2, while the polymer DM1P

HTM layer was the smoothest at 2.3 nm. These smooth surfaces are expected to provide perfect HTM/perovskite contact and thus to facilitate the hole transfer from perovskite to HTM and protect the perovskite absorber layer against moisture.

Optical and electrochemical properties

The optical properties of the three materials were examined in dilute solution and in the thin solid film state. The spectra are disclosed in Figure 4 and the related data are gathered in Table 1. In dilute solution, all materials have identical absorption behaviour with a strong absorption band located in the UV spectral domain centred at ca. 300-320 nm. Another weaker absorption band is observed at longer wavelength, in the UV-near visible region (350-450 nm). Obviously, all compounds, including the polymer DM1P, have the same optical behaviour as the model compound without double bonds in the side chain at the N-position of the carbazone, which is typical for bisdiphenylaminylcarbazole derivatives. [START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF] In our previous study, DFT calculations showed that the side chain at the N-position of the carbazole ring had no effect on the optical properties of these derivatives. [START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF] The polymerization did not affect the optical properties of the target compound as it created a saturated polystyrene polymer backbone. The spectra recorded for thin films are identical to that recorded in dilute solution. It indicates no aggregation in the solid state. This could be attributed to the tridimentional bulky structures of 3,6bisdiarylamine substituted carbazole derivatives. The optical bandgap was determined from the onset of the absorption spectrum in the solid state and the data is given in Table 1. These values are somewhat smaller than that of the Spiro-OMeTAD (Eg = 3.0 eV). [START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF] The redox behavior of carbazole materials were examined by cylic voltammetry. In combination with optical bandgap, the redox potentials were used to estimate their EHOMO and ELUMO. EHOMO were estimated from the oxidation potential by using the ferrocene (Fc) ionization potential value as reference (4.8 eV vs. vacuum [START_REF] Pommerehne | Efficient Two Layer Leds on a Polymer Blend Basis[END_REF][START_REF] Cardona | Electrochemical Considerations for Determining Absolute Frontier Orbital Energy Levels of Conjugated Polymers for Solar Cell Applications[END_REF] ). As expected, the nature of the side chain at the N-position of the carbazole and the polymerization did not significantly influence the electrochemical redox properties of final materials.

In Figure 5, the cyclic voltammogram of all materials are identical, showing multiple oxidation processes while no reduction process was observed in the cathodic scan. The first reversible oxidation potential could be attributed to the formation of radical cation of the carbazole moiety. DM1 exhibits an oxidation potential onset of -0.20 V/Fc/Fc + , which is nearly identical to that of the DM2 (-0.16 V) and DM1P (-0.15 V). These oxidation potentials are very close to that reported for similar molecules without polymeriable double and that of Spiro-OMeTAD measured under the same conditions. [START_REF] Bui | Carbazole-Based Molecular Glasses as Hole-Transporting Materials in Solid State Dye-Sensitized Solar Cells[END_REF] This confirms the trend that the side chains at N-position of carbazole have a negligible effect on the electrochemical properties of the targeted compounds. EHOMO was estimated from the first oxidation potential onset while ELUMO was estimated from the EHOMO and the optical bandgap. The EHOMO values of DM1, DM2 and DM1P are -4.60, -4.64, and -4.65 eV respectively. From the EHOMO and the optical bandgap, corresponding ELUMO were estimated and reported in Table 1.

These EHOMO and ELUMO suggest that all compound could be suitable for use as HTM in the PSC, ensuring good hole transfer and the electron blocking at the interface. PSCs employing these HTMs have been fabricated and their photovoltaic characteristics are discussed in the next section.

Photovoltaic properties

The new carbazole compounds were tested for an application as HTM in PSC with methylammonium lead iodide (CH3NH3PbI3) as the hybrid perovskite compound. [START_REF] Chai | Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections[END_REF] We used the same optimized device structure as we reported earlier with the optimized thickness of the HTM layer of ca. 200 nm in comparison to 150 nm of Spiro-OMeTAD. [START_REF] Ulfa | Impact of Organic Hole Transporting Material and Doping on the Electrical Response of Perovskite Solar Cells[END_REF] The general device structure used in this work and the energy level alignment of different device components are shown in Figure 6. The two-step perovskite deposition method was employed as fully described in our previous paper. [START_REF] Wang | Effects of Perovskite Monovalent Cation Composition on the High and Low Frequency Impedance Response of Efficient Solar Cells[END_REF] The photovoltaic characteristics of the best cells together with average values for each HTM are gathered in Table 2. They are compared to benchmark Spiro-OMeTAD cells produced during the same batches.

The J-V curves of the best devices are given in Figure 7. The maximum efficiency was measured for the DM1 HTM at 12.2% and it was approaching 12% for DM2. DM1P cells had a slightly higher PCE with a maximum at 13.24%. However, the averaged PCE values of all HTM were nearly identical with the Voc of all devices being close. The best result with DM1 compared to DM2 is due to the higher FF parameter. It is also noteworthy that the DM2 cells exhibited the highest Jsc, with values close to that of the Spiro-OMeTAD cells (Table 2). The cells with the polymer HTM had slightly ameliorated performances. Without dopant, the efficiency of DM1P cells was very low (below 0.7%); their J-V curves were typical of highly resistive devices and problems of charge injection and transport occurred. For DM1P HTM, without doping, the device suffered from a very low FF that can be attributed to the low HTM conductivity and therefore high charge transport resistances. For all the carbazole HTMs, the Voc of the cells were lower than the Spiro-OMeTAD ones. It suggests the occurrence of more recombinations at the perovskite/HTM interface. Interestingly, the obtaining of the maximum PCE with the molecular carbazole HTMs (DM1 and DM2) required a rather long light soaking. The first measurements gave low PCEs. However, with time, this parameter continuously increased before to stabilize as illustrated in Figure 8. The photovoltaic parameters provided in Table 2 are the stabilized ones. This light soaking effect was not observed for the polymer DM1P. We explain it by the low Tg of pristine molecular HTMs (Tg bellow 75°C). The presence of additives in the HTM (tert butylpyrridine and LiTFSI) even lower these Tg. [START_REF] Zhao | Effect of Selective Contacts on the Thermal Stability of Perovskite Solar Cells[END_REF] Under continuous light irradiation, the cells were heated. It conducted to morphological change of the HTL, to an improvement of the HTM/Perovskite interface and thus ameliorated the photovoltaic efficiency. In contrast, DM1P has a much higher Tg (186°C) and is morphologically stable at the operating temperature. Its amorphous nature ensures a good HTL/perovskite interface during the continuous light irradiation.

Conclusion

We have reported three novel carbazole-based molecular and polymeric materials. These compounds have thermal, morphological, optical and electrochemical properties suitable for an application in perovskite based photovoltaic devices. While their optoelectrochemical properties are nearly identical, the thermal properties differ in a large extent between the monomer and its corresponding polymer. These materials have the same photovoltaic performances, the polymeric HTM has led to more stable devices.

This illustrates the higher thermal stability of polymeric HTM versus its corresponding monomer.

29

 29

Scheme 2 .

 2 Scheme 2. Synthesis of DM1 and DM2

Figure 1 .

 1 Figure 1. DM1P 120° "triangular" conformation (left) and 180° conformation (right)

Figure 2 .Figure 3 .

 23 Figure 2. TGA curves of DM1, DM2 and DM1P

Figure 4 .

 4 Figure 4. UV/Vis absorption spectra of DM1, DM2 and DM1P in dilute CH2Cl2 solutions and in thin film.

Figure 5 .

 5 Figure 5. CV of DM1, DM2 and DM1P in 0.01M TBAPF6/CH2Cl2 solution and scan rate was 100 mV s -1 .

Figure 6 .

 6 Figure 6. Schematic presentation of PSC device structure (left) and the energy alignment of its different components (right)

Figure 7 .

 7 Figure 7. J-V reverse scan curves of the best PSC devices reported in this work

Figure 8 .

 8 Figure 8. Effect of light soaking of DM1 and DM2 PSCs

Table 1 .

 1 Thermal, optical and electrochemical properties of the investigated HTM's DM1, DM2 and DM1P.

	Compound	Td [a]	Tg [b]	λmax [c]	λonset [c]	Egap [d]	Eox onset [e]	EHOMO [f]	ELUMO [g]
		°C	°C		nm	nm	eV	V/Fc/Fc +	eV	eV
	DM1	415	61	305, 375	443	2.80	-0.20	-4.60	-1.80
	DM2	414	50	307, 372	444	2.80	-0.16	-4.64	-1.84
	DM1P	405	186	309, 382	454	2.73	-0.15	-4.65	-1.92
	[a]Obtained from TGA measurements. [b]Obtained from DSC measurements. [c] Measured in CH2Cl2 solution at
	concentration of 1 mM. [d] E g opt =hc/λonset =1240/λonset (eV). [e] in CH 2 Cl 2 solution containing TBAPF 6 as supporting
	electrolyte (0.1 M), scan rate 100 mV.s -1	. [f] EHOMO = -(4.8 + Eox onset ) (eV). [g] ELUMO = EHOMO + Egap (eV)

Table 2 .

 2 Photovoltaic parameters of the best device of each HTM. The values in parenthesis are average values.

	HTM	Scan direction	Voc (V)	Jsc (mA.cm -2 )	FF (%)	η(%)
	DM1	Reverse Forward	0.99 (0.97) 0.99 (0.97)	17.99 (20.69) 18.6 (20.92)	68.4 (58.18) 59.65 (49.89)	12.26 (11.63) 11.08 (10.07)
	DM2	Reverse Forward	0.97 (0.97) 0.97 (0.96)	21.31 (21.42) 21.02 (21.16)	56.18 (54.12) 48.00 (43.98)	11.73 (11.32) 9.84 (8.97)
	DM1P	Reverse Forward	0.92 (0.91) 0.95 (0.93)	18.98 (15.95) 16.94 (17.98)	75.53 (76.04) 48.01 (42.62)	13.24 (11.09) 7.70 (7.18)
	Spiro-OMeTAD	Reverse Forward	1.01 (1.01) 1.00 (1.01)	22.15 (21.32) 22.72 (22.14)	78.75 (75.63) 63.02 (62.81)	17.68 (16.33) 14.45 (14.05)

Acknowledgements

M.U. and D.M. acknowledge the Indonesia Endowment Fund for Education (LPDP) and the Mexican government for their PhD scholarship, respectively.


V Russier
email: russier@icmpe.cnrs.fr

Juan J Alonso
email: jjalonso@uma.es

I Lisiecki

A T Ngo

C Salzemann

S Nakamae

C Raepsaet

Phase diagram of a three-dimensional dipolar model on an fcc lattice

The magnetic phase diagram at zero external field of an ensemble of dipoles with uniaxial anisotropy on a FCC lattice is investigated from tempered Monte Carlo simulations. The uniaxial anisotropy is characterized by a random distribution of easy axes and its magnitude λ u is the driving force of disorder and consequently frustration. The phase diagram, separating the paramagnetic, ferromagnetic, quasi long range ordered ferromagnetic and spin-glass regions is thus considered in the temperature, λ u plane. This system is aimed at modeling the magnetic phase diagram of supracrystals of magnetic nanoparticles.

I. INTRODUCTION

Assembly of magnetic nanoparticles (MNP) in dense packed structures receives a growing interest as a way to build new composite materials in a bottom-up strategy. Besides promising applications, for instance in nanomedicine The modeling of the magnetic properties of dense packed structures of MNP is in principle a multi scale problem but is fortunately highly simplified for particles whose size falls in the single domain range if the non collinearity of the surface spins is neglected (the so-called spin canting). Then the effective one spin model (EOSM) can be used although this level of approximation has been questioned and extended for situations where the core-surface MNP morphology is expected to play a significant role 18,19 . In the framework of the EOSM and considering the MNP as frozen in position in the system one is led to model an ensemble of MNP free of super exchange interaction owing to the non magnetic coating layer and thus interacting only through the dipole dipole interaction (DDI) between the moments of the macro spins and undergoing the magneto-crystalline anisotropy energy (MAE). It is well known that strongly coupled dipolar systems present a low temperature phase dependent on the underlying structure and the amount of disorder [20][21][22][23][24] which originates either from the structure, the dilution or the MAE through the distribution of easy axes. Hence, the so-called super ferromagnetic (SFM) or super spin-glass (SSG) phases can be observed. Experimental evidence of the SSG state has been given either in randomly distributed MNP or on well ordered MNP supra-crystals with randomly distributed anisotropy easy axes. In the case of well ordered supra-crystals for which a FCC lattice is a representative situation, the low temperature phase results from the competition between the dipolar induced ferromagnetic phase of the purely dipolar system (free of MAE), and the frustration introduced by the disorder. In the case of uniaxial MAE and for a wholly occupied lattice, the latter comes only from either the amplitude of the MAE when the easy axes distribution is random, or the degree of alignment of the easy axes in the strong MAE limit. The former situation corresponds to a supra-crystal synthesized in the absence of external field. This strong MAE limit has been investigated recently in the random close packed and for MNP on a perfect FCC lattice 25,26 .

The influence of the random anisotropy such as introduced by the uniaxial MAE with random distribution of easy axes has already been widely studied in the framework of the random anisotropy model (RAM), first introduced by Harris et al. 27 , where the spin-spin interaction is quite generally the short-range Heisenberg or Ising exchange one, the latter being the strong anisotropy limit of the former. The relevant disorder control parameter is the anisotropy to exchange ratio (D/J) which governs the crossover from the FM state (D/J = 0) to the spin-glass state obtained in the D/J → ∞ limit 28 , where D and J are the anisotropy and exchange coupling constants respectively. When D/J increases the long-range (LR) FM state first disorders to a quasi long-range order (QLRO) state with a ferromagnetic character 29 . The Monte Carlo simulations of Nguyen and Hsiao 30,31 show clearly, on the basis of the dynamic behavior, that the PM/FM like magnetic phase transition at weak anisotropy disappears with the increase of D/J and a glassy phase transition takes place beyond a threshold value. This scenario is also in agreement with the simulation results for the spontaneous magnetization in terms of (D/J) 32 .

In the present work following related contributions on the dipolar Ising model 25,26 devoted to the effect of the MAE easy axes texturation in the infinite MAE limit, we investigate the magnetic phase diagram from Monte Carlo simulations of FCC supra-crystals of MNP interacting via DDI and undergoing a uniaxial MAE whose easy axes are randomly distributed. Here the disorder control parameter is the MAE to the DDI strengths ratio as is the case in the RAM. An important difference with the dipolar Ising model studied in Refs. 25,26 is the 3 dimensional nature of the model variables, as is the case in the Heisenberg model with important consequences essentially in the spin-glass phase and with the occurrence of a transverse spin-glass state associated with a QLRO in the FM region of the phase diagram. The choice of the FCC structure is mainly dictated by the MNP organizations obtained experimentally as outlined above. To investigate the phase diagram we calculate the ferromagnetic and overlap spin-glass order parameters, from which a finite size analysis of the corresponding Binder cumulants, spin-glass and transverse spinglass correlation lengths is performed. We also focus on the finite size behavior of the heat capacity. The paper is organized as follows. We introduce the model in section II, the simulation details and the observables we use are then presented in section II A and II B respectively. Section III is devoted to the analysis of the results and a conclusion is given in section IV.

II. MODEL

We model an assembly of MNP free of super exchange interactions, characterized by a uniaxial magneto-crystalline anisotropy (MAE) and self organized on a supra-crystal of face centered cubic (FCC) structure. To this aim we place ourselves in the framework of the effective one spin model where each single domain MNP is assumed to be uniformly magnetized with a temperature independent saturation magnetization M s . Hence we consider a system of N dipolar hard spheres of moment ⃗ µ i = µμ i located on the sites of a FCC lattice occupying a total volume V , interacting through the usual dipole dipole interaction (DDI) and subjected to a one-body anisotropy energy,

-K i v i (n i . μi ) 2 .
μi and ni are three dimensional unit vectors and in the following, hatted letters denote unit vectors. K i , ni and v i are the anisotropy constant, the easy axis and the volume of the particle i respectively. We consider a random distribution of easy axes on the unit sphere, namely the azimuthal angles are randomly chosen while the polar angle distribution follow the probability density sin(Θ)/2. The MNP ensemble is monodisperse with MNP diameter d. The particle moment is related to the material saturation magnetization

M s through µ = v(d)M s (v(d) = πd 3 /6
). The Hamiltonian of the system is given by

βH = 1 2 β d i≠j μi . μj -3(μ i . rij )(μ j . rij ) (r ij /d) 3 -βKv(d) i (n i . μi ) 2 with d = µ 0 4π µ 2 d 3 (1)
where rij is the unit vector carried by the vector joining sites i and j, r ij its length and β = 1/(k B T ) is the inverse temperature. The Hamiltonian (1) is the same as that used in 24 while in 25 , where the infinite anisotropy limit was considered, the second term of equation ( 1) was not included and the directions of the moments, μi were set equal to the easy axes ni representing then the Ising axes. Concerning the reduced temperature, instead of the natural choice T k B / d , we take advantage of the 1/r 3 dependence of the DDI to introduce the more convenient reduced temperature

T * = T k B /( d (d/d r )
3 ), with d r = d(Φ r /Φ) 1/3 where Φ and Φ r are respectively the volume fraction, (N /V )πd 3 /6, and a reference value (here, the maximum value for hard spheres on a FCC lattice). This choice of T * equivalent to measuring the temperature in terms of a dipole dipole energy weighted by the volume fraction instead of its maximum value at contact 26,[START_REF]Introducing d r = d(Φ r /Φ) involving the volume fraction Φ of the system ((N /V )πd 3 /6) and a reference value, Φ r (for instance the maximum value for hard spheres on a FCC lattice)[END_REF] . Equation ( 1) is then rewritten as

βH = 1 T * 1 2 i≠j μi . μj -3(μ i . rij )(μ j . rij ) (r ij /d r ) 3 -λ u i (n i . μi ) 2 λ u = Kv(d) d (d r /d) 3 (2a) ≡ 1 T * 1 2 i≠j μi Tij μj -λ u i (n i . μi ) 2 (2b)
which introduces the MAE coupling constant λ u . The simulation box is a cube with edge length L s = 2Ld r and the total number of dipoles is N = 4L 3 . (In the following without loss of generality we consider the case of a FCC lattice with d r = d). We consider periodic boundary conditions by repeating the simulation cubic box identically in the 3 dimensions. The long range DDI interaction is treated through the Ewald summation technique 21,[START_REF] Allen | Computer simulation of liquids[END_REF] , with a cut-off k c = 10 k m , k m = (2π/L s ), in the sum of reciprocal space and the α parameter of the direct sum chosen is α = 5.80, a value which permits to limit the sum in direct space to the first image term 21,[START_REF] Kusalik | [END_REF] . The Ewald sums are performed with the so-called conductive external conditions 21,[START_REF] Allen | Computer simulation of liquids[END_REF] , i.e. the system is embedded in a medium with infinite permeability, µ s = ∞, which is a way to avoid the demagnetizing effect and thus to simulate the intrinsic bulk material properties regardless of the external surface and system shape effects. For the different values of λ u considered in the following, the simulations are performed for system sizes up to either L = 7 (N = 1372) or L = 8 (N = 2048).

A. Simulation method

In order to thermalize in an efficient way our system presenting strongly frustrated states, we use parallel tempering algorithm 36,37 (also called tempered Monte Carlo) for our Monte Carlo simulations. Such a scheme is widely used in similar systems 25,26,38 , and we do not provide the details. The method is based on the simultaneous simulation runs of identical replica for a set of temperatures {T * n } with exchange trials of the configurations pertaining to different temperatures each N M Metropolis steps according to an exchange rule satisfying the detailed balance condition. The set of temperatures is chosen in such a way that it brackets the transition temperature while ensuring a satisfying rate of exchange between adjacent temperature configurations. Our set {T * n } is either an arithmetic distribution or an optimized one in order to make the exchange rate between adjacent paths as constant as possible in the whole range of {T * n } according to the efficient constant entropy increase method 39 . In the present work we take N M = 10, the number of temperatures is in between 36 and 48 according to the value of L and the amplitude of temperatures in the set {T * n } for L ≤ 7 and up to 96 for L = 8. When necessary, precise interpolation for temperatures between the points actually simulated are done through reweighting methods 40 . We use t 0 Monte Carlo steps (MCS) for the thermalization and the averaging is performed over the (t 0 , 2t 0 ) interval following MCS with t 0 = 5 10 5 .

We deal with frozen disorder situations where each realization of the easy axes distribution {n i } defines a sample. Accordingly, a double averaging process is performed first relative to the thermal activation, i.e. the Monte Carlo step, and second on the whole set of N s samples. Consequently, the mean value of an observable A, results from a double averaging denoted in the following as [< A >] where < . > corresponds to the thermal average on the MC sampling for a fixed realization of the axes distribution and [.] to the average over the set of samples considered. The number of samples necessary to get an accurate result depends strongly on the value of λ u . Obviously, for λ u = 0, N s = 1 should be sufficient for a very long MC run in order to get a satisfying average. In practice, we use N s = 10 for λ u = 0. We find that N s of the order of 100 to 200 is sufficient up to λ u ≤ 1.5 whereas accurate results for λ u = 4 require the use of at least N s = 400 realizations. However, it is worth mentioning that accurate results for the heat capacity C v are much less demanding than for the overlap order parameter and about 50 to 100 realizations allow to get very well averaged C v up to λ u = 100. In any case, in the present work we limit ourselves to λ u = 4 in the computations of the spin-glass overlap order parameter q 2 and related properties. The error bars of the averaged quantities are deduced from the mean squared deviations of the sample to sample fluctuations. The simulation code is massively parallelized, all temperatures and the 2 replicas running together, and the typical CPU time is c.a. 40h for the complete run of one sample of L = 7 (N = 1372) on the Intel Xeon E5 processors at the CINES center.

B. Observables

Our main purpose is the determination of the transition temperature between the paramagnetic and the ordered phase and on the nature of the latter, namely ferromagnetic or spin-glass, in terms of λ u . For the PM/FM transition, we consider the spontaneous magnetization

m = 1 N i μi (3)
computing its moments, m n = [< m n >], n = 1,2 and 4. We compute also the nematic order parameter P 2 together with the instantaneous nematic direction, d which are the largest eigenvalue and the corresponding eigenvector respectively of the tensor Q 21 , where Ī is the unit tensor of cartesian components I αβ = δ αβ . The spontaneous magnetization can also be studied in the ordered phase from the mean value projected total magnetization on the nematic direction, which defines

= 1 N ∑ i (3μ i μi -Ī)/2
m d = 1 N i μi . d (4)
We compute the mean value

m d1 = [< |m d | >] and the moments m dn = [< m n d >],
with n = 2, 4. To locate the transition temperature, T * c , as usually done, we will use the finite size scaling (FSS) analysis of the Binder cumulant 41 which is defined either from the moments m n or m dn characterized by 3 or one degree of freedom respectively

B m = 1 2 5 -3 m 4 m 2 2 , B md = 1 2 3 - m d4 m 2 d2 (5)
From these normalizations, B m , B md → 1 in the long range FM phase and B m , B md → 0 in the limit L → ∞ in the disordered PM phase. For the PM/SG transition, we consider the usual overlap order parameter 42,43

q 2 = αβ |q αβ | 2 , with q αβ = 1 N i µ (1) iα µ (2) iβ , (6)
where the superscripts (1) and (2) stand for two independent replicas of an identical sample and α, β denote the cartesian coordinates. From q 2 we calculate the mean value q 2 = [< |q 2 | >] and q 4 = [< q 4 >] and the corresponding Binder cumulant,

B q = 1 2 (11 -9 q 4 q 2 2 ) (7)
In order to determine the PM/SG transition temperature it may be more convenient to use the spin-glass correlation length, usually defined from a fit of the spin-glass overlap parameter correlation function on its Ornstein-Zernike form, assuming that it follows from a so called Φ 4 theory 44,45

ξ = 1 2sin(k m /2) q 2 (0)) q 2 ( ⃗ k m ) -1 1/2 with q 2 ( ⃗ k) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⟨ α,β q αβ ( ⃗ k)q * αβ ( ⃗ k)⟩ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ q αβ ( ⃗ k) = 1 N i µ (1)
iα µ

(2)

iβ e i ⃗ k ⃗ r i , (8)
with ⃗ k m = (k m , 0, 0). By exploiting the non dimensional character of the ratio ξ/L one is led to deduce T * c as the crossing point of the ξ/L curves in terms of T * for different values of L. Finally, the magnetic susceptibility, χ m and the heat capacity , C v are calculated from the magnetization and the energy fluctuations respectively

χ m = N T * < m 2 > -< m > 2 , C v = 1 N T * 2 < H 2 > -< H > 2 (9)
III. RESULTS

The phase diagram, which summarizes the present work, is displayed in figure (1). It separates the paramagnetic (PM), ferromagnetic (FM) and the spin-glass (SG) phases. The important feature of T * c (λ u ) along both the PM/FM and the PM/SG lines is its very weak dependency with λ u , the overall variation being limited to 0.6 ≤ T * c ≤ 1.0. Before going in the details of its determination, a qualitative overview of the evolution with respect to the MAE coupling constant λ u combined with the finite size effect of both the ferromagnetic and nematic order parameters, namely the magnetization m and the eigenvalue P 2 is useful. This is displayed in figure (2). The rising of the nematic order parameter below a threshold temperature, close to the transition temperature T * c is indicative of the occurrence of a collective spontaneous direction in the system induced by the DDI which competes with the random anisotropy contribution. It is worth mentioning that the finite value of the order parameters in the PM phase, well above T * c , results from the the short range correlations in finite size systems 29 and is a quite general result 22,23,46,47 . The curves of figure (2) suggest clearly an evolution of the low temperature phase from a well ordered FM phase in figure (2 a) corresponding to the pure dipolar case, λ u = 0 where the PM/FM transition is well established 20,21,24,47 to a gradual crossover to the absence of FM order, figure (2 c) with the complete disappearance of the latter at λ u = 2. Between these two regimes, figure (2 b), the similarity of m and P 2 in terms of T * with the pure dipolar case remains only on a qualitative level since then both order parameters P 2 and m decrease with the system size L conversely to the case at λ u = 0. The behavior of both P 2 and m at λ u = 0.2 is characterized, as is the case at λ u = 0 by the merging at low T * of the curves corresponding to different values of L, indicating a finite value of P 2 and m in the thermodynamic limit (L → ∞) and we find the onset of the decrease with L of the ferromagnetic order parameters at λ u ≃ 0.4, as is clarified in figure (3c) where we compare the behavior of m 2 (T * ) with the system size between λ u = 0.4 and 0.8. For small values of λ u , starting from λ u = 0, we check the ferromagnetic nature of the low temperature phase and determine the transition temperature T * c by following the FSS and thus we look for the crossing point of the Binder cumulant, B m (L) curves for different values of L, consequence of the scaling behavior

B m = b m ((T * -T * c )L 1/ν ) . (10)
However, the system sizes considered in the present work are too small to determine accurate values of the exponent ν. At λ u = 0 we get T * c = 0.614 ± 0.004 in agreement with Refs. 20,24,47 . The spontaneous magnetization orientates along the <111> directions in agreement with results of the literature 20,47 and as Alkadour et al.. 47 we do not find any spin reorientation at low temperature. Moreover we have checked that equation ( 10) is very well reproduced using ν = 0.692 48 , the value of ν obtained for dipoles on lattice of cubic symmetry from renormalization group developments, When λ u increases in the low MAE regime, the behavior of m 2 with respect to the system size changes as we observe a decrease of the magnetization with the increase of L at low temperature (see figure ( 6). From the behavior of q 2t in terms of T * , given in figure (6 a) for increasing values of λ u , we see that no transverse spin-glass state is expected for λ u < 0.4. Indeed, q 2t not only remains very small on the whole range of T * but presents no noticeable rising up at low temperature and is moreover a decreasing function of L down to very low temperature up to λ u = 0.35 as shown in figure (6 b,c) where we also compare figure (6 d) to the case λ u = 0.8 well inside the FM-QLRO region of the phase diagram. We emphasize that this is related to the onset at λ u = 0.4 of the the L-dependence of m 2 (T * ) at low T * as mentioned above. Given the very small values taken by q 2t for λ u = 0.35, we cannot interpret the crossing point of the q 2t (T * ) curves for different values of L as the onset of the transverse spin-glass state. We emphasize that q 2t (T * ) when T * < 0.6 remains smaller or comparable to its value at the hump ( q 2t ≤ 0.04 ≪ 1). This is no longer the case for instance at λ u = 0.8, figure (6a,d), in the FM-QLRO region of the phase diagram. The corresponding transition temperature, T * xy below which the FM and transverse spin-glass orders coexist is determined by the crossing point of the ξ t /L curves where ξ t is the transverse spin-glass correlation length as obtained from equation ( 8) with the overlap order parameter replaced by its transverse equivalent. The result we get is that T * xy is very close to and seems slightly larger than the PM/FM transition temperature T * c , see figures (7) and (1) Increasing λ u , we still find a PM/FM transition up to λ u = 1.3 as indicated by the behavior of B m as shown in figure (4 c,d). Indeed, we still get a crossing point in temperature, allowing the determination of T * c . Moreover, the orientation of the magnetization in the FM phase remains very close to the <111> directions (see figure (8b)). Of course the distribution of the moments around the mean polarization direction depends largely on λ u and this can be easily visualized from an instantaneous moment configuration as is done in figure (8 a,b) and estimated from the variation of the nematic order parameter mean value, P 2 in terms of λ u at constant temperature. Using a simple model to represent the distribution of the moment around the nematic direction, as for instance

P (Θ) ∝ sin(Θ)[exp(-Θ 2 /(2σ 2 ))], 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 0.3 0.5 0.7 0.9 1.1 C v T
* λ u = 0 a) L = 4 5 6 7 8
0.3 0.5 0.7 0.9 we can get the variance σ of this distribution for a given value of P 2 . For instance at T * = 0.1, we get σ = 0.226, 0.458 and 0.656 for λ u = 0, 0.8 and 1.2 where P 2 = 0.859, 0.543 and 0.304 respectively. Then, at λ u = 1.2, the value taken by B m at low temperature is slightly smaller than 1, and the clear crossing observed for T * < T * c at lower values of λ u becomes merely a merging at least for sufficiently large values of L. At λ u = 1.3 the different B m (L) curves merge in a very limited range of temperature, and from λ u = 1.4 and beyond B m (L) is a decreasing function of L whatever the value of T * indicating the absence of PM/FM transition as displayed in figure (4 e,f). This is an indication of the change in the nature of the transition from the PM/FM to the PM/SG one or at least to a PM/glassy phase. This is corroborated by the reduction of the finite size dependence of the heat capacity, as can be seen on figure (5b), since no singular finite size behavior of C v is expected in a PM/SG transition 52 . This latter feature is more evident on the dipolar component of C v in the strong MAE coupling case (discussed in more detail in section III C and in figure (11)). We conclude that the PM/FM line extends from λ u = 0 to λ T * FIG. 6. Spin-glass transverse overlap order parameter defined in the plane normal to the nematic direction for a) λ u = 0.1 to 1.2 as indicated and L = 7. Inset: detail of the hump in the vicinity of the PM/FM transition temperature. q 2t for λ u = 0.1 and 0.2 are only slightly distinguishable at the scale of the inset. Finite size effect on q 2t , for L = 4 to 7 as indicated and b) λ u = 0.2, c) λ u = 0.35 showing the very beginning of the low temperature rising up in the latter case. d) Case λ u = 0.8, well inside the QLRO region of the phase diagram.

B. Spin-glass transition

Beyond λ u = 1.4 we are clearly in the PM/SG transition region. As is the case for the PM / chiral glass transition of the 3D Heisenberg spin-glass model the Binder cumulant (7) does not present a well defined crossing point as in the FM/PM transition. It is instead characterized by a dip in the B q negative region, which deepens with the increase of L. The location of this dip must converge towards T * c 43,53,54 . Our simulations confirm such a behavior for B q as can be seen in figure (9) where B q is displayed for λ u = 4. Therefore, for λ u > 1.5, we determine the value of T * c as the crossing point of the ξ/L curves for different values of L on the one hand and the interpolation to 1/L = 0 of the linear fit in 1/L of the value of T * at the dip in B q on the other hand. Because of the heaviness of the computations, we consider only λ u = 1.6, 2 and 4. The two above mentioned determinations of T * c are in quite good agreement, as can be seen in figure (10) corresponding to λ u = 4. Furthermore the value of T * c in the limiting case, λ u → ∞, which is the so-called random axes dipoles (RAD) model 55,56 is known from the literature, T * c (∞) = 0.95 ± 0.05 26 .

The FM/SG transition line below T * c is not easy to locate. We have nevertheless estimated λ u (T * c ) along the FM/SG line from the condition that B m (L) at T * constant is an increasing (decreasing) function of the system size L in the FM (SG) side of this line, with however a rather large uncertainty (see figure (1)).

C. Large λ u limit

In the strong λ u limit, the MAE contribution to the total Hamiltonian (1) tends to align the moments μi on the anisotropy axes ni making the model coincide with the dipolar Ising model 38,46,55,56 , the anisotropy axes playing the role of the Ising axes. In this limit as is the the case for the D/J → ∞ of the Heisenberg RAM 28 , one is left with a much simpler system since the continuous 3D unit vector per site transforms in the two valued scalar variable of the Ising model. Therefore, it is important to determine the value of λ u beyond which the system behaves as a dipolar Ising system. To do this we note that when λ u ≫ 1 the moments directions μi get close to the axes ni leading to μi ≃ s i ni with s i = ±1 and thus the Hamiltonian from equation (2b) can be written as

βH = 1 T * (1/2) i≠j s i J ij s j -λ u i (μ i .n i ) 2 (11)
Where J ij are the coupling constants ni Tij nj , independent of the configuration of the moments {μ i } in the framework of the frozen easy axes distribution. Equation ( 11) reflects the fact that in this limit the MAE depends on the fluctuations ((μ i .n i ) 2 -1) conversely to the DDI which depends on the configuration of the {s i }. We then have two independent sets of variables per site namely {s i } and (μ i .n i )

2 . The second one, describing the MAE are totally uncoupled and represent two degrees of freedom per site (μ i are 3-dimensional unit vectors) and thus contribute, according to the Dulong and Petit law, or equivalently to the equipartition of energy, as N k B to the total heat capacity, or in other words

(C = C v /(N k B )) C an (T ) = 1 N k B d[< E an >] dT → 1 (12)
Then equation ( 12) is a clear criterion to determine whether the system of interacting dipoles plus MAE is Ising like since it is an indication that the variables {s i } and (μ i .n i ) 2 become uncoupled and equation (11) holds. In such a case, we expect the dipolar component of C, C dip to be close to the C of the corresponding dipolar Ising model which does not include the component C an , as for instance in the RAD case. From the heat capacity curves displayed in figure (11) we conclude that the Ising dipolar behavior is reached beyond λ u = 30. It is worth mentioning that a similar conclusion was obtained in the case of the totally textured distribution of easy axes

IV. CONCLUSION

In this work, we have determined a significant part of the magnetic phase diagram of an ensemble of dipoles with uniaxial anisotropy located on the nodes of a FCC lattice from tempered Monte Carlo simulations. This is motivated first by the search for the conditions under which a super-ferromagnetic phase induced by DDI can be reached in the supra-crystals of MNP synthesized experimentally, and more generally by the determination of the nature of the ordered low temperature phase in these systems. The nature of the low temperature ordered phase is found successively FM and SG with increasing the MAE strength. From the behavior of m 2 with respect to L at low T * beyond λ u = 0.4 we interpret the FM phase in this region of the phase diagram as a FM-QLRO phase. Moreover the FM-QLRO phase is related to a transverse spin-glass state.

The crucial points to relate our findings to actual experimental situations is to determine the corresponding value of the λ u parameter and the reference temperature T r defined as T * = T /T r (see section II). λ u and T r are related to the physical characteristics of the MNP, the uniaxial anisotropy constant K and the saturation magnetization M s by λ u = (24/µ 0 )(K/M 2 s )(Φ r /Φ) and T r = µ 0 /(4πk B )(π/6 M s ) 2 d 3 (Φ/Φ r ) respectively from the definitions given after equations (1,2a). As typical examples, in the case of maghemite or cobalt MNP we get a lower bound λ u > 5 which means that in the absence of texturation of the easy axes, the corresponding supra-crystals are in the SG region of the phase diagram at low temperature. A precise determination of T r and thus of the transition temperature T c = T * c T r necessitates the precise knowledge of both the MNP size d and M s which because of finite size and surface chemistry effects deviates from its bulk value, making thus T r specific to a given system, and not only to a given material. One may also wonder on the effect of the temperature dependence of the relevant experimental parameters, M s and K on the validity of such a model. Concerning M s one has to keep in mind that the ordering temperature (either the Néel or the Curie temperature) of the usual MNP materials (c.a. ∼ 850 -1000 K) is much larger than the room temperature which is somewhat an upper bond for the MNP properties of interest and thus only a very small variation of M s (∼ M s (0) ( Finally, we also have obtained the limiting value beyond which the DDI plus uniaxial MAE model is assimilable to a 1D dipolar Ising model. This value , λ u ≃ 30, seems larger that the experimental expectation of λ u . As a result, modeling the MNP ensemble by a dipolar Ising model may be not strictly speaking justified.

V. ACKNOWLEDGEMENTS.

This work was granted an access to the HPC resources of CINES under the allocations 2018-A0040906180 and 2019-A0060906180 made by GENCI, CINES, France. We thank the SCBI at University of Málaga and IC1 at University of Granada for generous allocations of computer time. J.-J. A. thanks for financial support from grant FIS2017-84256-P (FEDER funds) from the Agencia Española de Investigación.

FIG. 1 .

 1 FIG. 1. Phase diagram in the (T * , λ u ) plane separating the PM, FM, SG regions. We also indicate the estimated range in λ u corresponding to the FM-QLRO and transverse spin-glass phase. Open upward triangles: PM/FM line; open squares: PM/SG line; open downward triangles : transverse spin-glass transition temperature, T xy ; open circles: FM/SG line obtained from the behavior of B m in terms of L at constant T * .

FIG. 2 .

 2 FIG. 2. Polarization m d1 , open symbols and nematic order parameter, [< P 2 >] solid symbols in terms of T * for a) λ u = 0 , b) λ u = 0.8 and c) λ u = 1.6 for different system sizes L = 4 to 8 (N = 256 to 2048) for λ u = 0 and 0.8, and L = 4 to 7 for λ u = 1.6, as indicated. In c) the dashed line represents P 2 for λ u = 2 and L = 7 to illustrate the gradual vanishing of the nematic order at low T * with the increase of λ u . The lines result from the interpolation performed by the reweighting method 40 .

8 FIG. 3 .

 83 FIG. 3. Log-log plot of m 2 in terms of L for : a) λ u = 0.4 and T * = 0.2028, 0.3050, 0.4481, 0.5402, 0.5913, 0.6833, 1.910 from top to bottom; b) λ u = 0.8 and T * = 0.1983, 0.3084, 0.4429, 0.5651, 0.6263, 0.6385, 0.6874, 1.9100. The dotted lines are guides to the eyes for m 2 = cst. at the smallest temperature; the solid lines are linear interpolation determined from L = 4 and L = 6. excepted at bottom where they indicates the 1/L 3 behavior in the PM phase. The arrow localizes T

(

  see figure (4 a,b)). Another important feature of the PM/FM transition is the strong finite size dependence of both the heat capacity C v and the magnetic susceptibility χ m with a peak located at T * L (A), (A = C v , χ m ), and T * L (A) → T * c in the L → ∞ limit 26,49 . This behavior, leading to a divergence of χ m (not shown) at T c is displayed on figure (5a) for C v 50 .

  2 b)). Although we cannot rule out a crossover for much larger values of L leading to a convergence of m(T * ) at fixed and low T * we interpret this behavior as a FM state with quasi long range order with m(L → ∞) = 0 well below T * c , instead of a long range FM one. In the absence of MAE, the different m(L, T * ) merge at low T * (see figure (2 a) and Ref. 24 ) which is no longer the case beyond λ u ∼ 0.4 as can be seen in figures (2 b) and (3 c). This FM-QLRO is observed beyond λ u ≃ 0.4 as shown in figure (3) where m 2 in terms of L in log scale (a and b) and m 2 in terms of T * for different values of L (c) are displayed for λ u = 0.4 and 0.8. In any case the paramagnetic behavior, namely m 2 ∼ 1/N , is reached sufficiently far from T * c . The FM-QLRO is related to the onset of the transverse spin-glass state 51 where the components of the moments normal to the nematic direction freeze in a spin-glass like state and corresponds to the mixed state of the phase diagram outlined by Ayton et al. 22 . The transverse spin-glass sate is defined in the FM region of the phase diagram and characterized by the transverse overlap order parameter, q 2t , obtained by using the transverse components μit = (μ i -( d. μi ) d)/ (μ i -( d. μi ) d) of the moments instead of μi in equation (

8 FIG. 4 .

 84 FIG. 4. Magnetization Binder cumulant for different ssytem sizes as indicated. In the pure dipolar case (a) and b)). a) in terms of T * ; b) rescaled according to equation (10) with t * = (T * -T * c )L 1s/ν , T * c = 0.614 and ν = 0.692 48 .For different values of λ u , as indicated (c) to f)); including the vicinity, λ u = 1.4, e) and f) beyond, λ u = 1.5, the onset of the SG phase of the magnetic phase diagram. The lines are the interpolations according to the reweighting method40 . The dotted line in d) to f) is a guide to the B m = 1 limiting value, which is no longer reached for λ u ≥ 1.2.

FIG. 5 .

 5 FIG. 5. Evolution of the heat capacity C v with the system size L for a): the pure dipolar case (λ u = 0) and L = 4 to 8 and b): λ u = 1.4 and L = 4 to 7. The error bars are smaller than the size of the symbols. The lines result from the interpolation performed by the reweighting method 40 .

u = 1 . 3

 13 and we expect the low temperature phase to present a SG character beyond this value.

FIG. 7 .FIG. 8 .

 78 FIG. 7. Transverse spin-glass correlation length reduced by L for λ u = 0.4, bottom and 0.8, top, showing the crossing point behavior in terms of T * for different values of L. The vertical lines visualize the value of the PM/FM transition temperature T * c obtained from B m .

FIG. 9 .FIG. 10 .

 910 FIG.9. Spin glass Binder cumulant B q , a) and second moment of the spin-glass overlap order parameter, q 2 , open symbols, and of the magnetization, m d2 , solid symbols, b), for λ u = 4 in terms of the system size L.

FIG. 11 .

 11 FIG. 11. Anisotropy, a) and dipolar, b) components of the heat capacity. In the RAD case the total capacity includes only the C dip component. Symbol: L = 7 for λ u ≤ 4 and L = 4 otherwise. Solid blue lines for λ u ≤ 4: L = 4.

1 -

 1 (T /T c )) b , b ∼ 3/2) is expected. Concerning K, since it is in general an effective anisotropy constant including different contributions (magnetocrystalline, shape, surface chemistry), a precise evaluation of K(T ) is not a simple task. However, our results show only a very weak dependence of the transition temperature with respect to λ u . On an other hand an important bound is the T b /T r ratio, where T b is the blocking temperature related to the MAE through T b = Kv(d)/(ck B ), the constant c being specific to the measurement considered (c ≃ 30 for magnetometry). Below T b the moments are blocked in the MAE potential well during the measurement and thus the DDI induced transition can be observed only if T b < T r T * c . From the definition of T b we easily get T b /T r = λ u /c and consequently the DDI induced transition will be experimentally observable only if λ u < cT * c , the sufficient condition being in any case λ u < c/2 according to our results.