# Modelling the magnetic properties of 1D arrays of FePc molecules

R. E. Pico  Instituto de Física de Rosario, Consejo Nacional de
Investigaciones Científicas y Técnicas (CONICET) and Universidad Nacional de
Rosario, Bv. 27 de Febrero 210 Bis (2000) Rosario, Argentina  pico@ifir-
conicet.gov.ar A. F. Rebola  Instituto de Física de Rosario, Consejo Nacional
de Investigaciones Científicas y Técnicas (CONICET) and Universidad Nacional
de Rosario, Bv. 27 de Febrero 210 Bis (2000) Rosario, Argentina  J. Lasave
Instituto de Física de Rosario, Consejo Nacional de Investigaciones
Científicas y Técnicas (CONICET) and Universidad Nacional de Rosario, Bv. 27
de Febrero 210 Bis (2000) Rosario, Argentina  P. Abufager  Instituto de Física
de Rosario, Consejo Nacional de Investigaciones Científicas y Técnicas
(CONICET) and Universidad Nacional de Rosario, Bv. 27 de Febrero 210 Bis
(2000) Rosario, Argentina  I. J. Hamad  Instituto de Física de Rosario,
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) and
Universidad Nacional de Rosario, Bv. 27 de Febrero 210 Bis (2000) Rosario,
Argentina<EMAIL_ADDRESS>
###### Abstract

We investigate the magnetic properties of Fe Phthalocyanines (FePc) that are
experimentally arranged in quasi one-dimensional chains when they are grown in
thin films or powders. By means of DFT calculations we reproduce the
structural parameters found in experiments, and then we build a generalized
Heisenberg magnetic model with single ion anisotropy, and calculate its
parameters. The results show a anisotropic exchange interaction $J$ between
FePc molecules, and an easy plane single ion anisotropy $D$. By means of Monte
Carlo simulations, with this model, we found an explanation to the non-
saturation of the magnetization found at high fields, which we interpret is
due to the anisotropic exchange interaction $J$. Finally, we also investigate
the presence of magnetic solitons versus temperature and magnetic field. This
results provide additional evidence that FePc is a soliton bearing molecular
compound, with solitons easily excited mainly in the molecular $xy$ plane.

##  1 Introduction

Over the past decade, significant progress has been made in the field of
organic materials due to their potential applications as nanoscale spintronics
1, organic photovoltaics 2, 3, transistors 4, 5, 6, light-emitting diodes 7,
8, and electrochemical devices 9, 10. Organic materials constructed with one-
dimensional building blocks, in particular, have become a subject of great
interest because they serve as excellent model systems for studying and
applying low dimensional magnetic phenomena. It is important, then, to
understand the underlying principles of their magnetic properties. Among the
building blocks for organic materials, metal phthalocyanines (MPcs) are one of
the most extensively used 11.

The MPc molecule is composed of a conjugated planar macrocycle, MN8C3H16, with
M being a central metal atom that is 4-coordinated to the N isoindole or
pyrrole N atoms, Npy. The macrocycle is made up of four isoindole units linked
by a ring of four N atoms (Nbridge), as shown in Fig.1(a). Depending on the
growth conditions, MPc crystallizes into a range of polymorphs, that have been
observed in powders, thin films (TF), and nanowires systems 11. The building
blocks of these crystalline forms are one dimensional (1D) columns, composed
of MPc molecules, with the metal atom M at the center serving as the column
axis, or the stack axis (see Fig.1(d)). The stacking columnar structures
differ in the tilting angles ($\varphi$) between the stacking axis and the
molecular plane, and the interplanar distance ($d$) between two neighboring
molecular planes along the column, resulting in different Pc crystal
structures (see Fig.1(c)) 12, 11, 13, 14, 15, 16, 17. Modifying the magnetic
center of MPc introduces an additional degree of freedom to the properties of
its polymorphs, since it impacts the magnetic properties of the molecule and
hence of the crystals formed by them.

Two different arrangements have been proposed to describe the stacking
geometry of MPc molecules. The first one is referred to as the herringbone
structure, where the molecules stack with opposite angles between adjacent
columns 18. In the second one, the MPc molecules stack in the same orientation
in each column, forming the so-called brickstack structure 13.

Particularly, iron phthalocyanines (FePc) are a class of MPc molecules that
have garnered significant attention. In bulk arrangements, FePc exhibits two
polymorphs, $\alpha$ and $\beta$, with stacking angles of 26.5∘ and 44.8∘,
respectively. Previous studies of the $\beta$ phase have shown that it behaves
paramagnetically 19. However, in the $\alpha$ phase found in powder systems,
FePc shows ferromagnetic correlations below approximately 10 K 19. The
evidence in powders points to an arrangement of the herringbone type 20. On
the other hand, in thin films (TF), FePc molecules adopt a brickstack
structure when grown on a variety of substrates 17, 12, 14, 21. The
polymorphism found in TF has been called $\alpha_{+}$ 13 (but it has also been
directly called $\alpha$ 14 or TF 17), since, besides the brickstack
arrangement, it is different from the bulk $\alpha$ phase from the way each
molecule is shifted with respect to the one in the upper/lower molecular plane
(see section Stacking Geometry below for more detail). In this work we will
refer to this polymorphism as $\alpha_{+}$. FePc thin films have been reported
to become ferromagnetically correlated at temperatures below 20 K 14 and
exhibit a ferromagnetic transition at $4.5$ K 22, very similar to the
transition temperature in powders at $5$ K 12.

In order to explain the magnetic phenomena observed in experiments, both in
powders and in TF, different model Hamiltonians have been proposed. Moreover,
even the same parameters have very different proposed values in different
works (see section Magnetic Model below for more details). On the experimental
side, magnetization measurements made in powders and also also in TF show that
at relatively high magnetic fields, the magnetization does not saturate 19,
14. In the case of powders, this non saturation has been attributed to the
herringbone arrangement, due to the canting of the adjacent columns of FePc
molecules 19. However, in thin films the arrangement is of the brickstack
type20, so that other mechanism must be also contributing to the non
saturation of the magnetization as well. Finally, in a recent work,
experimental evidence of magnetic solitons has been found 14. Since FePc has
gained attention partly because of their electrical and magnetic properties,
it is important to study how magnetic solitons are developed in this system.
In order to shed light on these issues, a theoretical study beginning from
first principles calculations ending in a magnetic model was carried out, to
help to understand this and related systems further.

Figure 1:  (a) FePc molecule view from top. The four Nitrogen closer and
farthest to the Fe atom are the pyrrole (Npy) and the bridging (Nbridge)
nitrogen atoms, respectively. Fe, Npy, Nbridge, C and H are indicated by red,
cyan, blue, gray and white balls, respectively. (b) and (c) Stacking geometry
of FePc-film. $d$ is the distance between the molecular planes and $\varphi$
is the angle between the stacking axis $\bar{u}_{3}$ and a vector
perpendicular to the molecular plane. (d) The cell parameters for a supercell
containing 1 molecule are $|\overline{u}_{1}|=$13.24 Å,
$|\overline{u}_{2}|=$13.43 Å, $|\overline{u}_{3}|=$3.60 Å,
$(\overline{u}_{1},\overline{u}_{2})$=85∘,
$(\overline{u}_{1},\overline{u}_{3})$=63.6∘,
$(\overline{u}_{2},\overline{u}_{3})=$74.1∘. It is important to take into
account that for the magnetic model, we employed a cell size that is doubled
in the $\overline{u}_{3}$ direction and contains two molecules per cell.

In this work we theoretically study the $\alpha_{+}$ phase of FePc chains, by
means of density functional theory (DFT) and classical Monte Carlo. With DFT
calculations, we first obtain structural parameters, comparing with available
experimental data, and then proceed to calculate the parameters of a
generalized Heisenberg model with single ion anisotropy. We obtain a
anisotropic Heisenberg model, with an easy plane single ion magnetic
anisotropy $D$ and negligible $E$. Afterwards, we study this model in an array
of weakly coupled FePc chains, by means of classical Monte Carlo simulations.
We were able to explain the lack of saturation observed in magnetization
experiments. Finally, we investigate the magnetic excitations of this system,
focusing on exploring the existence of solitons in the phase diagram.

##  2 Methods

###  2.1 Computational Details

DFT calculations were performed with the VASP code 23, 24, 25, 26, 27, 28
within the slab-supercell approach and using the projector augmented-wave
(PAW) method 27.Wave functions were expanded using a plane wave basis set with
an energy cut-off of 550 eV. Different schemes were tested in order to improve
the description of vdW interaction: the PBE+D3 29 approach and the
optB86b-vdW, optB88-vdW, optPBE-vdW, vdW-DF2 non-local correlation functionals
30. From test calculations for the bulk-like FePc structure, we conclude that
the optb86b scheme is the one that gives the best interplanar distance with
respect to the experimental value, therefore, it is used for the structural
determination of bulk structures. However, when calculating the parameters of
the magnetic model, in particular the anisotropy parameters $D$ and $E$, which
require spin-orbit interaction turned on, opt86b is not suitable, and the
PBE+D3 scheme is used instead (see section magnetic model). Exchange coupling
parameters $J_{x}$, $J_{y}$ and $J_{z}$ are also calculated within this PBE+D3
scheme. For testing purposes, in an isotropic Heisenberg model the exchange
coupling $J$ has no appreciable difference when calculated with PBE$+$D3 or
with optb86b.

The Hubbard U correction is considered with the DFT+U approximation 31 to deal
at DFT level with the Fe d-electrons. $U_{eff}=U-J=3eV$ was chosen, since this
value has also been used in previous studies in FePc molecules32, 33.

The FePc $\alpha_{+}$ polymorphism was studied. For testing purposes, results
for the gas phase of FePc are also presented here. Gas phase calculations were
carried out at $\Gamma$ point. For the $\alpha_{+}$ phase, the k-point grid
was 2x2x15. Calculation without and with spin-orbit coupling (SOC) have been
performed. Calculations with SOC were carried out at fixed geometry. In all
cases, $k$ point convergences were achieved. All geometry optimizations were
carried out until the forces on every mobile atom were smaller than 0.02 eV/Å,
and all calculations were spin polarized.

###  2.2 Gas Phase calculations

Despite the numerous theoretical and experimental studies there is still no
definite consensus about the ground state of the molecule. However, there is a
broad consensus that the FePc molecule in gas phase has spin $S=1$. DFT
studies have predicted ${}^{3}A_{2g}$ [$(d_{z^{2}}^{\uparrow\downarrow})\
(d_{xz,yz}^{\uparrow\uparrow})\ (d_{xy}^{\uparrow\downarrow})$] 34 ,
${}^{3}E_{g}(a)$ [$(d_{z^{2}}^{\uparrow})\
(d_{xz,yz}^{\uparrow\downarrow\uparrow})\ (d_{xy}^{\uparrow\downarrow})$] 35,
36 and ${}^{3}B_{2g}$ [$(d_{z^{2}}^{\uparrow})\
(d_{xz,yz}^{\uparrow\downarrow\uparrow\downarrow})\ (d_{xy}^{\uparrow})$] 37
as possible ground states and the prediction actually depends on the choice of
the exchange-correlation (XC) functionals 38. The same conclusion has been
established by a recent work based on a diffusion Monte Carlo study where they
obtain mainly ${}^{3}A_{2g}$, but also ${}^{3}B_{2g}$ and even ${}^{3}E_{g}$
depending on XC 39. A recent theoretical work employing numerous and
sofisticated methodologies has concluded that the ground state of the molecule
is ${}^{3}A_{2g}$ 40.

Within our calculation, we obtained ${}^{3}E_{g}(a)$ as the ground state (GS),
a ${}^{3}A_{2g}$ state at 36 meV and ${}^{3}B_{2g}$ at 279 meV higher energy
in line with the results obtained by Ichiba et al. 39 for DFT+U calculations
with U$=$4eV, and also with Ref 41. However, as mentioned earlier, in the
isolated molecule, numerous configurations appear to be closely competitive.
For instance, in our case, the energy difference between ${}^{3}E_{g}(a)$ and
${}^{3}A_{2g}$ is merely 36 meV. Sophisticated methods like those of Refs. 39,
40 need to be employed to obtain a definite answer to the GS of the isolated
molecule. Our focus will be on FePc chains on TF, where the situation seems to
be different (see below). Finally, for all electronic configurations, the spin
magnetic moment of the molecule is close to 2$\mu_{b}$ concentrated in the Fe
center.

Figure 2: Directions were the total energy was computed in DFT calculations
including SOC. Perpendicular to the molecular plane (001) and two in plane
orthogonal directions (100 and 010), respectively.

###  2.3 $\alpha_{+}$ phase. Stacking Geometry

In FePc TF grown on different substrates, molecules stack layer by layer
almost parallel (with a tilting angle of $7$ or $9$ degrees) to the substrate
surface 17, 14, 21. In this ”flat-lying” configuration, two adjacent molecular
layers are shifted a distance $\Delta$ (see Fig. 1(c)) apart from one another
along the Fe-Npy bond direction indicated by the $a$, or alternatively the
$b$, axis in Fig.1(a). As a result, the Fe atom in a molecule is close to the
center of a Npy atom in the molecule beneath it 17, 14 (see Fig.1(b)). In this
stacking geometry Fe atoms form parallel unidimensional chains that deviate
from the normal to the molecular plane by an angle $\varphi$ (see. Fig.1(c)).
In the case of the herringbone arrangement found in powders, apart from the
two-sublattice canted structure the other difference is that the shift between
FePc molecules occurs along the Fe-Nbridge bond 17. In the case of a FePc TF
grown on Au(111) 17, a structural refinement based on XLPA simulations gives
$\Delta$=1.48 Å, which is very close to the estimated value from STM images,
$\Delta$=(1.0 $\pm$ 0.3) Å. The distance between adjacent molecular planes is
$d_{z}$ = 3.25 Å, and molecules within the molecular plane form a square
lattice with a lattice parameter of 13.0 $\pm$ 0.2 Å17 (see Fig.1(d)).

The estimated experimental values for FePc TF grown on PTCDA/Si(100) 14 are
$d_{z}$ = 3.42 Å and a distance $d_{Fe-Fe}=$ 3.77 Å between the nearest Fe
atoms within a chain. In the case of FePc TF grown Au-coated sapphire,
$d_{z}$= 3 Å and the lateral separation between chains is about 13 Å21.
Similar values, with slightly varying, thickness-dependent, stacking angles,
have been found in a recent study of FePc films in different substrates 42.
The similarities in structural data reported for FePc-TF in the ”flat-lying”
configuration suggest that once layer-by-layer growth is initiated due to the
templating effect of the substrate 15, intermolecular interactions seem to
drive the film structure. This leads to performing DFT structural calculations
without including the substrate.

We have performed DFT calculations for bulk-like structures with one molecule
per unit cell, simultaneously optimizing the unit cell and the molecular
degrees of freedom, with the aim of reproducing to the most FePc-TF structure.
For the obtained lowest energy structure, the shift of two adjacent layers
runs along Fe-Npy bond axis111The deviation is less than 2∘, with a distance
of $\Delta=1.61$ Å. The distance between molecular planes, the stacking angle
of the Fe chain and Fe-Fe distance within a chain are, respectively,
$d_{z}=3.22$ Å, $\varphi=26.6^{\circ}$ and $d_{Fe-Fe}=3.60$ $\AA$.
Furthermore, the dimensions of the unit cell vectors u1 and u2 and the angle
among them are 13.2 Å, 13.4 Å and 85∘, respectively. Because the molecular
plane is close to the defined plane by u1 and u2 222The angular deviation is
$\sim$4∘ their dimensions can be compared to the lateral distance between
chains experimentally measured. Overall, all parameters agree well with the
previously described experimental data, lending support to the DFT
calculations for the $\alpha_{+}$ phase. For this stacking geometry, the
calculations show that each molecule is in a $E_{g}(a)$ configuration
($(d_{z^{2}}^{\uparrow})\ (d_{xz,yz}^{\uparrow\downarrow\uparrow})\
(d_{xy}^{\uparrow\downarrow})$). This result coincides with previous
constrained DFT calculations 41 and with experimental evidence 43, 20, altough
other works predict a ground state in TF that is a superposition of
${}^{3}A_{2g}$,${}^{3}B_{2g}$ and ${}^{3}E_{g}$ 44.

Finally, the calculated spin magnetic moment in the Fe atom is
$m_{s}(Fe)=2.021\mu_{b}$, corresponding to $S=1$, the same than in the gas
phase.

##  3 Results and discussion

###  3.1 Magnetic model

On the basis of experimental evidence, several magnetic models have been
proposed for the bulk phase $\alpha$ and for the TF $\alpha_{+}$ of FePc.
There is consensus that the molecule has $S=1$, and also that there is a
positive single ion magnetic anisotropy $D$, aside from the ferromagnetic
exchange $J$. Working with powders of FePc molecules, Evangelisti et. al19
observed that the magnetization did not saturate even at high magnetic fields,
and this was interpreted in terms of the above mentioned herringbone
arrangement, a canted structure between two different magnetic sublattices
originally proposed for $MnPc$. An effective Ising $S=1/2$ model was proposed.
This Ising model was later backed up in Ref. 20. However, it was not
considered as appropriate to account for the observed behaviour in TF in Ref.
14, and an isotropic Heisenberg model with single ion anisotropy was proposed
instead. As for the mechanism for the magnetic exchange, in Ref. 19 the direct
Fe-Fe interaction was proposed, which renders the interaction ferromagnetic in
character. However, the value proposed for the exchange parameter $J/k_{B}$ or
$J_{z}/k_{B}$ varies significantly. In a Ising model written as
$-2J_{z}\sum_{i}S_{z,i}S_{z,i+1}$ for powder samples, it was estimated to have
the value of $25.7K$ in Ref. 19 and $76K$ in Ref. 20 . In isotropic Heisenberg
models $-J\sum_{i}S_{i}S_{i+1}$, $J/k_{B}$ was estimated to be $20K$ in Ref.
14 and $13K$ in Ref. 42, these last two working on TF. On the other hand,
while there is a consensus that the magnetic anisotropy $D$ is positive,
indicative of an easy plane situation, its value has been estimated from
experiments, varying from $D=92$ K in Ref. 45 to $D=53K$ in Ref. 14.

Spin orbit interaction energy is important in FePc molecules, since they
present a high unquenched orbital moment 43, 20. In this case, it is well
known that an effective spin hamiltonian can be derived, which contains terms
of the form $DS_{z}^{2}+E(S_{x}^{2}+S_{y}^{2})$, being the parameters $D$ and
$E$ matrix elements proportional to the orbital magnetic moment. The $D_{4h}$
symmetry present in isolated FePc molecules is broken when the chains are
formed in the bulk or thin film structure, since the stacking axis is not
perpendicular to the molecular plane. This suggests the idea that the
Heisenberg exchange term may in principle be anisotropic. An extreme case of
this is precisely the Ising model considered in Ref. 19. However, most recent
works in FePc bulk structures or TF have considered the isotropic Heisenberg
model only 14, 21, 46. On the other hand, the interchain exchange is expected
to be relatively small, due to relatively large lateral distances between Fe
centers, as confirmed by calculations done in CrPc chains 47. These
considerations render FePc thin films to be treated, as DFT calculations are
concerned, as one dimensional chains of FePc molecules, geometrically arranged
as described in section Stacking Geometry, with the possibility of anisotropic
magnetic exchange and with other spin-orbit derived interactions present, such
as single ion magnetic anisotropy terms. On the other hand, regarding the
Dzyaloshinskii-Moriya (DM) anisotropic exchange, it’s important to consider
that its order of magnitude is given by $D_{DM}\sim(\Delta(g)/g)J$ 48. With
the most recent estimation of $g$ in TF, $g=2.16$ 14 the magnitude can be
estimated in $D_{DM}\sim 0.08J$. Consequently, this interaction has not been
considered significant in this work.

With all this considered, the proposed model Hamiltonian for the one
dimensional chains of FePc molecules in the $\alpha_{+}$ polymorphism is:

$\begin{split}H&=\sum_{i}D(S_{i}^{z})^{2}+\sum_{i}E(S_{ix}^{2}-S_{iy}^{2})\\\
&-\frac{1}{2}\sum_{<i,j>}J_{z}(S_{i}^{z}\cdot
S_{j}^{z})-\frac{1}{2}\sum_{<i,j>}J_{x}(S_{i}^{x}\cdot S_{j}^{x})\\\
&-\frac{1}{2}\sum_{<i,j>}J_{y}(S_{i}^{y}\cdot S_{j}^{y})\hskip 56.9055pt\\\
\end{split}$ (1)

where the sum over i runs over all Fe atoms in the supercell and the notation
$<i,j>$ represents first neighbor atoms. The factor $1/2$ is introduced to
avoid double counting. The first term in the Hamiltonian describes the
uniaxial anisotropy as we choose z as the off-plane direction. The second term
describes the in-plane anisotropy. The anisotropic exchange interaction is
considered in the last three terms.

We now proceed to compute the magnitude of the parameters of the model
Hamiltonian with DFT calculations, by means of the widely used method of
energy mapping 49, 50. To do so, the energy difference between classical
configurations, described in Ec. 2, has been computed.

$\begin{split}&E_{x}^{FM}=2ES^{2}-2J_{x}S^{2}+E_{0}\\\
&E_{x}^{AF}=2ES^{2}+2J_{x}S^{2}+E_{0}\\\
&E_{y}^{FM}=-2ES^{2}-2J_{y}S^{2}+E_{0}\\\
&E_{y}^{AF}=-2ES^{2}+2J_{y}S^{2}+E_{0}\\\
&E_{z}^{FM}=2DS^{2}-2J_{z}S^{2}+E_{0}\\\
&E_{z}^{AF}=2DS^{2}+2J_{z}S^{2}+E_{0}\\\ \end{split}$ (2)

Such calculations were performed for a chain defined by two molecules per unit
cell. For each configuration the directions x, y, z are defined as in Figure
2. In each line of Eq. 2, $E_{direction}^{FM/AFM}$ is defined, where $FM$
means that the magnetic moments of both molecules are parallel and $AFM$,
antiparallel. For determining the parameters $J_{x}$, $J_{y}$, $J_{z}$, $D$
and $E$ we used the ground state energy obtained from our DFT calculations for
each configuration in equation 2. Then, the parameters expressed as function
of the ground states energies were

$J_{i}=\frac{E_{i}^{AF}-E_{i}^{FM}}{4S^{2}}$ (3)
$D=\frac{1}{2}\big{[}E_{z}^{FM}-\frac{E_{x}^{FM}+E_{y}^{FM}}{2}+2(J_{z}-\frac{J_{x}+J_{y}}{2})\big{]}$
(4) $E=\frac{1}{4}\big{[}E_{x}^{FM}-E_{y}^{FM}+2(J_{x}-J_{y})\big{]}$ (5)

where the index $i$ runs for $x,y,z$.

The obtained results for the Heisenberg exchange were $J_{x}=1.77$ meV,
$J_{y}=2.31$ meV and $J_{z}=2.07$ meV. The precision of the DFT calculations
in VASP was more than enough to assure the last digit. While the mean value of
the exchange is very similar to the proposed $J$ values in isotropic models,
our calculations show that it is anisotropic, with an important amount of
exchange anisotropy $J_{y}/J_{x}=1.31$. It is important to note that the
positive value of the Heisenberg exchange means a ferromagnetic coupling, in
line with previous experimental results 19, 14. The type and value of the
exchange coupling in MPCs is strongly influeced by the stacking angle, hence
by the $\Delta$ value47.

The calculated single ion anisotropy is $D=0.65$ meV, which is equivalent to
$\sim 7.5K$, a value much lower than the estimations from experiments
previously mentioned. Nevertheless, it is important to note that the
calculated value agrees with earlier theoretical findings in which the
preferred in plane magnetization direction was also identified 41. The easy
plane anisotropy obtained is due to the spin-orbit interaction trying to align
the spin with the orbital magnetic moment, this being the result of the
orbital motions of the Fe electrons. Our results show an orbital magnetic
moment greater in the molecular plane, with values of 0.118, 0.140, and 0.039
$\mu_{B}$ along the $x$, $y$, and $z$ directions, respectively. However, these
values are likely underestimated, as suggested by previous research43. DFT
calculations tend to underestimate the orbital magnetic moment when electron
correlation effects are important 51, 52, and it is expected that the actual
value of $D$ is higher than in our calculations. Finally, the calculations
showed than $E$ is very small compared to the other parameters of the model,
and consequently it will be disregarded from now on.

To this one-dimensional model for FePc chains we add interchain interactions
$J_{inter}$ that are expected to be one order of magnitude smaller than the
intrachain ones, due to the large interchain Fe-Fe distance compared with the
intrachain one 14, 47. Although the DFT calculations were made for the
particular geometry of the $\alpha+$ phase found in TF, we expect the model to
account for the physics of other polymorphisms as well. Our results are
summarized in table 1.

$\alpha_{+}$ | $J_{x}$ | $J_{y}$ | $J_{z}$ | D | E
---|---|---|---|---|---
(meV) |  1.77 |  2.31 |  2.07 |  0.65 |  0.04
(K) |  20.5 |  26.8 |  24.0 |  7.5 |  0.5
Table 1: DFT calculated parameters for the magnetic model of FePc TF
$\alpha_{+}$ phase.

###  3.2 Monte Carlo calculations.

Taking into account the ferromagnetic ground state of weakly coupled FePc
chains, with individual molecules having spin $S=1$, classical Monte Carlo is
a suitable method to study its magnetic properties. All the Monte Carlo
simulations presented in this work were performed on periodic FePc chains
within a $3$x$3$x$1000$ supercell. In such an arrange each spin, representing
the Fe magnetic moment, interacts with 2 nearest neighbors within its chain
and with 4 nearest neighbors belonging to different chains located along the
directions of the x and y axes. The model Hamiltonian is the one described in
Eq. 1, with the derived parameters: $J_{1x}=1.76$ meV, $J_{1y}=2.31$ meV,
$J_{1z}=2.07$ meV (equivalent to $20.5,27$ and $24$ Kelvin, respectively),
magnetic anisotropy $D=0.65$ meV with the addition of interchain interactions
of $J_{inter}$, which were shown to be relatively small but not negligible 14,
47. The chosen value is $J_{inter}=0.01$ meV, since for this value the
transition temperature (i.e. the temperature at which there is a peak in
susceptibility or specific heat curves, not shown) is of around 5 Kelvin, a
value very similar to that reported in experiments.

###  3.3 Magnetization.

In magnetization experiments both in powders and thin films, the results
showed that even at high fields of several Tesla, the magnetization per
molecule did not saturate to the expected value for $S=1$, that is
$2g\mu_{B}$. In FePc powders, this was interpreted in terms of the herringbone
structure (see above), due to the canting of the sublattices 19. On the other
hand, in thin films, there is evidence that the structure is of the brickstack
type 14, so that no canted two sublattice structure is present. All the same,
the magnetization in thin films does not saturate at high fields. This
suggests that another mechanism could be playing a role in this observed
behaviour.

Concretely in Ref. 14, working with TF, it is shown that even at the lowest
temperatures reached, of $2K$:

  * •

For the field in the molecular plane (Fig. 4 d, e, templated films), the
magnetization did not saturate for fields as high as 7 Tesla, reaching
$1.8\mu_{B}$ at this value of the field.

  * •

For the field perpendicular to the molecular plane (Fig. 4 a, b, non-templated
films), the magnetization is lower than that corresponding to the same value
of $B$ for $B$ in the molecular plane ($M<1\mu_{B}$ for $B<6$ Tesla).

In the experiments in thin films, grains are formed, so that when the field is
in the molecular plane, it can in principle form any angle with the molecular
’y’ axis defined above. To simulate this, we have performed Monte Carlo
simulations with the magnetic field $B$ at varying angles with the $y$ axis of
the molecule, and we have averaged the magnetization, in the direction of the
field, for all the directions of $B$ between $0$ and $90$ degrees. This gives
the results showed in table 2, all for a $T=2$ Kelvin, the lowest reported in
the experiments.

Table 2: Average magnetization (in Bohr magnetons) for different values of $B$, being the field in the molecular plane always. The Monte Carlo calculations correspond to $T=2\ K$. B (T) | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
$M_{avg}(\mu_{B})$ | 1.30 | 1.39 | 1.47 | 1.54 | 1.61 | 1.66 | 1.73
$M_{avg}/M_{sat}$ | 0.65 | 0.70 | 0.73 | 0.77 | 0.80 | 0.83 | 0.86

The results are in very good agreement with those presented in Ref. 14. For
fields below $3$ Tesla, the magnetization is below $1,5\ \mu_{B}$, i.e. $75\%$
of the saturation vale, and for a field of magnitude $7$ Tesla, we obtain a
magnetization of $1,73\ \mu_{B}\equiv 86\%M_{sat}$, close to the $1,8\
\mu_{B}$ reported in experiments that represented $83\%$ of the saturation
value. Thus, in the frame of our model, the magnetization does not saturate
even at fields of 7 Tesla, due to the anisotropy of the exchange $J$.
333Within our model, saturation is indeed reached at much larger fields, and
the non-saturation of the magnetization in FePc powders at $B=20$ Tesla
reported in Ref. 19 has to be explained adding the canting of the FePc
sublattices found in the herringbone phase.

With our model, the non saturation of magnetization within the molecular plane
can be explained as a consequence of the spatial anisotropy in the exchange
$J$, since the difference between $J_{y}$ and $J_{x}$ is of $\sim 0.5$ meV.
This can explain why the magnetization does not saturate with a field of
several Tesla, since a field of 1 Tesla involves an energy of $g\mu_{B}B\sim
0.1$ meV for $S=1$.

For a field perpendicular to the molecular plane, that is, in the $z$
direction, the Monte Carlo results show that the magnetization is always
lower, for any value of $B$, than the average magnetization for $B$ in the
molecular plane. Although the results show that $|J_{x}|<|J_{z}|<|J_{y}|$,
with $J_{x}=$1.77 meV, $J_{z}=$2.07 meV and $J_{y}=$2.31 meV, which would
indicate that it is more favourable for the magnetic moment to move in the yz
plane, if positive value of $D=0.65$ meV is taken into account, the final
result is that it is energetically more favorable for the magnetic moment to
remain in the molecular plane. This explains the result mentioned above.

The previous results are in line with magnetization measurements in the $xy$
and $z$ direction 43, and can also explain why in the non-templated films 14,
in which the molecules are perpendicular to the substrate, when a field
parallel to the substrate is applied, the magnetization is always lower than
for the templated films, in which the field is always in the molecular plane.
This is because in non-templated films, for some grains, the field is, or
almost is, perpendicular to the molecular plane, while for some others it is
in the molecular plane, with intermediate situations for other grains, so that
the average magnetization can be simulated as an average of the in plane
magnetization with the out of plane one.

All in all, our model gives an explanation that can, at least partly, account
for the reported results for the non-saturation of the magnetization in FePc
thin films. It can be expected that this results also contribute to understand
the experiments in FePc powders.

###  3.4 Magnetic Solitons

Ferromagnetic one dimensional systems have magnetic solitons as natural
excitations above the ground state 53, 14. Solitons possess special properties
like stability, mobility without loosing their shape, and energy efficiency,
making them promising candidates for various applications in the field of
spintronics or information processing. For example, magnetic solitons have
been proposed as a mechanism for remotely manipulating single qbits 54.
Besides, the density of magnetic solitons influences experimental measurements
of specific heat 55, neutron scattering experiments 56, Mossbauer spectra 57
and dynamic structure factor 58.Several studies in similar model Hamiltonians
53, 59, 60 have found this type of excitations. They have also been reported
in experimental studies of cuasi 1D compounds, being $CsNiF_{3}$ 56, 55 the
paradigmatic ferromagnetic soliton bearing compound, and also others like
quasi 1D antiferromagnets $Li_{2}Mn_{0.98}Fe_{0.02}F_{5}$ and
$Na_{2}Mn_{0.98}Fe_{0.02}F_{5}$ 57. In relation to FePc chains, in Ref. 20
Mossbauer and a.c. susceptibility measurements, and in Ref. 14 susceptibility
measurements, were interpreted in terms of magnetic solitons, but with
different models, or by means of an equation for the energy of the solitons
derived from a renormalized Sine-Gordon theory 53. Here, in the context of our
first-principles derived model Hamiltonian, we study their existence, type and
density, in particular as a function of temperature and magnetic field.

In figure 3 the three components of the spins along a segment of 200 sites of
our model is plotted, for $T=8K$ and a magnetic field of $0.25$ Tesla in the
’y’ direction. $2\pi$ solitons can be visualized. The observed solitons are
predominantly of the xy type, i.e. the spins remain mostly in the plane of the
molecule, but a certain amount of $z$ component is also present. We have also
made calculations with $D$ 10 times higher, and checked that the number of
solitons changes little but their $z$ component decreases much further. In
ferromagnetic isotropic Heisenberg models with single ion anisotropy, solitons
are expected to be purely of the xy type for temperatures much lower than
$T=(JD)^{1/2}/k_{B}$, which amounts to between $12$ and $14$ K for the
parameters of this model, but $\sim 40$ K for $D$ 10 times higher, and this is
the reason why we observe a certain amount of z component in the solitons for
the parameters of our model, but almost none for $D$ enhanced by a factor of
10.

In order to calculate the density of magnetic solitons, the following
quantitative criterion is established, following similar criteria established
in previous works 59, 61:

Figure 3: Snapshot of the spins along a selected segment of 200 sites (out of
a total of 1000 sites per chain), representing the 1D arrangement of FePC
molecules, for $T=8$ K and a magnetic field of $B=0,25$ Tesla. Several
magnetic solitons, marked with circles, can be observed.

  1. 1.

For a given temperature T we average the y components of all the sites and and
obtain $\bar{S_{y}}$.

  2. 2.

A soliton is counted if $\bar{S_{y}}>0.4\times S$ (in our case S=1) and at a
certain site $i$ $S_{i}^{y}<-0.4\times\bar{S_{y}}$, and this remains valid for
at least 3 more sites, i.e, the soliton width is larger than 3 sites. Thus we
obtain $n_{s}$.

  3. 3.

We then average the $n_{s}$ values by repeating 2) every $10^{6}$ MC steps out
of a total of $10^{9}$ steps.

Besides, for low fields, like in Fig. 3, we can roughly count the number of
solitons per 100 sites visually, as is sketched in that figure. The criterion
above gives the same number as the visual count. With this criterion, the
density of magnetic solitons is plotted in Fig. 4, together with a fit by an
Arrhenius equation 59. We have tested that across the full range of
temperatures and magnetic field strengths investigated, the statistical error
associated with our results remained below $1.10^{-2}$ and decreases as
$n_{s}$ decreases.

$n_{s}=(A/T)\exp(-E_{s}/T)$ (6)

Figure 4: Density of magnetic solitons as a function of temperature for
several magnetic fields. The lines are fittings with eq. 6. The error bars are
covered by the symbols.

We can see that the density of solitons increases exponentially as a function
of temperature up to the transition temperature. From this temperature on,
magnetic solitons decay abruptly to zero (not shown), and this occurs
naturally because the magnetic order disappears. The fitting with the
Arrhenius equation gives excellent results. Soliton energies are shown in Fig.
5 as a function of the magnetic field. For low magnetic fields, the energy is
around $73$ K, and grows linearly with the field. This energy is much higher
than any of the $J$, and this is in line with what was found in
antiferromagnetic Heisenberg models, where it was shown that a small amount of
anisotropic exchange leads to the formation of solitons 59, of energy much
higher than $J$. The value obtained for the soliton energy has a very good
agreement with what was obtained from Mossbauer/a.c. susceptibility
experiments 20, 17, i.e. $72$ K with $B=0.08$ Tesla. It is also in between the
energy values from the classical domain wall theory $E_{s}=2\ \pi\
\sqrt{JD}\simeq 90K$ and the Sine-Gordon theory $E_{s}=4\ \sqrt{JD}\simeq 60K$
57 for an isotropic Heisenberg model with single ion anisotropy $D$, adopting
for $J$ the value of $J_{y}$ for our model.

Figure 5: Energy of magnetic solitons as a function of magnetic field
according to the fitting done with eq. 6. A linear regression is depicted by a
dashed line.

Finally, in Fig. 6 we plot a colour map for the density of mangetic solitons
as a function of magnetic field and temperature. It can be clearly observed
how magnetic solitons are excited with temperature and magnetic field up to
the transition to the paramagnetic phase. The magnetic field increases the
transition temperature (for example at $B=1$ Tesla the transition temperature
is $\sim 12K$ and for $B=4$ Tesla is $\sim 20K$) and consequently solitons are
excited up to higher temperatures as the field increases, but at the same time
the soliton density reaches lower values. The highest solitons densities are
reached for fields below $2$ Tesla at temperatures between $8$ and $12$ K.

Figure 6: Density of magnetic solitons as a function of temperature for
several magnetic fields.

All in all, Monte Carlo simulations show the existence of magnetic solitons in
the model proposed for $\alpha_{+}$ FePc in thin films, which are thermally
activated in the ferromagnetic phase and increase their energy with the
application of a magnetic field, which renders them more difficult to be
thermally excitated.

##  4 Conclusions.

In this work we have theoretically studied $\alpha_{+}$ polymorphism of FePc
Phthalocyanines that is found in thin films, consisting of quasi one
dimensional chains of FePc molecules. By means of DFT calculations, after
reproducing the structural properties of the phase, we have calculated the
parameters of a Heisenberg model with single ion anisotropy, finding that due
to the spin-orbit interaction the magnetic exchange interaction is
anisotropic, i.e. that $J_{x}\neq J_{y}\neq J_{z}$, with $J_{x}=1.77$ meV,
$J_{y}=2.31$ meV, and $J_{z}=2.07$ meV. Calculations also show a single ion
anisotropy $D$ of $0.65$ meV and negligible $E$. Adding to this one
dimensional model a small interchain hopping, as estimated from experiments
and theory 14, 47, with classical Monte Carlo calculations we calculated the
magnetization at low temperatures, finding that the non saturation found in
thin films at fields up to 7 Tesla can be explained by the model. However, in
FePc powders, the non-saturation of the magnetization extends to much higher
fields, as a consequence of the canted structure of the herringbone
arrangement. We have also found that our model is consistent with the
existence of magnetic solitons in FePc TF. The solitons observed in MC
calculations are of the $2\pi$ or double kink type. It is found that there is
a non-negligible density of solitons even at zero magnetic field, and that the
soliton energy for zero magnetic field is around $73$ K, in line with
Mossbauer and a.c. susceptibility experiments. Their energy increases linearly
with magnetic field, and the thermal energy window in which they are exited is
at higher temperatures as the field increases, but at the same time their
density also decreases.

##  5 Conflicts of interest

There are no conflicts to declare.

##  6 Acknowledgements.

R. P. , P. A and I. H. acknowledge financial support from CONICET PIP0883 and
the computational time provided by the CCT-Rosario computational center and
Computational Simulation Center (CSC) for Technological Applications, members
of the High Performance Computing National System (SNCAD-MincyT Argentina). We
acknowledge fruitful conversations with N. O. Vargas, L. O. Manuel, C. J.
Gazza and A. E. Trumper.

## References

  * Guo et al. 2019 Guo, L.; Gu, X.; Zhu, X.; Sun, X. Recent Advances in Molecular Spintronics: Multifunctional Spintronic Devices. _Advanced Materials_ 2019, _31_ , 1805355
  * Inganäs 2018 Inganäs, O. Organic photovoltaics over three decades. _Advanced materials_ 2018, _30_ , 1800388
  * Grant et al. 2019 Grant, T. M.; Josey, D. S.; Sampson, K. L.; Mudigonda, T.; Bender, T. P.; Lessard, B. H. Boron Subphthalocyanines and Silicon Phthalocyanines for Use as Active Materials in Organic Photovoltaics. _The Chemical Record_ 2019, _19_ , 1093–1112
  * Fleet et al. 2017 Fleet, L. R.; Stott, J.; Villis, B.; Din, S.; Serri, M.; Aeppli, G.; Heutz, S.; Nathan, A. Self-Assembled Molecular Nanowires for High-Performance Organic Transistors. _ACS Applied Materials & Interfaces_ 2017, _9_ , 20686–20695
  * Klauk 2010 Klauk, H. Organic thin-film transistors. _Chemical Society Reviews_ 2010, _39_ , 2643–2666
  * Zhou et al. 2021 Zhou, W.; Yutronkie, N. J.; Lessard, B. H.; Brusso, J. L. From chemical curiosity to versatile building blocks: unmasking the hidden potential of main-group phthalocyanines in organic field-effect transistors. _Mater. Adv._ 2021, _2_ , 165–185
  * Kalyani and Dhoble 2012 Kalyani, N. T.; Dhoble, S. Organic light emitting diodes: Energy saving lighting technology—A review. _Renewable and Sustainable Energy Reviews_ 2012, _16_ , 2696–2723
  * Plint et al. 2016 Plint, T.; Lessard, B. H.; Bender, T. P. Assessing the potential of group 13 and 14 metal/metalloid phthalocyanines as hole transport layers in organic light emitting diodes. _Journal of Applied Physics_ 2016, _119_ , 145502
  * Rivnay et al. 2018 Rivnay, J.; Inal, S.; Salleo, A.; Owens, R. M.; Berggren, M.; Malliaras, G. G. Organic electrochemical transistors. _Nature Reviews Materials_ 2018, _3_ , 1–14
  * Basova and Ray 2020 Basova, T. V.; Ray, A. K. Review—Hybrid Materials Based on Phthalocyanines and Metal Nanoparticles for Chemiresistive and Electrochemical Sensors: A Mini-Review. _ECS Journal of Solid State Science and Technology_ 2020, _9_ , 061001
  * Cranston and Lessard 2021 Cranston, R. R.; Lessard, B. H. Metal phthalocyanines: thin-film formation, microstructure, and physical properties. _RSC Adv._ 2021, _11_ , 21716–21737
  * Bartolomé et al. 2014 Bartolomé, J.; Monton, C.; Schuller, I. K. In _Molecular Magnets: Physics and Applications_ ; Bartolomé, J., Luis, F., Fernández, J. F., Eds.; Springer Berlin Heidelberg: Berlin, Heidelberg, 2014; pp 221–245
  * Hoshino et al. 2003 Hoshino, A.; Takenaka, Y.; Miyaji, H. Redetermination of the crystal structure of $\alpha$-copper phthalocyanine grown on KCl. _Acta Cryst. B_ 2003, _59_ , 393
  * Wu et al. 2019 Wu, Z.; Robaschik, P.; Fleet, L. R.; Felton, S.; Aeppli, G.; Heutz, S. Controlling Ferromagnetic Ground States and Solitons in Thin Films and Nanowires Built from Iron Phthalocyanine Chains. _Advanced Functional Materials_ 2019, _29_ , 1902550
  * Yim et al. 2003 Yim, S.; Heutz, S.; Jones, T. S. Influence of intermolecular interactions on the structure of phthalocyanine layers in molecular thin film heterostructures. _Phys. Rev. B_ 2003, _67_ , 165308
  * Wang et al. 2010 Wang, H.; Mauthoor, S.; Din, S.; Gardener, J. A.; Chang, R.; Warner, M.; Aeppli, G.; McComb, D. W.; Ryan, M. P.; Wu, W. et al. Ultralong Copper Phthalocyanine Nanowires with New Crystal Structure and Broad Optical Absorption. _ACS Nano_ 2010, _4_ , 3921–3926, PMID: 20527798
  * Bartolome et al. 2015 Bartolome, F.; Bunau, O.; Garcia, L. M.; Piantek, C. R. N. M.; Pascual, J. I.; Schuller, I. K.; Gredig, T.; Wilhelm, F.; Rogalev, A.; Bartolome, J. Molecular tilting and columnar stacking of Fe phthalocyanine thin films on Au(111). _Journal of Applied Physics_ 2015, _117_ , 17A735
  * Ashida et al. 1966 Ashida, M.; Uyeda, N.; Suito, E. Unit Cell Metastable-form Constants of Various Phthalocyanines. _Bulletin of the Chemical Society of Japan_ 1966, _39_ , 2616–2624
  * Evangelisti et al. 2002 Evangelisti, M.; Bartolomé, J.; de Jongh, L. J.; Filoti, G. Magnetic properties of $\alpha$-iron(II) phthalocyanine. _Phys. Rev. B_ 2002, _66_ , 144410
  * Filoti et al. 2006 Filoti, G.; Kuzmin, M. D.; Bartolomé, J. Mössbauer study of the hyperfine interactions and spin dynamics in $\alpha$-iron(II) phthalocyanine. _Phys. Rev. B_ 2006, _74_ , 134420
  * Vargas et al. 2020 Vargas, N. M.; Torres, F.; Baker, A. A.; Lee, J. R. I.; Kiwi, M.; Willey, T. M.; Monton, C.; Schuller, I. K. Helical spin structure in iron chains with hybridized boundaries. _Applied Physics Letters_ 2020, _117_ , 213105
  * Gredig et al. 2012 Gredig, T.; Werber, M.; Guerra, J. L.; Silverstein, E. A.; Byrne, M. P.; Cacha, B. G. Coercivity Control of Variable-Length Iron Chains in Phthalocyanine Thin Films. _Journal of Superconductivity and Novel Magnetism_ 2012, _25_ , 2199–2203
  * Kresse and Hafner 1993 Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. _Phys. Rev. B_ 1993, _47_ , 558
  * Kresse and Hafner 1993 Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Open-Shell Transition Metals. _Phys. Rev. B_ 1993, _48_ , 13115
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. _Comput. Mater. Sci._ 1996, _6_ , 15
  * Kresse and Furthmüller 1996 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
  * Kresse and Joubert 1999 Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. _Phys. Rev. B_ 1999, _59_ , 1758
  * Hafner 2008 Hafner, J. Ab-Initio Simulations of Materials Using VASP: Density-Functional Theory and Beyond. _J. Comput. Chem._ 2008, _29_ , 2044–2078
  * Grimme et al. 2010 Grimme, S.; Ehrlich, J. A. S.; ; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. _J. Chem. Phys._ 2010, _132_ , 154104
  * J.Klimes et al. 2010 J.Klimes; Bowler, D. R.; Michaelides, A. Van der Waals density functionals applied to solids. _Phys. Rev. B_ 2010, _83_ , 195133
  * Dudarev et al. 1998 Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron energy loss. _Phys. Rev. B_ 1998, _57_ , 1505
  * Mugarza et al. 2012 Mugarza, A.; Robles, R.; Krull, C.; Korytár, R.; Lorente, N.; Gambardella, P. Electronic and magnetic properties of molecule-metal interfaces: Transition-metal phthalocyanines adsorbed on Ag(100). _Phys. Rev. B_ 2012, _85_ , 155437
  * Wang et al. 2018 Wang, Y.; Li, X.; Zheng, X.; Yang, J. Manipulation of spin and magnetic anisotropy in bilayer magnetic molecular junctions. _Phys. Chem. Chem. Phys._ 2018, _20_ , 26396–26404
  * Marom and Kronik 2009 Marom, N.; Kronik, L. Density functional theory of transition metal phthalocyanines, II: electronic structure of MnPc and FePc—symmetry and symmetry breaking. _Applied Physics A_ 2009, _95_ , 165–172
  * Białek et al. 2003 Białek, B.; Kim, I. G.; Lee, J. I. First-principles study on the electronic structures of iron phthalocyanine monolayer. _Surface Science_ 2003, _526_ , 367–374
  * Kuz’min et al. 2009 Kuz’min, M. D.; Hayn, R.; Oison, V. Ab initio calculated XANES and XMCD spectra of Fe(II) phthalocyanine. _Phys. Rev. B_ 2009, _79_ , 024413
  * Brena et al. 2011 Brena, B.; Puglia, C.; de Simone, M.; Coreno, M.; Tarafder, K.; Feyer, V.; Banerjee, R.; Göthelid, E.; Sanyal, B.; Oppeneer, P. M. et al. Valence-band electronic structure of iron phthalocyanine: An experimental and theoretical photoelectron spectroscopy study. _The Journal of Chemical Physics_ 2011, _134_ , 074312
  * Marom and Kronik 2009 Marom, N.; Kronik, L. L. Density functional theory of transition metal phthalocyanines, II: electronic structure of MnPc and FePc- symmetry and symmetry breaking. _Appl. Phys._ 2009, _95_ , 165–172
  * Ichibha et al. 2017 Ichibha, T.; Hou, Z.; Hongo, K.; Maezono1, R. New Insight into the Ground State of FePc: A Diffusion Monte Carlo Study. _Sci. Rep._ 2017, _7_ , 2011
  * Phung et al. 2023 Phung, Q. M.; Nam, H. N.; Saitow, M. Unraveling the Spin-State Energetics of FeN4 Complexes with Ab Initio Methods. _The Journal of Physical Chemistry A_ 2023, _127_ , 7544–7556
  * Nakamura et al. 2012 Nakamura, K.; Kitaoka, Y.; Akiyama, T.; Ito, T. Constraint density functional calculations for multiplets in a ligand-field applied to Fe-phthalocyanine. _Phys. Rev. B_ 2012, _85_ , 235129
  * Xia et al. 2023 Xia, H.; Li, L.; Fang, M.; Jones, T. S.; Yang, J. Molecular-orientation-dependent magnetic properties of iron phthalocyanine (FePc) thin films and microwires. _Organic Electronics_ 2023, _121_ , 106870
  * Bartolomé et al. 2010 Bartolomé, J.; Bartolomé, F.; García, L. M.; Filoti, G.; Gredig, T.; Colesniuc, C. N.; Schuller, I. K.; Cezar, J. C. Highly unquenched orbital moment in textured Fe-phthalocyanine thin films. _Phys. Rev. B_ 2010, _81_ , 195405
  * Natoli et al. 2018 Natoli, C. R.; Krüger, P.; Bartolomé, J.; Bartolomé, F. Determination of the ground state of an Au-supported FePc film based on the interpretation of Fe $K$\- and $L$-edge x-ray magnetic circular dichroism measurements. _Phys. Rev. B_ 2018, _97_ , 155139
  * Barraclough et al. 1970 Barraclough, C. G.; Martin, R. L.; Mitra, S.; Sherwood, R. C. Paramagnetic Anisotropy, Low Temperature Magnetization, and Electronic Structure of Iron(II) Phthalocyanine. _The Journal of Chemical Physics_ 1970, _53_ , 1643–1648
  * Torres et al. 2020 Torres, F.; Kiwi, M.; Vargas, N. M.; Monton, C.; Schuller, I. K. Chiral symmetry and scale invariance breaking in spin chains. _AIP Advances_ 2020, _10_ , 025215
  * Wu et al. 2013 Wu, W.; Harrison, N. M.; Fisher, A. J. Suitability of chromium phthalocyanines to test Haldane’s conjecture: First-principles calculations. _Phys. Rev. B_ 2013, _88_ , 224417
  * Moriya 1960 Moriya, T. Anisotropic Superexchange Interaction and Weak Ferromagnetism. _Phys. Rev._ 1960, _120_ , 91–98
  * Guterding et al. 2016 Guterding, D.; Valentí, R.; Jeschke, H. O. Reduction of magnetic interlayer coupling in barlowite through isoelectronic substitution. _Phys. Rev. B_ 2016, _94_ , 125136
  * Whangbo and Xiang 2017 Whangbo, M.-H.; Xiang, H. _Handbook of Solid State Chemistry_ ; John Wiley and Sons, Ltd, 2017; Chapter 10, pp 285–343
  * Wäckerlin et al. 2022 Wäckerlin, C.; Cahlík, A.; Goikoetxea, J.; Stetsovych, O.; Medvedeva, D.; Redondo, J.; Švec, M.; Delley, B.; Ondráček, M.; Pinar, A. et al. Role of the Magnetic Anisotropy in Atomic-Spin Sensing of 1D Molecular Chains. _ACS Nano_ 2022, _16_ , 16402–16413
  * Gallardo et al. 2019 Gallardo, I.; Arnau, A.; Delgado, F.; Baltic, R.; Singha, A.; Donati, F.; Wäckerlin, C.; Dreiser, J.; Rusponi, S.; Brune, H. Large effect of metal substrate on magnetic anisotropy of Co on hexagonal boron nitride. _New Journal of Physics_ 2019, _21_ , 073053
  * Samalam and Kumar 1982 Samalam, V. K.; Kumar, P. Statistical mechanics of a magnetic chain. _Phys. Rev. B_ 1982, _26_ , 5146–5152
  * Cuccoli et al. 2016 Cuccoli, A.; Nuzzi, D.; Vaia, R.; Verrucchi, P. Single-qubit remote manipulation by magnetic solitons. _Journal of Magnetism and Magnetic Materials_ 2016, _400_ , 149–153, Proceedings of the 20th International Conference on Magnetism (Barcelona) 5-10 July 2015
  * Ramirez and Wolf 1982 Ramirez, A. P.; Wolf, W. P. Specific Heat of ${\mathrm{CsNiF}}_{3}$: Evidence for Spin Solitons. _Phys. Rev. Lett._ 1982, _49_ , 227–230
  * Kjems and Steiner 1978 Kjems, J. K.; Steiner, M. Evidence for Soliton Modes in the One-Dimensional Ferromagnet CsNi${\mathrm{F}}_{3}$. _Phys. Rev. Lett._ 1978, _41_ , 1137–1140
  * Frommen et al. 1996 Frommen, C.; Mangold, M.; Pebler, J. Magnetic Solitons in the 1-D Antiferromagnetic Chains of Li2Mn0.98Fe0.02F5 and Na2Mn0.98Fe0.02F5. _Zeitschrift für Naturforschung A_ 1996, _51_ , 939–949
  * Leung and Huber 1979 Leung, K.; Huber, D. Soliton dynamic structure factors in a planar ferromagnetic chain. _Solid State Communications_ 1979, _32_ , 127–130
  * Gaulin 1987 Gaulin, B. D. Soliton spin configurations along the classical anisotropic Heisenberg chain. _Journal of Applied Physics_ 1987, _61_ , 4435–4437
  * Etrich et al. 1985 Etrich, C.; Mikeska, H. J.; Magyari, E.; Thomas, H.; Weber, R. Solitons on a discrete ferromagnetic spin chain. _Zeitschrift für Physik B Condensed Matter_ 1985, _62_ , 97–111
  * Gerling and Landau 1984 Gerling, R.; Landau, D. Spin solitons in the classical xy-chain. _Journal of Magnetism and Magnetic Materials_ 1984, _45_ , 267–271