Monte-Carlo phase diagram of a Hubbard-Peierls model in the search for spin crossover transition in \pi-conjugated polymers

# Monte-Carlo phase diagram of a Hubbard-Peierls model in the search for spin crossover transition in π-conjugated polymers

## Abstract

We present a Monte Carlo study of the finite temperature properties of an extended Hubbard-Peierls model describing one dimensional -conjugated polymers. The model incorporates electron-phonon and hyperfine interaction and it is solved at the mean field level for half filling. In particular we explore the model as a function of the strength of electron-electron and electron-phonon interactions. At low temperature the system presents a diamagnetic to antiferromagnetic transition as the electron-electron interaction strength increases. At the same time by increasing the electron-phonon coupling there is a transition from a homogeneous to a Peierls dimerized geometry. As expected such a Peierls dimerized phase disappears at finite temperature as a result of thermal vibrations. More intriguing is the interplay between the electron-phonon and the electron-electron interactions at finite temperature. In particular we demonstrate that for a certain region of the parameter space there is a spin-crossover, where the system transits from a low-spin to a high-spin state as the temperature increases. In close analogy to standard spin-crossover in divalent magnetic molecules such a transition is entropy driven. Finally we discuss the rôle played by the hyperfine interaction over the phase diagram.

## I Introduction

Over the past few years the field of organic spintronics has received interest from a growing and diverse scientific community Alex (); Stefano (); Stefano2 (); dediu1 (); WW1 (). This is primarily driven by its potential for opening new avenues to cheap, low-weight, mechanically flexible, chemically inert and bottom up fabricated spin-devices. The most important advantage of organic materials Greg (); harris (); VIK (); sbandy () over their inorganic counterparts das () is the strong suppression of any efficient spin flipping mechanisms, resulting in extremely long spin lifetimes. Spin orbit coupling is very weak in organic semiconductors, as Carbon has a low atomic number () and the strength of spin orbit interaction is proportional to . Furthermore also the hyperfine interaction is weak in molecules, since there are very few nuclei with non-vanishing magnetic moment in the upper part of the periodic table and the molecular orbitals typically responsible for the electron transport are extended -type Alex2 (). A recent experimental work fert () on (La,Sr)MnO(LSMO)/tris[8-hydroxyquinoline] Aluminum(Alq)/Co tunnel junctions reported a giant tunnel magneto resistance of up to 300% at low temperatures. This is a value that compares well with the best inorganic tunnel junctionsikeda () and it is attributed to a favorable interaction between the organic molecules and the magnetic metallic surface fert (); spinterface ().

Amongst the various possible materials for organic spintronics, -conjugated organic semiconductors appears very appealing dediu2 (); vardeny (). This is because of their very long and relatively temperature independent spin relaxation time and their ability to form good interfaces with metal electrodes when incorporated in spin-valve-like devices. The most relevant structural feature of such polymers is their planar shape. The -electron wave function distends from the molecular plane and it can easily interact with the wave function of adjacent molecules. Therefore a face-to-face molecular configuration is usually stabilized through strong - bonding and diffusive Van der Waals interaction. This leads to a molecular stacking arrangement resulting in the formation of a low dimensional lattice. Because of this peculiar structural conformation many of the -conjugated polymers can be described by simplified one dimensional (1D) model Hamiltonians picozzi ().

An important aspect for the successful integration of organic semiconductors in magnetic memories and in magnetic switching devices is the feasibility of manipulating the spin orientation in the organic media. This is difficult to achieve in non-magnetic molecules by the tiny non-equilibrium spin population originating from spin injection. In fact the standard techniques for manipulation used in inorganic semiconductors, for example optical methods, are ineffective because of the weak spin-orbit interaction. A more promising option is that of manipulating the internal spin degrees of freedom of the organic medium, when this is magnetic. Intriguingly there is a vast class of molecules, generally known as spin crossover compounds, whose spin state can be changed from low spin to high spin by an external perturbation SC (). Since the crossover transition is entropy driven, it is most typically achieved by increasing the temperature, although also light, pressure and electro-chemical redox reactions can all produce it. Most recently the possibility of spin crossover driven by static electric fields has been proposed theoretically Nadjib (); Kim (); Andrea (). Notably all these processes occur in compounds incorporating a transition metal [usually Fe(II)], which is responsible for the magnetic moment. Thus, it is interesting to explore the possibility of obtaining spin crossover in organic -conjugated polymers.

In this work we use a combination of energy minimization techniques and Monte Carlo (MC) simulations to investigate the temperature dependent phase diagrams for a model Hamiltonian describing -conjugated polymers. In particular, we explore the parameter space of the model in the search of a spin crossover transition. Intriguingly we find that the interplay between on-site Coulomb repulsion and electron-phonon (el-ph) coupling, leading to Peierls distortion, can be responsible for spin crossover.

## Ii Computational Methods

We consider an extended single-site Hubbard-Peierls model Hubbard () for a 1D lattice described by the following Hamiltonian

 H = ∑iσϵiσc†iσciσ+ (1) + ∑ij[tij+α(qi−qj)](c†iσcjσ+c†jσciσ)+ + U∑i^ni↑^ni↓+JH∑iαβ→Si⋅[c†iα(→σαβ)ciβ]+ + 12k∑i,j(qi−qj)2

where () denotes the creation (annihilation) operator for an electron at site with spin , is the onsite energy and is the transfer integral for an uniform undimerized lattice. Here we consider only nearest neighbour hopping, i.e. for and for any other pair. The other microscopic parameters of the model are the el-ph coupling parameter, , the Hubbard repulsion strength, , the hyperfine exchange , and the elastic constant, . Thus, the second term of the Hamiltonian, in addition to electron hopping, describes the el-ph coupling, with being the atomic displacement of site (we consider 1D longitudinal motion only). The third term is the standard on-site Hubbard repulsion while the fourth one describes a Heisenberg-like interaction between the electron spins, , and a set of classical vectors, , representing the nuclear spins. Finally the last term is classical elastic energy.

In all our calculations we consider only the half-filling case (one electron per site) so that the model exhibits an insulating behaviour at for an undistorted infinite chain. The data presented here are for 10 site long chains as further tests for longer chains (20 sites) gave rather similar results. Our approach consists in replacing the Hubbard term with its unrestricted mean-field approximation meanfield () and then in solving the Hamiltonian self-consistently for different lattice displacements . Energy minimization is performed by conjugate gradient over and further verified by additional simulated annealing SA (). The bond lengths, , are calculated from the ground state displacements as , where is the equillibrium bond distance. The main observables calculated are the dimerization parameter, , and the local magnetic moment per site, . These are defined respectively in Eq. (2) and Eq. (3)

 D=N−2∑i=1|xi−xi+1|+|xi+1−xi+2|, (2)
 mi=|n↑i−n↓i|. (3)

Note that in Eq. (3) is the electron spin density at site , which is calculated as the diagonal matrix elements of the density matrix , where and are respectively -th electronic eigenvalues and eigenvectors. Our results are then plotted in a phase-diagram like form, where the different phases are presented as a function of and .

For finite temperature simulations we consider the system described by the Hamiltonian of Eq. (1) and evolve the classical dynamical variables and by using the standard Metropolis algorithm metropolis (). Note that since the hyperfine interaction has little effect on the - phase diagram, so that in what follows we will neglect it unless otherwise indicated. As such the only classical dynamical variables are the atomic displacements . In the Metropolis algorithm the acceptance probablity of a new state is unity if the new configuration has an energy lower than that of the old configuration. Otherwise it is given by the Boltzmann factor , where is the difference in the Gibbs free energy between the old and new configuration. By using the grand-canonical ensemble the Gibbs free energy can be written as

 G({qi})=−1βNo∑n=1ln(1+e−β[En({qi})−μ]), (4)

where the chemical potential is obtained from

 N=∑n1eβ[En({qi})−μ]+1, (5)

with being the inverse temperature and the total number of electrons ( here). For every value of (, ) and each temperature the system is allowed to reach equilibrium. Then both and are calculated over several million MC steps.

In what follows we will express all energy related quantities (including the temperature) as a function of , which sets the energy scale of the problem. The on-site energy is taken to be zero and is 5 . Note that for 2.5 eV this corresponds to  eV/Å, which is in between the value for the H molecule and that of Au monoatomic chains Will (). Also it is important to note that Å is a value commonly used in recent literature about transport in organic polymers picozzi (); saxena ().

## Iii Results and Discussion

### iii.1 Ground State

Let us begin our discussion by investigating the - phase diagram at , which is presented in Fig. 1. This is populated by four difference phases characterized by the different combined structural and magnetic properties of the chain. In particular there are two magnetic states and two geometrical configurations. For small and the chains are undistorted (the atomic spacing is approximately uniform throughout the chain) and in a non-spin-polarized diamagnetic (DM) state. We denote this phase as DM-undistorted. As increases eventually a spin-polarized solution develops. This is the one expected from the mean field Hubbard model at half-filling, i.e. it is a antiferromagnetic (AFM) phase, where local moments form at each atomic site, but their orientation alternates along the chain. Such an AFM phase may or may not be accompanied by a structural distortion, depending on the value of . In general however there is a vast region of the model parameter space, where no significant distortion appears for the AFM spin state. We denote this phase as AFM-undistorted.

As increases for moderate the system progressively developed a Peierls instability and makes a transition to a geometry where long bond distances alternate to short ones. This dimerized phase, expected for the non-interacting case, remains diamagnetic (DM-dimerized) for small but can coexist with a AFM solution for a significant range of parameter. In summary the phase diagram is characterized by a competition between the on-site repulsion, driving the magnetic instability, and the el-ph coupling, driving the Peierls distortion.

In the discussion we have assigned the phase boundary of the magnetic transition to the condition . In contrast assigning the phase boundary to the Peierls transition is more complicated since changes continuously upon increasing . Thus we have used the operational definition of placing the phase boundaries at , which is interpreted as representative of strong dimerization. The complete evolution of as a function of and is presented in the right hand side panel of Fig. 1. The figure clearly reveals the interplay between el-ph coupling and Coulomb repulsion. In fact grows almost linearly with for small but then is drastically reduced as grows.

### iii.2 Finite temperature phase diagram

We now move to study the finite temperature properties of our model. These are summarized in Fig. 2, where we present the - phase diagram for two representative temperatures, respectively  1/eV (these correspond respectively to or, for  eV, to  K).

The most important feature of the finite temperature plots is the complete absence of a structural phase transition. This means that in general the system does not dimerize any longer as the temperature is increased. The dimerization is instead replaced by a general increase in bond length and by a random distribution of the various bonds along the chain.

More details about the structure of the chains at finite temperature can be found in the right-hand side panes of Fig. 2, which show the quantity [defined in Eq. (2)] as a function of and for and . Note that at finite temperature the quantity plotted is a measure of the disorder in the bond distances of the chain as thermal vibrations onset. As such is maximum at high values of the el-ph coupling and low but falls rapidly as increases.

The second most striking feature of the finite temperature - phase diagram is the movement of the DM-AFM phase boundary towards lower as the temperature increases. This essentially means that as the temperature grows it takes less on-site Coulomb interaction to drive the system towards a magnetic instability. We further explore this finding in figure 3, where we present the critical value, , at which the magnetic phase develops. This effectively represents the position of the DM-AFM phase boundary. In particular is plotted as a function of the temperature, , and for three different el-ph strengths .

In general we find that the DM-AFM phase boundary moves in response to the disappearing of the distorted phase. Thus for the lower value of (0.12), for which there is no distorted phase even at , the phase boundary does not change as the temperature is increased and remains constant at . For the larger values of the el-ph interaction strengths investigated the phase diagram presents both a distorted and a homogeneous structural phase depending on . In this case the DM-AFM phase boundary (i.e. ) decreases fast at low temperatures in response to the melting of the distorted phase and then becomes essentially constant.

An important consequence of the movement of the DM-AFM phase boundary as the temperature increases is the fact that there is a vast region in the - diagram where the system makes a DM to AFM transition as the temperature increases. Such a region is the one enclosed between the two vertical lines marking respectively the phase boundary at finite temperature and at in figure 2. For the particular values of and found in such a region (called the spin crossover -SC- region) there is a temperature driven spin crossover. This is analyzed next.

### iii.3 Spin crossover

We now explain the spin crossover transition by using the standard framework of spin crossover usually employed for magnetic molecules incorporating divalent transition metals SC (). In general the thermodynamically stable phase at finite temperature of a system that can assume different competing configurations is the one with the lowest Gibbs free energy, . For the present case, where the competition is among the DM and the AFM phases, the relevant quantity is the difference, , between their Gibbs energies,

 ΔG=ΔH−TΔS (6)

where and are respective the enthalpy and entropy differences.

For standard spin crossover , so that the most stable phase at low temperature in low spin (DM here). However, since the crossover transition is associated to the softening of the phonon modes of the first coordination shell around the transition metal and to the formation of a local magnetic moment, we also have . Hence as the temperature increases the entropic contribution to the Gibbs energy may eventually dominate over the enthalpic one and drive a phase transition. We now want to establish that the same mechanism holds for the spin crossover region of the - diagram of Fig. 2.

We have already demonstrated (see figure 1) that for the spin crossover region is occupied by the DM phase, meaning that . Therefore one has only left to show that also . In general the entropy comprises of two main contributions, an electronic, , and a vibrational one, . Since the AFM phase is characterized by local spins, which are absent for the DM phase, we can immediately conclude that . A more precise evaluation of can be obtained by computing

 Sel=−kBTr[^ρln^ρ], (7)

where is the system density matrix. The calculated as a function of for the representative value of , for which the spin crossover region is quite large, is presented in figure 4. The electron densities entering the evaluation of have been calculated as follows. For the low temperature phase (DM) is calculated by fixing the occupation to for every site and the geometry is that obtained from the diagram. In contrast the entropy of the high-temperature phase is computed from a density matrix in which the occupation is fixed to the proper antiferromagnetic state (the temperature is ) and the geometry is again that of the solution. We have checked that the finite temperature geometry is rather similar to that obtained at for such a density-constraint solution.

The general trend can be understood as follows. For small values of there is no local magnetic moment formation, regardless of the chain geometry, so that . Then, as gets larger a local magnetic moment gradually appears at each site, producing a linear increase of . Such an increases continues until the local moment reaches the maximum value compatible with the chosen electron filling, which is 1  ( is the Bohr magneton) for half-filling. At this point there is no further change in the electronic entropy and saturates to a positive value.

Similarly we can also estimate the vibrational contribution to the entropy. This is obtained from the molecule phonon spectrum, , as

 Svib=kBlnNv∑i1eβℏωi−1+kBβNv∑iℏωi, (8)

where is the vibration frequency of the -th mode and is the total number of modes. We calculate the phonon spectra of the AFM and DM configurations by diagonalizing the associated dynamical matrices. These are constructed by finite difference, i.e. by displacing the atomic sites by a small fraction of the equilibrium bond length (0.1 %) and then by numerically evaluating the energy gradient (the force) associated to such a displacement. The density matrices and the initial geometries used to construct the finite difference dynamical matrices are the same used for calculating the electronic contribution to the entropy. Also in this case as a function of for is presented in figure 4.

In general the vibrational contribution to the entropy difference shows only a small dependence on the Coulomb on-site repulsion and approximately . However we also report a relatively sharp decreases as approaches 1.7. This is value close to for , i.e. corresponds to a region in the parameter space where our Monte Carlo analysis does not find a DM-AFM transition and therefore system remains in the AFM state at any of the temperatures investigated.

In summary the picture emerging from figure 4 is that of a region in which is always positive. This is the only region of the parameter space where the entropy can drive the spin crossover and substantially matches the spin crossover region observed in our Monte Carlo simulations (see Fig. 2). We then conclude that also in this case, where the magnetic moment is not associated to the shell of a transition metal, the spin crossover is entropy driven.

Finally before summarizing, we wish to comment on the rôle played by the hyperfine interaction. In general we expect very little changes to the phase diagram obtained by neglecting the hyperfine contribution to the total energy, since this is rather small for realistic hyperfine coupling strengths. In particular we have verified that minor modifications to the phase diagram start to appear for is the region of , which corresponds to local magnetic fields of 10 T (considering and of the order of 1 eV).

## Iv Conclusion

In summary we have discussed the phase diagram of an extended Hubbard-Peierls model chosen to emulate the electronic structure of conjugated polymer chains. The model has been explored both in its ground state and at finite temperature by Monte Carlo methods combined with a mean field treatment of the Hubbard many-body interaction. At the model presents four different phases depending on the relative strength of the Coulomb on-site repulsion and the el-ph coupling strength . The four phases are characterized either by a diamagnetic or an antiferromagnetic magnetic state and by the possible presence of a dimerized geometrical configuration. By increasing the temperature the structural distortion disappears and only the phase boundary between the DM and the AFM solution remains. Intriguingly the position in the - parameter space of such a phase boundary changes with temperature so that there is a region of parameters where a temperature driven DM-AFM spin crossover transition can be found. We have investigated the nature of the phase transition by calculating the relative entropy of the different magnetic phases and found that this is indeed entropy driven as in the most conventional case of divalent magnetic molecules. This suggests that -type magnetism can be achieved in organic polymers and that this can be tuned by temperature.

## Acknowledgments

This work is sponsored by Science Foundation of Ireland CSET grant underpinning CRANN. Computational resources have been provided by the HEA IITAC project managed by the Trinity Centre for High Performance Computing.

### References

1. A.R. Rocha, V. M. Garcìa-Suarez, S.W. Bailey, CJ. Lambert, J. Ferrer and S. Sanvito, Nature Materials 4, 335 (2005).
2. S. Sanvito, J. Mat. Chem. 17, 4455 (2007).
3. S. Sanvito, Nature Materials 6, 803 (2007).
4. V.A. Dediu , L.E. Hueso, I.  Bergenti and C. Taliani, Nature Materials 8, 707 (2009).
5. L. Bogani and W. Wernsdorfer, Nature Materials 7, 179 (2008).
6. G. Szulczewski, S. Sanvito and J.M.D. Coey, Nature Materials 8, 693 (2009).
7. C.B.  Harris, R.L. Schlupp and H. Schuch, Phys. Rev. Lett. 30, 1019 (1973).
8. V.I. Krinichnyi , S.D. Chemerisov and Y.S. Lebedev, Phys. Rev. B 55, 16233 (1997).
9. S. Pramanik, C.G. Stefanita, S. Patibandla, S. Bandyopadhyay K. Garre, N. Harth and M. Cahay, Nature Nanotech. 2, 216 (2007).
10. I.  Zutic, J. Fabian and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).
11. S. Sanvito and A.R. Rocha, J. Comp. Theo. Nano. 3, 624 (2006).
12. C. Barraud, P. Seneor, R. Mattana, S. Fusil, K. Bouzehouane, C. Deranlot, P. Graziosi, L. Hueso, I. Bergenti, V.A. Dedie, F. Pertroff and A. Fert, Nature Phys. , (2010).
13. S. Ikeda, J. Hayakawa, Y. Ashizawa, Y. M. Lee, K. Miura, H. Hasegawa, M. Tsunoda, F. Matsukura, and H. Ohno, Appl.Phys.Lett. 95, 082508 (2008).
14. S. Sanvito, Nature Physics 6, 562 (2010).
15. V. Dediu, M. Murgia, F.C. Matacotta, C. Taliani and S. Barbanera, Sol. State Commun. 122, 181 (2002).
16. Z. H. Xiong, D. Wu, Z. V. Vardeny, and J. Shi, Nature (London) 427, 821 (2004).
17. G. Giovannetti, S. Kumar, A. Stroppa, J.V.D. Brink and S. Picozzi, Phys. Rev. Lett. 103, 266401 (2009).
18. P. Gütlich, H.A. Goodwin, in Spin Crossover in Transition Metal Compounds (eds. P. Gütlich, H.A. Goodwin) (Springer, 2004).
19. N. Baadji, M. Piacenza, T. Tugsuz, F. Della Sala, G. Maruccio, S. Sanvito, Nature Materials, 8, 813 (2009).
20. M. Diefenbach, K.S. Kim, Angew. Chem. 119, 7784 (2007).
21. A. Droghetti and S. Sanvito, in preparation.
22. J. Hubbard, Proc. Roy. Soc. A 276, 238 (1963).
23. P.W. Anderson in Fizika Dielectrikon edited by G.I. Skanov and K.V. Filippov (Academy of Sciences, Moscow, 1960) p. 290.
24. S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220, 4598 (1983).
25. N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, J. Chem. Phys. 21, 1087 (1953).
26. W. Lee, N. Jean and S. Sanvito, Phys. Rev. B B 79, 085120 (2009).
27. S.J. Xie, K.H. Ahn, D.L. Smith, A.R. Bishop and A. Saxena, Phys. Rev. B 67, 125202 (2003).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters