text:
# Magnetic ordering, electronic structure and magnetic anisotropy energy in the high-spin Mn10 single molecule magnet

## Abstract

We report the electronic structure and magnetic ordering of the single molecule magnet \[Mn10O4(2,2’-biphenoxide)4Br12\]4- based on first-principles all-electron density-functional calculations. We find that two of the ten core Mn atoms are coupled antiferromagnetically to the remaining eight, resulting in a ferrimagnetic ground state with total spin $`S=13`$. The calculated magnetic anisotropy barrier is found to be 9 K in good agreement with experiment. The presence of the Br anions impact the electronic structure and therefore the magnetic properties of the 10 Mn atoms. However, the electric field due to the negative charges has no significant effect on the magnetic anisotropy.

The interest in magnetic molecular clusters of transition metal ions has been continuously growing since the observation of magnetic bi-stability of a purely molecular origin in the so-called Mn12-ac Mn12 , which shows a magnetic hysteresis cycle below 4 K similar to that observed for bulk magnetic materials. The magnetic bi-stability associated with the hysteresis cycle has created an interest in these clusters for information storage, although at low temperature quantum effects affect the reversal of the magnetization, resulting in steps in the hysteresis Mn12-steps . This phenomenon of quantum tunneling of magnetization is governed by the magnetic anisotropy energy (MAE)vFleck  barrier which is due to directional dependencies of the spin-orbit-coupling operator.

Recently, Pederson and Khanna have developed a method for accounting for second-order anisotropy energiesped-so . This method relies on an exact and simple expression for the spin-orbit coupling operator used in a second-order perturbative treatment to determine the dependence of the total energy on spin projection. Initial applications to the uniaxial Mn12-ac lead to a density-functional-based second-order anisotropy energy ped1  of 55.7 K, in agreement with the experimentally deduced values Mertes ; barra-Mn12  of 54.8(3) K or 55.6 K.

Because the second-order anisotropy energy scales with the square of the magnetization it was generally believed that a high-spin ground state $`S`$ would be beneficial for a large barrier. The \[Mn10O4(2,2’-biphenoxide)4Br12\]4- cluster has been reported to have a $`S=12`$ high-spin ground state Mn10  but only a small energy barrier of about 7.7 K. In this work we investigate the electronic and magnetic properties and the magnetic anisotropy energy of this high-spin single molecule magnet. The information obtained here may be useful in the search for single molecule magnets with a greater magnetic anisotropy.

Fig. 1 shows the structure of the \[Mn10O4(2,2’-biphenoxide)4Br12\]4- molecular magnet. The 10 Mn atoms form a tetrahedron like structure with Mn atoms at the corners and at the middle of the tetrahedron edges. Two of the Mn atoms, the top and the bottom spheres in Fig. 1, are coupled antiferromagnetically to the rest of the Mn atoms. The Mn atoms are bridged by O atoms. The magnetic core is further stabilized by organic rings that are also connected to the O atoms. The negatively charged cluster is compensated by \[(CH3CH2)3NH\]2\[Mn(CH3CN)4(H2O)2\] in the molecular crystal, but experimental results suggest that the magnetic anisotropy is due to the localized valences of the 10 Mn atomsMn10 . In order to use the higher symmetry and make the problem computationally feasible, we carried out our calculation on the negatively charged \[Mn10O4(2,2’-biphenoxide)4Br12\]4- cluster which contains 114 atoms. The eight symmetry operations reduce the complete cluster to 18 inequivalent atoms.

The theoretical studies were carried out using a linear combination of atomic-orbitals molecular-orbital (LCAO-MO) approach within the density-functional framework hohenberg ; kohn . The molecular orbitals were expanded as linear combinations of Gaussian functions centered at the atomic sites. The calculations were carried out at the all-electron level and the multicenter integrals required in the solution of the Kohn-Sham equation were calculated by integrating numerically over a mesh of points pederson1990 .

Our density functional-based calculations were performed with the all-electron Gaussian-orbital based Naval Research Laboratory Molecular Orbital Library (NRLMOL) program pederson1990 ; jackson1990 ; briley ; porezag1996 ; pore99 ; pss2000 , using the Perdew-Burke-Ernzerhof (PBE) generalized-gradient approximation for the exchange and correlation functionalpbe . NRLMOL combines large Gaussian orbital basis sets, numerically precise variational integration and an analytic solution of Poisson’s equation to accurately determine the self-consistent potentials, secular matrix, total energies and Hellmann-Feynman-Pulay forces. The exponents for the single Gaussian have been fully optimized for DFT calculations pore99 . The basis set for the Mn10 cluster consisted of a total of 3756 contracted orbitals. The minimum and maximum exponent of the bare Gaussians, the number of bare Gaussians, and the number of contracted $`s`$-, $`p`$\- and $`d`$-like basis functions are given in Table 1 for each atomic species. The contraction coefficients for atomic orbitals were obtained by performing an SCF-LDA calculation on the spherical unpolarized atom where the total energy of the atom was converged to within 10 meV. The basis functions that do not correspond to atomic wavefunction were constructed from the longest range bare Gaussians in the basis set.

Here we repeat some of the formulas needed for discussion of the magnetic anisotropy energy. The same definitions and notation are used as in Ref. ped-so, . In the absence of a magnetic field the second-order MAE $`\mathrm{\Delta }_2`$ resulting from the spin-orbit coupling, for an arbitrary symmetry, reduces to

$$\mathrm{\Delta }_2=\underset{\sigma \sigma ^{}}{}\underset{ij}{}M_{ij}^{\sigma \sigma ^{}}S_i^{\sigma \sigma ^{}}S_j^{\sigma ^{}\sigma },$$
(1)
which is a generalization of Eq. (19) of Ref. ped-so, . The matrix elements $`S_i^{\sigma \sigma ^{}}=\chi ^\sigma |S_i|\chi ^\sigma ^{}`$ implicitly depend on two angles ($`\theta ,\beta `$) defining the axis of quantization. The matrix elements $`M_{ij}^{\sigma \sigma ^{}}`$, which are related to of the induced orbital moment, are given by

$$M_{ij}^{\sigma \sigma ^{}}=\underset{kl}{}\frac{\varphi _{l\sigma }|V_i|\varphi _{k\sigma ^{}}\varphi _{k\sigma ^{}}|V_j|\varphi _{l\sigma }}{\epsilon _{l\sigma }\epsilon _{k\sigma ^{}}},$$
(2)
where $`\varphi _{l\sigma }`$, $`\varphi _{k\sigma }`$ and $`\epsilon _{l\sigma }`$ are, respectively, the occupied, unoccupied and the corresponding energies of states. $`V_i`$ is same as defined in Eq. (7) of Ref. ped-so,  and is related to derivatives of the Coulomb potential. The matrix elements can be evaluated by integrating products of the Coulomb potential with partial derivatives of the basis functions. This procedure avoids the time consuming task of calculating the gradient of the Coulomb potential directly.

In addition to the magnetically interesting complex, the crystal also contains single Mn complexes to balance the charges. Using high-field EPR spectroscopy Barra et al. Mn10  found that this \[(CH3CH2)3NH\]2\[Mn(CH3CN)4(H2O)2\] unit is paramagnetic. We also find that this unit is paramagnetic, with the Mn atom in a $`+2`$ charge state and a spin of $`S`$=5/2. The complex exhibits easy-plane behavior with an energy well of 0.1 K. We therefore focus the remainder of our work on the Mn10 unit only.

Calculations on a $`S`$=12 high spin state revealed that this spin state would not be magnetically stable because the Fermi levels in the majority and minority spin channel would not be aligned. The Fermi level misalignment indicated further transport of electrons of the minority to the majority spin channel. As a result we obtained a $`S`$=13 high spin state as the magnetic ground state instead of the $`S`$=12 state obtained from high-field EPR spectroscopy. Our result is consistent with experiment since it is difficult to differentiate experimentally between the two possibilities roberta .

A plot of the spin density shown in Fig. 1 clearly confirms the antiferromagnetic coupling of two of the Mn atoms to the remaining Mn atoms. The spin density around the minority spin Mn and the four majority spin Mn at the corners of the tetrahedron-like magnetic core show a nearly spherical spin density which one would expect for a closed $`d^5`$ shell. The other four Mn atoms which are on the edges of the tetrahedron (large dark spheres) show a less spherical spin density, indicating another charge state for these atoms. In order to analyze the magnetic ordering we calculated the spin density in spheres around the atoms. This gives a measure of the localized moment at the atom, but will generally underestimate the exact value. Using spheres with a radius of 2.23 $`a_B`$ around the Mn atoms we obtained for the three non-equivalent Mn atoms (majority spin tetrahedron edge, majority spin tetrahedron vertex, and minority spin tetrahedron edge) a local moment of 3.61, 4.33 and -4.33 $`\mu _B`$, respectively. This result suggests an ionic picture that the first Mn has an Mn3+ ($`S=2`$) state, whereas the other two are Mn2+ ($`S=5/2`$). This picture is fully in accord with the spin density plot in Fig. 1. Due to the symmetry of the cluster, the two types of majority spin Mn atoms have a multiplicity of 4 whereas the minority spin Mn atom has a multiplicity of 2, resulting in the previously mentioned $`S=4\times 2+4\times 5/22\times 5/2=13`$ magnetic ground state.

The electronic density of states (DOS) for the majority and minority spin channels is shown in Fig. 2. For each spin, the DOS is further decomposed into the $`3d`$ contributions of all Mn atoms. It is evident from the plot that the states around the Fermi level are clearly connected with $`3d`$ states of the Mn atoms. This result also agrees well with the experimental picture that the states near the Fermi level are well localized and do not show strong hybridization with other atoms, although we find some oxygen, bromine and nitrogen contributions for the occupied states.

Starting from the experimental geometry roberta  we carried out about 30 steps of a conjugate-gradient algorithm using the Hellmann-Feynman-Pulay forces for optimization of the geometry. For each new geometry we calculated the complete Hamiltonian of the magnetic anisotropy and the second-order contribution to magnetic anisotropy barrier $`DS_z^2`$. In accord with experimental data, we find that the Mn10 single molecule magnet is an easy-axis system. The barrier showed no strong dependence on geometry varying between 8.8 K and 10.4 K with a value of 9.5 K for the lowest energy geometry. Expressing the barrier in the form of a simple spin Hamiltonian $`H=DS_z^2`$ we obtained a value of $`D=0.056`$ K. For the calculation of the spin orbit matrix elements, we included all valence electrons and all unoccupied states in an energy window of 13.6 eV above the highest occupied state. The difference between the second-order treatment and exact diagonalization of the Hamiltonian including the spin-orbit matrix elements, which includes some higher order effects too, was of the order of 0.1 K.

Eq. 2 shows that the barrier is related to matrix elements between occupied and unoccupied orbitals in the majority and minority spin channels. In order to give a deeper insight into which states are forming the barrier, we analyze these contributions in more detail. First, we focus on the contributions of the different spin channels. Table 2 summarizes the result in form of the $`D`$ parameter allowing only a given spin channel, for example including only matrix elements $`M_{ij}^{\sigma \sigma ^{}}`$ between occupied majority states and unoccupied minority states, in the calculation of the barrier. All matrix elements from the occupied majority electrons prefer an easy-axis system, whereas the matrix elements from the occupied minority spin channel would result in an easy-plane system. Only due to the larger values of the contributions of the occupied majority spin channel the system ends up as an easy-axis system. This destructive interference between the different spin channels seems to be the reason for the small barrier compared with Mn12. In Mn12 constructive interference between the different spin channels was observed.

Besides the spin channel contribution, we can analyze which electronic states contribute most to the matrix elements $`M_{ij}^{\sigma \sigma ^{}}`$. In Fig. 3 we display plots of the square of the wavefunctions of the occupied majority state and the unoccupied minority state which contribute to the matrix element $`M_{ij}^{\sigma \sigma ^{}}`$ with the largest absolute value. The view is from the top with respect to Fig. 1. It is clearly visible that the states of interest are $`d`$-states localized at the same Mn atom. In this case, the states are localized at the minority spin Mn atoms (light spheres in Fig. 1). In order to emphasize the $`d`$ character of the wavefunctions, we have chosen the top view, although the wavefunctions of the other minority Mn atom are just below the top ones and are not visible.

While Mn is the only magnetically active species in the complex, the remaining atoms affect the magnetic properties of the molecule. In particular the electric field of the twelve Br- ions can affect the MAE through its effects on the electronic structure and on the spin-orbit coupling. Effects of the Br ions on the electronic structure are of a chemical nature, and a detailed analysis is beyond the scope of this work. However, direct effects on the spin-orbit coupling energy could raise the possibility that even small variations in the positions of these ions could affect the magnetic properties of the molecule.

To measure the effects of the Br- ions on the electronic structure we redid the calculations with various subsets of the Br- ions. For each removed Br we also removed an extra electron, keeping the remaining molecule isoelectronic with the original complex. For these systems with either zero, four, or eight Br atoms, we observed a range of behaviors. In some cases the electronic structure near the Fermi level was similar to the original molecule, although it never showed truly rigid-band-like behavior. In other cases, however, the electronic structure changed significantly, sometimes completely closing the HOMO-LUMO gap. Associated with these electronic structure changes were large changes in the MAE, including changes in the magnitude of the anisotropy barrier, as well as instances of changes to an easy-plane system.

To measure the direct effect of the electric field of the Br ions on the spin-orbit coupling, external Coulomb potentials which acted to cancel the long range affects of the Br anions were added and the spin-orbit interaction and magnetic anisotropy were recomputed. This neutralized the electric field due to the Br- ions near the Mn sites without changing the electronic structure of the molecule. We tested the effects with neutralizing charge distributions of various widths, and by neutralizing four, eight, or all twelve Br anions. The MAE changed by less than 1 K in all of these calculations. We therefore conclude that the electric fields created by the Br- ions do not have a significant effect on the magnetic properties of the molecule.

In conclusion, we present a study of the electronic and magnetic properties of the Mn10 single molecule magnet. We confirm the experimentally suggested magnetic ordering, although we find that a state with $`S=13`$ is the magnetic ground state in contrast to the $`S=12`$ state suggested from high field EPR measurementsMn10 . In agreement with experiment we find the Mn10 unit is an easy-axis system with a small barrier of 9.5 K and the compensating cluster in the molecular crystal, which has one Mn atom with $`S=5/2`$, is an easy plane system with a MAE of 0.1 K, negligible compared with the Mn10 unit. We show that the magnetic anisotropy is determined by a competition between different spin channels involved. The electric field caused by the negative charges of Br anions has no significant direct effect on the spin-orbit coupling or the MAE, although their chemical interactions do have significant effect on the electronic structure and therefore on the MAE. The states most important for the magnetic anisotropy energy involve transitions between occupied majority and unoccupied minority $`d`$-states at the same Mn-atom.

J.K. would like to thank the Schloeßmann Foundation for financial support and R. Sessoli for helpful discussions. The work was partially supported by the Office of Naval Research and the DOD HPCMO.


format: text

text:
# A Bell telegraph

## 1 Introduction

“There are no Bell telephones”, says an old joke in the community of foundations of quantum mechamics, meaning that one cannot (should not) be able to use Bell’s nonlocality (collapse of entangled quantum states at a distance due to local measurements, of the kind of the EPR argument ) for transmitting information. There is no general proof of this statement, and that’s fortunate since Bell himself claimed once that all an impossibility theorem proves is lack of imagination. However, if quantum mechanics admitted such information transmission to be performed, it would enter in direct conflict with very firmly established consequences of special relativity. In such a case one of these theories (or both) would be in serious trouble, and the main suspect would be the Copenhagen interpretation of quantum mechanics.

What we present here is an apparent “Bell telegraph” and the reasons why it does not work. It is a replacement of a “youth sin” version from 2002, in which I posed it in a “paradox format”, expecting to spark discusions and reactions (I did, indeed). But, to my surprise, the manuscript attracted also some attention from the media, where it was regarded not as a theoretical puzzle to be solved but as an actual “invention”, a proposal to actually communicate instantly at arbitrary distances 111see for instance http://www.freerepublic.com/focus/news/735915/posts, last accessed 10/4/2017. I hope this version will frame the discussion correctly as an effort in thinking about the interpretations of quantum mechanics.

This manuscript is organized as follows: the telegraph is first presented exactly as in the 2002 version. In the next section (also exactly that of the 2002 version) the conflict with special relativity is discussed. Then, in section (4) (which is new) an explanation of why the telegraph does not work is shown (you might challenge yourself trying to understand it before reading that section). We also briefly discuss to what extent the ‘no Bell telephone” conjecture can be considered a feature of quantum mechanics under the Copenhagen interpretation. Finally, we discuss the description of collapse in a relativistic realm and pose some questions.

## 2 The device

A sketch of the device is shown in the figure above. It is inspired in Feynman’s discussion on quantum interference . A beam of excited atoms (which will emit a two photon Bell state) is divided in two by a potential barrier in such a way that a given atom has the same probability of following each path, which is itself a Bell state. At the decay, the resulting photons are collected by some kind of pipes (for instance, two optic fibers), so every pair of photons is in any of those, forming a kind of doubly entangled state. If no one collapses this system, the photons at one end of the pipes can be set to interfere at a screen. But one could ’see’ where each pair is just looking at the second photon at the other end, destroying the interference pattern. So, turning on and off detectors at one end sets the other end at two observably different states (random arrival and interference respectively), which can obviously be used as a telegraph.

Still we have to show that this can indeed be used to send information in a time as short as necessary to be considered instantaneous. We need several photon pairs (say M) to conclude that we are before an interference or not, with a given confidence interval. One solution is to use many uncorrelated pairs simultaneously, and invoke superposition principle. But in order not to obscure the conclusions by adding hypothesis, we will give an alternative proof using single pairs. Suppose that for a given telegraph we need a time separation T for producing each pair being sure that we deal with a single pair at a given instant. So we need a time MT for sending a ”byte”. Now take N such telegraphs, in which the first pair in each is produced in an instant at random between 0 and T. With this ensemble we build an improved telegraph, turning on and off all detectors simultaneously for sending signals. Now the time needed reduces to MT/N, and can be made as little as we wish. Thus we can send information instantly (at the large N limit) using (lots of) single pair states.

## 3 Troubles

Let’s suppose for the moment that we haven’t make any mistakes; if so, quantum mechanics imply we can send information instantly. What do we mean by ”instantly” in a relativistic context? Or, put differently, instantly in what reference frame? This is a difficulty with the very concept of quantum collapse. Even if it turns that no information could be transmitted via Bell states, this is still a conceptual puzzle. Whatever the answer is, it fits into either of the following two posibilities:

* Deny relativity: There is a privileged frame F in which quantum collapse occurs instantly. This possibility is not inconsistent but deeply nasty.
* Preserve relativity: The frame F in which the collapse is instantaneous depends on the state (for instance, it could be the one in which the center of mass is at rest). This possibility seems more acceptable, but a closer look shows it leads to inconsistencies, which can be seen as paradoxes. Consider the setting shown in the figure below: two identical telegraphs A and B with F(A) and F(B) moving in opposite directions results in sending a signal to the past. As F(A) moves to the left, A reception happens before A emmision in this frame. The message is instantly retransmitted trough B (B emmision) and again, as F(B) moves to the right, B reception happens before B emission, and at the same space location than A emmision. In the setting, the automaton located there transmit $`m_1`$ if and only if it transmit $`m_2`$ and vice versa, which is absolutely inconsistent.

Of course, there remains the possibility that we have made some mistake and quantum mechanics really forbids such a device to work. Hopefully this work will call someone’s attention in order to find it out.

## 4 The resolution and the “astuteness” of the Copenhagen interpretation

When calculations are carried carefully, the result is that there is never interference in the screen, irrespective of what is done in the detectors. Call the “pipes” $`1`$ and $`2`$, call $`L`$ the side of the screen and $`R`$ the side of the detectors. Observe that prior to any measurement, the quantum state of the system is

$$|\mathrm{\Psi }_0=\frac{1}{\sqrt{2}}\left(|1L|1R+|2L|2R\right)$$
(1)
where $`|iL`$, $`|jR`$ are supposed to be normalized, localized, time dependent and mutually orthogonal state vectors. If we call $`|x`$ the eigenstates of position on the screen, the amplitude on it will read

$$x|\mathrm{\Psi }_0=\frac{1}{\sqrt{2}}\left(x|1L|1R+x|2L|2R\right)$$
(2)
which is still a state vector. The probability density of finding a photon in a position on the screen is thus its modulus:

$`p(x)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(|x|1L|^21R|1R+|x|2L|^22R|2R\right)`$ (3)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(|x|1L|^2+|x|2L|^2\right)`$ (4)
where we used the orthonormality of the vectors $`|1R`$ and $`|2R`$. Thus, even if we do not collapse the wavefunctions with the detectors, the pattern is just a statistical average showing no interference. Of course, if the detectors on the right operate, the state at the moment of the screening will not be (2); assume the detectors measure the observable O, and that there is a complete ortonormal basis for the observables $`|O_k^{(1)}`$ and $`|O_k^{}^{(2)}`$. If the outcome is for instance the eigenvalue $`O_l^{(1)}`$, then the state will be

$$\frac{1}{\sqrt{2}}x|1LO_l^{(1)}|1R|O_l^{(1)},$$

which will hapen with (the conditional) probability

$$P\left(x|O_l^{(1)}\right)=\frac{1}{2}x|1L|^2O_l^{(1)}|1R^2.$$

But since at the screen there is not information about what happened at the detectors, the probability of finding the photon at $`x`$ will be the sum on all outcomes on the detector, and using the completeness of the base and its orthonormality the probability (4) is recovered.

It should be noticed that the “no Bell telegraph” feature is far from being built-in the formalism, which is indeed highly non-local, but at the same time very astute. It comes as somewhat of a miracle that the entanglement, which provides the possibility of influencing at a distance, provides also an orthogonality preventing the system of making observability in one end dependent on what is done in the other. To appreciate the extent of this astucy, let us consider the following “upgraded telegraph”: instead of trying to destroy interference with the detectors, we try to recover the interference by re-joining the path of the $`R`$-side photon with a system of mirrors, with the hope that ending with the physical separation we will make the $`R`$-states non-orthogonal. This does not work, provided the evolution operator is unitary, thus preserving the inner products in the state space. This, in turn, is a consequence of the hermiticity of the Hamiltonian, but an actual system of mirrors will be slighty absorptive, and thus non-Hermitian. ¿Can we be absolutely positive that any conceivable system that recombinates the states $`|1R`$ and $`|2R`$ will preserve their orthogonality?

## 5 So, is collapse compatible with special relativity?

Let us suppose that we take for granted that there is no Bell telegraph. We will not find then any observational inconsistency between quantum mechanics and special relativity. But the description of the states will change radically in different reference systems. There is no question that, in a given reference frame, the collapse is to be described as instantaneous. But then what measurement collapses the state will depend on the observer, and it seems impossible to avoid this ugly non-covariance at the level of the description of quantum collapses. We might wonder if we aren’t mixing theoretical frames by placing entanglement and collapse in a relativistic limit, but we can see easily that such a caveat is baseless: relativistic quantum field theory is a standard quantum theory whose state space can be construed as a system of an unlimited number of identical particles, and so entangled states are part of the state space of the relativistic field. Measurement in the context of relativistic fields has been pursued in the heydays of quantum theory, but to our knowledge it has not been developed in order to incorporate the developements involving entanglement.

Let us consider two remote detectors measuring a system of entangled particles in two laboratories $`L`$ and $`R`$ separated by a space-like four-vector, as shown in figure (3). A historian $`H`$ is placed in a common region of the future light cones of $`R`$ and $`L`$.

An entangled state is prepared in the state $`|L|R`$. At $`L`$ a detector will measure one from the complete set of states $`|O^{(L)}_k`$, while at $`R`$ another detector will measure one from $`|O^{(R)}_k^{}`$. In terms of these states, the entangled state can be cast in the form

$$\left(\mathrm{\Sigma }_kO_k^{(L)}|L|O_k^{(L)}\right)\left(\mathrm{\Sigma }_k^{}O_k^{}^{(R)}|R|O_k^{}^{(R)}\right)$$

Let us suppose the historian determines that at $`L`$ the value $`O_A^{(L)}`$ has been obtained, while on the left the outcome is $`O_B^{(R)}`$. He will thus conclude that the state after both laboratories performed their measurements is $`|O_A^{(L)}|O_B^{(R)}`$. What was the state at time zero in the figure?

In the figure of the left the movement of the observer is such the measurement at $`L`$ happen at $`t<0`$, while measurement at $`R`$ is at $`t>0`$. Thus, at time zero the state is:

$$|O_A^{(L)}\left(\mathrm{\Sigma }_k^{}O_k^{}^{(R)}|R|O_k^{}^{(R)}\right).$$

On the other hand, for an observer which moves such that the order is the one shown on the right, the state at time zero is

$$\left(\mathrm{\Sigma }_kO_k^{(L)}|L|O_k^{(L)}\right)|O_B^{(R)}.$$

So, can these be actually the states of the system? How can we understand that such state will depend on the movement of the observer? Even if this highly non-covariant description shows no observable non-covariance, it seems difficult to state that the vectors shown correspond to the states of the system under study. It might be more proper to think that it is impossible to separate system from observer, even while no measurement is being performed.

## 6 Concluding remarks: why should we care about all this?

We could separate the caveats about collapsing entangled states in “hard” (observable) possible covariance problems and “soft” (interpretational) worries. If one takes a positivist point of view, only the hard problems should bother us. As we have seen, quantum theory plus Copenhagen interpretation seems astute and, in spite of permitting conspicuous non-local and non-covariant descriptions at the unobservable level of the states, preparing a system that violates covariance at the level of the observables is seemingly impossible. We could try to prove stringent “no-go” theorems, but history warns us that it is impossible to have control over hidden lemmas: conditions assumed implicitly in the proof of such theorems. It is in this sense, I think, that the celebrated quote “all an impossibility theorem proves is lack of imagination” should be interpreted. If the theory actually had such “hard” non-covariance, it should be interpreted, I think, as a formal difficulty of quantum mechanics, since relativistic causality is so firmly established, and not as a prediction of the technical possibility of superluminal communication. It should be said that the possibility of such “hard” noncovariance seems extremely remote.

However, the “metaphysical” soft non-covariance should not be underestimated. Recall that the Viena circle failed in rooting all meaning on “empiric facts” . It is important to bear in mind that quantum collapse remains the most difficult problem in the interpretation of quantum mechanics, and understanding it is unavoidable if the potential practical applications of entanglement (quantum computing, teletransportation, etc) are to work at a scale large enough to be more than expensive toys. In this sense, wondering “what is the state of a system” when the sequence order of collapse is frame-dependent is indeed an intriguing problem. I feel that holding the belief that the state is strictly a property of the system is untenable. Our examples of remote collapse seemingly show that the state vector is a description of the state of the system-observer complex, becoming impossible to separate subject and object even while measurements are not being performed.

## 7 Aknowledgements

I am gratefull to Alejandro Daleo and Monica Manceñido for very useful comments and discussions. I am also grateful to the many that wrote after the 2002 version was released. Regrettably, I did not preserve those e-mails, not even the original e-mail box. I deeply regret to have waited fifteen years to make this replacement. The delay has to do with many personal and national events here in my country.


text:
# 1 Introduction

## 1 Introduction

Several many particle systems are characterized by a steady mass transport. Examples for this kind of systems can be found in biological transport  or vehicular traffic . From a theoretical point of view these processes are of particular interest, because they exhibit generic non-equilibrium behavior. Due to the large number of important applications, many microscopic models for particle transport have been suggested in recent years .

Among these, the most important microscopic model for non-equilibrium particle transport is the so called asymmetric exclusion process (ASEP) . In this model, particles jump on a one-dimensional lattice, either to the right (with probability $`pdt`$) or to the left (with probability $`qdt`$), if the corresponding sites are empty. The model shows a number of generic effects  that are characteristic for non-equilibrium particle transport and maintain for the more specialized variants of the model . At the same time the ASEP is simple enough to obtain several exact results for the system, which is of great importance, because the general theoretical framework of non-equilibrium physics is less developed. Exact results exist, e.g. for the stationary state of the system with periodic  and open boundary conditions . The case of open boundary conditions is of special interest because one observes boundary induced phase transitions . In this paper, we restrict ourselves to the totally asymmetric exclusion process (TASEP), that is $`q=0`$. This case includes the most important phenomena but simplifies considerably the discussion of the model. Open boundary conditions are implemented by two particle reservoirs that are coupled to the chain. The capacities of the reservoirs determine, together with the bulk hopping rates, the actual state of the system.

The complexity of the open system is also reflected in the mathematical structure of the stationary solution. While the steady state of the periodic system is given by a simple product measure , it is highly non-trivial for the open system. Nevertheless, it can be calculated and it is possible to obtain several non-trivial quantities, e.g. current- or density fluctuations in the stationary state , or large deviation functions .

The dynamic properties of the open system are, however, more puzzling. Exact analytical results for the largest relaxation time $`\tau `$ and the corresponding dynamical exponent $`z`$ are so far only possible for the periodic chain by applying the Bethe ansatz . For the open chain estimates for $`\tau `$ can be obtained from a phenomenological approach, that models directly the dynamics of the boundary layer separating the high and low density domains imposed by the particle reservoirs . It has been shown, that this approach gives asymptotically correct results in a certain parameter regime . For finite systems, as well as for general in- and output rates, relaxation times have to be calculated numerically.

This has been done by U. Bilstein and B. Wehefritz , and later by M. Dudzińsky and G.M. Schütz  who calculated the relaxation times by using exact diagonalization. This method is, however, restricted to very short chains (less than twenty), that are not in the asymptotic limit. In this work we used the density matrix renormalization group (DMRG) technique , that enables us to calculate the relaxation times for much larger system sizes, compared to . By treating large chains we have obtained more conclusive results for the dynamical exponent $`z`$ in the maximal current phase and could also achieve convergence between the exact numerical values of the relaxation times and the estimates of the phenomenological approach.

The paper is organized as follows. In the next section we discuss briefly the relevant physical concepts and the applied numerical techniques. In the third section we show the comparison of the domain wall predictions for finite systems with our results. Section 4 is devoted to the special case of the disorder line $`\alpha +\beta =1`$. This section is followed by a discussion of the dynamical behavior when approaching the phase boundaries as well as in the maximum current phase.

## 2 Analytic predictions and the DMRG method for non-equilibrium systems

### 2.1 The TASEP with open boundary conditions

For self-containedness we will repeat the definition of the model. The TASEP is defined on a one-dimensional lattice with $`L`$ sites. The boundary sites of the chain are coupled to two particle reservoirs, one reservoir on the left that controls the particle input and a second on the right that governs the output of particles.

We regard the process in continuous time (see  for a comparison of the different update procedures), which corresponds to a random sequential update in computer simulations. If a link between sites $`i`$ and $`i+1`$ is selected, a particle located at $`i`$ moves to site $`i+1`$ if site $`i+1`$ is empty (for convenience we set the hopping rate to one). In case of choosing the link $`(0,1)`$ one introduces a particle with probability $`\alpha `$ if the first site is empty. Finally a particle may leave the system with probability $`\beta `$, if the link $`(L,L+1)`$ is chosen.

Each configuration $`\sigma `$ can be written in terms of boolean lattice gas variables $`\sigma _i`$, i.e. $`\sigma _i=0(1)`$ if the site is empty (occupied). If we introduce an orthonormal basis $`|\sigma =|\sigma _1,\mathrm{},\sigma _L`$ in the $`2^L`$-dimensional configuration space, we can define the probability vector $`|P(t)`$ as $`|P(t)=_{\{\sigma \}}P(\{\sigma \},t)|\sigma `$. The time evolution of $`|P(t)`$ is determined by means of the master equation, that can be written as a Schrödinger equation in imaginary time :

$$\frac{}{t}|P(t)=|P(t),$$
(1)
where $``$ denotes the stochastic Hamiltonian. The matrix elements of $``$ are the rates $`w(\sigma \sigma ^{})`$ for a transition $`\sigma \sigma ^{}`$. Explicitly $``$ is given by $`\sigma ||\sigma ^{}=w(\sigma ^{}\sigma )`$ for the off-diagonal elements ($`\sigma \sigma ^{}`$) and by
$`\sigma ||\sigma =_{\{\sigma \}\{\sigma ^{}\}}w(\sigma \sigma ^{})`$ for the diagonal elements.

The exact stationary solution of the master equation can, for general $`\alpha ,\beta `$, be written as a product of infinite dimensional matrices . This solution allowed to calculate the $`\alpha ,\beta `$ dependence of the stationary quantities, e.g. the average flux. The results are summarized in the phase diagram of the system shown in fig. 1. Three phases can be distinguished by means of a different functional behavior of the flow. In the low density phases AI and AII, i.e. for $`\alpha <\beta `$, $`\alpha <1/2`$, the flow $`J=J(\alpha )`$ is given by $`J=\alpha (1\alpha )`$ and analogously in the high density phases BI and BII by $`J=\beta (1\beta )`$. If both $`\alpha >1/2`$ and $`\beta >1/2`$ (phase C) the capacity of the particle reservoirs exceeds the capacity of the chain. Then the flux is independent of $`\alpha `$ and $`\beta `$ and given by $`J=1/4`$. This phase is called maximal current phase.

Both the high and low density phases are divided into two subphases. In phase AI and BI, the capacities of both reservoirs are below the chain capacity. In phase AII (BII), only the capacity of the exit (entrance) exceeds the chain capacity. This has, e.g. consequences for the asymptotics of the density profile , and, a question we address in this article, possibly also for the dynamics of the chain. Another important line is given for $`\alpha +\beta =1`$. On this one-dimensional line the stationary solution is much simplified, i.e. it is given as product measure.

Now we shall summarize known results for the relaxation times. For the periodic system, exact results for the energy-gap, i.e. the inverse of the largest relaxation time $`\tau _1`$, have been obtained by Bethe-ansatz techniques . The largest relaxation time scales for any density of the system asymptotically as $`\tau _1L^z`$, where $`z=3/2`$.

As translational invariance is lost in the open system, the technique does not apply. Besides, the dynamics is profoundly modified by the presence of the two particle reservoirs, which impose the coexistence of two domains into the system. In case of the phases $`AI,BI`$, it is known that the two domains imposed by the reservoirs have a simple factorized structure. In this parameter regime the so-called domain wall (DW) theory can be applied . The DW theory uses a coarse grained description of the dynamics of the process: Each particle reservoir that is coupled to the chain imposes independently a domain of a given constant density $`\rho _L`$ and $`\rho _R`$. The two domains are separated by a localized domain wall, which performs a biased random walk. The bias is due to the different capacity of the two reservoirs and can be calculated simply by using the conservation of mass. In a finite system the position of the domain wall is confined between two reflecting walls. The description of the process allows to calculate the stationary as well as the fully time dependent probability distribution of the domain wall positions . Then, it is in particular possible to estimate the largest relaxation times of the system, that are given by

$$\tau _n=\left[D^++D^{}2\sqrt{D^+D^{}}\mathrm{cos}(\pi n/(L+1))\right]^1n=1,2,\mathrm{}$$
(2)
with

$$D^+=\frac{\beta (1\beta )}{1\alpha \beta };D^{}=\frac{\alpha (1\alpha )}{1\alpha \beta }.$$
(3)
These results are valid for $`\alpha ,\beta <0.5`$. The remaining parameter space has to be explored numerically.

### 2.2 The DMRG method for stochastic models

In the previous section we have shown that the master equation can be rewritten as a Schrödinger equation in imaginary time. Now, in order to calculate the longest relaxation time of the system, one has to calculate the two lowest eigenvalues of eq. (1) (the first eigenvalue is trivially zero, with an eigenvector corresponding to the stationary state). The eigenvalue calculations have to be done with very efficient diagonalization methods, in order to reach sufficiently large system sizes. This is possible by applying the DMRG method .

This method was first developed to study the properties of strongly correlated electrons . It has recently been generalized in order to treat stochastic many particle systems .

The idea of the DMRG is the following. One starts with a small system (here 12 sites) that can be diagonalized with standard numerical methods. Then one performs a large number of renormalization cycles in order to increase the system size. At each renormalization cycle, first the system is enlarged by adding 2 sites in the bulk of the system, i.e. in a place where it is less likely that the eigenmodes will be perturbed. The largest eigenvalues of the corresponding enlarged Hamiltonian are computed. These eigenvectors are used in order to construct the density matrix of the system, which will be used in the next stage. Second, in order to avoid an exponential growth of the hamiltonian, a projection onto the “most important” modes has to be done. It has been shown that the best choice is to keep the $`m`$ leading eigenvectors of the density matrix . The parameter $`m`$ has to be chosen small enough in order to allow a fast calculation of the eigenvectors, but large enough in order to obtain a high numerical precision. The error due to each truncation is well controlled .

Although no fundamental differences between the DMRG method for stochastic and quantum mechanical systems exist, one has to make a certain effort in order to overcome the numerical difficulties. This has been done following the suggestions of ref. . We now discuss some details of our implementation.

The main difference between the quantum mechanical and stochastic problem is that the stochastic hamiltonian is not hermitian. Therefore, the calculation of the eigenvectors of the hamiltonian required at each renormalization step cannot be done by standard diagonalization techniques (as the basic Lánczos or Davidsson algorithms) but one has to apply, e.g. the Arnoldi method or Lánczos for non-symmetric matrices. These methods are quite efficient but, compared to their analogues for symmetric matrices, less stable . We use the Arnoldi method, that turned out to be numerically more stable than the Lánczos method for non-symmetric matrices. Nevertheless the method is numerically well controlled, because the accuracy of the calculated eigenvector can be obtained via the residual norm, i.e. we check that it is indeed an eigenvector. Actually, we have used the Arnoldi method in an iterative way, where the initial guess is the result of the previous Arnoldi run, until the desired precision is reached.

Apart from being non-hermitian, there is another difference between stochastic and quantum mechanical problems. It is well known that the lowest eigenvalue of stochastic hamiltonians is always zero, corresponding to the existence of a stationary state, and that the corresponding left eigenvector $`|0_l`$ has all its coordinates equal to 1. Therefore the first task will be to calculate the right eigenvector $`|0_r`$. As we are interested in the asymptotic dynamics of the system, we have to calculate the energy gap G(L) as well, i.e. the next smallest eigenvalue and the corresponding left and right eigenvectors $`|1_l`$ and $`|1_r`$. Then the longest relaxation time $`\tau =\tau _1`$ will be given by

$$\tau =G(L)^1.$$
(4)
After the calculation of these eigenmodes, one has to eliminate the least important degrees of freedom of the system, in order to keep the system at a manageable size. It has been shown that the relevant degrees of freedom are the eigenmodes associated to the largest eigenvalues of the density matrix $`\widehat{\rho }`$. The latter is defined depending on how many eigenmodes we need to calculate. Here, for the calculation of the gap, we take

$$\widehat{\rho }=\frac{1}{2}\widehat{Tr}\left\{|0_l0_l|+|0_r0_r|+|1_l1_l|+|1_r1_r|\right\}.$$
(5)
Note that it is always possible to construct a symmetric form for the density matrix. Therefore, standard algorithms can be applied in order to calculate all eigenvalues and eigenvectors of $`\widehat{\rho }`$ as for quantum systems.

Another trick that allowed to gain accuracy on the gap calculation is to compute first the right eigenmode $`|0_r`$ for the fundamental state (the left is known), and then to define a new hamiltonian

$$H^{}(\mathrm{\Delta })=H+\mathrm{\Delta }|0_r0_l|.$$
(6)
This hamiltonian $`H^{}`$ has, apart from the groundstate, the same spectrum as the original hamiltonian $`H`$. The gap is now given by the fundamental state of $`H^{}`$ .

Note that one requirement of the DMRG method is that the requested eigenvalue is well separated from the next one. This may not be true when the system size increases, and then the calculation becomes instable. The system size limitations come from these instabilities rather than from runtime or memory requirements. So, the DMRG method either gives very accurate results for the eigenvalues, or does not converge at all.

Finally we want to mention another particularity of the system. As rules are different at each end of the system (input or output), we cannot use the left/right symmetry of the system as people do for closed systems, and the left and right parts of the system have to be computed alternatively.

## 3 Test of the Domain Wall Theory at Small System Sizes

The comparison between DW predictions and the exact analytical for the stationary state shows that they agree in the limit of large system sizes and $`\alpha ,\beta <1/2`$ . For finite systems, however, the coarse grained picture slightly deviates from the exact results.

In case of the dynamical properties, so far no exact analytical results exist. Therefore the DW predictions have to be checked with numerical methods.

In a previous work, we have compared time-dependent density profiles in a non-stationary regime, and found good agreement with the domain wall predictions . Here, we calculate explicitly the relaxation times of the process, in order to compare with the predictions of the DW method. The same has been done already by using the non-symmetric Arnoldi method , but the system sizes that one can treat by the non-symmetric Arnoldi method are too small ($`L16`$ in ) in order to obtain convergence with the DW predictions. Our aim is to improve this convergence thanks to the ability of DMRG to treat larger systems.

First we consider cases with both $`\alpha `$ and $`\beta <1/2`$, where the domain wall theory is expected to work well, though this was not fully verified by the exact numerical calculation, due to size limitations. Fig. 2 shows the comparison of the inverse relaxation time with the domain wall predictions in the high density phase ($`\alpha =0.2`$, $`\beta =0.1`$). For comparison we also included the result for $`L=16`$ obtained in , which perfectly coincides with the DMRG calculations. Obviously the theoretical predictions and the numerical results are in good agreement, if the length of the chain reaches thirty sites.

For $`\alpha =0.4`$ and $`\beta =0.2`$, the DW predictions for small system sizes are slightly less accurate, as shown in figure 3. This could be expected, as the density difference between the two domains $`\mathrm{\Delta }\rho =1\beta \alpha `$ is smaller and therefore the width of the shock larger than in the previous case. If, however, the shock touches the boundaries the approximation of microscopic states by $`\mathrm{\Theta }`$-functions is less appropriate.

For $`\alpha `$ or $`\beta `$ larger than 0.5, the coexistence is not between a high- and low-density domain, but, e.g. for $`\alpha >1/2`$, between a maximal current and high-density domain. This case is much more complicated than the coexistence of domains in the phase $`AI`$ or $`BI`$. It is known from the exact results, that long ranged correlations exist in the maximal current regime, which implies that the maximal current domain has no simple structure, and cannot be easily described.

A naive attempt to describe the maximal current domain, which uses the fact that the flux is a constant $`J=1/4`$ and would lead to exponentially decaying density profiles, is to build the density profile iteratively from the boundary with a mean field formula like $`j=\rho _i(1\rho _{i+1})`$. But this yields a wrong density profile. It was also suggested  to use the exact results of the density profile. However, this approach has intrinsic shortcomings. First the DW is no longer self-contained and, second, it is no longer possible to specify the relaxation times in a closed form, because the hopping rates of the random walk are now site-dependent. But above all, the description of the system states as the juxtaposition of two phases with an independent boundary layer between them breaks down. Indeed, the wall dynamics is coupled to the whole structure of the maximal current phase, due to the long ranged correlations.

So we explore this parameter regime guided by two questions: $`(i)`$ Do we obtain a finite relaxation time in the phases $`AII/BII`$ in the limit of large system sizes? $`(ii)`$ And, if it is the case, is it possible to model the maximum current phase as a flat domain with an effective density $`\rho _{eff}`$?

First we checked the domain wall predictions close to the transition line at $`\alpha =0.5`$, i.e. for $`\alpha =0.51`$ and $`0.55`$ (see figure 4). In this parameter regime one already has long ranged correlations in the maximal current domain, but at the same time the magnitude of the deviations from a flat domain is rather small. Our results indicate, that within this parameter regime the DW formulas (2-3) still lead to satisfying results for quite large system sizes. On the other hand, the DW prediction with $`\alpha =1/2`$, which was conjectured at first to be the proper one for any $`\alpha >1/2`$ , clearly underestimates the gap for any system size.

For larger values of $`\alpha `$ (e.g. $`\alpha =0.65`$ in figure 4), the algebraic corrections of the density profile become non-negligible. Equations (2-3) do not give anymore a good estimate.

Next, we checked whether our results are compatible with a finite relaxation time in the limit $`L\mathrm{}`$. Therefore we interpolated our results for $`G(L)`$ using the functional form $`f(L)=aL^b+c`$. We find a very good agreement between the fits and our numerical data, with a finite value of $`c`$ (see the inset of fig. 4), which indicates that the relaxation times remain finite. Our estimates for the relaxation times are $`\tau (\mathrm{})=116.(6)`$ for $`\alpha =0.51`$, $`\tau (\mathrm{})=101.(7)`$ for $`\alpha =0.55`$, and $`\tau (\mathrm{})=96.(1)`$ for $`\alpha =0.65`$.

A finite relaxation time for $`L\mathrm{}`$ implies at the same time, that the dynamical exponent vanishes. This can be checked systematically by calculating the size dependent dynamical exponent given as:

$$z(L)=\frac{\mathrm{ln}\left[G(L+2)\right]\mathrm{ln}\left[G(L)\right]}{\mathrm{ln}(L+2)\mathrm{ln}(L)}.$$
(7)
Figure 5 shows that our results are compatible with $`z=0`$ for all values of $`\alpha `$ we took into account. So these results for $`\alpha >0.5`$ further support that the relaxation times are finite within these subphases.

Now, one can simply take the values $`\tau (\mathrm{})`$ in order to calculate an effective density of the maximal current domain. This procedure leads to DW predictions for small systems, which do not agree very well with our numerical findings. We thus believe that the non-trivial structure of the maximal current domain must be considered in modeling the domain wall motion.

## 4 Disorder line

For general $`\alpha ,\beta `$ it has been shown that the stationary weights are given as products of infinite dimensional matrices. Nevertheless a one-dimensional line in the parameter space exist, i.e. $`\alpha +\beta =1`$, where the stationary weights are simple products as for the periodic systems. These product states are in analogy with disordered states of quantum spin models. At the disordered line the system is homogeneous, i.e. one can not identify two different domains. This implies that the relaxation at this particular line is not governed by the motion of the domain wall.

Then it is an open question to know whether dynamical properties are changed on this line and how the dynamical properties compare to the periodic system. For the periodic system, Bethe ansatz predicts the dynamical exponent $`z=3/2`$ . Besides, the first non-zero eigenvalue has an imaginary part except for the density $`1/2`$. The divergent relaxation times of the periodic system are, however, related to its translational invariance .

By contrast, for the open system it is expected that the relaxation times are finite if $`\alpha 0.5`$, because the density fluctuations spread with a nonzero drift velocity .

First, we have applied the DMRG method to the special point $`\alpha =\beta =0.5`$, where the drift velocity of an excess density vanishes.

Fig. 6 shows that our numerical results for $`z(L)`$, obtained up to $`L=244`$, agree with the exponent $`z=3/2`$ predicted for the periodic system.

We have also explored the remaining part of the disorder line. We present our simulation results in fig. 7. They seem to exclude the case $`z=3/2`$, and indicate that the relaxation time would be finite in the large $`L`$ limit ($`z(\mathrm{})=0`$).

Although these results are consistent with the expected relaxation behavior, we cannot exclude that a branch corresponding to complex eigenvalues could cross our solution when $`L`$ becomes large. The DMRG method does not allow to distinguish whether the numerical instabilities are due to the intersection or the convergence of two eigenvalue branches.

Nevertheless, we state that our numerical results are in accordance with a finite relaxation time and, therefore, consistent with the physical picture which was presented in . Further support for this scenario comes from its position in the phase diagram. The disordered line touches a phase boundary only in a single point, i.e. for $`\alpha =\beta =1/2`$. All other points of the disordered line are located inside a phase where the relaxation times of the system are finite. Therefore, it would be rather counter-intuitive to observe a diverging relaxation time for these parameters.

## 5 Phase transitions and the maximal current phase

If the capacity of both particle reservoirs exceeds the capacity of the chain, the maximal current phase is realized. In the maximal current phase both localization lengths are infinite. This divergence of length scales leads, e.g. to an algebraic slope of the density profiles.

In this phase, previous numerical calculations by Bilstein et al  for $`L20`$ extrapolated to large systems seemed to indicate an exponent $`3/2`$ in the whole maximal current phase.

We have computed the dynamical exponent in various points in the maximal current phase. Results are presented in figure 8, for a constant $`\beta =0.8`$, $`\alpha `$ being varied through the second order transition ($`\alpha =0.5`$) and inside the maximal current phase (we have checked that of course, results are the same when $`\alpha `$ and $`\beta `$ are interchanged). We clearly see a transition when the maximal current phase is entered. On the transition line $`\alpha =0.5`$, we recover $`z=3/2`$ in the large system size limit. Besides, from our numerical results, the infinite size dynamical exponent seems to be equal to $`3/2`$ in the whole maximal current phase, confirming the extrapolation in .

## 6 Conclusion

The use of a DMRG approach to compute the largest relaxation time in the TASEP model has allowed to give new strong evidence of the validity of the domain wall picture when reservoirs control the flow.

For $`\alpha `$ or $`\beta `$ larger than 1/2, the capacity of one of the reservoirs becomes larger than the chain capacity, and one of the coexisting domains is then a maximal current phase. Our numerical results indicate that the relaxation times $`\tau (\mathrm{})`$ are finite for $`L\mathrm{}`$ as in the phases $`AI/BI`$. As the bulk density in a maximal current phase is $`1/2`$ independently of the boundaries, one could have expected that the relaxation time $`\tau (\mathrm{})`$ would be independent of the most efficient reservoir (i.e. independent of $`\beta `$ in AII and $`\alpha `$ in BII). This is not the case, as boundary layers in the maximal current domain decay only algebraically, and lead to a site-dependence of the density. Besides, strong correlations that cannot be easily characterized exist in this phase. This means that the maximal current phase cannot be modeled as a flat domain with density $`1/2`$.

But even if we model the maximal current domain as a domain of a constant effective density $`\rho _{eff}`$, which leads to the estimated value of $`\tau (\mathrm{})`$, we observe large deviations between DW predictions and relaxation times for finite system sizes - again due to the non trivial structure of the maximal current domain. The DW predictions would probably be much improved if one considers site dependent hopping rates, but they are difficult to obtain since the structure of an isolated maximal current domain is not known.

Anyhow, as in phase AI and BI, the domain wall theory still offers a simple explanation for a finite relaxation time of the system in phases AII and BII.

On the disorder line, this picture cannot be applied, as no domain wall can be identified anymore, and a qualitatively different dynamical behavior cannot be excluded. However, though a wall cannot be defined anymore, we could still trace the dynamics of the density fluctuations by introducing a second class particle. The second class particle performs a biased motion to one of the boundaries for any density different from $`1/2`$. Therefore one expects that the density fluctuations are also driven out of the system with a finite velocity, which should lead to finite relaxation times, if $`\alpha 1/2`$ .

Contrary, in the periodic system, which has the same simple structure of stationary state, the relaxation times diverge as $`LL^{3/2}`$ for arbitrary densities. The divergent relaxation times have been related to the translational invariance of the periodic system  and are, therefore, not expected in case of the open system, except if $`\alpha =\beta =1/2`$.

Indeed, our numerical results support this scenario. Besides, we have shown that the dynamical exponent $`z=3/2`$ is also recovered in the whole maximal current phase, for which no theoretical prediction exists.

So we could summarize the expected behavior of the relaxation time as follows :

(i) In the whole phases AI, AII, BI, and BII, except on the transition line $`\alpha =\beta `$, our results are consistent with a converging relaxation time. Its limit value as $`L\mathrm{}`$ is known from the DW theory in AI and BI.

(ii) On the line $`\alpha =\beta `$, and in the whole maximal current phase, the relaxation time diverges when $`L\mathrm{}`$. The associated dynamical exponent $`z`$ is $`2`$ on the line $`\alpha =\beta `$ and $`3/2`$ in the maximal current phase and at the point $`\alpha =\beta =1/2`$.

Further improvement could be obtained from the use of Finite Size algorithm for DMRG (FSM), as described by Carlon et al , in order to gain some precision on our calculations - and thus to reach larger system sizes - in the whole phase diagram, and especially near critical lines.

It would also be interesting to apply the DMRG method to other models, which have a non-trivial but not necessarily site dependent domain structure, in order to check whether the DW theory still describes correctly the relaxation behavior of the system.

Acknowledgments: We acknowledge fruitful discussions with B. Derrida, J. Krug and E. Carlon, and precious remarks from G. Schütz. L.S. acknowledges kind hospitality at the L.P.S. and financial support by C.N.R.S.


text:
# Comment on “Multi-phonon Raman scattering in semiconductor nanocrystals: Importance of non-adiabatic transitions”

## Abstract

In a recent paper \[E. P. Pokatilov et al, Phys. Rev. B 65, 075316 (2002)\] the Raman selection rules in spherical nanocrystals with a degenerate valence band are analyzed. Some precisions are given here.

The authors state that the one-phonon Raman spectrum of spherical CdSe nanocrystals in crossed polarization configuration is determined only by d-phonons with angular momentum projection $`m=\pm 1`$. This is not a general result, as the number $`m`$ is rather conventional. First, to define the angular momentum projection, we need to set a quantization axis, respect to which the phonons are defined. Let the quantization axis be the $`Z`$-axis of our coordinate system. If we choose another quantization axis, let say $`Z^{}`$, we shall find new $`m^{}=\pm 1`$ phonons, which are linear combinations of $`Z`$-quantized phonons with all values of $`m`$.

Now, consider a Raman experiment with incident an scattered light linearly and crossly polarized. We may choose the $`Z`$\- and $`Y`$-axis of our reference frame to coincide with the incident and scattered polarization vectors, respectively. According to the theory of angular momentum and the dipole approximation, the absorption of an incident $`Z`$-polarized photon creates a $`P`$-type electronic excitation with angular momentum projection $`M_1=0`$. The nature of the electronic excitation is not important for this analysis. Furthermore, the emission of a $`Y`$-polarized photon is possible only by annihilation of a $`P`$-type electronic excitation with $`M_2=\pm 1`$. The transition between these electronic excitations in the case of CdSe nanocrystals is carried out through the creation or annihilation of phonons with $`m=\pm 1`$, as stated by the authors. Now, let rotate the reference frame, so that the new $`X^{}`$\- and $`Y^{}`$-axis of coincide with the incident and scattered polarization vectors, respectively. Then, the electronic excitations created when the incident photon is absorbed can have $`M_1^{}=\pm 1`$. The transitions to $`M_2^{}=\pm 1`$ electronic states can be mediated only by phonons with $`m^{}=0,\pm 2`$ ($`M_2^{}M_1^{}=0,\pm 2`$).

Both descriptions describe the same reality: the calculated Raman cross section, adding up the contributions of all the phonons, has the same value in both cases. The distinction between phonons with different angular projections is purely conventional if the system has spherical symmetry.

Due to the same reasons, in parallel polarization configuration, the Raman spectrum has contributions from phonons with $`m0`$ if the quantization axis does not coincide with the light polarization vectors.


text:
# Absence of correlation between built-in electric dipole moment and quantum Stark effect in InAs/GaAs self-assembled quantum dots

## Abstract

We report significant deviations from the usual quadratic dependence of the ground state interband transition energy on applied electric fields in InAs/GaAs self-assembled quantum dots. In particular, we show that conventional second-order perturbation theory fails to correctly describe the Stark shift for electric field below $`F=10`$ kV/cm in high dots. Eight-band $`𝐤𝐩`$ calculations demonstrate this effect is predominantly due to the three-dimensional strain field distribution which for various dot shapes and stoichiometric compositions drastically affects the hole ground state. Our conclusions are supported by two independent experiments.
PACS numbers: 78.67.Hc, 73.21.La, 31.15.Md

preprint:  For Physical Review Letters

Self-assembled InAs/GaAs quantum dots (SADs) are three-dimensional (3D) semiconductor nanostructures in which electrons and holes are completely confined along the three dimensions of space by the band gap difference between InAs and GaAs materials bim . In these nanoscale systems, the determination of the electronic spectra of both particles represents a major challenge because of the low symmetry of the 3D confined nanostructures that take (truncated) pyramidal or lens shapes, and are affected by the strong influence of strain due to the lattice mismatch between the InAs and GaAs crystals keg . Experimentally, the difficulty of ascertaining the exact dot shape, and the non-uniformity in size distribution resulting from the growth process are also a handicap ebi . In this context, the knowledge of the respective positions of electrons and holes can provide information on the confining potential experienced by both particles, which is of primary importance for fundamental gam  as well as practical reasons dep .

It has been recently argued that because of the non-homogenous stoichiometric composition of SADs, resulting from atomic InAs-GaAs inter-diffusion, the center of mass of the ground state electrons and holes are displaced from one another with the hole above the electron, thereby inducing a built-in electric dipole oriented from the top to the base of the dot. This inverted electron-hole alignment fry  has been derived experimentally by assuming the usual linear relation between the electric dipole moment and the Stark shift in Stark effect spectroscopy on SADs in p-i-n structures, i.e.,

$$E(F)=E_0+pF+\beta F^2$$
(1)
where $`E`$ is the electron-hole ground state transition energy, $`p`$ is the built-in dipole moment, and $`\beta `$ measures the polarization of the electron and hole states.

By estimating the electric field at which the maximum transition energy occurs, one can determine whether the structure has a positive (if that happens at positive fields) or negative (if otherwise) built-in dipole moment. This quadratic dependence is well-known for quantum-well systems and has been confirmed in nanocrystallite quantum dots emp . Recent theoretical works have made use of this linear relation to support the universality of the inverted alignment in SADs bak .

In this letter, we demonstrate that the linear relation between electric dipole moment and Stark shift as expressed in Eq. 1 is usually violated in SADs because of the unique strain field distribution in the dot, which strongly influences the rapidly varying confining potential for the holes. We show that this peculiar effect invalidates the conventional perturbation approach to establish the electric field dependence of the optical transition in SADs.

The system of SAD structures, containing both electrons and holes, can be well described by an eight-band strain-dependent $`𝐤𝐩`$ Hamiltonian bah , which reads,

$$(𝐇_{kp}^0+𝐇_{kp}^s+|e|Fz)\psi =E\psi .$$
(2)
Here $`𝐇_{kp}^0`$ and $`𝐇_{kp}^s`$ is the kinetic and strain-dependent part of the eight-band Hamiltonian pry , respectively, and $`\psi =(\psi _1,\psi _2,\mathrm{},\psi _8)`$ is the envelop eigenvector. By using continuum elasticity theory, the strain tensor is calculated on a large grid of $`150\times 150\times 150`$ sites. The Hamiltonian is then solved by Lanczos algorithm. The same technique has been successfully applied to study few-particle effect lan  and optical transitions swd2  in quantum dots and very good agreement with experiments has been achieved.

Figure 1 shows the calculated dependence of interband transition energies on electric fields for two SAD structures. Fig. 1(a) illustrates the data for a $`\text{In}_{0.8}\text{Ga}_{0.2}\text{As}`$ dot with a constant composition throughout the structure. The pyramidal dot is $`18`$ nm wide and $`7.8`$ nm high, which is of similar size as that in the experiment of Fry et al. fry . The electric field is applied vertically to the structure, pointing from the base to the top. For comparison, the corresponding experimental data is shown in the top-right inset.

An clear deviation from the quadratic dependence of the transition energy on the electric field is obtained from our numerical simulations that show two distinct branches, one for roughly each field direction, merging at around zero field. It is interesting to notice that, if the data at negative electric fields are shifted upwards by $`11`$ meV, and if one ignores a few data points around zero-field, one can easily fit the remaining data with a single parabola characterized by a maximum at around $`F=40`$ kV/cm. In this case, according to Eq. 1, the displaced curve yields a negative built-in dipole ($`14`$ Å) in spite of the fact that the center of mass of the electron state is above that of the hole state by $`7`$ Å as shown on the lower inset, thereby achieving normal alignment.

This fitting procedure was adopted by Fry et al. fry  to show the overall quadratic field dependence of the experimental transition energies (see upper right inset). The points marked with solid dots are the data measured from experiments. We notice however the structure may be of different stoichiometric composition than our SAD in Fig. 1(a). The authors shifted their data to the position marked with open dots to match the curve fitting from their data at positive electric fields. With this procedure, they claimed negative built-in dipole in the quantum dot structure. This experimental procedure was justified by the fact that for p-i-n structures, the Stark effect can only be measured under reverse biases. Therefore, given the asymmetric SAD shape and the fact that the Stark effect was measured for both field directions, two sets of samples i.e p-i-n and n-i-p structures, were used in the experiment.

While the conventional perturbation theory fails for the pyramidal quantum dot, it makes a successful prediction of Stark shifts for another quantum dot which is a largely truncated pyramidal InAs quantum dot shown in Fig. 1(b), which has the same base size, but is only $`5.4`$ nm high. It has smaller transition energies than the pyramidal structure because of its uniform InAs composition. Unlike the dot shown in Fig. 1(a), it has nearly perfect quadratic dependence at electric fields as strong as $`F=\pm 100`$kV/cm. The parameter extracted from the fitting curve gives a positive built-in dipole of $`4.57`$ Å which agrees well with the actual value $`4.8`$ Å.

In Fig. 2, we show the electron and hole ground state energies for the two structures as a function of electric fields. Electron energies for both structures are seen to have a linear profile with a slight bowing that can be well described by the usual quadratic dependence on the electric field, i.e.,

$$E_e(F)=E_e(0)+p_eF+\beta _eF^2,$$
(3)
in the second-order perturbation theory. The hole state in the truncated structure exhibits a similar feature and its energies can be well fitted with a parabola, i.e.,

$$E_h(F)=E_h(0)+p_hF+\beta _hF^2.$$
(4)
Here $`p_{e(h)}/|e|=\psi _{e(h)}|z|\psi _{e(h)}`$ is the center of mass of the electron (hole) state $`\psi _e`$ ($`\psi _h`$) and $`\beta _{e(h)}`$ measures the polarization of the electron(hole) states. The correlation between the built-in dipole moment $`p=p_ep_h`$ and quantum Stark effect is then obtained by taking the difference between the electron and hole energies to reflect Eq. 1. In the pyramidal dot, however, a quadratic fitting for the whole field range is not possible, although the left and right branches of the spectrum could be fitted piecewise by two parabolas, respectively.

The different behavior of electron and hole energies is attributed to their respective band edge profiles shown in the insets of Fig. 2. In both SADs, the conduction band edge is almost flat inside the dot with hard walls on the sides, which provides electrons a constant, strong confinement. Consequently, under the influence of external electric fields, electrons experience a smoothly varying confining potential inside the dot, which can be treated as a small perturbation with electron energies following the quadratic dependence given by Eq. 3. This type of behavior is also observed for holes in the truncated SAD structure (Eq. 4).

In the pyramidal dot, the valence band edge exhibits a more complicated, double-triangular profile with the lower part extending from the SAD base to about two-third of the structure hight, which effectively confines the ground hole state close to the bottom of the dot (see the inset in Fig. 1). Therefore, the holes are localized in a rapidly varying potential while the valence band edge changes abruptly close to the top of the structure. Hence, as vertical electric fields force the holes to move along the structure, they experience significantly different local confinement at different fields, which invalidates the perturbation approach as described in Eq. 1. This behavior is seen to occur at electric fields as low as $`F=\pm 10`$ kV/cm.

The rapidly varying confining potential which is responsible for the breakdown of the perturbation theory in the pyramidal structure is induced by the unique 3D strain field distribution within the SAD. Because the effect of other strain components on the valence bands are inferior in importance, we concentrate on the biaxial component of the strain which affects mostly the valence bands bah2 . Fig. 3 shows a 3D plot of the magnitude of biaxial strain along a  plane across the center of the quantum dots. It is pointed out that the biaxial component of strain takes negative value inside the dot that is under biaxial compression and positive value in the barrier that is under biaxial tension.

Despite the fact that the structural difference between the two dots results only from the truncation, both structures exhibit substantially different strain distribution. The biaxial strain distribution in the pyramidal dot shows rapid variation inside the structure along the growth direction while the truncated dot exhibits a smoother profile. The strain field is so coherent that the tip of the pyramidal dot has a significant effect on the whole strain distribution.

In Fig. 4, we show the ground hole states at both negative and positive electric fields for the two structures. For the pyramidal structure of Fig. 1(a), the wave function at negative field $`F=20`$ kV/cm looks similar to that at zeros field except that it is localized closer to the dot bottom. At a positive field $`F=10`$ kV/cm, the hole wave function develops a pair of small ‘wings’ around the facet edges pry2 . Because of a strong local confinement imposed by the sharp triangular potential (see the band diagram inset in Fig. 2), the hole state could not extend into the upper half part of the structure as electron states do. At a stronger positive field $`F=20`$ kV/cm, these wings even dominate over the other portions of the wave function.

The obvious asymmetric behavior of the ground hole states implies that the corresponding wave functions cannot be expanded to the first order in the electric field as

$$\psi _h(F)=\psi _h(0)+\varphi _hF,$$
(5)
which would result into symmetric wave functions with respect to positive and negative electric fields for the whole range of electric fields. Instead, if we chose different values of the perturbation, $`\varphi _h`$, for negative and positive fields, respectively, it is possible to fit the hole energies with two parabolas (see Fig. 2).

In the truncated structure, the hole states are seen to have a different probability density distribution. The wave functions change very smoothly with the electric fields, because of a relatively flat valence band edge (see Fig. 2 inset). In addition, the lateral confinement close to the SAD top is also weaker compared with that in the pyramidal SAD as shown Fig. 1. Therefore, hole states behave coherently throughout the whole range of electric fields, exhibiting nearly perfect quadratic dependence of their eigen energies on electric fields.

Table I lists the actual and fitted values of dipole moment for several SADs. All structures have the same base dimensions, but different heights or composition profiles. It should be stressed that the fitted values are obtained from the model data only for positive electric fields, i.e., the right branch of the spectrum. Except for the last structure, all values exhibit large discrepancy between the fitted curve and the model (calculated) data at negative fields. The first structure has a negative built-in dipole induced by the inhomogeneous diffusion swd , responsible for two bumps in the transition energy for positive and negative fields, which makes the fitting over a broad field range, meaningless. This structure is similar to the SADs investigated recently by Chen et al. mad  where a negative dipole moment was extracted (the orientation of their field was opposite to ours) from a quadratic fitting over a smaller range of field ($`7.5`$ kV/cm $`F`$ $`10`$ kV/cm), in good agreement with our data. The next three structures have homogeneous composition, and all exhibit positive dipoles. However, the fitted values for these structures are all negative. Therefore, it is a general feature of high pyramidal quantum dots with less than $`25\%`$ truncation to exhibit erroneously a ’negative’ built-in dipole when fitted to the Stark shifts. The fourth structure has about $`27\%`$ truncation with a dipole larger than those in the two previous dots. The corresponding fitted value is positive but is $`30\%`$ smalled than the actual value. In the last structure, which is strongly truncated, the difference between the actual value and the fitted one becomes negligible.

In conclusion, we have shown that the unique 3D strain field distribution is responsible for the breakdown of the linear dependence of the built-in electric dipole moment in the Stark shift in SADs. This effect is mostly noticeable in high SADs where holes are mostly affected by rapidly varying components of the strain induced confining potential, which invalidates the second order perturbation theory. While perturbation theory is expected to fail at strong electric fields, we point out that significant deviations from the usual Stark effect already occur at fields lower than $`F=10`$ kV/cm. For largely truncated or flat SADs, the Stark shift retains its usual quadratic dependence on external fields.

###### Acknowledgements.

This work is supported by Army Research Office, grant #DAAD 10-99-10129 and National Computational Science Alliance, grant #ECS000002N.