Modeling the Neutral-Ionic Transition with Correlated Electrons Coupled to Soft Lattices and Molecules

Abstract

Neutral-ionic transitions (NITs) occur in organic charge-transfer (CT) crystals of planar $\pi$-electron donors (D) and acceptors (A) that form mixed stacks ... D^+ρA^−ρD^+ρA^−ρD^+ρA^−ρ ... with variable ionicity $0<\rho <1$ and electron transfer t along the stack. The microscopic NIT model presented here combines a modified Hubbard model for strongly correlated electrons delocalized along the stack with Coulomb intermolecular interactions treated in mean field. It also accounts for linear coupling of electrons to a harmonic molecular vibration and to the Peierls phonon. This simple framework captures the observed complexity of NITs with continuous and discontinuous ρ on cooling or under pressure, together with the stack’s instability to dimerization. The interplay of charge, molecular and lattice degrees of freedom at NIT amplifies the nonlinearity of responses, accounts for the dielectric anomaly, and generates strongly anharmonic potential energy surfaces (PES). Dynamics on the ground state PES address vibrational spectra using time correlation functions. When extended to the excited state PES, the NIT model describes the early (<1 ps) dynamics of transient NIT induced by optical CT excitation with a fs pulse. Although phenomenological, the model parameters are broadly consistent with density functional calculations.

1. Foreword

The physics of strongly-correlated electron systems is characterized by a complex phenomenology that includes multiple competing phases, divergent responses, and anomalous metallic states [1]. In these materials the interplay between correlated electrons and phonons further enhances the complexity both in terms of the appearance of new phases and of amplification of materials responses. Mixed-stack charge-transfer (ms-CT) crystals offer an interesting opportunity to study strongly correlated electrons coupled to molecular vibrations and lattice modes. In these crystals, planar π-conjugated molecules with strong electron donor (D) and acceptor (A) character alternatively form face-to-face stacks. Electrons are delocalized along the stack, leading to fractional charge: ... D^+ρA^−ρD^+ρA^−ρD^+ρA^−ρ ... The most impressive demonstration of strong correlations in these systems is offered by the so-called neutral-ionic phase transition (NIT) [2], a transition from a valence insulator ($\rho <0.5$) to a Mott insulator ($\rho >0.5$), that can be induced by temperature [3], pressure [4], or even light [5,6,7,8,9]. Both continuous and discontinuous NITs are known, and a precise control of chemical pressure drives the transition to very low temperatures, in the regime of quantum phase transitions [10,11,12].

NIT is always accompanied by lattice dimerization [2], eventually leading to ferroelectric behavior [13]. Large dielectric anomalies at NIT are due to the electronic charge fluxes induced by the dimerization mode [14,15]. Other evidence of the intimate entanglement of correlated electrons and phonons is found in several anomalies in infrared (IR) and Raman (R) spectra [2], and in an intense diffuse X-ray scattering signal [16].

Competing electronic ground states on a deformable lattice can result in the simultaneous stability of different phases. Multistability is a rare phenomenon in molecular materials that requires the precise balance between strong and competing interactions [17,18]. Multistability in the close proximity of the discontinuous pressure-induced NIT of TTF-CA has been experimentally observed [19] and extensively discussed in Ref. [20]. Light absorption in multistable systems can populate metastable states, that may be hidden to conventional investigation tools [21], in the so-called photoinduced NIT (PINIT) [5,6,7,8,9]. The possibility to switch a neutral crystal into an ionic ferroelectric metastable state by light absorption is intriguing both from a fundamental perspective and in view of innovative applications.

Ground state properties, including not just charge distribution and lattice structure [20,22], but also other specific and experimentally accessible quantities, like vibrational spectra [23], diffuse X-ray scattering data [16], the dielectric anomaly [14,15], have been satisfactorily described in recent years based on a modified Hubbard model accounting for the inequivalence of lattice sites, and introducing the coupling to a molecular vibration, modulating on-site energies, and to a lattice mode, modulating hopping integrals. The same model, extended to describe the first electronic excited state, captures the basic physics driving the early time dynamics following impulsive excitation at the PINIT [24].

This feature article has several objectives. One is to collect NIT results from papers over many years, and to present in a single, coherent picture our best current modeling of ms-CT stacks. We hope that this article will be a useful entry to newcomers to the field. The NIT literature is covered in depth in the companion articles of this issue. Specifically, the contribution from Masino et al. [25] discusses the vibrational spectroscopy and experimental characterization of NIT in ms-CT crystals. There is considerable interest nowadays in strongly correlated electronic systems, in low-dimensional systems, in frustrated interactions and quantum phases, in nonlinear responses and in electron-vibration coupling. The physics of NIT systems such as ms-CT crystals overlaps these topics and illustrates them.

2. Microscopic Models for NIT

2.1. An Introductory Model: The Mulliken dimer

The Mulliken dimer illustrates CT in a simple two-level system. The neutral state $|DA\rangle$ describes an electron donor D and acceptor A molecule, both with closed shells. When an electron is transferred from D to A, the open shell system can be in a singlet or triplet state. Only the singlet ionic state $|D^{+}{A}^{-}\rangle$ is mixed with the neutral state by $-\sqrt{2}t=\langle DA|\widehat{H}|D^{+}{A}^{-}\rangle$, where t is the hopping integral between D and A. The energy difference between the ionic (I) and the neutral (N) states, $2z=\mathcal{I}-\mathcal{A}+V$ is a model parameter accounting for the ionization potential of D ($\mathcal{I}$), the electron affinity of A ($\mathcal{A}$) and the Coulomb interaction V between the ions D⁺ and A⁻. The ground state is:

$$|G\rangle =\sqrt{1-\rho }|DA\rangle +\sqrt{\rho }|D^{+}{A}^{-}\rangle .$$

(1)

Accordingly, the ground state may be described as $D^{\rho +}{A}^{\rho -}$ with $\rho$ measuring the amount of CT. The ionicity is fully defined by the dimensionless parameter $\eta =z/t$:

$$\rho =\frac{1}{2}\left(1-\frac{\eta }{\sqrt{{\eta }^2+2}}\right)$$

(2)

The left panel of Figure 1 shows $\rho \left(\eta \right)$: as expected, $\rho =1/2$ for $z=0$. The ground state is neutral with $\rho <0.5$ for $z>0$ and becomes ionic, $\rho >0.5$, for negative z.

As an instructive toy model of NIT, first introduced in Ref. [26], we consider a crystal with one DA dimer per unit cell and electrostatic (Coulomb) interactions between dimers. Delocalization (t) within the dimer favors $\rho \sim$ 0.5 while attractive Coulomb interactions are maximized at $\rho \sim$ 1. Competing interactions are a recurrent theme in the context of NIT. In a mean field approach we can reduce the problem of interacting dimers to that of an isolated dimer in the electrostatic potential of the surrounding dimers, each with the average ionicity $\rho$. The mean field Hamiltonian for an embedded dimer is the same as for an isolated dimer. The parameter z for electron transfer in an isolated pair is simply renormalized to $z-{ϵ}_c\rho$, where ${ϵ}_c=V-\mathcal{M}$ and $\mathcal{M}$ is the Madelung (electrostatic) energy of the fully ionized lattice. Positive values of ${ϵ}_c$ correspond to attractive interactions, favoring charge separation. The self-consistent problem can be easily solved extracting from the universal $\rho \left({\eta }_{eff}\right)$ a new $\rho \left(\eta \right)$ curve for each ${ϵ}_c$ value, where ${\eta }_{eff}=(z-{ϵ}_c\rho )/t$, as shown in the left panel of Figure 1 for a few ${ϵ}_c$ values.

For increasing attractive electrostatic interactions, the $\rho \left(\eta \right)$ curves in Figure 1 become steeper and steeper until $\rho \left(\eta \right)$ becomes a multivalued, S shaped curve for ${ϵ}_c/t>4/\sqrt{2}$. Much as in the gas-liquid transition (van der Waals equation), the NIT of non-overlapping DA dimers becomes first order (discontinuous) beyond the critical point ${ϵ}_c/t=4/\sqrt{2}$. In this case, the system presents two stable solutions up to the spinodals and the region of the $\rho \left(\eta \right)$ curve with positive slope corresponds to unstable, thermodynamically forbidden solutions. The Maxwell construction excludes the unphysical region by drawing a vertical line crossing the curve at $\rho =0.5$, discriminating between the thermodynamically stable phase and metastable solutions.

The mean field treatment of a crystal of Mulliken dimers captures a central feature of NIT, the possibility of continuous or discontinuous $\rho$ depending on ${ϵ}_c/t$. In this respect the electronic properties of ms-CT stacks mainly require a numerical analysis of the dependence of $\rho$ on model parameters. On the other hand, DA dimers or dimer approximations to ms-CT stacks do not have inversion symmetry at D and A sites, which is a crucial property for the dimerization instability and for the spectroscopy of ms-CT crystals.

2.2. The Modified Hubbard Model

The modified Hubbard Model (MHM) is a 1D generalization of the Mulliken dimer, although its actual origin was a limiting case of the Hubbard model with site energies $\pm \mathsf{\Delta }$ at A and D sites, and two electrons per DA pair (half-filling) [27]. Hence we consider a linear chain of equally-spaced molecules (regular stack, as opposed to dimerized stacks that will be introduced in Section 2.4) with intermolecular hopping t and Hubbard on-site repulsion U. The Hamiltonian reads:

$$H=\mathsf{\Delta }\sum _i{(-1)}^i{\widehat{n}}_i+U\sum _i{\widehat{n}}_{i,\alpha }{\widehat{n}}_{i,\beta }-t\sum _{i,\sigma }\left({\widehat{a}}_{i,\sigma }^{†}{\widehat{a}}_{i+1,\sigma }+H.c.\right)$$

(3)

where H.c. is Hermitian conjugate, i runs on the N sites of the stack, ${\widehat{a}}_{i,\sigma }$ and ${\widehat{a}}_{i,\sigma }^{†}$ are the annihilation and creation operators for an electron with spin $\sigma$ on site i, ${\widehat{n}}_{i,\sigma }={\widehat{a}}_{i,\sigma }^{†}{\widehat{a}}_{i,\sigma }$ counts the electrons with spin $\sigma$ on site i. The computational effort needed to attack the problem can be reduced by adopting the so-called reduced basis with three rather than four states per site on excluding high-energy doubly ionized sites, D²⁺ or A²⁻ [28,29]. The reduced basis approximation corresponds to letting U and $\mathsf{\Delta }$ go to infinity while maintaining a finite difference $2\mathsf{\Gamma }=2\mathsf{\Delta }-U$ between the ionization potential of D and the electron affinity of A. The modified Hubbard Hamiltonian in the reduced basis reads:

$$H_{MHM}=\mathsf{\Gamma }\sum _i{(-1)}^i{\widehat{n}}_i-t\sum _{i,\sigma }\left({\widehat{a}}_{i,\sigma }^{†}{\widehat{a}}_{i+1,\sigma }+H.c.\right).$$

(4)

Two remarks are in order here: first, Hamiltonian 4 apparently contains only one-electron terms but it actually accounts exactly for strong correlations, according to the $U\to \infty$ limit. In this limit, D sites (odd i) can only have ${\widehat{n}}_i$ = 1 or 2, the A sites (even i) have ${\widehat{n}}_i$ = 0 or 1. Second, the reduced-basis leads to fewer model parameters, allowing to explore the phase diagram of ms-CT salts as function of the $\mathsf{\Gamma }/t$ ratio, while preserving the rich physics of electronic correlations.

The fully neutral state $|N\rangle =|...DADADA...\rangle$, obtained for $\mathsf{\Gamma }/t\to \infty$, is a singlet state with a pair of electrons residing on the HOMO of D, the LUMO of A being empty. The fully ionic state $|I\rangle =|...D^{+}{A}^{-}{D}^{+}{A}^{-}{D}^{+}{A}^{-}...\rangle$ ($\mathsf{\Gamma }/t\to -\infty$) is vastly more complicated due to the twofold spin-degeneracy of radical ions. There are even more states containing both neutral molecules and radical ions. An analytical solution is of course not possible, but valence bond methods [30] can be used to get the numerically exact ground state for finite site systems up to ∼20 sites. Figure 2 illustrates finite size effects on the ionicity

$$\rho =\langle (1/N)\sum _j{\widehat{\rho }}_j\rangle ,$$

(5)

where the operators ${\widehat{\rho }}_j$ measure the absolute value of the charge on D (${\widehat{\rho }}_j=2-{\widehat{n}}_j$, odd j) and A (${\widehat{\rho }}_j={\widehat{n}}_j$, even j) sites. The $\rho (\mathsf{\Gamma }/t)$ curves in Figure 2 are calculated for linear chains with periodic boundary conditions and increasing number of sites from 6 to 16. Systems with $4n$ and $4n+2$ sites converge to the same limit from opposite directions [31], setting narrow uncertainties on the $N\to \infty$ extrapolation of $\rho$. Apart from a small region at intermediate ionicity, 14 and 16 sites results are superimposed in the scale of the figure. The black curve that represents our best estimate of the infinite chain is just the average of $N=14$ and 16 results.

We underline that $4n$ chains have discontinuous $\rho$ at the point where the ground state wavefunction changes from symmetric to antisymmetric with respect to the inversion center at each molecular site of the regular chain. $4n+2$ chains do not show any symmetry crossover nor discontinuity. Indeed the discontinuity is washed out in the thermodynamic limit. A remnant of the discontinuity remains in the slope variation of the $\rho \left(\mathsf{\Gamma }\right)$ curve, that points to a crossover that moves away from $\rho =0.5$ of the dimer towards ${\rho }_c=0.686$ at ${\mathsf{\Gamma }}_c/t=-0.666$. At this critical point there is a vanishing gap between symmetric and antisymmetric states with respect to inversion symmetry and the gap stays closed in the I regime of the extended stack [31,32].

Additional important information on the nature of the ground state can be gained thanks to the modern theory of polarization of extended systems as proposed by Resta [33]. The complex quantity

$$Z=\langle \psi |\exp (i\frac{2\pi \widehat{M}}{N})|\psi \rangle ,$$

(6)

is central to Resta’s approach, with $|\psi \rangle$ representing the ground state wavefunction, i the imaginary unit and $\widehat{M}={\sum }_j{(-1)}^{j+1}j{\widehat{\rho }}_j$ the dipole moment operator of the regular open-boundary chain. Equation (6) allows the definition of the electric polarization, which is proportional to the phase of Z and vanishes by symmetry in regular stacks. The polarizability in systems with periodic boundary conditions is discussed in Section 2.6. Figure 2b shows the polarizability calculated for systems of increasing size and periodic boundary condition. A very sharp peak in the polarizability is observed at the crossover point, marking the presence of delocalized charges. To better understand the nature of the very special crossover point, panel c in the same figure shows the absolute value of Z. As discussed by Resta and Sorella [34], the inverse of $\left|Z\right|$ measures the wavefunction delocalization, so that a vanishing $\left|Z\right|$ corresponds to a metallic state of the extended system. Data in Figure 2 clearly point to a metallic state at the critical point at ${\mathsf{\Gamma }}_c=-0.666t$.

The critical point ${\mathsf{\Gamma }}_c$ is associated with the charge degrees of freedom of the MHM and has major consequences for NIT. The MHM also has a magnetic critical point ${\mathsf{\Gamma }}_m$ due to spin degrees of freedom [35], a feature that was not recognized previously. The infinite chain is diamagnetic with a finite energy gap between the singlet ground state and the lowest triplet state for $\mathsf{\Gamma }>{\mathsf{\Gamma }}_m=-0.91$. The singlet-triplet gap is zero for $\mathsf{\Gamma }\le {\mathsf{\Gamma }}_m$, and the infinite chain is paramagnetic with finite magnetic susceptibility at 0 K. Quantum phase transitions related to the opening of a singlet-triplet gap have been extensively treated theoretically in spin-1/2 chains with frustrated exchange interactions [36,37,38] and in half-filled extended Hubbard models [39,40,41]. The magnetic critical point has no role so far in NIT models, mainly due to structural instabilities that are discussed below.

2.3. Electrostatic Interactions

Since NIT is induced by intermolecular interactions that change due to volume variations upon cooling or under pressure, several extensions of the MHM are called for. The first is to introduce intermolecular Coulomb interactions and supplement Hamiltonian (4) with:

$$H_V=\frac{1}{2}\sum _{i,j}{V}_{ij}{\widehat{\rho }}_i{\widehat{\rho }}_j$$

(7)

where $i,j$ indices run on all sites of the 3D crystal and $V_{ij}$ measures the electrostatic interaction between fully charged sites at $i,j$. Attacking a problem of correlated electrons in 3D is a formidable task, so that simplified approaches are needed. As originally proposed by Soos and Mazumdar [29] and extensively validated by Painelli and Girlando [22], electrostatic intermolecular interactions can be dealt within the mean field approximation. As shown above for a crystal of Mulliken dimers, this reduces the 3D problem to an effective 1D problem. Indeed, as for intermolecular electrostatic interactions, we view the 1D stack as an array of DA dimers. For each dimer, $2z=2\mathsf{\Gamma }+V$ measures the energy required to transfer an electron from D to A. All other intermolecular electrostatic interactions within the stack and all interstack interactions are treated in the mean field approximation, leading to the following Hamiltonian [20]:

$$\begin{array}{ccc}\hfill H_{mf}& =& z_{eff}\sum _i{(-1)}^i{\widehat{n}}_i-t\sum _{i,\sigma }\left({\widehat{a}}_{i,\sigma }^{†}{\widehat{a}}_{i+1,\sigma }+H.c.\right)\hfill \\ \hfill & +& \frac{N}{2}{ϵ}_c{\rho }^2-N{ϵ}_c\rho \hfill \end{array}$$

(8)

where the effective on-site energy alternation $z_{eff}=z-{ϵ}_v\rho$ self-consistently depends on the ionicity $\rho =\langle {\widehat{\rho }}_i\rangle$ through a Coulomb correlation strength defined by ${ϵ}_c=V-2\mathcal{M}$, where $\mathcal{M}=\frac{1}{2N}{\sum }_{ij}{V}_{ij}$ measures the Madelung energy of a 3D lattice of D⁺ and A⁻ ions.

We can now follow a similar strategy as described in Section 2.1 for the analytical DA dimer: the $\rho (\mathsf{\Gamma }/t)$ curve extrapolated for the infinite chain coincides with the universal $\rho (z_{eff}/t)$ curve relevant to the mean field Hamiltonian. From this curve, it is possible to extract, for each ${ϵ}_c$ value, the relevant $z=z_{eff}+{ϵ}_c\rho$ to get the $\rho \left(z\right)$. The $\rho \left(z\right)$ curve shown in the right panel of Figure 1 becomes steeper and steeper with increasing ${ϵ}_c$, and for ${ϵ}_c$ larger than a critical value, estimated as 1.8t, it becomes S-shaped, with a region of positive slope corresponding to unstable states. In other terms, electrostatic intermolecular interactions and, as we will see in the following, the on-site vibrational coupling cooperate towards bistability in ms-CT crystals, making states of intermediate ionicity unattainable, and opening a bistability window where both quasi-neutral and quasi-ionic states are accessible.

The optical CT transition is the lowest electronic excitation of ms-CT crystals in the N phase, experimentally located around $E_{CT}\sim 1$ eV and polarized along the stack axis. The Mulliken dimer has $E_{CT}=2\sqrt{2}t$ at $\rho =1/2$, while the MHM has $E_{CT}=0$ at ${\rho }_c=0.686$. The optical gap opens up again in the I lattice leading to an almost symmetrical $E_{CT}\left(\rho \right)$ behavior across the critical ${\rho }_c$ [42]. The critical point, characterized by a vanishing charge gap, is however never attained experimentally, not even in case of a smooth $\rho$ behavior: as it will be discussed below, soft lattices always dimerize before reaching the critical point.

2.4. Dimerization Instability of the Soft Lattice

The Peierls instability of 1D systems has been extensively discussed in connection with polyacetylene, (CH)_x [43], as has the spin-Peierls instability of chains with exchange interactions between sites with $s=1/2$ spins [44,45]. Most NIT systems, including the best characterized ones, have planar D and A molecules and molecular ions with high symmetry (typically $D_{2h}$). Molecular inversion centers are retained in crystals with regular stacks, but are lost upon lattice dimerization. Since the ionic side has radical ions, a generalized Peierls transition is expected and found in ionic chains.

The Peierls mode $\delta$, describing the relative shift of D and A sublattices, is the only lattice mode we consider. At variance with (CH)_x, this optical mode is located at the Brillouin zone center in mixed stacks with two molecules, D and A, in the unit cell. This has important implications for vibrational spectroscopy. Following the Su-Schrieffer-Heeger model for (CH)_x [46], we allow the molecular sites to displace around the equilibrium position of the regular stack and retain only the linear term in the Taylor expansion of t with intermolecular distance. Accordingly, a new term enters the Hamiltonian $H_{MHM}+H_V+H_d$ with

$$H_d=-t\sum _{i,\sigma }{(-1)}^i\delta \left({\widehat{a}}_{i,\sigma }^{†}{\widehat{a}}_{i+1,\sigma }+H.c.\right)+\frac{N{t}^2}{2{ϵ}_d}{\delta }^2$$

(9)

where the lattice relaxation energy ${ϵ}_d$ measures the strength of the electron-phonon coupling and the last term accounts for a purely harmonic potential for lattice vibrations. We work in the adiabatic (Born-Oppenheimer) approximation to neglect the kinetic energy of the lattice mode.

The ground-state wavefunction of $4n$ regular rings belongs to the antisymmetric subspace with respect to inversion, whereas $4n+2$ rings have a symmetric ground state: this suggests that in the thermodynamic limit a vanishing gap has to be found between the two subspaces. As it was recognized in the early days of NIT, this leads to an unconditional instability of the ionic regular stack towards dimerization [31,47,48,49,50]. The resulting instability is often called a spin-Peierls instability, but a generalized Peierls instability is more appropriate. TTF-BA and M₂P-TCNQF₄ have highly ionic ($\rho \sim$ 0.9) regular ms-CT stacks at ambient temperature and dimerize at $T_{SP}=53$ K [51,52] and 122 K, [53,54] respectively. The charge degrees of freedom participate to a small extent [55]. The magnetic susceptibility in either case indicates a substantial thermal population of excited states that drives the transition. Generalized Peierls transitions typically have variable intermediate $\rho$ in systems with large singlet-triplet gap and weak paramagnetism even at room temperature [55].

Mixed stack CT crystals show a rich physics governed by the interplay of the charge-instability, as driven by electrostatic interactions, and the generalized Peierls instability. Some considerations are in order about the nature of possible transitions in experimental systems. Starting from the neutral side and a regular stack at ambient conditions, decreasing temperature or increasing pressure leads to lattice compression and hence to an increased electrostatic energy and ionicity. Systems with t small with respect to ${ϵ}_c$ will be driven towards a discontinuous phase transition with an abrupt ionicity increase. Since the I-phase is unconditionally unstable towards dimerization at 0 K, if the transition temperature is sufficiently low, we expect a discontinuous jump from a neutral regular phase to an I-dimerized phase, as observed in TTF-CA and TTF-BrCl₃Q (see the companion paper [25]). A continuous NIT is expected in systems with a larger $t/{ϵ}_c$. Upon lattice compression, a soft regular stack with small ρ will dimerize in a continuous transition before reaching the metallic point. The ionicity may increase after the transition, possibly driving the system towards large ρ, but without any additional phase transition. This is the behavior observed for DMTTF-CA and DMeTTF-2,5Br₂Cl₂Q [25]. Below we will discuss in some detail TTF-CA and DMTTF-CA as prototypical examples of discontinuous and continuous NIT.

2.5. Soft Molecules and the Complete NIT Model

The molecular sites have internal degrees of freedom related to molecular vibrations. Totally-symmetric modes in the molecular point group modulate on-site energies, and then couple to delocalized electrons (e-mv or Holstein coupling). Again we account for a single vibration of either D or A, and introduce the corresponding Brillouin zone-center mode q that describes the in-phase oscillations of all D (or of all A) molecules. The coordinate q linearly modulates the difference of site energy between D and A sites and is the final ingredient of the microscopic NIT model. For the sake of clarity and for future reference we write here the full Hamiltonian:

$$\begin{array}{ccc}\hfill H_{msCT}& =& (\mathsf{\Gamma }+q)\sum _i{(-1)}^i{\widehat{n}}_i-t\sum _{i,\sigma }\left[1+{(-1)}^i\delta \right]({\widehat{a}}_{i,\sigma }^{†}{\widehat{a}}_{i+1,\sigma }+H.c.)+\frac{1}{2}\sum _{i,j}{V}_{i,j}{\widehat{\rho }}_i{\widehat{\rho }}_j\hfill \\ \hfill & +& \frac{N{t}^2}{2{ϵ}_d}{\delta }^2+\frac{N}{2{ϵ}_v}{q}^2.\hfill \end{array}$$

(10)

The last term accounts for the harmonic elastic energy related to q, with ${ϵ}_v$ measuring the small-polaron binding energy. In other contexts, this energy is also called the Franck-Condon relaxation energy, and, divided by the vibrational frequency, it corresponds to the Huang-Rhys factor.

As illustrated in the next sections, the above Hamiltonian represents the most general model that we will discuss. It can be numerically diagonalized for finite systems with different $\delta$ and q values to obtain the potential energy surface, as discussed in Section 3. However the problem is numerically tractable only for 1D chains, while intermolecular electrostatic interactions, the last term in the first line of Equation (10), lead to a 3D problem. In this respect, the mean field approximation, discussed above, is very useful, as it reduces the 3D problem to an effective 1D problem with renormalized on-site energies (cf Equation (8)).

Before closing this section, we illustrate a simple approach to e-mv coupling to elucidate its role on the nature of NIT. In line with the adiabatic approximation, we may write the electronic Hamiltonian substituting q in the above expression with its equilibrium value: $q_{eq}={ϵ}_v(1-\rho )$. If we consider a regular chain, $\delta =0$, the electronic Hamiltonian in the mean field approximation can be written exactly as in Equation (8) but with $z_{eff}\to z_{eff}+q_{eq}=z-({ϵ}_c+{ϵ}_v)\rho +{ϵ}_v$. Apart from a rigid shift of the effective z parameter by ${ϵ}_v$, the most important result here is that e-mv coupling adds to electrostatic interactions in favoring discontinuous NIT.

2.6. Polarization and Polarizability of the Modified Hubbard Model

Polarization and polarizability play a central role in the phenomenology of ms-CT crystals as these are key quantities for the dielectric anomaly at NIT, vibrational spectra and the ferroelectric polarization in dimerized phases. We recall here some general results on the dielectric properties of the MHM that will be useful in the following.

The modern theory of polarization in solids [33,56,57,58] provides a rigorous definition of the electrical polarization in terms of a Berry phase, overcoming the very fundamental issues arising from the ill-defined dipole moment operator in a system with periodic boundary conditions. The electronic polarization of a many-body quantum system is defined as [59]

$$P_{el}=\frac{1}{2\pi }\mathrm{Im}[ln\langle \psi |\exp (i\frac{2\pi \widehat{M}}{N})|\psi \rangle ]=\frac{1}{2\pi }\mathrm{Im}[ln\left(Z\right)],$$

(11)

where the quantity Z and the dipole operator $\widehat{M}$ of an open-boundary 1D chain have been introduced in Equation (6). $P_{el}$ is expressed in units of $ea$ per site, being e the elementary charge and a the lattice constant. It is worth noticing that $P_{el}$, being related to the phase of the wavefunction $\psi$ and not to its amplitude, is defined modulo one in current units, the so-called polarization quantum. Indeed only the variations of $P_{el}$ are physically meaningful, and in the following we will be concerned with its derivatives with respect to vibrational coordinates, related to IR absorption intensity of molecular and lattice modes. We will also discuss $P_{el}$ derivatives with respect to an applied electric field, i.e., the polarizability, as well as the $\mathsf{\Delta }{P}_{el}\left(\delta \right)=P_{el}\left(\delta \right)-P_{el}\left(0\right)$ difference, corresponding to the spontaneous polarization of a dimerized stack.

To illustrate the dielectric properties of the MHM we step back to a minimal NIT Hamiltonian $H=H_{MHM}+H_d$ (see Equations (4) and (9)) that only depends on the two parameters $\mathsf{\Gamma }$ and $\delta$, with $t=1$ setting the energy scale. $\mathsf{\Gamma }$ can implicitly accounts for the effect of electrostatic interactions and intramolecular vibrations, as discussed in Section 2.3 and Section 2.5. Figure 3 shows $P_{el}$, calculated with Equation (11), and the $\mathsf{\Delta }{P}_{el}(\delta ,\mathsf{\Gamma })=P_{el}(\delta ,\mathsf{\Gamma })-P_{el}(0,\mathsf{\Gamma })$ difference, computed for every fixed value of $\mathsf{\Gamma }$, in panels a and b, respectively. In the regular stack $P_{el}$ has a discontinuity from 0 (N phase) to 0.5 (I phase) at ${\mathsf{\Gamma }}_c\approx -0.666t$. $P_{el}$ and $\mathsf{\Delta }{P}_{el}$ are therefore equivalent only in the N system, while the discontinuity in $P_{el}$ across $\delta =0$ for $\mathsf{\Gamma }<{\mathsf{\Gamma }}_c$ is only apparent as it exactly amounts to one quantum of polarization.

$\mathsf{\Delta }{P}_{el}$, the physically meaningful quantity, is an odd function of the dimerization $\delta$ that vanishes in the regular chain, as expected on the basis of elementary symmetry arguments.The magnitude of $\mathsf{\Delta }{P}_{el}$ increases upon approaching the metallic N-I boundary (located at $\mathsf{\Gamma }/t=0.666$ for the regular chain, see Figure 2). $\mathsf{\Delta }{P}_{el}$ changes its sign around $\delta =0$ with a smooth and almost linear behavior far from NIT, while the variation becomes sharp with a marked sigmoid character close to the transition (see Figure 3c). The derivatives ${(\partial \mathsf{\Delta }{P}_{el}/\partial \delta )}_0$ diverge at the NIT of the rigid chain conferring giant IR intensity to the Peierls mode [14,15]. Charge fluxes induced by the lattice distortion are amplified at the N-I crossover as a result of the enhanced electron delocalization and of the instability towards the formation of a bond charge density wave.

It is worth remarking that ${(\partial \mathsf{\Delta }{P}_{el}/\partial \delta )}_0$ derivatives diverge with opposite sign on the N and I side, and $\mathsf{\Delta }{P}_{el}$ undergoes a discontinuous change of sign at finite $\delta$ (see Figure 3d). These are the signatures of the different nature of the two insulating phases and the sign reversal of dynamical charges (electronic charge flow) has been discussed as a general consequence of strong electronic correlations in the Mott-insulator phase in systems of partial covalent-ionic character [34].

The definition of $P_{el}$ in Equation (11) does not take into account the physical dimerization of the lattice, i.e., the unequal molecular spacing. The definition of the open-chain dipole operator $\widehat{M}$ can be generalized to account for this effect, allowing the quantification of the total polarization that includes also a contribution from frozen charges at molecular ions in a dimerized lattice $\mathsf{\Delta }{P}_{ion}$. It has been shown that this term is much smaller in magnitude than that of dynamical electron charges [15] and that $\mathsf{\Delta }{P}_{el}$ and $\mathsf{\Delta }{P}_{ion}$ have opposite sign in the strongly-correlated I phase [34]. These theoretical results found an elegant confirmation in the experiment by Kobayashi et al., who showed the spontaneous polarization of the I-dimerized phase of TTF-CA points in the opposite direction and is ∼23 times larger in magnitude than what expected for frozen charge $\pm \rho$ at the molecular sites of the dimerized stack [60]. TTF-CA represents a notable example of a colossal Born atomic charges of opposite sign with respect to the respective molecular valence.

Data in Figure 3 can be compared to the spontaneous polarization measured in the ferroelectric phases of TTF-CA [60] and TTF-BA [52]. The polar phase of TTF-CA is located in Figure 3b at $\delta =0.28$ and $\mathsf{\Gamma }=-0.5$, which correspond to the experimental ionicity in the I phase, $\rho =0.6$, and to the dimerization amplitude obtained from density functional theory calculations performed on the 40 K crystal structure (see Section 6) [61]. For these values, we obtain $\mathsf{\Delta }{P}_{el}=8.6 \mathsf{\mu }$C/cm² which compares well with the experimental value at 59 K of 6.3 μC/cm² [60]. TTF-BA is an almost fully-I system with $\rho =0.95$ at all temperatures that dimerizes below $T_c=53$ K. For the crystal structure at 20 K we estimate $\delta =0.14$, leading to $\mathsf{\Delta }{P}_{el}=0.7\mathsf{\mu }$C/cm², a value that is in line with the experimental 0.2 μC/cm² [52]. The relative magnitude of the polarization in TTF-CA and TTF-BA testify to the greater electron delocalization in the system undergoing NIT, where the lattice dimerization drives large charge flows.

The electronic polarizability is the first derivative of the polarization with respect to a uniform external electric field. Although a static electric field F is not compatible with periodic boundary conditions, Nunes and Gonze [62] showed how to combine P and F by introducing the energy of a polarized system as the functional

$$E(F,{\psi }_F)=\langle {\psi }_F|H_0|{\psi }_F\rangle -NF{P}_{el}\left({\psi }_F\right),$$

(12)

where $H_0$ is the Hamiltonian at $F=0$ and ${\psi }_F$ is the polarized wavefunction at finite field. The minimization of $E(F,{\psi }_F)$ through a self consistent iterative scheme allows the calculation of ${\psi }_F$ that can be then used to compute the induced polarization $P_{el}\left(F\right)$ via Equation (11) [15,58,62]. For the calculation of the linear susceptibility ${\alpha }_{el}$, it is sufficient to correct ${\psi }_F$ at the first order, obtained on iterating Equation (12) only once [15].

The electronic polarizability of the MHM, shown in Figure 4, diverges at the metallic N-I boundary of the infinite regular chain (notice the logarithmic scale; see also Figure 2b), signaling the vanishing charge gap and the enhanced electronic delocalization. The derivatives of ${\alpha }_{el}$ over the vibrational coordinates give the cross section of non-resonant Raman scattering. Vibrational dynamics exploring this singularity are therefore expected to show nonlinear amplification of the Raman signal, as experimentally observed in DMTTF-CA [23] (see Section 4). The electronic polarizability also contributes to the large peak observed in the dielectric constant in close proximity of either continuous or discontinuous NIT [13,63]. However this contribution has been shown to be a marginal one in both continuous and discontinuous transitions, where dimerization shifts the system away from the unstable metallic point, largely suppressing ${\alpha }_{el}$ (see Figure 4) [14,15]. To rationalize the dielectric anomaly observed at either continuous or discontinuous NIT, it is actually crucial to include the vibrational contribution, that is related to the IR intensity of the Peierls mode and to the electronic charge flux induced by the lattice distortion [14,15].

3. Potential Energy Surfaces

3.1. Discontinuous NITs

Now that the elements for the modeling of ms-CT crystals have been summarized, we are in the position to discuss the NIT on the basis of the complete model Hamiltonian in Equation (10), accounting for inter-site electrostatic interactions at the mean field level, and for the coupling to one effective intramolecular mode q and the Peierls lattice phonon $\delta$. The following discussion will be centered around the adiabatic potential energy surface (PES) of the ground state $E(q,\delta )$. PESs bring important information on the competing valence and structural instabilities at NIT that result in either continuous or discontinuous transitions. We recall that molecular and lattice vibrations enter the Hamiltonian as strictly harmonic modes, and anticipate that coupling to electrons leads to strongly anharmonic PES, especially for model parameters close to NIT [20].

We begin with the discontinuous NIT, similar to that occurring with temperature in TTF-CA, that we induce by acting on electrostatic interactions by increasing V, mimicking the effect of lattice contraction. Reliable estimates of model parameters are known from previous studies [22,42,64]. The left panels of Figure 5 report the dimerization ${\delta }_{eq}$ and ionicity ${\rho }_{eq}$, calculated at the equilibrium geometry of the system. The transition is clearly discontinuous: ${\rho }_{eq}$ smoothly increases from ∼0.25 to ∼0.4 with increasing V and then abruptly jumps toward I values (${\rho }_{eq}>0.6$). The stack remains regular (${\delta }_{eq}=0$) on the N side and dimerizes in the I regime (${\delta }_{eq}\sim$ 0.12). The dotted lines in Figure 5a,b mark the region of bistability close to the transition where two non-equivalent minima are found. Bistability is a typical feature of discontinuous transitions.

Left panels of Figure 5 show the ground state PESs $E(\delta ,q)$ calculated for three different V values marked by arrows in Figure 5a. Figure 5c ($V=-2.26t$), corresponds to a N regular system, whose equilibrium is found at ${\delta }_{eq}=0$ (blue circle). The PES in panel d ($V=-2.42t$) is instead characteristic of an I dimerized stack, which, unconditionally unstable to dimerization, presents two equivalent minima at finite $\pm {\delta }_{eq}$ (red circles), corresponding to I dimerized phases with opposite polarization. The PES of a system within the bistability region ($V=-2.34t$), shown in Figure 5e,f, is instead characterized by three minima, corresponding to coexisting N-regular (blue circle) and I-dimerized (red circles) phases.

3.2. Continuous NIT

The features of NIT are strongly affected by slight changes in model parameters [20]. Figure 6 shows the results for a continuous transition, obtained for the same parameters as used in Figure 5, except for $\mathcal{M}=-1.1t$. Panels a and b of Figure 6 show that now dimerization occurs on the N side at ${\rho }_{eq}\sim 0.4$, where a steep but continuous increase of the ionicity occurs. This behavior is similar to that of DMTTF-CA with temperature, where however the ionicity remains always lower than ∼0.45, while we calculate ${\rho }_{eq}\sim$ 0.7 at $V=-3.2t$. This disagreement is not too surprising considering the simplicity of the model.

Figure 6c–f show the PESs across the continuous NIT, calculated for the V values marked by arrows in panel a. The PES in Figure 6c ($V=-2.6t$) corresponds to a N-regular system, with a single minimum at ${\delta }_{eq}=0$, while that in panel d ($V=-3.1t$) describes an I dimerized stack with ${\delta }_{eq}\sim 0.12$. At variance with the discontinuous transition discussed above, in this case the PES smoothly evolves from a single minimum in the N phase to a double minimum in the I phase. The PES in Figure 6e,f, calculated for $V=-2.94t$, represents a system just at the continuous N-I interface.

Continuous NITs such as that of DMTTF-CA are driven by the Peierls dimerization instability that develops through the softening of the Peierls mode, whose harmonic frequency goes to to zero at the critical point. The corresponding PES is flat along $\delta$ (at the equilibrium q) and is characterized by banana-shaped isolines (see Figure 6e) as a result of the mutual coupling between the lattice and the molecular mode mediated by electrons. The softening of the Peierls mode is only incipient in discontinuous NITs driven by the valence instability, as for instance in TTF-CA. In these cases a full softening is prevented by the sudden ionicity jump bringing the system to an I phase of finite dimerization (see Figure 5) [15,20].

4. Vibrational Spectra at the Continuous NIT of DMTTF-CA

Correlated and mobile electrons in soft lattices manifest in large anharmonicity and electronic charge flows driven by vibrational motion in ms-CT crystals undergoing NIT. It is, therefore, not surprising that the interpretation of vibrational spectra is rather complex but also that it offers detailed information about the entanglement of electronic and vibrational degrees of freedom. The theoretical modeling of vibrational spectra can follow different approaches. In the most common harmonic approximation, vibrational spectra are calculated from the electronic properties at the PES equilibrium positions. The harmonic frequency is indeed related to the PES curvature, while IR and Raman intensities are proportional to the square of the first derivatives with respect to relevant coordinates of the polarization and polarizability, respectively. While results based on the harmonic approximation provided key insight on the vibrational properties at NIT (e.g., the softening of the Peierls mode and its huge IR intensity [14,15]), the strongly anharmonic PESs in Figure 5 and Figure 6 suggest relaxing the harmonic approximation. In the following, we will rely on a time-correlation function (TCF) approach that makes use of the anharmonic PESs, the explicit dependence of $P_{el}$ and ${\alpha }_{el}$ on vibrational coordinates (see Figure 3 and Figure 4), and finite temperature [23].

To this aim, we exploit the fluctuation-dissipation theorem of linear response theory to express optical spectra in terms of the Fourier transform of relevant time-correlation functions (TCFs). For IR absorption spectroscopy the imaginary (dissipative) part of the dielectric constant of the material is obtained from the Fourier transform of the TCF of the dielectric polarization $P_{el}$ [65]:

$${\epsilon }^{″}\left(\omega \right)=\frac{1-e^{-\frac{\hslash \omega }{k_{B}T}}}{4{\epsilon }_{0}V\hslash }{\int }_{-\infty }^{\infty }{e}^{-i\omega t}\langle P_{el}\left(0\right)P_{el}\left(t\right)\rangle dt$$

(13)

where ℏ is the reduced Planck constant, $k_B$ is the Boltzmann constant, ${\epsilon }_0$ is the vacuum permittivity, V is the volume. Similarly, the TCF of the polarizability ${\alpha }_{el}$ is related to the intensity of the non-resonant Raman scattering [65]:

$$I_R\left(\omega \right)={\int }_{-\infty }^{\infty }{e}^{-i\omega t}\langle {\alpha }_{el}\left(0\right){\alpha }_{el}\left(t\right)\rangle dt.$$

(14)

While Equations (13) and (14) are exact, resulting from a quantum mechanical treatment of the interaction between light and matter in the weak coupling regime, approximations are introduced in the calculation of TCFs. We adopt a semiclassical approach to TCFs, with the vibrational motion governed by classical equations of motion in the quantum mechanical potential obtained from the model Hamiltonian. The two effective vibrational coordinates q and $\delta$, defined in Section 2.4 and Section 2.5, are considered and polarization and polarizability are evaluated along many classical trajectories to compute the respective correlation functions [23]. Since our electronic model is strictly 1D, the components of the electronic polarization and polarizability that are probed correspond to those along the stack axis $\overrightarrow{a}$ (see Section 2.6). Equation (13) therefore targets the IR absorption for incident light polarized along $\overrightarrow{a}$, while Raman spectra calculated via Equation (14) correspond to those with both incident and scattered light along the same direction. Raman spectra in the molecular vibrations region are usually taken with incident and scattered light polarized perpendicular to the stack axis, probing instead the molecular polarizability. To simulate these spectra we approximate the molecular polarizability by a linear expansion on q.

We explicitly focus on the vibrational spectra across the continuous NIT of DMTTF-CA. To allow for a direct comparison with experimental spectra, the unperturbed harmonic frequencies of the Peierls and molecular modes are set to ${\omega }_P^0=90{\mathrm{cm}}^{-1}$ and ${\omega }_H^0=1540{\mathrm{cm}}^{-1}$, respectively, with coupling constants ${ϵ}_d\approx 0.4t$ and ${ϵ}_v=0.7t$. The small value of ${ϵ}_v$ is due to the fact here we focus on a specific mode, namely the one responsible for the band observed at about 1450 ${\mathrm{cm}}^{-1}$ in the room temperature Raman spectra of DMTTF-CA [66]. The contribution from other modes to the charge instability can be accounted for by increasing the strength of electrostatic interactions. The continuous NIT of DMTTF-CA is induced by increasing the electrostatic interaction V according to an empirical $V\left(T\right)$ relationship derived by reproducing the evolution of vibrational spectra in the molecular vibrations region (see inset of Figure 7 and discussion below). Other details of the calculation can be found in the original paper [23].

Figure 7 shows the calculated IR (polarized parallel to the stack axis) and Raman (perpendicular polarization) spectra in the spectral region of molecular vibrations. The Raman spectra (black lines) are characterized by the softening of approximately 100 cm⁻¹ of the totally symmetric molecular vibration band upon approaching NIT at $T_c\approx 65$ K. In the N phase (see $V=-0.5t$, $T=250$ K), the IR spectra (red lines) are characterized by the presence of the combination bands ${\omega }_H\pm {\omega }_P$. Combination bands of appreciable intensity naturally result from the strong anharmonicity of the PES, which in turn arises from the coupling between vibrations and delocalized electron (see Section 3). Upon approaching the transition (see $V=-1.3t$, $T=110$ K) combination bands broaden and, because of the softening of the Peierls mode, their spacing from the central Raman band at ${\omega }_H\approx 1400$ cm⁻¹ is reduced. We observe that the hot band (${\omega }_H-{\omega }_P$) is generally broader and more intense than the cold band (${\omega }_H+{\omega }_P$). For $V=-1.3t$ the PES is locally flat, and while the harmonic frequency of the Peierls mode vanishes, our finite temperature simulation ($T=110$ K) yields ${\omega }_P\approx 45$ cm⁻¹, measured as half the difference of the position of the two bands. The PESs for $V<-1.5t$ present the double minima characteristic of dimerized system, causing the progressive disappearance of the side bands and the observation of the intramolecular mode at ${\omega }_H\approx 1340$ cm⁻¹ in both IR and Raman spectra. For instance, for $V=-1.9t$ and $T=60$ K the IR spectrum shows the coexistence of the peak of the molecular mode, signaling stack dimerization, with combination bands, that are instead characteristic of regular stacks. This results from our ensemble average of classical dynamics that includes both trajectories crossing the barrier separating the two minima at finite δ and trajectories that are instead trapped in one of the two potential wells.

The comparison of calculations in Figure 7 with experimental spectra of DMTTF-CA [23,66] in the region of molecular vibrations is very good. While this result holds for any reasonable parameterization of a continuous NIT, a truly quantitative agreement on the softening of the Raman active molecular vibration as well as the evolution of the related IR combination bands upon approaching NIT is obtained by fine tuning the V as shown in the inset of Figure 7). In effect, $V\left(T\right)$ is optimized for Raman and IR spectra at temperature T.

We now turn our attention to the region of lattice phonons. Panel a of Figure 8 shows the IR spectra, largely dominated by the absorption band of the Peierls mode. The Peierls mode softens upon approaching NIT from about 60 cm⁻¹ ($V=-0.5t$, $T=250$ K) to 40 cm⁻¹ ($V=-1.9t$, $T=60$ K). This softening is consistent with the half difference between combination bands in Figure 7. The intensity of the Peierls mode grows upon approaching NIT as can be better appreciated by the integrated intensity in the 0–100 cm⁻¹ interval shown in Figure 8b. The spectral weight grows by approximately a factor 4 upon approaching the transition at $T_c=65$ K (dashed line) and decreases in the I phase.

These results were recently experimentally confirmed for TTF-CA [67]. Although the model is optimized to describe the continuous transition of DMTTF-CA, it is interesting to note that the relative increase of the intensity of the Peierls mode is consistent with the measurements on TTF-CA, where the 0–100 cm⁻¹ spectral weight grows from 600 ω⁻¹cm⁻² at 300 K to 2000 ω⁻¹cm⁻² at the transition temperature of 81 K (approximate values, from Figure 7 of Ref. [67]).

The growing intensity of the Peierls mode at NIT can be understood from a joint inspection the PES and of the $P_{el}(q,\delta )$ surface (i.e., the dependence of $P_{el}$ on the vibrational coordinates q and δ) and their evolution as a function of V. In fact, while the PES governs the classical dynamics of the systems, the IR absorption intensity depends on the amplitude of fluctuations of the $P_{el}$ correlation function, according to Equation (13). Panel c of Figure 8 superimposes the $P_{el}(q,\delta )$ surface (color map) and the PES (isolines) for a N system ($V=-1.3t$). We notice that the $P_{el}(q,\delta )$ surface closely resembles $P_{el}(\mathsf{\Gamma },\delta )$ in Figure 3a, as expected considering that q enters Hamiltonian (10) as an additive term to Γ. Since DMTTF-CA does not undergo a true NIT (ρ is always smaller than 0.5 [66]) it is equivalent to consider $P_{el}(q,\delta )$ or $\mathsf{\Delta }{P}_{el}(q,\delta )$ (see Section 2.6)). For this system with a largely N ground state, the single minimum of the PES (cross) lies in the region where the $P_{el}(q,\delta )$ surface is smooth. In this case thermal vibrational motion cannot induce large $P_{el}$ fluctuations and the IR intensity remains modest. As the system approaches NIT, the PES minimum shifts towards the region of large $P_{el}$ variations along δ, and thermal lattice fluctuations induce larger and larger electronic charge flows along the stack that maximize at the dimerization transition. Large IR intensities still persist in a dimerized system such as that in Figure 8d ($V=-1.3t$). In this case the PES presents two shallow minima separated by a small energy barrier that can be overcome at $T=60$ K, causing large $P_{el}$ variations. Further decrease of V and/or T traps the system in one of the two minima, suppressing the IR intensity. The anharmonic simulation provided here confirms the results of the harmonic calculation [14,15]. Indeed the harmonic approximation catches the main physics of the phenomenon, simply relating the IR intensity to the squared δ derivative of the polarization at equilibrium ${\left(\frac{\partial P}{\partial \delta }\right)}_0^2$.

IR spectra in the lattice phonon region can be understood to a first approximation, in the harmonic limit. However, in the same spectral region Raman spectra are largely dominated by anharmonicity effects. Panel a in Figure 9 shows Raman spectra calculated for light polarized along the stack axis, hence probing the polarizability due to delocalized electrons. In the high-T N phase, the Peierls mode is silent in Raman spectra by symmetry, hence in the harmonic approximation no signal is expected in the region of lattice phonons. At $T=250$ K ($V=-0.5t$, black line) the spectrum is characterized by the presence of a single band around 120 ${\mathrm{cm}}^{-1}$. The comparison with IR spectra in the region of lattice phonons (see Figure 8a) and of molecular vibrations (see Figure 7) allows identifying the Raman band as the overtone of the Peierls mode $2{\omega }_P$. The overtone band is symmetry-allowed in the regular stack and its presence is a manifestation of the PES anharmonicity. This band softens upon approaching NIT and acquires an asymmetric shape that can be also ascribed to the anharmonicity of the potential. Most importantly, the intensity of the Raman peak greatly increases approaching the transition and maximizes at $T=55$ K (blue line), where the appearance of the higher overtone $4{\omega }_P$ is predicted by the calculations. As can be seen by the integrated intensity plot in Figure 9b, the maximum Raman signal occurs just below $T_c=65$ K. In the dimerized phase at $T=20$ K, the Peierls mode is directly observable both in Raman and IR (not shown) spectra at ${\omega }_P\approx 80$ cm⁻¹ [23].

These theoretical results are confirmed by the experimental observation of a broad and intense Raman signal appearing below 100 cm⁻¹ in proximity of the T-induced transition of DMTTF-CA [23] and of the similar compound DMTTF-QBr₂Cl₂ [68]. More interestingly, simulations provide a key to understand the anomalous Raman signal that is generated by $\rho$ vibrations and electronic polarization at finite dimerization (see Figure 4). Figure 9c,d show the ${\alpha }_{el}(q,\delta )$ surface (log-scale color map) with PES isolines superimposed. As the system evolves from a N regular phase (Figure 9c, $V=-1.3t$) towards the transition, the PES minimum approaches the singularity of ${\alpha }_{el}(q,\delta )$ and the system dynamics explores more and more regions that greatly contribute to the correlation function in Equation (14). In the just-dimerized phase at $V=-2.2t$ (Figure 9d), the PES minima are located just across the singularity that is however accessible by thermal fluctuations at $T=55$ K, conferring huge Raman signal to the overtones of the Peierls mode.

More generally, this careful theoretical analysis sheds light on the entanglement of charge and vibrational degrees of freedom that confers a strongly anharmonic character to systems close to the NIT. The effect is twofold: on one hand there is the so-called mechanical anharmonicity, i.e., the PES dependence on vibrational coordinate largely deviates from quadratic behavior, leading to the appearence of combinations and overtone bands in vibrational spectra. On the other hand, the dependence of the electronic polarization and polarizability is strikingly different from lowest-order approximations. This electrical anharmonicity is responsible for the large intensity of anharmonic features that are prominent in the vibrational spectra of ms-CT crystals [23,66,69]. Finally, the observation of an anomalous Raman scattering in the lattice phonon region at NIT [23] provided an indirect experimental evidence for the divergent electronic polarizability of the regular lattice.

5. Multistability and Photoinduced NIT

Multistability is the characteristic feature of discontinuous phase transitions, opening the way to metastable states, coexisting phases and photoinduced phase transitions. The energy of metastable domains depends on the domain length and on the boundaries, as discussed in [20]. The T-induced NIT of TTF-CA occurs at too low temperature to allow for the population of metastable domains, but the pressure-induced NIT at room temperature shows a sizable population of metastable domains, leading to coexisting N and I features in vibrational spectra as shown in Figure 20 of Masino et al. [25]. While the thermal population of metastable domains is an interesting feature, the possibility to populate metastable states via photoexcitation in the so-called photoinduced NIT is even more interesting.

Light induced transitions have been observed in TTF-CA crystals at temperature close to the discontinuous NIT. Starting from the I phase a TTF-CA crystal is driven by light absorption to a metastable state similar to the N high-temperature phase [5,6]; and vice versa for a N crystal to a I-like state [8]. Due to the long lifetime of metastable states, usually in the range of nanoseconds, this process is referred to as a photoinduced NIT (PINIT). PINIT has been demonstrated in TTF-CA through spectroscopic measurements (time resolved pump-probe experiments) [5,6,8] . Moreover, in time resolved X-ray experiments it was possible to observe the PINIT through the evolution of specific Bragg reflections after photo excitation [9], confirming the bulk nature of the observed phenomenon.

Theoretical modeling of PINIT is fairly recent and requires the addition of a pump pulse to systems of correlated electrons coupled to soft lattice and molecules. The Peierls-Holstein-Hubbard model for photoinduced transitions reviewed by Yonemitsu [70,71] is broadly similar to the full Hamiltonian in Equation (10) with linear coupling to harmonic molecular and lattice modes. The model is applied to TTF-CA as well as to charge ordering transitions and Mott insulator to metal transitions in BEDT-TTF₂ salts that correspond to 3/4 filled 2D Hubbard models of D = TTF. Yonemitsu keeps the full Hubbard basis, explicitly treats the pump pulse of finite duration using the vector potential, and finds the time evolution of lattice and molecular vibrations in finite systems during and following an excitation pulse. We start instead with an impulsive, delta-function pulse at $t=0$ and follow the time evolution of the system on anharmonic PES derived from other TTF-CA spectra. Accordingly, we explicitly and quite naturally address phenomena that results from the large anharmonicity of PES close to NIT and specifically the large coupling between the dimerization and the molecular-vibration motions, as discussed below.

In spectroscopic pump-probe experiments the crystal reflectivity in the visible region is measured after photoexcitation, to give information on the molecular properties [6,8]. This experiment will be discussed in the following. A threshold effect on the incident light intensity was observed for the PINIT along with a large non-linearity. Through a single photon absorption, multiple ionized DA pairs are created; the reported numbers vary from ten [72] to a few hundred [6] pairs. This behavior is very atypical, since, according to the standard exciton theory, transitions involving CT between more then one DA pair would be forbidden [73,74,75].

PINIT in TTF-CA occurs in two steps. In the first few picoseconds (ps) after photoexcitation the CT mechanism is triggered and the first ionized chains are created; then from ∼10 ps to the first nanosecond (ns) the ionized chains reorganize and long-range ferroelectric order is created. This comparatively slow dynamics can be conveniently followed by time-resolved X-rays experiments [9]. Optical pump-probe experiments, instead, offer a very detailed description of the early time dynamics of the system (0–10 ps) that we aim at modelling.

In particular we address very interesting experiments by Okamoto et al. [72,76], where the variation of reflectivity ($\mathsf{\Delta }R/R$) of TTF-CA is measured in the visible region following an ultrafast (15 fs) laser pump pulse. The data were collected at different temperatures, different excitation photon density and at different probe wavelength. The pump beam is polarized parallel to the stacking axis $\overrightarrow{a}$ of TTF-CA with an energy interval of $0.65÷1.2$ eV that matches the broad CT band seen in the polarized absorption spectra of TTF-CA. The probe beam is polarized perpendicularly to $\overrightarrow{a}$ with an energy of 2.2 eV in the region of localized electronic excitations that are sensitive to the ionicity of the crystal. $\mathsf{\Delta }R/R$ spectra collected at different time delays between the pump and the probe are compared with the difference between spectra collected at 4 K (I phase) and at 90 K (N phase): the strong similarity between the shape of the pump-probe and the steady state spectra suggests that an I-like state is generated by the pump signal. Accordingly, the $\mathsf{\Delta }R/R$ variations were interpreted [72,76] as ionicity variations.

Since the pump photon irradiation density per DA pair was measured during the experiment, Okamoto et al. [72,76] estimated the number of ionized pairs created by a single photon absorption as $\sim 10$ from the $\mathsf{\Delta }R/R$ jump at $t=0$. A detailed analysis of the temporal evolution of $\mathsf{\Delta }R/R$ shows an oscillatory component with frequencies related to molecular or lattice vibrations. In particular the $\sim 52$ ${\mathrm{cm}}^{-1}$ component was interpreted as related to the Peierls dimerization mode. Additionally, a wavelet analysis was performed on the oscillating signal to study the temporal evolution of the frequency components: some of them showed a significant shift in time. For example, the ∼957 cm⁻ component, which was assigned to an e-mv coupled vibration of CA ($a_g{\nu }_3$), oscillates in time at almost the same frequency as the Peierls mode. We will demonstrate that this interesting result is again a consequence of the strong anharmonicity of PES close to the phase transition.

To model PINIT we need not only the ground state PES (G-PES) that has been used to discuss vibrational spectra, but the excited state PES (E-PES) as well. We report results obtained for standard TTF-CA model parameters, as discussed in Section 3.1, but relax the mean field approximation, and explicitly account for nearest neighbor interactions. Figure 10, left panel, shows the adiabatic G-PES and E-PES calculated for $V=-1.55t$. The G-PES shows a stable neutral regular phase and two metastable dimerized minima at $\pm \delta$. The two surfaces touch at a conical intersection, the red dot in the middle panel, that corresponds to the critical point ${\mathsf{\Gamma }}_c$ of the MHM discussed in Section 2. The vertical arrow in the left panel marks the strong, optically allowed absorption polarized along the stack that initiates PINIT by generating a lattice-unrelaxed exciton of multielectron transfer character.

After a fs pump pulse, the photoexcited system evolves on the E-PES. The conical intersection makes it possible to tunnel towards one of the metastable minima. To model this complex dynamics we take advantage of the nature of the G-PES close to a discontinuous phase transition, where multistability indicates nearby surfaces. More specifically, the diabatic PES relevant to the neutral regular state and to the ionic dimerized state cross and generate the G-PES and E-PES when they are adiabatically mixed. Therefore, starting from the adiabatic PES for the ground and first excited state in Figure 10 we construct the diabatic PES relevant to the ionic dimerized (metastable) state, as shown in Figure 10, middle and left panels. This PES practically coincides with the E-PES in the region of the vertical excitation and with the G-PES in the dimerized metastable region. The middle panel shows the PES with minimum gap between E-PES and G-PES. The right panel is the diabatic PES obtained by smoothing out the gap. While the smoothing procedure is somewhat arbitrary, we emphasize that the dynamics is marginally affected by the details of the smoothing. In fact this critical region is only crossed once at the beginning of the photoexcited dynamics.

A detailed analysis of the early-time dynamics of the photoexcited system is reported in Ref. [24]. Here in Figure 11 we show the trajectory of the system after instantaneous excitation and the Fourier transforms of the $\delta \left(t\right)$ and $q\left(t\right)$ signals that show oscillating components at $\nu \sim$ 55 and 945 cm⁻¹, respectively. The wavelet analysis of the same signal, shown as a color map on panel b in Figure 11, clearly shows that the molecular vibration oscillates at a frequency that is not constant in time, but that in turn oscillates about a mean values (∼945 cm⁻¹) at a frequency that corresponds to that of the dimerization mode (∼55 cm⁻¹). In the equilibrium state, before photoexcitation, the system is in a neutral regular state. After vertical photoexcitation q is wildly out of equilibrium and a very fast initial dynamics brings it towards the equilibrium position relevant to the photoexcited ionic states. This fast initial dynamics accounts for the observed initial rise of the reflectance variation. After the initial increase, $q\left(t\right)$ oscillates about the new equilibrium position with an amplitude that progressively decreases due to friction. Indeed, due to anharmonicity, the motion along the two coordinates is mixed and while δ oscillates, different regions of the PES are explored, characterized by a different curvature along the q coordinate, then explaining the modulation of the observed frequency for the q mode by the δ frequency, as shown in the inset of Figure 11c, and in line with experimental observation.

Figure 11d shows the temporal evolution of the system ionicity, that closely resembles $q\left(t\right)$. The variation of ionicity of the system is proportional to the number of D⁺A⁻ pairs created upon photoexcitation, as shown on the right axis of Figure 11d. Results in the figure tell us that upon vertical excitation (the dot located at $t=0$) almost 1.5 electrons are transferred from D to A sites, a result that calls for sizable correlation of the motion along the chain. However, ultrafast lattice relaxation, and more specifically the relaxation along the q coordinate, further increases the number of transferred electrons, that oscillates with q around an average value of about 4. This intriguing result is in line with the experimental observation that photoexcitation generates a domain (D⁺A⁻)_n with $n\sim 10$ within 20 fs.

Vertical excitation at 260 K, far from NIT for TTF-CA parameters, also generates ionic domains [76]. In the simplest approximation, coherent oscillations of coupled molecular vibrations are again modeled by q, but there are no δ-oscillations since δ = 0 is the minimum on both surfaces. Now the thermal population of the Peierls mode with ${\omega }_P\sim$ 50 cm⁻¹ must be taken into account. We do so along the same lines discussed in the previous Section and in Ref. [23]. We choose a symmetric δ-distribution in $E_G(q,\delta )$, the G-PES, follow the coherent dynamics after a vertical excitation, and weight the results of many runs by a Boltzmann factor in $(E_G(q,\delta )-E_G(q,0))$. Coherent δ-oscillations are generated with greatly reduced amplitude, as found experimentally [76]. The most interesting feature is the doubled frequency $2{\omega }_P\sim$ 100 cm⁻¹ in the Fourier transform of ρ that corresponds to Figure 11e at 90 K. The frequency of the oscillations around $\delta =0$ is doubled since ρ is symmetric on δ. For completeness, we note that the 90 K results near the NIT in Figure 11 hardly change upon accounting for thermal population of the PES along δ.

6. First Principles Model Parameters

The MHM provides a robust framework to understand the valence and structural instabilities in ms-CT crystals. Input parameters are required for MHM, as is typical in a semiempirical model. Parameters have traditionally been obtained and progressively refined for prototypical systems such as TTF-CA or DMTTF-CA on the basis of optical, vibrational and other experimental data. In this section we cover a recent alternative approach that instead relies on DFT electronic structure calculations [61,77]. This has the advantage that the crystal structure is the only experimental input needed, possibly at different temperature or pressure.

The most refined model for ms-CT crystals accounts for several interactions and hence several parameters enter $H_{msCT}$, Equation (10). The parameter set includes those describing the electronic structure (t, z), the coupling to molecular and lattice vibrations (${ϵ}_v$ and ${ϵ}_d$), as well as intermolecular electrostatic interactions (V, $\mathcal{M}$ or the derived ${ϵ}_c$). Reliable estimates of these quantities can be obtained from atomistic calculations, with the notable exception of ${\epsilon }_d$.

As discussed in Section 2, the MHM can be considered a 1D extension of the Mulliken dimer model, whose mean-field electronic properties are fully determined by t and $z_{eff}$. To fix these model parameters, we extract DA dimers from the experimental crystal structure of specific ms-CT crystals and run DFT calculations to map relevant states into the Mulliken two-state model. The mapping proceeds through the identification of two properties of choice that can be accessed by both approaches. In our case, we opt for quantities obtained from ground-state calculations in the singlet and in the triplet subspace, namely the ionicity ${\rho }_{dim}$ and the singlet-triplet gap ${\mathsf{\Delta }}_{ST}$, leading to the following closed expressions:

$$\begin{array}{ccc}\hfill t& =& {\mathsf{\Delta }}_{ST}\sqrt{\frac{{\rho }_{dim}}{2(1-{\rho }_{dim})}}\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill z_{eff}& =& t\frac{1-2{\rho }_{dim}}{\sqrt{2{\rho }_{dim}(1-{\rho }_{dim})}}.\hfill \end{array}$$

(16)

The choice of ground-state properties is considered to be safer with respect to excited state calculations, which are known to be troublesome in CT systems within the framework of time-dependent DFT. Range-separated hybrid functionals, such as CAM-B3LYP or ωB97XD provide reliable estimates for the sought quantities, in line with earlier empirical estimates [61,77]. The favorable comparison with experimental data for the ionicity ρ in the crystal validates this choice a posteriori.

In the adiabatic approximations underlying Hamiltonian (10), the strength of the coupling to effective vibrational modes is measured by the relaxation energies ${ϵ}_v$ and ${ϵ}_d$. ${ϵ}_v$ is defined as the energy gained upon intramolecular structural relaxation associated with the electron transfer from D to A and includes the contribution from several vibrational modes of the molecule, yet q in Equation (10) is an effective coordinate that accounts for their combined effect. It is, therefore, preferable to compute intramolecular reorganization through DFT calculations performed at the geometries optimized for isolated neutral and charged molecules as ${ϵ}_v={\lambda }_D^{+}+{\lambda }_A^{-}$, where ${\lambda }_D^{+}$ (${\lambda }_A^{-}$) is the hole (electron) relaxation energies for the D (A) defined in standard electron transfer theories.

On the other hand, a first principles estimate of ${ϵ}_d$ is computationally challenging as it would require the calculation of the energy for the formation of a bond-order charge density wave in the crystal (excluding dispersion and steric interactions that are modeled by an elastic term in the harmonic approximation, see Equation (10)) along an effective dimerization coordinate that enforces rigid molecules. Given these practical difficulties, in this case a different approach with respect to that of Section 3 has been proposed, and calculations are performed at the dimerization obtained from experimental crystal structure (${\delta }_{exp}$) and not at the optimal lattice geometry for a given ${ϵ}_d$, or scanning the PES near equilibrium. This is straightforward in the case of regular stacks where ${\delta }_{exp}=0$ by definition. In dimerized systems, two independent DFT calculations are required for nonequivalent dimers in symmetry-broken stack, obtaining $t^{\pm }=t(1\pm {\delta }_{exp})$.

The evaluation of long range intermolecular electrostatic interactions goes beyond the reach of first principles calculations, calling for multiscale approaches that can be parametrized ab initio. Electrostatic models of atomistic resolution represent here an ideal compromise between accuracy and computational cost [78]. In these approaches point atomic charges approximating the charge density of neutral or charged molecules and experimental atomic positions are employed to compute intermolecular interactions via direct electrostatic summations, possibly accounting also for the screening provided by molecular polarizabilities [61]. Since V is readily computed for a dimer and $\mathcal{M}$ for a crystal, we obtain the quantity ${ϵ}_c=V-2\mathcal{M}$ that enters the mean field approximation to (10).

TTF-CA is by far the most studied ms-CT crystal and its crystal structure has been experimentally determined at temperatures ranging from 15 to 300 K, i.e., in a large interval around $T_c=81$ K. In Ref. [61] the parametrization described above has been performed for the different crystal structures, allowing the theoretical description of the NIT as a quantum phase transition where temperature only enters through the structural change. Apart from the dimerization occurring below T_c, the most remarkable variation with T of the parameters is the 45 meV increase of the Madelung energy shown in Figure 12a. The effect of the increasing strength of electrostatic interactions on the MHM ground state energy and ionicity are shown in Figure 12b,c, respectively. At room T TTF-CA is largely neutral and its ground state is progressively destabilized by the increasing $\mathcal{M}$ due to the thermal contraction of the lattice. Upon decreasing temperature, the system enters a bistable region with coexisting N and I phases, the latter stabilized by the $\mathcal{M}$ increase. This bistability window is a characteristic of discontinuous transitions and relevant energies (free energies at finite T) discriminate between the thermodynamically stable and the metastable phase (see also Section 2.3 and Section 3). The calculated evolution with temperature of the ionicity $\rho$, shown in Figure 12c, favorably compares to experiment, although quantitative deviations are evident. In fact, TTF-CA turned out to be always too neutral or too ionic in the calculation, a flaw that can be attributed to the mean field treatment of electrostatics, which may suppress intermediate $\rho$ states, or by finite temperature effects such as thermal population or disorder. Nevertheless, capturing the discontinuous NIT of TTF-CA in the correct temperature range and without introducing adjustable parameters is a remarkable result that validates the MHM parametrization through comparison with experiment. This result further confirms the nature of the NIT as a quantum phase transition driven by the increase of the Madelung energy with the lattice contraction.

In a more recent study [77], some of us extended this approach to eleven ms-CT crystals, demonstrating the wide applicability of the first principles parametrization. With the only exception of just one among the 11 investigated systems, TMPD-TCNQ [79] for which the method was inapplicable, in all other cases it has been possible to find good correlations between the experimental and the calculated ionicity $\rho$ [77]. Moreover, the different ms-CT crystals have been classified in terms of their bistable character and their propensity towards a discontinuous transition. The key parameter for the classification was identified as the ratio $({ϵ}_v+{ϵ}_c)/t$ (see Section 2.3 and Section 2.5), which returns a good correlation between parameters obtained from atomistic calculations for each system and their experimental behavior.

Theoretical calculations on specific experimental systems, explicitly accounting for the strong correlations characterizing ms-CT crystals, have been recently proven to be possible and reliable. The approach discussed here collapses the complexity of several systems into a few physically meaningful parameters, allowing to draw insightful connections between molecular chemistry, supramolecular organization in the solid state and the electronic properties of materials. All this can be then rationalized within the framework of the universal picture of NIT in ms-CT crystals provided by the MHM accounting for electron-phonon coupling.

7. Summary

The general situation discussed here with respect to modeling NIT is typical of the evolution of understanding of condensed phase systems. The initial focus on ionicity $\rho$ was soon expanded to the structural instability to dimerization and the $\rho$-dependence of molecular vibration frequencies. Mixed stack CT crystals are a subgroup of the larger family of quasi-1D CT crystals that also include segregated stacks. The optical, magnetic, electrical, vibrational and structural properties of crystals with quasi-1D stacks have posed interesting experimental and theoretical challenges for decades. The present consensus is in terms of correlated electrons in extended Hubbard-like models, coupled to the Peierls phonon and molecular vibrations, and 2 or 3D considerations as suggested by experiment. We have not gone beyond 1D modeling of a single stack, but ferroelectricity in ms-CT crystals clearly requires that all stacks dimerize in the same sense.

The microscopic NIT model presented here is minimally required to account for the electronic properties of ms-CT stacks. The model has one lattice mode ($\delta$), one molecular mode (q) and usually a mean field approximation for Coulomb interactions. The experimental situation is far more nuanced and complex. The IR and Raman spectra of TTF-CA or of DMTTF-CA presented in the companion work by Masino et al. [25] require additional lattice phonons and more sophisticated vibrational analysis. Depending on one’s perspective, more quantitative treatment either underlines the approximate nature of the microscopic model or merely the need to introduce additional parameters. Although we have taken as few parameters as possible, there is a handful: t and $\mathsf{\Gamma }$, the basic parameters of the modified Hubbard model in the reduced basis, ${ϵ}_c$ measuring the strength of electrostatic intermolecular interactions in the ionic lattice, V the electrostatic interaction between adjacent D⁺ and A⁻ sites and ${ϵ}_v$ and ${ϵ}_d$, the strength of the coupling of electrons to molecular and lattice modes, respectively. The model may be minimal but it is not simple and a consistent set of parameters is based on many experiments. While ${ϵ}_c$, V, and perhaps $\mathsf{\Gamma }$ are in the eV scale and can be evaluated with a reasonable accuracy, NIT inherently features competition (differences) between large energies. The state of the system depends sensitively on energies (t, ${ϵ}_v$ and ${ϵ}_d$) that are an order of magnitude smaller, and still quite challenging to evaluate reliably to 10% accuracy. Nevertheless, progress with parameter evaluation is an important aspect of modeling organic solids in general that makes credible small additional adjustments for comparison with experiment.

The MHM focuses on the HOMO of D and LUMO of A in extended correlated-electron stacks. The optical CT excitation around 1 eV is polarized along the stack and is the lowest electronic excitation, but not by much; $\pi \to {\pi }^{*}$ excitations polarized in the molecular plane may start around 1.5 eV. Modern treatments of systems with a few D or A moieties sometimes distinguish among electronic excitations but often do not. A more general approach to the electronic structure of an extended stack and eventually to many stacks may well be possible and even required in the future.

The NIT model that we have discussed is spectroscopic in nature. The ground state energy $E_0(\delta ,q)$ is a PES with a lattice and molecular coordinate, both strongly coupled to correlated electrons. The PES has a single minimum on the neutral side, two minima on the ionic side, and metastable states for very special parameter values in systems with discontinuous NIT. The PES accounts semi-quantitatively for vibrational spectra discussed in Section 4 and by Masino et al. [25]. The emergence of highly anharmonic surfaces from a purely harmonic lattice and molecular vibration is an interesting consequence of coupling to correlated electrons that are delocalized along the stack. As emphasized throughout, the model provides a unified treatment of continuous and discontinuous NIT with small changes of parameters. That is our principal result. We have invoked the excited-state PES to model the early dynamics of photoinduced NIT in TTF-CA, which speaks to the model’s robustness.