Incommensurate quantum-size oscillations in acene-based molecular wires - effects of quantum fluctuations

TL;DR

This paper investigates the quantum-size oscillations in acene-based molecular wires, demonstrating that these oscillations persist despite electron interactions, using a combination of theoretical models and numerical methods.

Contribution

It provides a detailed analysis showing that incommensurate gap oscillations remain robust in acene ladders even when electron-electron interactions are considered.

Findings

Incommensurate gap oscillations persist for a range of interaction strengths.

Interactions can potentially open a large gap, but do not eliminate oscillations.

The study combines Hartree-Fock and DMRG methods for analysis.

Abstract

Molecular wires of the acene-family can be viewed as a physical realization of a two-rung ladder Hamiltonian. For acene-ladders, closed-shell ab-initio calculations and elementary zone-folding arguments predict incommensurate gap oscillations as a function of the number of repetitive ring units, $N_{\text{R}}$, exhibiting a period of about ten rings. %% Results employing open-shell calculations and a mean-field treatment of interactions suggest anti-ferromagnetic correlations that could potentially open a large gap and wash out the gap oscillations. % Within the framework of a Hubbard model with repulsive on-site interaction, $U$, we employ a Hartree-Fock analysis and the density matrix renormalization group to investigate the interplay of gap oscillations and interactions. % We confirm the persistence of incommensurate oscillations in acene-type ladder systems for a significant fraction of parameter space spanned by $U$ and $N_{\text{R}}$.

Figures

Equations

$$\hat{H}=\sum _{i\ne j}{t}_{ij}{\hat{c}}_{i\sigma }^{†}{\hat{c}}_{j\sigma }^{}+U\sum _i\left({\hat{n}}_{i\uparrow }-\frac{1}{2}\right)\left({\hat{n}}_{i\downarrow }-\frac{1}{2}\right)$$

$${ϵ}_{sb}(k)=2t^{\prime }\cos (k)+\frac{1}{2}s\left(t+t_{\perp }^{\prime \prime }+2t_{\times }^{\prime \prime }\cos (k)\right)+\frac{sb}{2}\sqrt{{\left(t-t_{\perp }^{\prime \prime }+2t_{\times }^{\prime \prime }\cos (k)\right)}^2+8(1+\cos (k))(t+s{t}^{\prime }{)}^2},$$

$${\Delta }_{1,N_{\mathrm{R}}}^{\mathrm{HF}}=\sqrt{{\Delta }_0^2+(M^{\mathrm{HF}}U{)}^2},$$

$${\Delta }_{1,\infty }^{\text{HF}}(U)\approx 0.5U{e}^{\frac{-D}{U^2+D_{1}U}},$$

$${\frac{{\Delta }_{\infty }^{\text{HF}}}{{\Delta }_{\infty }^{\text{DMRG}}}|}_{U\approx D}\sim \frac{D^2}{t^2};$$

$$\cos (\pi -k_{\mathrm{D}})=\frac{t^2-\frac{1}{2}t{t}_{\perp }^{\prime \prime }}{t^2-t_{\perp }^{\prime \prime }{t}_{\times }^{\prime \prime }}.$$

$$\pi -k_{\mathrm{D}}$$

$$v_{\mathrm{F}}$$

$$(1-k_{\mathrm{D}}/\pi {)}^{-1}=\pi \sqrt{t{t}_{\perp }^{\prime \prime }}\left(1+\frac{t_{\times }^{\prime \prime }}{t}+\mathcal{O}(t^{\prime \prime }/t{)}^2\right).$$

$${\Delta }^{\text{DMRG}}(N_{\mathrm{R}})$$

$$\Delta k(N_{\mathrm{R}})$$

$${\Delta }_1(L)$$

$$g(x)=\frac{{\sum }_{\ell =1,3,5,\cdots ,}^n\cos (\ell x)/{\ell }^2}{{\sum }_{\ell =1,3,5,\cdots }^n{\ell }^{-2}}$$

$${\mathcal{H}}^{\mathrm{HF}}={\hat{H}}^{U=0}+U\sum _{x,\sigma }\left({\bar{n}}_{\sigma ,x}-0.5\right){\hat{n}}_{-\sigma ,x},$$

$$E^{\mathrm{HF}}=E^{\mathrm{SC}}-\frac{UM}{4}+U\sum _j{\bar{n}}_{\sigma ,x}{\bar{n}}_{-\sigma ,x}$$

$$m_j$$

$$M^{\mathrm{HF}}$$

Full text

I Introduction

The question how properties of a macroscopic system emerge when more and more atoms or molecules accumulate exhibits many different facets and, for that reason, reappears every once in a while in different contexts. On the simplest level of “emergence”, one could consider the appearance of band-structures in crystal growth. A first qualitative description of this phenomenon can already be given within a picture of non-interacting quasi-particles. A conceptually analogous situation arises in mesoscopic physics: when a quantum dot couples to an electrode, the dot states hybridize with contact states, thereby acquiring the level broadenings, $\Gamma$, that signal finite lifetime effects. Qualitatively new physics emerges upon accounting for interaction effects. For instance, in quantum-dots, the Kondo phenomenon results from the interplay of hybridization and the Coulomb-energy. Another example displaying the appearance of strong correlation effects with growing system size is superconductivity. vonDelft ; EfetovReview There, one is facing a crossover scenario, where the super-conducting gap, $\Delta_{\text{sc}}$, competes with the single particle level spacing, $\Delta_{0}\gg\Gamma$.

In this paper, we investigate an incarnation of a similar motif as it appears in molecular electronics: we ask how the electronic spectral properties of a molecular wire built out of $N_{\text{R}}$ repetitive units (“rings”) evolve with increasing $N_{\text{R}}$. Motivated by earlier experiments kiguchi2008highly ; yelin2013 ; yelin2015 , our specific example will be the acene-family such as benzene, naphthalene, anthracene, and others. Due to their peculiar electronic structure, oligo-acenes have received a considerable amount of attention in organic chemistry over the past 50 years.bettinger10 ; bettinger16 ; yang16 A recent technological interest is based on suggestions to use them as light harvesters in organic photo-voltaic devices. dabo13 Besides being a simple realization of nano-wires, acenes can be viewed as the narrowest graphene nano-ribbons with zig-zag terminated edges. Due to recent progress in synthesis methodology, oligoacenes up to length $N_{\text{R}}=9$ (nonacene) have already been produced.toenshoff10 ; bettinger15

Previously, we have found that linear (oligo)acenes exhibit oscillations in the $N_{\text{R}}$ (or length) dependence of optical excitation gaps. korytar2014 Oscillations of electronic properties with the geometric size of a nano-system are per se frequently encountered. An established example is given by the threefold gap-oscillation of zigzag carbon nanotubes (armchair-ribbons) as a function of tube diameter (width). reichBook ; brey06 ; dasgupta12 In this context, the three-fold periodicity reflects the symmetry of the hexagonal Brillouin zone, with the Dirac points located at the corners. A remarkable aspect of the acene oscillations is that they are incommensurate in the sense that the oscillation period is not dictated by lattice symmetry; it can reach periods of ten times the length of the unit cell. The oscillations come with an important consequence, namely, a close-to closing gap for molecules of certain (periodically repeating) finite length.

We have demonstrated that similar to the case of the nanotubes, the acene-oscillations can likewise be described (nearly) quantitatively by a selection rule (zone-folding) for the wave-vectors applied to the band structure of the underlying extended tight-binding model (Fig. 1). The procedure is analogous to the threefold periodicity in band-gaps of nanotubesreichBook and nanoribbonsbrey06 ; dasgupta12 ; in the acene case, however, incommensurate periodicity is reached because of the non-universal position of $k_{\text{D}}$ labelling the longitudinal momentum at the acene-analogue of the Dirac point. Based on the computational analysis of the acene-like Hubbard model, we find that the zone-folding argument is not invalidated by interactions beyond the independent quasi-particle picture; the oscillations persist given the interaction strength $U$ does not exceed a threshold of the order of the bandwidth - at least as long as interactions are screened, i.e., short-ranged.

Building upon our previous work, we intend to give a comprehensive study of interactions and incommensurate oscillations (IO) of acenes by contrasting Hartree-Fock (HF) mean field treatments against microscopic numerical studies. We model acene-type molecular wires by finite length Hubbard ladders. A principal physical incentive is given by the potentially rich phase diagram of interacting acene wires karakonstantakis13 . In this context, one could hope to manipulate acenes via chemical synthesis, or engineering of the environment so as to explore some part of this phase diagram. (Admittedly before this is experimentally feasible, major technological difficulties have to be overcome that are realated to chemical stability. Momentarily, long acene wires need to be kept within a noble-gas matrix to prevent further reactionsbettinger15 .)

Further methodological incentive is provided by recent open-shell $\text{sc}GW$-calculations for the acene-serieswilhelm16 ; they indicated that in short wires the IO might disappear in the presence of anti-ferromagnetic order in the ground state. At first sight, this result may look very plausible. However, effects of quantum-fluctuations are neglected in its derivation and in quasi-onedimensional systems their effect typically is very strong. Specifically, in the long-wire limit, quantum fluctuations tend to destabilize mean-field ordered phases that break a continuous symmetry PhysRevLett.17.1133 . In this (quantum-disordered) phase, the order-parameter exhibits significant short-range correlations that leave a characteristic signature in local probe measurements.

In this work, we clarify the fate of the anti-ferromagnetic correlations in acenes-type ladders in the presence of quantum fluctuations, i.e., by providing a joint analysis from Hartree-Fock (HF) mean field and density matrix renormalization group (DMRG). Within HF mean field analysis, we reproduce the observation made in the ab-initio study Ref. wilhelm16, for acene-molecules with Coulomb interactions: acene-ladders with onsite Hubbard interactions exhibit anti-ferromagnetic order on the mean field level along with a large charge gap, thereby diminishing IO. By microscopic DMRG analysis, however, we show that also in short ladders quantum fluctuations destroy the long-range order. In particular, they significantly reduce the HF excitation gap, restoring the IO signature and up to intermediate interaction strengths of the order of the bandwidth.

II Model and Methods

II.1 Model definition

We study the acene Hubbard-ladder

[TABLE]

at half filling, where the first term describes the tight-binding dynamics, and the second term in Eq. (1) the onsite Hubbard repulsion with strength $U$. In the case of periodic boundary conditions (PBC), we have a linear length of $L=2N_{\mathrm{R}}$ sites, where $N_{\mathrm{R}}$ gives the number of hexagon (benzene) rings. For hard wall boundary conditions (HWBC), we complete the last ring leading to $L=N_{\mathrm{R}}+1$. The total number of sites is $2L$. The nearest neighbor hopping is characterized by an amplitude $t_{\langle ij\rangle}=-t$, and the longer-range hoppings accordingly by $t^{\prime}$, $t^{\prime\prime}_{\perp}$, and $t^{\prime\prime}_{\times}$ (Fig. 1). If not stated otherwise, $t$ will be taken as the unit of energy and $t^{\prime\prime}=t^{\prime\prime}_{\perp}$ = $t^{\prime\prime}_{\times}$. Note that $t^{\prime\prime}_{\perp}$ and $t^{\prime\prime}_{\times}$ give rise to a non-avoided crossing of conduction- and valence band at non-commensurate wave-vectors $k_{\text{D}}$ predicted for the case $U{=}0$.Kivelson For $t^{\prime\prime}_{\perp},t^{\prime\prime}_{\times}=0$, the system features a quadratic band touching karakonstantakis13 . For $U=0$, the energy spectrum of (1) reduces to

[TABLE]

where $b,s=\pm 1$ denote a band and parity index, respectively, yielding four bands resulting from the four-site unit cell.

In the following, we exclude next-nearest neighbor hopping, $t^{\prime}=0$, so there is no hopping between sites sharing the same sublattice; the model becomes bipartite and a chiral symmetry is imposed connected to the particle-hole symmetry of the spectrum (Fig. 1). This choice can be made without loss of generality regarding IO, as the main effect of breaking the chiral symmetry in more general models with non-vanishing $t^{\prime}$ would be to shift the position of the Dirac point $k_{\text{D}}$. The corresponding analytical calculations are given in the Appendix A.

General comments on Hubbard-models. We briefly comment about the applicability of Hubbard models to real systems. Indeed, models along the lines of (1) have been investigated intensively over the last fifty years in physics and in chemistry. They provide substantial simplifications by (a) ignoring the coupling to molecular vibrations; (b) including only one orbital per atom and (c) by including all other degrees of freedom in terms of effective parameters, only, such as a short range interaction that mimics screening. In fact, information about the molecular geometry (e.g. the relative postion of the atoms) is absorbed in hopping parameters and interaction matrix elements. These parameters span a multi-dimensional space which has a hyperplane that features those particular combination of parameters that can be reached in realistic situations.

In return, because Hubbard models are in a certain sense minimal, they allow for a numerically exact treatment of the many-body problem, with a size of the computational Hilbert space that is much larger as compared to most other many-body treatments that attempt to be more “first principles” on the Hamiltonian level.

Keeping this in mind, Hubbard models like (1) are by no means as oversimplistic as they may appear at first sight. First, if properly parametrized such models have a good chance to give a realistic, sometimes even quantative description of low-energy excitations also in those situations where quantum-fluctuations are dominant, so that most competing methods fail. Second, one often is interested in the system behavior under varying external conditions, i.e. in the phase diagram of the model; one would like to figure out what phases exist, what are their properties and, ideally, also where are the phase-boundaries located at. Since the qualitative aspects of phase diagrams tend to be very robust under a change of the Hamiltonian, one often adopts convenient choices for the model parametrization, even though they are sometimes far from realistic.

Once the phase-diagram has been mapped out, it may be possible to predict the qualitative behavior for a given experiment, even though the actual parameters representative for this situation may not be known with high accuracy, at least in those situations where very stable phases have been identified.

II.2 Computational details

In this spirit, we continue our analysis and treat the model Hamiltonian (1) numerically within DMRG White:PRL92 ; White:PRB93 . Similar to previous authorsJiang ; Chan , we also observe that the ground state is a spin singlet independent of the wire length. Thereby, we track the many-body Hamiltonian at fixed particle number (equal to the number of lattice sites) and in the $S^{z}_{\text{tot}}{=}0$ spin sector. By targeting not only the ground state but also excited states, we obtain the optical gap by taking the energy difference between the first excited and the ground state. For non-interacting systems, the optical gap corresponds to the HoMO/LuMO gap, as it is given by the lowest particle-hole excitation. Due to spin-rotational invariance, the excited states of the model in (1) come in spin-multiplets according to irreducible representations of SU(2). By performing reference calculations in the spin-sectors $S^{z}_{\text{tot}}{=}1,2$, we have ascertained that the lowest excitation is a triplet state, and that the ground state is indeed a singlet state. Details of the DMRG implementation are given in Appendix C. To investigate the effect of quantum fluctuations, and to link with previous mean field treatments, it is useful to compare the DMRG results to a HF description. We have also implemented a self-consistent HF cycle to obtain mean field estimates for the energy gaps, ground state, and magnetization (Appendix C). In HF theory, we look for the lowest lying (charge-neutral) excitation energy to determine the energy gap, for which we compute the difference in energy between the lowest unoccupied and the highest occupied eigenvalue of the Fock operator.

III Results and Discussion

III.1 Hartree-Fock analysis

Before investigating finite size features of the spectrum, we analyze the infinite wire limit. In Fig. 2 (top), we display the evolution of the HF gap for repulsive interaction strength $U$. We extract the bulk gap from finite-size calculations by extrapolating the finite size gap vs. the inverse linear length of the system (Appendix B). While the non-interacting wire exhibits a band crossing at the Fermi-energy and therefore is metallic, the repulsive interaction $U$ opens up a correlation gap.

Within HF theory, the ground state is found to be unstable against antiferromagnetic order; the corresponding staggered magnetization is shown in Fig. 2 (bottom). As we will see in our DMRG-calculations below, this order will be destroyed by quantum fluctuations. We here witness the established fact that mean-field theories overestimate the tendency for symmetry breaking, at least if the interactions are short ranged. In passing, we mention that broken symmetries refer to the long-distance properties of a quantum state, which are not well described on the HF-level. In contrast, the total energy is predominantly dictated by the pair-correlation function at short-distances. The latter is largely insensitive to the global symmetries. Therefore, with respect to total energy calculations, HF can often give useful answers even though the long-range physics is described qualitatively incorrect. (Löwdin’s “symmetry breaking dilemma” of HFlykos63 )

In the perturbative limit of weak repulsion, $U{\ll}1$, we can restrict ourselves to states with AF ordering on top of the non-interacting ground state. In this situation the gap opens in long (but finite-length) wires as

[TABLE]

where $\Delta_{0}$ is the finite size gap of the non-interacting wire, $\Delta_{0}\equiv\Delta^{U=0}_{1,N_{\mathrm{R}}}$; $M_{\mathrm{HF}}$ denotes the staggered magnetization amplitude. We recall that at $U{\ll}1$ the response of energies to the staggered field is linear in $U$, implying a quadratic gap opening $\Delta^{\text{HF}}_{1,\infty}(U)\sim U^{2}$ in the thermodynamic limit (where $\Delta_{0}\rightarrow 0$).

At intermediate interactions, the HF values for the gap opening allow for a phenomenological two-parameter fit,

[TABLE]

where $D$ is of the order of the non-interacting bandwidth; we obtain $D\approx 5.5\pm 0.4$ and $D_{1}\approx 1.0\pm 0.1$, provided $U$ is not too small. As expected, in the large $U$ limit, the HF particle-hole gap is set by $U$.

Its evolution follows closely the one of the anti-ferromagnetic order parameter $M^{\text{HF}}(U)$, Fig. 2 (bottom).

We now turn to the HF analysis of acene wires with finite length. In this case, one expects that the opening of the HF gap $\Delta^{\text{HF}}_{1,N_{\text{R}}}$ decreases the ground-state energy only if $\Delta^{\text{HF}}_{1,N_{\text{R}}}$ exceeds the non-interacting (“native”) level-spacing $\Delta_{0}$ at the Fermi energy significantly. Hence, for shorter wires, the gap opening is most effective at particular values of $N_{\mathrm{R}}$ that have an anomalously small single-particle gap – and ineffective at the others. In other words, only at $N_{\mathrm{R}}$-values situated close to a minimum of the IO an AF-order develops and at the others the non-interacting ground-state prevails. In this situation IO remain visible even in $\Delta^{\text{HF}}_{1,N_{\mathrm{R}}}$.

For longer wires, however, the IO fade away, because the native band-gap is suppressed as $\sim 1/N_{\mathrm{R}}$ and eventually no longer overcomes $\Delta^{\text{HF}}_{1,N_{\text{R}}}$. As soon as $\Delta^{\text{HF}}_{1,N_{\text{R}}}$ exceeds the maxima of the native band gaps, the IO disappear from the HF spectra. Indeed, the IO in $\Delta^{\text{HF}}_{1,N_{\text{R}}}$ are seen to become weaker with increasing $N_{\mathrm{R}}$ (and $U$) (Fig. 3). This explains why large systems sizes $N_{\mathrm{R}}{\gg}1$ are needed in order to obtain the HF gap in the small $U$ regime: the finite size induced gap has to fall below the HF gap.

As the magnetization closely follows $\Delta^{\text{HF}}_{1,N_{\mathrm{R}}}$, it emerges preferably in acene wires that exhibit a relatively small (native) band-gap (Fig. 3).

Vice versa, the HF-ground state remains non-magnetized at wire-lengths, such as $N_{\mathrm{R}}=11,19,28,\ldots$, where the native gap is particularly large; there mean-field magnetism only appears at larger interaction values.

In contrast, for large interaction values the HF gap surpasses the non-interacting gap already for small system size (Fig. 4), and the IO disappear in that regime.

III.2 DMRG calculations

We reignite the discussion from the viewpoint of DMRG by the infinite wire limit. Within our DMRG-implementation, this limit can only be achieved through finite size extrapolation: for every value of $U$, the extrapolation is performed to yield one data point of the blue curve for the excitation gap in Fig. 2 (top), where it is compared against the HF result. We mention that at small values of $U$ numerical calculations are more challenging to converge in DMRG, so that the low $U$ asymptotics is not fully resolved in Fig. 2(top). The difficulty is that due to the metallic character of the wire, the entanglement content is spread over a large manifold of modes that have non-negligible weight in the reduced density matrix.

By comparing to the HF result, one directly observes the effect of strong quantum fluctuations in the DMRG-trace shown in Fig. 2 (top). Fluctuation effects are visible by a reduction of the excitation gap, but also by a modified gap evolution with $U$ as compared to the HF trajectory. The asymptotic power-law at large $U$ seen in Fig. 2 is readily understood in terms of the mapping to the Heisenberg spin-chain. It proceeds via the effective coupling constant $J(U)=4t^{2}/U$ that sets the scale for $\Delta^{\text{DMRG}}_{1,\infty}$ in this limit. Anderson:PR50 ; Anderson:PR59 ; Bulaevskii:ZETF66 ; Takahashi:JPC77 As is seen in Fig. 2 (top), the Heisenberg limit is reached at values of $U$ exceeding the bandwidth, $D$, by about a factor of two. 111We adopt the ratio $\Delta^{\text{HF}}_{\infty}/\Delta^{\text{DMRG}}_{\infty}$ taken at a cross-over scale, $D$, to quantify the strength of quantum-fluctuations in a parametric fashion. A natural choice for $D$ is the (non-interacting) bandwidth that according to Fig. 1 (right) is given by $D\approx 6t$. Furthermore, $\Delta^{\text{DMRG}}_{\infty}(U{\approx}D)\sim t^{2}/D$ while $\Delta^{\text{HF}}_{\infty}(U{\approx}D)\sim D$ so that we expect a large ratio

$\left.\frac{\Delta^{\text{HF}}_{\infty}}{\Delta^{\text{DMRG}}_{\infty}}\right|_{U\approx D}\sim\frac{D^{2}}{t^{2}};$

it takes values about $36$ broadly consistent with the factor $35\pm 10$ seen in Fig. 2(top).

Our final Fig. 5 displays the most important result of our work. It confronts the mean-field data for the excitation gaps with the results obtained from the DMRG-calculations for finite-length wires. Specifically, the evolution of their DMRG-excitation gaps, $\Delta^{\text{DMRG}}_{1,N_{\mathrm{R}}}$, is shown for increasing length of the wire $N_{\mathrm{R}}$ at different $U$. It is clearly seen there, that the DMRG-traces exhibit pronounced IO even for those values of $U$ at which the residual oscillations in $\Delta^{\text{HF}}_{1,N_{\mathrm{R}}}$ have almost faded away, already. Very much in the spirit of the long-wire limit, we understand this result as a manifestation of quantum-fluctuations in finite-size wires. They destroy the tendency of mean-field theories to open up excitation gaps, thereby restoring the IO as a hallmark of the symmetry-unbroken phase.

IV Conclusions

We have studied the acene Hubbard model for a large range of on-site repulsion $U$. The gap oscillations, that are characteristic of the weakly-correlated phase, survive as a periodic modulation also in the strongly-correlated phase. The latter still persist for sizable interaction strength, as quantum fluctuations significantly reduce the formation of the correlation gap.

With an eye on acene-molecules, we expect that our results have the following implications. We have performed model-calculations employing short-range interactions. They may apply to acenes in the presence of a screening environment, such as substrates or electrochemical solutions. With screened Coulomb-interactions, the effective value of $U$ should be significantly smaller than the molecular bandwidth, implying that our weakly-correlated scenario should apply. In contrast, for gas-phase molecules screening is very weak and the applicability of our model-calculation is not immediately guaranteed. This is the situation investigated in Refs. wilhelm16, ; bettinger16, ; yang16, ; we conclude with a brief overview.

Concerning Ref. wilhelm16, , it remains to be seen if quantum-fluctuations can restore the spin-rotational symmetry also in the presence of a long-range Coulomb interaction. Very recent quantum-chemistry calculations on the level of density-functional multi-reference theory (DFT/MRCI) for oligoacenes (up to $N_{\mathrm{R}}{=}9$) and pp-RPA (up to $N_{\mathrm{R}}{=}12$) may give an indication that this indeed could be the case.bettinger16 ; yang16

The authors of Ref. bettinger16, confirm earlier findingsmondal09 that best agreement between DFT/MRCI results and experimental IR-excitation energies is achieved for closed-shell reference states (RB3LYP) even though spin-unrestricted calculations may produce a lower total energy. It is thus suggested that quantum-fluctuations are indeed strong in long acene-wires so that single-determinant approximations of the ground-state are not representative. In agreement with this interpretation, these authors observe that the weight of the reference-state in the DFT/MRCI ground state is progressively decreasing with growing wire length, so that more and more Slater determinants have to be included to accurately represent this state.

The authors of Ref. yang16, come to a qualitatively similar conclusion. The pp-RPA method is an interesting variant of RPA-type approaches, since it incorporates vertex-corrections that are missed, e.g., in standard RPA bubble-summations. Thus, it allows to include a certain amount of quantum fluctuations when treating unscreened wires (i.e. with Coulomb-interactions) that would otherwise be too long for more accurate approaches. A certain drawback of the method is is that it is uncontrolled: while it can singalize qualitatively that quantum fluctuations are likely to be important, there is no way of telling with confidence what they actually do in that case, quantitatively.

On the qualitative level, the conclusions of Ref. bettinger16, and Ref. yang16, are consistent. In particular, also the authors of yang16, emphasize that the multi-determinant character of the groundstate increases rapidly with growing wire length. It hinders the formation of long-range magnetic order and implies that mean-field results (i. e. unrestricted open-shell calculations) are to be interpreted with a grain of salt.

It remains to be seen, however, if restoring the spin-symmetry in the ground state is enough to also restore the IO. Namely, even assuming that long-range and short-range models of molecular wires share the same phases, it is still unclear whether the effective value of $U$ in the gas-phase calculations is small or large as compared to the (non-interacting) band-width.

Interestingly, the authors of Ref. yang16, propose an exponential decay of the lowest lying spin-triplet gap with increasing $N_{\mathrm{R}}$ saturating near $N_{\mathrm{R}}{=}11$ at roughly 0.13 eV without a trace of IO. If indeed correct, this finding could still be interpreted within the short-range model. It could be taken as a hint that the effective interaction $U$ as it appears in the study Ref. yang16, should be considered large. However, the authors’ data underlying their exponential fit (Table 1 in Ref. yang16, ) appears to be inconsistent with an exponential asymptotics. In fact, the data for the spin-restricted cases exhibits even non-monotonic behaviour (upturn at large wire length), which is qualitatively consistent with the onset of IO as predicted by us.

Acknowledgements.

We express our gratitude to Holger Bettinger, Reinhold Fink, Steve Kivelson and Jan Wilhelm for enspiring discussions. Also, we thank Jan Wilhelm for sharing a manuscript wilhelm16 prior to publication. PS and RT are supported by DFG-SFB 1170 and the ERC starting grant TOPOLECTRICS (ERC-StG-336012). FE acknowledges support of the DFG under projects EV30/8-1 and EV30/11-1. Most of the calculations were performed on the compute cluster of the YIG group of Peter Orth.

Appendix A Non-interacting band-structure and an effective model

We consider a non-interacting ($U=0$) polyacene wire with the lattice and single particle Hamiltonian of (1) explained in Fig. 1. For this model, we derive analytic expressions for the wave number of the nodal point $k_{\mathrm{D}}$, Fermi velocity $v_{\mathrm{F}}$, and the oscillation period of gaps of finite chains. The ladder-like lattice features a symmetry with a line group D*∞h*. An important member of the group is a mirror symmetry with respect to the interchange of stringers which we call parity. The third-neighbor hoppings can be divided in two groups: a hopping that connects atoms related by the parity and the two hoppings that connect different rungs (Fig. 1). We denote the former by $t^{\prime\prime}_{\perp}$ and the latter by $t^{\prime\prime}_{\times}$. We label the 4 Bloch bands $\epsilon_{sb}(k)$ by the wave-number in units of $a^{-1}$ (the inverse lattice spacing), a parity $s=\pm 1$ and a remaining band index $b=\pm 1$. The explicit expression for the band dispersions is given in (2). The bands are particle-hole symmetric only if $t^{\prime}=0$, whence $\epsilon_{sb}(k)=-\epsilon_{-sb}(k)$.

The two bands with $b=-1$ give rise to a linear dispersion $bv_{\mathrm{F}}(k-k_{\mathrm{D}})$ at the Fermi level if the condition $tt^{\prime\prime}_{\perp}>0$ is met, as we will show below. We derive analytic exressions for $k_{\mathrm{D}}$ and $v_{\mathrm{F}}$ under the condition $t^{\prime}=0$. The nodal (“Dirac”) point $k_{\mathrm{D}}$ for $t^{\prime}=0$ is the solution of the equation

[TABLE]

In the limit of $(\pi-k_{\mathrm{D}})$ small the nodal (Dirac) point and the Fermi velocity are expressed by

[TABLE]

For oligoacenes within the tight-binding approximation parameterized as above, the band crossing implies gap oscillations with the unit cell period of

[TABLE]

In the equations (6-8) the symbol $\mathcal{O}(t^{\prime\prime}/t)^{2}$ denotes all terms that are of second order in $t^{\prime\prime}_{\perp}$ and $t^{\prime\prime}_{\times}$. We see that the band crossing sets in with $\sqrt{t^{\prime\prime}_{\perp}}$, provided that the signs of $t$ and $t^{\prime\prime}_{\perp}$ are equal.

The hopping $t^{\prime}$ preserves the band crossing and weakly renormalizes $k_{\mathrm{D}}$ and $v_{\mathrm{F}}$. We remark that the low-energy behavior of the non-interacting polyacene model used in the present paper can be mapped onto a simpler effective model with $t^{\prime}=t^{\prime\prime}_{\times}=0$ and $t$ and $t^{\prime\prime}_{\perp}$ non-zero, i.e. a ladder model with alternating rungs. This model has been introduced originally by Kivelson and Chapman Kivelson .

Appendix B Extracting the bulk gap from finite-size wires

Assuming that the gap opening only affects the band structure close to the Dirac point, one would expect that the gaps are given by a minimal model of a type

[TABLE]

that follows the zone-folding argument of Ref. korytar2014, . The sqrt term in Eq. (9) accounts for the level crossing with parameters $\Delta(N_{\mathrm{R}})$, $v_{\text{F}}$, and $\phi$; $\Delta(N_{\mathrm{R}})$ and $\phi$ are fit-parameters that accommodate additional finite-size effects.

Somewhat unexpectedly, our calculations indicate that the band-structure underlying Eq. (9) does not provide a faithful description of the true bandstructure of the infinite wire with the range of $U$-values that we have been able to investigate. At values $U\sim 1$ the band-structure deviates significantly from the massive Dirac-shape and Eq. (9) no longer applies. (For instance we find a stronger damping of the oscillations as expected from Eq. (9) indicating an interaction induced flattening of the band.)

In the limit of large interaction $U$ the finite size oscillations are completely washed out and the thermodynamic limit of the interaction induced gap can be extracted from a linear fit in the inverse system size $1/L$, see Fig. 6. However note, that for $U=2.5$ an oscillatory part is already visible for not too large system sizes. Therefore we take an heuristic approach modeling the decay of the oscillations with a power law decay

[TABLE]

as shown in Fig. 7. Here, $a/L$ describes the decay of the non-oscillatory part, while the cosine term describes a decay of the oscillation faster than $1/L$, i.e. $\alpha>0$. $b$ gives the amplitude of the oscillations, $k_{\mathrm{D}}$ gives periodicity and $\eta$ allows for the shifting the phase of the oscillation. $\alpha$ is an heuristic exponent modelling the damping of the oscillations. In our fits we always find $\alpha>0$ and $\alpha$ approaches zero for $U\rightarrow 0$.

The gaps shown in Fig. 2 are obtained from the linear fits for $U>2.0$, see Fig. 6, and from fits of Eq. (11) for $U\leq 2.0$. We mention that the fit for $U=1.0$ is still sensitive to the fit range, and the result should be taken with care. For this reason we refrained from including results for $U<1.0$. In addition, in that regime it is not sufficient to fit a single cosine, as the results are getting close to the saw tooth like behavior as in the non-interacting case.

In order to fit the results for small $U$ we introduce a (truncated) saw tooth function

[TABLE]

and use it instead of the cosine in Eq. (9). Here we imposed a limit of n=17. The results for $U=0.2$, $0.5$ are displayed in Fig. 8. While the fits do capture the oscillation quiet well, the finite size oscillations are orders of magnitudes larger than the extracted gaps. In return we do not trust the extracted values and did not include them in Fig. 2.

Appendix C Methodological details

C.1 Numerical Hartree-Fock

Due to SU(2) symmetry combined with a local on-site interaction, we can restrict our HF equations to the Hartree terms

[TABLE]

where $\hat{H}^{U=0}$ denotes the non-interacting system according to Eq. (1), $\sigma$ denotes up and down spins, $-\sigma$ the opposite spin to $\sigma$, and $\bar{n}_{\sigma,x}$ the expectation value of the local density operator $\hat{n}_{\sigma,x}$ at site $x$ with spin $\sigma$. We have checked that by breaking the remaining $S_{\text{tot}}^{z}$ conservation, therefore allowing for Fock terms of the form $\hat{c}^{\dagger}_{\sigma,x}\hat{c}_{-\sigma,x}$, we only find solutions where the spin quantization axis is rotated with respect to the Hartree solution. Note that in this case the solution depends on the initialization of the self consistency loop.

Since up and down spins do not mix in (13), we can represent the up and down spin sectors by independent matrices. We start the HF calculations by first diagonalizing the non-interacting system in order to obtain the ground state energy $E_{0}$ as a reference. Note that in this case we obtain a homogeneous solution $\bar{n}_{\sigma,x}=0.5$. We then explicitly introduce a staggered magnetic order by setting $\bar{n}_{\sigma,x}=0.5\pm 0.01\sigma$, where we take ’$+$’ on one sublattice of the bipartite system and ’$-$’ for the other sublattice (Fig. 1). With these initial values we diagonalize ${\cal{H}}^{\mathrm{HF}}$ and calculate the locate densities $\bar{n}_{\sigma,x}$ which we then iteratively insert into (13). We perform at least five iterations and continue until the ground state energy $E^{\mathrm{HF}}$ of Hamiltonian ${\cal{H}}^{\mathrm{HF}}$ changes by less then $10^{-10}t$. In order to avoid getting stuck in an oscillation between two solutions, we damp the HF self consistency loop by taking an average of 0.7 times the new density and 0.3 times the densities of the preceding iteration.

Finally, we check whether the ground state energy $E^{\mathrm{HF}}$ including the double counting corrections

[TABLE]

is smaller than the ground state energy $E_{0}$ of the non-interacting system. Indeed we always find $E^{\mathrm{HF}}\leq E_{0}$ The spin-averaged local density $\bar{n}_{\uparrow,j}+\bar{n}_{\downarrow,j}$ of our HF solution is homogenous and the total staggered magnetisation $M^{\mathrm{HF}}$ is given by the difference between the magnetization of the two sub lattices

[TABLE]

C.2 Hartree-Fock: commensurate case

The main problem of obtaining the $U\ll 1$ HF solutions in Fig. 2 roots in the finite size oscillations being incommensurable. Therefore we need very large system sizes in order to start with an $U=0$ finite size gap $\Delta_{0}$ being smaller than the interaction induced gap in the HF calculations. This extraction of the HF gap gets simplified if we switch to the case of commensurate oscillations by choosing a fine-tuned $t^{{}^{\prime\prime}}$ from our solution Eq. (5). In Fig. 9 we display the results in the case of a periodicity of twenty rings. We obtained the corresponding $t^{{}^{\prime\prime}}=0.13$ by solving Eq. (5) numerically.222For reference, the precise value used in the numerics is $t^{\prime\prime}=0.130064666287555$. Since we now always hit a minimum of the gap oscillations precisely, i.e. a zero in the non-interacting case, we can extract the HF gap without difficulties. As expected, the HF gap indeed opens $\Delta^{\text{HF}}_{1,N_{\mathrm{R}}}\sim U^{2}$. In addition, we see the same exponential behavior for intermediate to large $U$ as displayed in Fig. 2.

C.3 DMRG

Within our DMRG calculations we keep up to 8000 states per block; the dimension of the target space grows up to $4.6\cdot 10^{7}$, the discarded entropies are considerably below $10^{-4}$ for $U>1.0$, while for $U\leq 1.0$ the discarded entropies reach $4.2\cdot 10^{-2}$ for systems with $N_{\mathrm{R}}>20$ in each DMRG step, typically significantly smaller. We employ at least 7 finite lattice sweeps with an optimized infinite lattice warm up for the smallest system sizes. In order to reach large system sizes we iteratively restart a $N_{\mathrm{R}}$ simulation in order to investigate a $N_{\mathrm{R}}+1$ ring. Note that for this we only have to deactivate the wave function prediction in the initial step of a restart. In order to initialize the simulations we apply a sliding environment block approach, that we already employed successfully to fractional quantum hall systems.Hu:PLA12 ; Johri:PRB14 ; Johri:NJP16 Instead of using a reflection or a growing right block during the DMRG infinite lattice warm up sweeps, we employ a right block of up to 5 sites, which we can always build exactly in such a way that we always have $4N_{\mathrm{R}}$ sites for PBCs and $4N_{\mathrm{R}}+2$ sites for HWBCs. With this initial warm up, we avoid running into excited states in the initial sweep.