###### Abstract

The magnetic anisotropy and exchange coupling between spins localized at the positions of 3d transition metal atoms forming two-dimensional metal-organic coordination networks (MOCNs) grown on the Au(111) metal surface are studied. In particular, we consider MOCNs made of Ni or Mn metal centers linked by TCNQ (7,7,8,8-tetracyanoquinodimethane) organic ligands, which form rectangular networks with 1:1 stoichiometry. Based on the analysis of X-ray magnetic circular dichroism (XMCD) data taken at T= 2.5 K, we find that Ni atoms in the Ni-TCNQ MOCNs are coupled ferromagnetically and do not show any significant magnetic anisotropy, while Mn atoms in the Mn-TCNQ MOCNs are coupled antiferromagnetically and do show a weak magnetic anisotropy with in-plane magnetization. We explain these observations using both a model Hamiltonian based on mean-field Weiss theory and density functional theory calculations that include spin-orbit coupling. Our main conclusion is that the antiferromagnetic coupling between Mn spins and the in-plane magnetization of the Mn spins can be explained neglecting effects due to the presence of the Au(111) surface, while for Ni-TCNQ the metal surface plays a role in determining the absence of magnetic anisotropy in the system.

Departamento de Física de Materiales UPV/EHU, 20018 Donostia-San Sebastián, Spain

Donostia International Physics Center (DIPC), 20018 Donostia-San Sebastián, Spain

Departamento Física Aplicada I, Universidad del País Vasco, 20018 Donostia-San Sebastián, Spain

Department of Materials, ETH Zürich, Hönggerbergring 64, 8093 Zürich, Switzerland

Swiss Light Source, Paul Scherrer Institute, 5232 Villigen PSI, Switzerland

Centro de Física de Materiales (CFM-MPC), Centro Mixto CSIC-UPV/EHU, 20018 Donostia-San Sebastián, Basque Country, Spain

IKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Basque Country, Spain

Tomsk State University, 634050 Tomsk, Russia

Saint Petersburg State University, 198504 Saint Petersburg, Russia

Academic Editor: Jordi Cirera

Version July 15, 2019 submitted to MDPI

Keywords: magnetism; metal-organic network; XMCD; Density Functional Theory

## 1 Introduction

There exists an exciting type of two-dimensional systems that can be grown on surfaces by self-assembly techniques with interest both from a fundamental point of view and also due to its potential application in the fabrication of electronic and spintronic devices. They are called metal-organic coordination networks (MOCNs) and consist of metal centers linked by organic ligands that permit, in principle, the design of overlayers with specific electronic and magnetic properties [1]. The synthesis and growth of a given MOCN with a given composition, essentially defined by its stoichiometry and coordination, depends on the relative strength of the interactions between the constituents (organic ligands and metal centers) and their interaction with the underlying surface [2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Indeed, the chemical state of the organic ligands and metal centers can be modified due to vertical electronic charge transfer from the surface [12]. Additionally, lateral charge transfer between the MOCN’s constituents is crucial for the bonding and equally important for the electronic and chemical properties of the overlayers. Particularly interesting is the role of the metal centers in the formation of the two-dimensional networks by favoring a given coordination and stoichiometry determining the charge and magnetic moment of the metal center and, occasionally, also of the organic ligand that can acquire spin polarization. The important point is that this very fact could be used to control the electronic and magnetic properties of the interface.

The case of 3d transition metal atoms as metal centers and molecules with large electronegativity, like TCNQ (7,7,8,8-tetracyanoquinodimethane) or F4TCNQ (2,3,5,6-tetrafluoro-7,7,8,8-tetracyanoquinodimethane), on metal surfaces is of special interest because they form well-ordered MOCNs with few defects [13, 4, 7] and different stoichiometry; this latter depending both on the underlying surface and preparation conditions. The experimental techniques typically used to characterize the geometric structure and chemical composition of MOCNs on surfaces are scanning tunneling microscopy (STM), low-energy electron difraction (LEED) and X-ray photoemission spectroscopy (XPS), while for the electronic and magnetic properties the standard techniques are angle resolved photoemission spectroscopy (ARPES) , scanning tunneling spectroscopy (STS) and X-ray magnetic circular dichroism (XMCD).

The theoretical description of this type of systems has been shown to be quite reasonable using standard electronic structure methods, like density functional theory (DFT) [13, 4, 11, 14], although one has to be aware of the limitations of each method when aiming at achieving quantitative agreement with the experimental data. However, the explanation of the observations at a qualitative level and its understanding with the help of a model Hamiltonian is the recipe that we follow in this work. In any case, it is worth to mention that an essential problem, which is hard to overcome, concerns the accuracy of the calculations when dealing with very low energy scales, as it is the case of the determination of exchange couplings or magnetic anisotropies in the sub-meV energy range. Apart from this limitation imposed by the methodology itself, it is also important to stress that exchange coupling and magnetic anisotropy energy are extremely sensitive to both slight geometrical distortions or band filling, i. e., electronic charge transfer. Therefore, it is important to balance the different effects that each and every approximation can have in the final results when trying to explain an observation that can be deduced, e. g., from XMCD data, regarding the strength and type (ferro- or antiferro-magnetic) of exchange coupling or the kind of magnetic anisotropy (easy axis or easy plane).

The exchange coupling between magnetic centers in two-dimensional MOCNs is affected to a major or lesser extent by the underlying substrate. The presence of the surface represents a difficulty to describe the full system (MOCN/surface) because the overlayer structure is not necessarily commensurable with the crystal surface or because the size of the commensurable supercell is too large. In the case of weak coupling between the overlayer and the surface, as it is the case of Au(111) surfaces, the essential features can be described by neglecting the role of surface electrons, in first approximation. As a rule of thumb, when the lateral bonds between the metal centers and the organic ligands are much stronger than the metal-surface or ligand-surface bonds, this approximation is expected to be a reasonable way to describe the magnetic coupling between metal atoms in the MOCN. However, in case of lateral (metal-molecule or intermolecular) and vertical (MOCN-surface) couplings of comparable strengths, an explicit inclusion of the surface is required to describe the system, something that could happen either due to strong coupling with the surface or weak lateral coupling. Next, one should consider the role of charge transfer between surface and overlayer, even in the case of weak coupling, as it can be important in determining magnetic moments, magnetic coupling or even magnetic anisotropy. Indeed, when intermolecular coupling is weak, the role of surface electrons can be relatively more important in determining the magnetic coupling between spins of the metal centers via RKKY interaction [15], which may appear not only on metal surfaces but also on the surface of topological insulators [16, 17]. Very recently, it has been proposed that the RKKY interaction is responsible for the long range ferrimagnetic order [18] in a two-dimensional Kondo lattice with underscreened spins by the conduction electrons in a FeFPc - MnPc mixture on the Au(111) surface.

In this work, we consider the case of MOCNs that consist of Mn or Ni magnetic atoms and TCNQ or F4TCNQ molecules grown on Au(111) surfaces, which show 1:1 stoichiometry with each metal center (Mn or Ni) coordinated with four organic ligands. For these systems, our XMCD data show a substantial difference between Mn and Ni networks both for TCNQ and F4TCNQ ligands that we explain using a model Hamiltonian approach and DFT calculations. The results of our calculations for the free-standing neutral Mn-TCNQ overlayers are consistent with both the antiferromagnetic (AFM) coupling between Mn centers and the weak magnetic anisotropy with in-plane magnetization, while for Ni-TCNQ overlayers we need to call for effects due to the presence of the underlying metal surface, like charge transfer and changes in coordination, to explain the absence of anisotropy in the system. Model calculations based on mean-field Weiss theory permit us to extract exchange coupling constants from the fits to XMCD curves, as well as to obtain additional information about the magnetic anisotropy and the different magnetic configurations that may appear in the networks. Finally, it is worth mentioning that, although the organic ligands TCNQ and F4TCNQ have a different electronegativity (higher in F4TCNQ than in TCNQ), based on the acquired XMCD data, there are not substantial differences in the magnetic properties of the corresponding Ni and Mn networks. Therefore, in the core of the paper we present the results for TCNQ networks and leave the F4TCNQ results for Section I of the Supplementary Material.

The paper is organized as follows. After describing the XMCD experiments and the technical details of the calculations in section 2, we present our XMCD data for Mn-TCNQ and Ni-TCNQ on Au(111) in Section 3.1 , together with fitting curves from model calculations that permit to explain the observations and extract information about the type of magnetic coupling and magnetic anisotropy (Section 3.2). Next, in Section 3.3 we present the results of our spin polarized DFT+U electronic structure calculations for Mn-TCNQ and Ni-TCNQ free-standing overlayers that confirm the observed behavior in the type of magnetic coupling between spins at the 3d metal centers. Then, in Section 3.4, we present the magneto-crystalline anisotropy analysis of the two considered systems under study based on calculations that include spin-orbit coupling (SOC). Finally, in Section 4, we present a discussion of our findings and establish the main conclusions that aim at explaining the XMCD observations and suggest that for Ni-TCNQ networks the Au(111) metal surface plays a role in determining the magnetic properties of the MOCN, while this is not the case for Mn-TCNQ.

## 2 Materials and Methods

### 2.1 XMCD experiments

The x-ray absorption spectroscopy (XAS) experiments were carried out at the X-Treme beamline of the Swiss Light Source. The samples were prepared in ultra-high vacuum chambers with a base pressure in the range of low 10 mbar. The pressure in the magnet-cryo-chamber was always better than 10 mbar. The Au(111) surface was cleaned by repeated cycles of Ar sputtering and subsequent annealing to 800 K. The molecules TCNQ (7,7,8,8-tetracyanoquinodimethane, 98% purity, Aldrich) and F4TCNQ (2,3,5,6-tetrafluoro-7,7,8,8-tetracyanoquinodimethane, 97% purity, Aldrich) were thoroughly degassed before evaporation. The organic adlayers were grown by organic molecular-beam epitaxy (OMBE) using a resistively heated quartz crucible at a sublimation temperature of 408 K onto the clean Au(111) surface that was kept at room temperature. The coverage of molecules was controlled to be below one monolayer. Ni or Mn was subsequently deposited using an electron-beam heating evaporator at a flux of about 0.01 ML/min on top of the molecular adlayers that were heated to 350-400 K to promote the network formation. The sample was checked in-situ by STM at the beamline and subsequently transferred to magnet chamber without breaking the vacuum.

The polarization-dependent XAS experiments were performed in total electron yield detection. Magnetic fields were applied collinear with the photon beam at sample temperatures between 2.5 and 300 K. The data were acquired by varying the photon energy at the L edges of Ni and Mn, as well as the K edges of O and N using circular and linear polarized light. The absorption spectra were normalized with respect to the total flux of the incoming x-rays and were further treated to be normalized to the absorption pre-edge due to TEY variations. The background obtained from clean or molecule-covered Au(111) was subtracted to allow comparison of the spectral features. The XMCD is obtained from the difference of the left and right circular polarized XAS spectra whereas the XAS is obtained from the average of the two circular polarizations. The sample was rotated between normal x-ray incidence with respect to the sample surface at and grazing incidence with . All shown spectra were acquired at K at external magnetic fields up to T. The magnetization curves were recorded by acquiring the maximum of the XMCD signal at the L edge as a function of the external magnetic field, normalized by the corresponding pre-edge of the XAS signal. To facilitate the extraction of the easy and hard magnetization axis the magnetization curves at different angles of the magnetic field were normalized to the same value at the highest magnetic field point.

### 2.2 Density Functional Theory calculations

DFT calculations have been carried out using the Vienna Ab Initio Simulation Package (VASP) [19, 20, 21]. For the description of electron-ion interactions the Projector Augmented-Wave (PAW) method has been employed, whereas the Perdew, Burke and Ernzerhof (PBE) [22] functional has been used to describe exchange and correlation within the generalized gradient approximation (GGA). A Hubbard-like Coulomb repulsion correction term (U= 4 eV) has been added to describe the 3d metal electron states, based on Dudarev’s approach [23] , as implemented in VASP. A previous study [11] has already corroborated that the results concerning magnetic moments and 3d level occupations do not change appreciably in the 3-5 eV range of the U parameter.

For the geometrical optimization of the free-standing Mn-(F4)TCNQ and Ni-(F4)TCNQ systems, periodic supercell boundary conditions have been imposed. The optimal cell dimensions and atomic positions have been obtained by an energy minimization procedure with a convergence criterion of eV for the energy and 0.02 eV/Å for the forces to assure that we reach sufficient accuracy in numerical values of the calculated magnitudes. The Kohn-Sham wave functions have been expanded in a plane wave basis set with a kinetic energy cutoff of 400 eV for all the systems considered. A Monkhorst-Pack k-point sampling [24] equivalent to in the surface unit cell and Methfessel-Paxton integration with smearing width 0.1 eV [25] have been used. Symmetry considerations have been switched off from the calculations and a preconverged charge density with a fixed value of the total spin for the unit cell has been used to relax all the networks. For the obtained relaxed geometries, where the layer is constrained to be flat, we have evaluated the magnetic anisotropy energies with adjusted parameters. Total energies have been converged with a tolerance of eV. A k-point sampling and the corrected tetrahedron method of integration [26] have been used instead of smearing methods.

Fig. 1 shows a top view visualization of the systems considered. The optimized geometrical parameters are included in Table 1, where and denote the lattice vectors while d and d denote the values of the Mn-N or Ni-N bond lengths indicated in Figure 1.

cell | Mn-TCNQ | Ni-TCNQ |
---|---|---|

11.52 | 11.32 | |

7.38 | 7.16 | |

d | 2.12 | 2.01 |

d | 2.12 | 1.95 |

## 3 Results

### 3.1 XMCD data

The XMCD intensity variation as a function of the applied magnetic field () defines a curve that is proportional to the system magnetization. Therefore, when the value of the spin magnetic moments at the metal centers () , the temperature () and the Landé -factor are known, one can use simple models to simulate the magnetization response. A good reference to be considered is the case of paramagnetic behavior (spins responding individually to the applied magnetic field) that can be represented by a Brillouin function. Whenever a preference for ferromagnetic (FM) or antiferromagnetic (AFM) coupling between spins appears, the corresponding magnetization curves will show higher or lower curvature, respectively, than the corresponding Brillouin function for the same , and -factor values. In this way, in principle, one can decide about the type of magnetic coupling between localized spins at the metal centers, as long as the value of the spin () is known. Note that, in the presence of strong magnetic anisotropies and high orbital angular moments, the analysis becomes more involved [27]. However, here, we can follow this simplified scheme, as shown below. According to our DFT calculations, described in section 3.3, Mn atoms in Mn-TCNQ have a localized spin magnetic moment close to , although somewhat lower, while Ni atoms in Ni-TCNQ have a much lower spin close to , although somewhat higher. Therefore, we use the values and for Mn and Ni, respectively, to do our XMCD analysis that includes fitting curves to XMCD data based on Weiss mean-field theory described in the next section, where and are defined, and also a comparison with the corresponding Brillouin functions.

The results are shown in Fig. 2 a) and b) for Mn-TCNQ and Ni-TCNQ, respectively. It is evident that in Mn-TCNQ the coupling between Mn spins is AFM, while in Ni-TCNQ it is FM. Additionally, the fitted values of the exchange coupling constants reveal a weaker coupling between Mn spins (meV) , as compared to the coupling between Ni spins ( meV), while the single ion anisotropy parameter meV corresponds to a weak anisotropy with in-plane magnetization for Mn-TNCQ and to the absence of anisotropy for Ni-TCNQ. In order to learn more about the magnetic anisotropy of these systems, in Fig. 3 we plot a comparison of XMCD data obtained for perpendicular and grazing incidence for Mn-TCNQ and Ni-TCNQ, the former showing a mild angular dependence with stronger intensity for grazing incidence, i. e., a fingerprint of magnetic anisotropy in the system with in-plane magnetization. Incidentally, this weak anisotropy is only observed at low temperatures. However, in the Ni-TCNQ XMCD data there is no significant angular dependence, which means a negligible magnetic anisotropy. A value of the Ni atom spin corresponds to the absence of single ion anisotropy [28].

### 3.2 Model for Mn-TCNQ and Ni-TCNQ

In Mn-TCNQ, the coupling between local moments is antiferromagnetic and occurs by means of the Anderson superexchange mechanism [29, 15]. In perturbation theory, the superexchange interaction was found [11] to be dominated by a virtual process in which two electrons hop from the doubly occupied LUMO of the TCNQ molecule onto two adjacent Mn atoms. Inclusion of additional molecular orbitals, such as the HOMO, leads to a generic superexchange interaction with coupling constants , , and , as shown in Fig. 4. The model Hamiltonian describing the magnetic properties of Mn-TCNQ, thus, reads

(1) |

where denotes the local moment of the Mn atom () on site , is the single-ion anisotropy energy, g is the Landé g-factor (), and is the magnetic field. The Heisenberg exchange constant is restricted to the nearest ( and ) and next-to-nearest () neighbors on the rectangular lattice. The summation in the Heisenberg interaction term accounts twice for each pair of interacting sites; hence the presence of a factor in Eq. (1).

A quick insight into the tendency to order the spins in this model is granted by the Fourier transform of the exchange coupling ,

(2) | |||||

where is the two-dimensional wave vector and is the position of the Mn atom on site . For ferromagnetic couplings (), the maximum of occurs at , which indicates that the spin order could be uniform from a mean-field point of view, not addressing the question about its stability against fluctuations in two dimensions. Additional terms, such as the single-ion anisotropy or the Zeeman interaction, may stabilize the uniform spin order.

In contrast, for antiferromagnetic couplings (), the maximum of occurs usually at the edge of the Brillouin zone, indicating that the magnetization is staggered in some way over the unit cell. When only nearest neighbors are coupled ( and ), the maxima lie at and its equivalent points, which results in the usual checkerboard-like antiferromagnetic order [see Fig. 4 (a)]. As the diagonal coupling is turned on (assuming an antiferromagnetic ), for a sufficiently large magnitude of there is a transition from the checkerboard pattern to a so-called superantiferromagnetic state of antiferromagnetically ordered rows or columns. For , by requiring at , we find at the transition point for antiferromagnetic column formation [see Fig. 4 (b)].

The effect of the diagonal coupling consists in introducing magnetic frustration [29, 30] in the spin lattice. We remark here that the special point is realized to a good approximation in our Mn-TCNQ lattice, because (i) the LUMO of the TCNQ molecule has overlap with the and orbitals of the Mn atom, thus, dominating the superexchange, and (ii) the direct coupling between the LUMOs of neighboring TCNQ molecules is rather weak. The latter makes it possible to consider two independent paths of superexchange for the nearest neighbors, with each path going separately via one of the two TCNQ molecules connecting the two neighboring Mn atoms. For the diagonal coupling, only one path is possible, which leads to a reduction of the diagonal coupling by a factor of as compared to the nearest-neighbor coupling. With the approximations (i) and (ii), the coupling constants obey (see Ref. [11] for further details).

Despite the fact that the Mn-TCNQ lattice may well be in a frustrated magnetic state consisting of a mixture of the two phases in Fig. 4, the XMCD data appear to be consistent with a much simpler description of the magnetization as a function of the -field, which is derived from the Weiss mean-field theory, and it captures faithfully weak deviations from the paramagnetic state. The superexchange couplings are rather weak, [11] of the order of , and the Zeeman term soon dominates. Additionally, there exists a fair amount of single-ion anisotropy, described by the term in Eq. (1).

We make the mean-field approximation for the model in Eq. (1),

(3) |

where gives the local description of the interacting system in terms of the Weiss fields . The spin averages can be regarded as variational parameters of the theory. The last term in the first line of Eq. (3) compensates for the double counting of interaction energy occurring in the local Hamiltonian and plays an important role when calculating the free energy of the interacting system. The minimization of the free energy allows us to determine the values of the order parameters . The procedure is described in the Appendix.

Next, we focus on the XMCD data taken at normal incidence (), for which the magnetic field is applied along the -axis, . For the (checkerboard) antiferromagnetic phase, we use two order parameters and , which represent the -components of the spins in the unit cell as shown in Fig. 4 (a), and minimize the upper bound to the free energy [] with respect to the order parameters and . Alternatively, one can require stationarity of free energy, and , which yields two coupled equations,

(4) |

where is the free energy of a single isolated spin. The mean-field solution is obtained from these self-consistent equations. As a rule, several solutions are found. The choice of the physical solution relies again on the least value of the free energy. For the superantiferromagnetic phase, we use again two order parameters, and , but now they are distributed in the unit cell as shown in Fig. 4 (b). The mean-field approximation takes into account only the connections (i.e. bonds) between the spins on a local scale, whereas the constrains related to the dimensionality of the systems go unaccounted. We can, therefore, adapt here all the results derived for the phase in Fig. 4 (a) by simultaneously replacing and in all expressions as

(5) |

The factors and appear here because each connector counts as half a bond in the unit cell, whereas each connector counts as a full bond.

We fit the experimental data for normal magnetic fields in Fig. 2 assuming the relation , which corresponds to the case when a single orbital of the ligand is dominating the superexchange. We reach a good fit to the experimental data for . Our working assumption was that the critical temperature () is sufficiently low as to allow application of the Weiss theory, i.e. . This means also that the order parameters and are never of opposite sign and are, in fact, equal to each other over the full range of applied magnetic fields. Therefore, the experimental data can equally well be fitted by a ferromagnetic mean-field theory with antiferromagnetic coupling constants. To simplify the matter even further, we consider a square lattice with a single coupling constant . Effectively, this coupling constant will be related to the previous coupling constants by equating to each other at for both models, which immediately yields . Using the above value, we arrive at and for .

The same effective model derived from a mean field Hamiltonian with , and parameters can be used for Ni-TCNQ, although its relation with the microscopic Hamiltonian described in Ref. [11] is different. In this case, we find a good fit with and for .

### 3.3 Spin polarized DFT+U calculations

We first consider a two-dimensional free-standing overlayer description for Mn-TCNQ and Ni-TCNQ networks. Both the lattice vectors and atomic positions have been optimized by using an energy minimization procedure within DFT, as described in the Materials and Methods section. The projected densities of states (PDOS) onto different atomic orbitals of the Mn and Ni atoms are shown in Figs. 5 and 6, respectively. The insets show the PDOS onto atomic orbitals of the C and N atoms of the organic ligand, as well as onto Mn and Ni 3d states without m number resolution, in a narrow energy range close to the Fermi level. A close inspection of Figs. 5 and 6 reveals important differences between the two systems under study. The most significant is the half-filling of the states with all the majority spin states occupied in Mn-TNCQ, which corresponds to a value of the spin approximately equal to localized at the Mn atoms, while in Ni-TCNQ only one minority spin state is fully unoccupied (), which corresponds to a value of the spin localized at the Ni atom approximately , although somewhat higher, as the minority spin states and are partially occupied. Additionally, in Ni-TCNQ, the and states are hybridized with TCNQ orbitals close to the Fermi level, in particular the LUMO, giving rise to a delocalized spin density [11]. This can be seen by comparing the PDOS onto atomic orbitals of the C and N atoms of the TCNQ organic ligand shown in the insets of Figs. 5 and 6 for Mn-TCNQ and Ni-TCNQ, respectively. In Ni-TCNQ the LUMO orbital is spin polarized but this is not the case in Mn-TCNQ, for which the TCNQ LUMO practically does not hybridize with Mn states and it is fully occupied.

Next, using these two optimized structures calculated with a surface unit cell within the DFT+U method with spin polarization as starting point, we proceed to double the size of the surface unit cell into a cell that contains two metal centers (Mn or Ni atoms) and two TCNQ molecules. In this way, we can decide which is the most favorable type of magnetic coupling ( ferro- or antiferro-magnetic) between spins localized at the Mn or Ni centers by comparing the values of the corresponding total energies. We consider a checkerboard configuration using oblique vectors in the surface unit cell and confirm that ferromagnetic coupling is favorable in Ni-TCNQ networks, while in Mn-TCNQ networks antiferromagnetic coupling is preferred in agreement with Ref. [11]. The corresponding spin densities are shown in Fig. 7 for Mn-TCNQ and Ni-TCNQ. In section II of the Supplementary Material we also include other configurations obtained by using a rectangular surface unit cell, in which other AFM configurations with spins aligned in rows or columns are considered as well [31], showing the importance of next to nearest neighbors (diagonal) couplings in the networks that has been discussed in the previous section.

### 3.4 Magnetocrystalline anisotropy

The magnetocrystalline anisotropy energies (MAE) can be obtained from DFT calculations that include SOC effects. The resulting total energies, thus, depend on the orientation of the magnetization density. For extended systems, where the transition metal atomic orbital momentum is expected to be partially or totally quenched, the MAE appears as a second order SOC effect. In systems where the PDOS is characterized by sharp peaks and devoid of degeneracies at the Fermi level, a second-order perturbative treatment of the SOC makes it possible to establish a few guidelines for the likelihood of an easy axis or plane. The perturbation couples states above and below the Fermi level and it is inversely proportional to the energy difference between states. When the spin-up -band is completely filled, it can be shown that the energy correction is proportional to the expected value of the orbital magnetic moment and that the spin-flip excitations are negligible [32, 33, 34].

The total energy variation as a function of the magnetization axis direction is very subtle, often in the sub-meV range per atom. When spin-orbit effects are not strong, it is common practice to use the so-called second variational method [35], where SOC is not treated self-consistently. First, a charge density is converged in a collinear spin-polarized calculation. Next, a new Hamiltonian that includes a SOC term is constructed and diagonalized for two different magnetization directions. Then, the MAE is calculated from the difference of the two band energies. Alternatively, a more precise MAE can be obtained from total energy calculations that include SOC self-consistently. Using the latter method, in this work we have calculated MAE values for free-standing Mn-TCNQ and Ni-TCNQ networks.

The small energies involved in the anisotropy are a challenge for DFT calculations. The MAE is highly sensitive to the geometry and electronic structure calculation details, such as the exchange and correlation functionals and basis set types. From a technical perspective, a reliable MAE is only achieved with demanding convergence criteria. For example, it has been observed that fine k-point samplings of the Brillouin zone are needed [36, 37, 38]. An account of the convergence details as well as MAE dependence on the parameter can be found in section III of the Supplementary Material.

Table 2 shows the obtained values for eV in cells (i.e. , only ferromagnetic ordering is considered in this Section). For the Mn-TCNQ rectangular network, we find in-plane magnetization with negligible azimuthal dependence, i.e., easy plane anisotropy. The MAE, calculated as the total energy difference between magnetic configurations with Mn magnetizations parallel to the OX and OZ axes, is 0.2 meV. In the Ni-TCNQ rectangular network the energetically preferred magnetization is out-of-plane and the MAE values vary significantly with the azimuthal direction. As shown in Table 2, the values change as much as 1.50 meV with azimuthal angle variations.

Mn | 0.20 | 0.19 | 0.20 | 0.20 |

Ni | -1.44 | -0.95 | -1.95 | -0.45 |

Ni (oblique) | -0.07 | -0.04 | 0.03 | -0.09 |

The different behavior of the MAE with the azimuthal angle in Mn and Ni networks can be understood in terms of the differences in the metal-molecule bonds, particularly the Mn-N and Ni-N bonds. In both networks the (with magnetic quantum number ), , and orbitals remain rather localized, whereas the and orbitals are spread over a wider energy range of a few eV below the Fermi level (see Figs. 5 and 6). The delocalization of electronic charge in these and orbitals is stronger in the Ni-TCNQ case, where the latter two sub-bands are partially occupied and form hybrid states at the Fermi level with the TCNQ LUMO. As these hybrid states lie at the Fermi level, they have a dominant role in the magnetic anisotropy and, since they yield markedly directional charge and spin density distributions along the Ni-N bonds, they are likely to produce azimuthal MAE variations. Conversely, the Mn -electrons hybridize weakly with the TCNQ orbitals close to the Fermi level, i. e., with the LUMO, and have essentially no weight at the Fermi level. The spatial extent of these relevant Ni-TCNQ hybrid states is manifested in the delocalized electron spin densities depicted in Fig. 7 b), as compared to the case of Mn-TCNQ shown in Fig. 7 a) with a spin density more localized at the Mn sites and its neighboring cyano groups.

The existence of an easy axis (plane) of magnetization for Ni (Mn) cannot be anticipated from the electronic structure details. In the Mn-TCNQ system, since the -band is half-filled, the MAE is led by spin-flip excitations and, therefore, the value of the exchange splitting is determinant. In the absence of same-spin excitations, the anisotropy would be associated to the anisotropic part of the spin distribution. More precisely, the MAE would be proportional to the anisotropy of the expected values of the magnetic dipole operator [33, 34]. However, Fig. 7 a) shows an anisotropic spin distribution extended towards the cyano groups of the organic ligand TCNQ in the network plane by the crystal field. The quadrupolar moment of this distribution should promote out-of-plane magnetization. This interpretation is at variance with the SOC-self-consistent DFT result. A more elaborated model has been proposed for systems with localized -orbitals. It states that the spin-flip excitations that keep the quantum number constant favor an in-plane magnetization [39]. The calculated PDOS of Fig. 5 shows that the two peaks () are those closer to the Fermi level, for majority and minority spin states, respectively. This situation is, in principle, compatible with an easy plane behavior. The conclusion we draw is that the basic qualitative feature of the magnetic anisotropy, namely the magnetization direction, cannot be accounted for by rules of general character, not even in a case like Mn-TCNQ, where the -electrons have a rather localized character that would make this system seem a priori a good playground for these models.

Next, we turn our attention back to the case of Ni-TCNQ, where the DFT calculations yield a relatively large value for the MAE with out-of-plane magnetization, as well as significant variations of the MAE in the network plane. This theoretical result contrasts with the experimental absence of anisotropy in this system and, thus, requires a further analysis oriented at finding an explanation. As we discuss below, the discrepancy could be explained by substrate effects, mostly due to electronic charge transfer from the metal Au(111) surface. However, if we tried to calculate MAE values from DFT calculations with SOC using the supported Ni-TCNQ/Au(111) model structures presented above, we would not obtain informative results, since it would be very difficult in practice to disentangle the anisotropy effects originated by different aspects of the system. The most significant of them is the unavoidable artificial strain introduced in the system by forcing a commensurable Ni-TCNQ overlayer on top of the Au(111) surface due to the use of periodic boundary conditions in a finite size system imposed by our DFT calculations. However, these limitations can be more conveniently understood using free-standing models.

In the rectangular Ni-TCNQ free-standing overlayer we can attribute the large MAE values to the partially occupied Ni() states. If these bands were completely filled by transfer of 0.5 electrons from the metallic substrate, their contribution to the MAE would be dramatically reduced. Additionally, the Ni atom spin would become close to , a case for which no single-ion anisotropy is possible [28]. However, it is hard to give a precise estimate of the amount of charge transfer and, on top of this, other sources of anisotropy reduction could be at play, like a reduction of Ni coordination due to a geometrical distortion. Indeed, the lowest-energy configuration of this rectangular unit cell is obtained upon a small symmetry-lowering distortion where the four Ni-N bonds are inequivalent: the bonds at degrees with the OX axis (d) have a length of 1.95 Å and the other pair at (d) of 2.01 Å (see Fig. 1). The former direction is that of the hardest magnetization axis. This symmetry breaking, though subtle from the geometry point of view, is nevertheless associated to a noticeable asymmetry in the electronic structure, which is in turn behind the strong azimuthal variability of the MAE. In a closer inspection of the PDOS we find that the Ni() peaks at the Fermi level hybridized with the molecule LUMO are contributed by -orbitals lying on the plane containing the short Ni-N bonds (d) and the surface normal (see section III of the Supplementary Material). The long bonds (d), to which states at the Fermi level do not contribute , correspond to a softer magnetization direction.

To understand the consequences of this distorted geometry on the magnetic anisotropy, we have constructed a free-standing flat Ni-TCNQ model in an oblique unit cell, in which the angle between the lattice vectors and is varied (the rectangular cell corresponds to ). As described in section III of the Supplementary Material, two cases have been considered: a weakly distorted case with and a larger distortion with . The unit cell angle has been reduced while uniformly scaling the lattice constants to keep the unit cell area equal to that of the rectangular equilibrium unit cell. Then, the atomic coordinates have been relaxed to satisfy the same convergence criteria as in other models of the present work. For a larger distortion of the rectangular cell with , one could force a commensurate supercell on Au(111) [4]. In the optimized structure the TCNQ is barely deformed, but one Ni-N bond at the azimuthal direction is broken because of the cell distortion and the pair of bonds at the direction have their lengths reduced to 1.85 Å (see section III of the Supplementary Material). The magnetic anisotropy is significantly reduced with respect to that of the rectangular cell, but the hardest direction is still the one along the shortest pair or Ni-N bonds (see Table 2). The main consequence of the Ni coordination reduction caused by the cell shape change is to partially quench its spin. We observe that the local magnetic moment is reduced by about , approaching the ideal state that would yield no anisotropy in the single-atom picture. We observe, nevertheless, that this distorted configuration still has partially filled Ni states at the Fermi level (see section III of the Supplementary Material). Thefore, we note that this mechanism of anisotropy reduction and the charge transfer effect proposed above are of a different nature, albeit both originate from the interaction with the substrate.

All in all, the observed lack of magnetic anisotropy in the Ni-TCNQ/Au(111) XMCD data is clearly a substrate effect, which reduces the Ni-TCNQ anisotropy by a combined effect of charge transfer and change of coordination. Nonetheless, other subtle substrate effects not considered here might also have a role, such as fluctuations in the Au-Ni charge transfer due to the incommensurability and corrugation of the network.

## 4 Discussion and Conclusions

Motivated by the XMCD data, we have performed a thorough analysis of the magnetic properties that characterized Mn and Ni metal-organic coordination networks, focusing on the magnetic coupling and anisotropy. By fitting the XCMD data using a model Hamiltonian based on mean-field Weiss theory and comparing with Brillouin functions, we find a completely different behavior for Mn and Ni networks: while in Mn networks the spins localized at the Mn centers are coupled antiferromagnetically with a mild preference to in-plane magnetization, in Ni networks the spins localized at the Ni atoms are coupled ferro-magnetically and do not show any sizable magnetic anisotropy.

These observations are also rationalized with the help of density functional theory calculations in two steps: first we focus on the magnetic coupling and next we address the subtle question of the magnetic anisotropy. Spin-polarized DFT calculations using a surface unit cell to describe the free-standing-overlayers reveal a very different electronic structure close to the Fermi level for the two systems under study. The system Mn-TCNQ is insulating and has weak hybridization between Mn and TCNQ states close to the Fermi level, while in Ni-TCNQ hybridization between Ni(3d) states and the TCNQ LUMO at the Fermi level is rather significant. This difference permits to explain the observed trends in XMCD data with antiferromagnetic (ferromagnetic) coupling for Mn (Ni) networks that is also confirmed by another set of DFT calculations using a surface unit cell.

We find that the basic qualitative feature of the magnetic anisotropy, namely the magnetization direction, cannot be accounted for by rules of general character. Actually, the magneto-crystalline anisotropy is contributed by many electron excitation channels and it clearly shows an intricate dependence on the fine electronic structure details of each particular system. While in Mn-TCNQ/Au(111) the observed magnetic anisotropy with in-plane magnetization agrees with the DFT calculations for the neutral Mn-TCNQ overlayer, the observed lack of magnetic anisotropy in Ni-TCNQ/Au(111) suggests the existence of a substrate effect, which reduces the Ni-TCNQ anisotropy due to a combination of electronic charge transfer and change of Ni-N coordination.

Supplementary Materials: The following are available online at www.mdpi.com/link; I. Supplementary XAS and XMCD data for TCNQ and F4TCNQ networks, II. Energetics of different magnetic configurations using a 2x2 unit cell, III. MAE convergence details, dependence with U and with lattice distortions.

Acknowledgments: We thank MINECO and the University of the Basque Country (UPV/EHU) for partial financial support, grant numbers FIS2016-75862-P and IT-756-13, respectively, the first one covering costs to publish in open access journals. MMO acknowledges the support by the Tomsk State University competitiveness improvement programme (project No. 8.1.01.2017) and the Saint Petersburg State University grant for scientific investigation (No. 15.61.202.2015).

Author Contributions: M.B.R., A.S., M.M.O, V.N.G. and A.A contribute with all the theoretical calculations, analysis of the data and dicussion. S. S., C. N., L. P., C. S., C. P. and P. G. contribute with all the experimental data, their analysis and interpretation. All the authors have participated in the discussions and revision of the manuscript written by M.B.R., A.S., S.S., V.N.G. and A.A.

Conflicts of Interest: The authors declare no conflict of interest.

## A

### a.1 Minimization procedure to obtain the self-consistent mean field equations

The general procedure for the minimization of the free energy of the metal-TCNQ model outlined in Section 3.2 is described here. Given an interacting Hamiltonian and a variational Hamiltonian , the Bogoliubov upper bound for the free energy reads

(A1) |

where, by definition,

(A2) |

Using the mean-field Hamiltonian of Eq. (3) in the place of , we obtain the upper bound for the free energy, which needs subsequently to be minimized with respect to the order parameters .

To analyze the XMCD data taken at normal incidence, we consider a magnetic field applied along the -axis, . The paramagnetic partition function of a single isolated spin reads

(A3) |

The corresponding spin average value can be found in this case by differentiating the free energy , where .

Two order parameters and are needed to describe the (checkerboard) antiferromagnetic phase. They represent the -components of the spins depicted in Fig. 4 (a). The upper bound to the free energy reads (per unit cell):

with

(A5) |

The terms in the last two lines of Eq. (LABEL:eqFreeAF) come from the average in Eq. (A1). These terms are required only when looking for the global minimum of , which is carried out over the domain and . The values of and at the global minimum correspond then to the mean-field solution. Alternatively, we can use the stationarity condition

(A6) |

to obtain the two coupled equations (4).

These self-consistency equations of the mean-field theory need to be solved for and by substitution of the expressions for and from Eq. (A5). Since several solutions can be found, we use the condition of least free energy value to select the physical solution. In practice, it is convenient to find a rough approximation for and by looking for the global minimum of on a discrete grid and then refine the obtained solution by iteratively substituting it into the self-consistency equations (4).

## References

- Dong et al. [2016] Dong, L.; Gao, Z.; Lin, N. Self-assembly of metal–organic coordination structures on surfaces. Progress in Surface Science 2016, 91, 101 – 135.
- Wegner et al. [2008] Wegner, D.; Yamachika, R.; Wang, Y.; Brar, V.W.; Bartlett, B.M.; Long, J.R.; Crommie, M.F. Single-Molecule Charge Transfer and Bonding at an Organic/Inorganic Interface: Tetracyanoethylene on Noble Metals. Nano Letters 2008, 8, 131–135.
- Bedwani et al. [2008] Bedwani, S.; Wegner, D.; Crommie, M.F.; Rochefort, A. Strongly Reshaped Organic-Metal Interfaces: Tetracyanoethylene on Cu(100). Phys. Rev. Lett. 2008, 101, 216105.
- Faraggi et al. [2012] Faraggi, M.N.; Jiang, N.; Gonzalez-Lakunza, N.; Langner, A.; Stepanow, S.; Kern, K.; Arnau, A. Bonding and Charge Transfer in Metal–Organic Coordination Networks on Au(111) with Strong Acceptor Molecules. The Journal of Physical Chemistry C 2012, 116, 24558–24565.
- Umbach et al. [2012] Umbach, T.R.; Bernien, M.; Hermanns, C.F.; Krüger, A.; Sessi, V.; Fernandez-Torrente, I.; Stoll, P.; Pascual, J.I.; Franke, K.J.; Kuch, W. Ferromagnetic Coupling of Mononuclear Fe Centers in a Self-Assembled Metal-Organic Network on Au(111). Phys. Rev. Lett. 2012, 109, 267207.
- Abdurakhmanova et al. [2012] Abdurakhmanova, N.; Floris, A.; Tseng, T.C.; Comisso, A.; Stepanow, S.; De Vita, A.; Kern, K. Stereoselectivity and electrostatics in charge-transfer Mn- and Cs-TCNQ4 networks on Ag(100). Nature Communications 2012, 3, 940 EP –.
- Abdurakhmanova et al. [2013] Abdurakhmanova, N.; Tseng, T.C.; Langner, A.; Kley, C.S.; Sessi, V.; Stepanow, S.; Kern, K. Superexchange-Mediated Ferromagnetic Coupling in Two-Dimensional Ni-TCNQ Networks on Metal Surfaces. Phys. Rev. Lett. 2013, 110, 027202.
- Giovanelli et al. [2014] Giovanelli, L.; Savoyant, A.; Abel, M.; Maccherozzi, F.; Ksari, Y.; Koudia, M.; Hayn, R.; Choueikani, F.; Otero, E.; Ohresser, P.; Themlin, J.M.; Dhesi, S.S.; Clair, S. Magnetic Coupling and Single-Ion Anisotropy in Surface-Supported Mn-Based Metal–Organic Networks. The Journal of Physical Chemistry C 2014, 118, 11738–11744.
- Bebensee et al. [2014] Bebensee, F.; Svane, K.; Bombis, C.; Masini, F.; Klyatskaya, S.; Besenbacher, F.; Ruben, M.; Hammer, B.; Linderoth, T.R. A Surface Coordination Network Based on Copper Adatom Trimers. Angewandte Chemie International Edition 2014, 53, 12955–12959.
- Rodríguez-Fernández et al. [2015] Rodríguez-Fernández, J.; Lauwaet, K.; Ángeles Herranz, M.; Martín, N.; Gallego, J.M.; Miranda, R.; Otero, R. Temperature-controlled metal/ligand stoichiometric ratio in Ag-TCNE coordination networks. The Journal of Chemical Physics 2015, 142, 101930.
- Faraggi et al. [2015] Faraggi, M.N.; Golovach, V.N.; Stepanow, S.; Tseng, T.C.; Abdurakhmanova, N.; Kley, C.S.; Langner, A.; Sessi, V.; Kern, K.; Arnau, A. Modeling Ferro- and Antiferromagnetic Interactions in Metal–Organic Coordination Networks. The Journal of Physical Chemistry C 2015, 119, 547–555.
- Otero et al. [2017] Otero, R.; de Parga, A.V.; Gallego, J. Electronic, structural and chemical effects of charge-transfer at organic/inorganic interfaces. Surface Science Reports 2017, 72, 105 – 145.
- Shi et al. [2010] Shi, X.Q.; Lin, C.; Minot, C.; Tseng, T.C.; Tait, S.L.; Lin, N.; Zhang, R.Q.; Kern, K.; Cerdá, J.I.; Van Hove, M.A. Structural Analysis and Electronic Properties of Negatively Charged TCNQ: 2D Networks of (TCNQ)2Mn Assembled on Cu(100). The Journal of Physical Chemistry C 2010, 114, 17197–17204.
- Mabrouk and Hayn [2015] Mabrouk, M.; Hayn, R. Magnetic moment formation in metal-organic monolayers. Phys. Rev. B 2015, 92, 184424.
- Yosida [1996] Yosida, K. Theory of Magnetism; Springer-Verlag Berlin Heidelberg, 1996.
- Abanin and Pesin [2011] Abanin, D.A.; Pesin, D.A. Ordering of Magnetic Impurities and Tunable Electronic Properties of Topological Insulators. Phys. Rev. Lett. 2011, 106, 136802.
- Caputo et al. [2016] Caputo, M.; Panighel, M.; Lisi, S.; Khalil, L.; Santo, G.D.; Papalazarou, E.; Hruban, A.; Konczykowski, M.; Krusin-Elbaum, L.; Aliev, Z.S.; others. Manipulating the topological interface by molecular adsorbates: adsorption of Co-phthalocyanine on Bi2Se3. Nano letters 2016, 16, 3409–3414.
- Girovsky et al. [2017] Girovsky, J.; Nowakowski, J.; Ali, M.E.; Baljozovic, M.; Rossmann, H.R.; Nijs, T.; Aeby, E.A.; Nowakowska, S.; Siewert, D.; Srivastava, G.; Wäckerlin, C.; Dreiser, J.; Decurtins, S.; Liu, S.X.; Oppeneer, P.M.; Jung, T.A.; Ballav, N. Long-range ferrimagnetic order in a two-dimensional supramolecular Kondo lattice. Nature Communications 2017, 8, 15388 EP –.
- Kresse and Hafner [1993] Kresse, G.; Hafner, J. Ab initio. Phys. Rev. B 1993, 47, 558–561.
- Kresse and Furthmüller [1996a] Kresse, G.; Furthmüller, J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B 1996, 54, 11169–11186.
- Kresse and Furthmüller [1996b] Kresse, G.; Furthmüller, J. Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set. Computational Materials Science 1996, 6, 15 – 50.
- Perdew et al. [1997] Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple [Phys. Rev. Lett. 77, 3865 (1996)]. Phys. Rev. Lett. 1997, 78, 1396–1396.
- Dudarev et al. [1998] Dudarev, S.L.; Botton, G.A.; Savrasov, S.Y.; Humphreys, C.J.; Sutton, A.P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. Phys. Rev. B 1998, 57, 1505–1509.
- Monkhorst and Pack [1976] Monkhorst, H.J.; Pack, J.D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188–5192.
- Methfessel and Paxton [1989] Methfessel, M.; Paxton, A.T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B 1989, 40, 3616–3621.
- Blöchl et al. [1994] Blöchl, P.E.; Jepsen, O.; Andersen, O.K. Improved tetrahedron method for Brillouin-zone integrations. Phys. Rev. B 1994, 49, 16223–16233.
- Gambardella et al. [2009] Gambardella, P.; Stepanow, S.; Dmitriev, A.; Honolka, J.; de Groot, F.M.F.; Lingenfelder, M.; Gupta, S.S.; Sarma, D.D.; Bencok, P.; Stanescu, S.; Clair, S.; Pons, S.; Lin, N.; Seitsonen, A.P.; Brune, H.; Barth, J.V.; Kern, K. Supramolecular control of the magnetic anisotropy in two-dimensional high-spin Fe arrays at a metal interface. Nature Materials 2009, 8, 189 EP –.
- Dai et al. [2008] Dai, D.; Xiang, H.; Whangbo, M.H. Effects of spin-orbit coupling on magnetic properties of discrete and extended magnetic systems. Journal of Computational Chemistry 2008, 29, 2187–2209.
- Fan and Wu [1969] Fan, C.; Wu, F.Y. Ising Model with Second-Neighbor Interaction. I. Some Exact Results and an Approximate Solution. Phys. Rev. 1969, 179, 560–569.
- Chakravarty et al. [1989] Chakravarty, S.; Halperin, B.I.; Nelson, D.R. Two-dimensional quantum Heisenberg antiferromagnet at low temperatures. Phys. Rev. B 1989, 39, 2344–2371.
- Otrokov et al. [2015] Otrokov, M.M.; Chulkov, E.V.; Arnau, A. Breaking time-reversal symmetry at the topological insulator surface by metal-organic coordination networks. Phys. Rev. B 2015, 92, 165309.
- Bruno [1989] Bruno, P. Tight-binding approach to the orbital magnetic moment and magnetocrystalline anisotropy of transition-metal monolayers. Phys. Rev. B 1989, 39, 865–868.
- van der Laan [1998] van der Laan, G. Microscopic origin of magnetocrystalline anisotropy in transition metal thin films. Journal of Physics: Condensed Matter 1998, 10, 3239.
- Stöhr [1999] Stöhr, J. Exploring the microscopic origin of magnetic anisotropies with X-ray magnetic circular dichroism (XMCD) spectroscopy. Journal of Magnetism and Magnetic Materials 1999, 200, 470 – 497.
- Li et al. [1990] Li, C.; Freeman, A.J.; Jansen, H.J.F.; Fu, C.L. Magnetic anisotropy in low-dimensional ferromagnetic systems: Fe monolayers on Ag(001), Au(001), and Pd(001) substrates. Phys. Rev. B 1990, 42, 5433–5442.
- Abdelouahed and Alouani [2009] Abdelouahed, S.; Alouani, M. Magnetic anisotropy in Gd, GdN, and tuned by the energy of gadolinium states. Phys. Rev. B 2009, 79, 054406.
- Ayaz Khan et al. [2016] Ayaz Khan, S.; Blaha, P.; Ebert, H.; Minár, J.; Šipr, O.c.v. Magnetocrystalline anisotropy of FePt: A detailed view. Phys. Rev. B 2016, 94, 144436.
- Otrokov et al. [2017] Otrokov, M.M.; Menshchikova, T.V.; Vergniory, M.G.; Rusinov, I.P.; Vyazovskaya, A.Y.; Koroteev, Y.M.; Bihlmayer, G.; Ernst, A.; Echenique, P.M.; Arnau, A.; Chulkov, E.V. Highly-ordered wide bandgap materials for quantized anomalous Hall and magnetoelectric effects. 2D Materials 2017, 4, 025082.
- Ke and van Schilfgaarde [2015] Ke, L.; van Schilfgaarde, M. Band-filling effect on magnetic anisotropy using a Green’s function method. Phys. Rev. B 2015, 92, 014423. © 2019 by the authors. Submitted to MDPI for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).