format: text

text:
# Long-time behavior of Granular Gases with impact-velocity dependent coefficient of restitution

## 1 Introduction

Gases of dissipatively colliding particles (Granular Gases) in the absence of external forces reveal a variety of interesting phenomena, such as characteristic deviations from the Maxwell distribution , overpopulation of the high energy tail of the distribution function , anomalous diffusion , and others (see  for an overview). The most striking phenomenon which makes Granular Gases different from molecular gases is the self-organized formation of spatio-temporal structures such as clusters  and vortices .

The loss of mechanical energy of dissipatively colliding particles $`i`$ and $`j`$ is characterized by the (normal) coefficient of restitution, which relates the normal component of the relative velocity before a collision, $`g`$, to that after the collision, $`g^{}`$:

$$\epsilon \frac{g^{}}{g}=\frac{\stackrel{}{v}_{ij}^{}\stackrel{}{e}_{ij}}{\stackrel{}{v}_{ij}\stackrel{}{e}_{ij}},\stackrel{}{e}_{ij}\frac{\stackrel{}{r}_i\stackrel{}{r}_j}{\left|\stackrel{}{r}_i\stackrel{}{r}_j\right|}$$
(1)
with $`\stackrel{}{v}_{ij}\stackrel{}{v}_i\stackrel{}{v}_j`$ and $`\stackrel{}{v}_{ij}^{}`$ being the corresponding after-collision value. Frequently it is assumed that the coefficient of restitution is a material constant which simplifies analytical calculations considerably. This assumption, however, contradicts experiments  and disagrees with a dimension analysis . Instead the coefficient of restitution is a function of the impact velocity.

Obviously, the coefficient of restitution can be obtained by integrating Newton’s equation of motion for the collision with appropriate choice of the contact force. The elastic component of the contact force of spheres of diameter $`\sigma `$ is given by Hertz’ law $`F^{(\mathrm{el})}=B\xi ^{3/2}`$, $`\xi (t)\sigma \left|\stackrel{}{r}_i\stackrel{}{r}_j\right|`$ with $`B(\sigma )`$ being the elastic material parameter . Assuming viscoelastic material properties the dissipative part of the contact force reads $`F^{(\mathrm{dis})}=A\sqrt{\xi }\dot{\xi }`$  with the dissipative material constant $`A(\sigma )`$. Integration of Newton’s equation of motion yields then

$$\epsilon _v(g)=1D_1g^{1/5}+D_2g^{2/5}\mathrm{},$$
(2)
where the coefficients $`D_{1/2}`$ depend on the elastic and dissipative particle material properties and on the radii of the particles (for details see ).

The temperature of Granular Gases in the absence of spatial correlations is defined as the mean kinetic energy of the particles, just as for molecular gases. Due to dissipative collisions the temperature decays persistently and, therefore, the thermal velocity $`v_T(t)\sqrt{2T(t)/m}`$ decays too. For vanishing temperature the particles collide elastically, i.e., $`lim_{T0}\epsilon =1`$. Since the formation of clusters in Granular Gases is a consequence of the dissipative particle collisions the question arises whether for gases of viscoelastic particles clusters may occur and persist. The linear stability analysis of the hydrodynamic equations of such gases supports the speculation that the cluster state is only a transient phenomenon for the case of viscoelastic particles .

For principal reasons (to be discussed below) the numerical investigation of the long time behavior of force free Granular Gases of viscoelastic particles is problematic. Therefore, in this paper we assume a coefficient of restitution with an extremely simplified impact velocity dependence: The particles collide with a constant coefficient of restitution $`\epsilon ^{}`$ if the impact velocity exceeds a certain value $`g^{}`$, otherwise they collide elastically:

$$\epsilon (g)=\{\begin{array}{cc}\epsilon ^{}\hfill & \text{for }g>g^{}\hfill \\ \text{1}\hfill & \text{for }gg^{}\text{ .}\hfill \end{array}$$
(3)
We wish to mention that the coefficient of restitution of the form (3) was used in  as a numerical trick to avoid the inelastic collapse in event driven MD simulations. In this context the transition velocity was chosen to be much smaller than the thermal velocity, $`gv_T`$. Before investigating the cluster formation process first we will derive the law of temperature decay in the homogeneous cooling state.

## 2 Decay of temperature in the homogeneous cooling state

In the absence of external forces a Granular Gas which is initialized at uniform number density $`n`$ stays homogeneous during the first stage of its evolution, called homogeneous cooling state. In this regime the kinetic energy of the gas may be characterized by the granular temperature $`T`$, defined by the second moment of the velocity distribution function $`f(v)`$

$$\frac{d}{2}nT\frac{mv^2}{2}f(v)𝑑\stackrel{}{v},$$
(4)
($`m`$ is the mass of particles and $`d`$ is the dimension) which decays persistently due to inelastic collisions. For $`\epsilon =\text{const.}`$ the decay is described by Haff’s law $`T_H(t)`$ , whereas for gases of viscoelastic particles the temperature is given by $`T_v(t)`$

$$T_H(t)=\frac{T_0}{\left(1+t/\tau _H\right)^2},T_v(t)=\frac{T_0}{\left(1+t/\tau _v\right)^{5/3}},$$
(5)
where $`\tau _H`$ and $`\tau _v`$ are the relaxation times .

Equations (5) have been first derived with the assumption of a Maxwell distribution for the particle velocities. It is known that the velocity distribution function of Granular Gases deviates slightly from a Maxwell distribution , however, the functional form of Eqs. (5) keeps conserved also when we take into account these deviations. Before investigating the temperature decay let us first check whether the assumption of an approximative Maxwell distribution is justified for our case too.

Starting at a certain temperature which corresponds to a thermal velocity $`v_Tg^{}`$ almost all collisions occur with a constant coefficient of restitution $`\epsilon ^{}`$, hence, in this range we expect to find in the homogeneous cooling state a velocity distribution close to a Maxwell distribution with the small deviations mentioned above. At late times when the thermal velocity is small with respect to the transition velocity, $`v_Tg^{}`$, the majority of the collision occurs elastically. In this range we expect, therefore, that the particle velocities obey a Maxwell distribution just as for molecular gases. In the intermediate range $`v_Tg^{}`$ a sizable part of the collisions, i.e. collisions between slow particles, occurs elastically, whereas fast particles mainly undergo dissipative collisions. Therefore, it is not à priori clear whether the distribution function is close to the Maxwell distribution.

We have simulated a two-dimensional force free cooling Granular Gas of $`N=10^5`$ particles of unit mass which collide according to the rule (3) by event-driven Molecular Dynamics. The parameters are: $`\epsilon ^{}=0.6`$, $`g^{}=0.1`$, $`\eta =0.1`$, and the initial temperature is $`T(0)=1`$ which corresponds to an initial thermal velocity $`v_T(0)=\sqrt{2}`$. Figure 1 shows the scaled velocity distribution function $`\stackrel{~}{f}(c)`$

$$f(\stackrel{}{v},t)=\frac{n}{v_T^d(t)}\stackrel{~}{f}\left(\stackrel{}{c}\right),\stackrel{}{c}\frac{\stackrel{}{v}}{v_T(t)},$$
(6)
for the three mentioned cases together with the Maxwell distribution. From the results we conclude that even in the transition region $`v_Tg^{}`$ the Maxwell distribution is a good approximation for the velocity distribution function.

Using the standard approach based on the Boltzmann equation (e.g. ) with the assumption of a Maxwell distribution the equation for the time evolution of temperature may be obtained

$$\frac{dT}{dt}=\frac{2}{d}g_2(\sigma )\sigma ^{d1}nv_TT\mu _2,$$
(7)
where $`\sigma `$ is the diameter of the particles, $`g_2(\sigma )`$ is the contact value of the pair distribution function ($`g_2(\sigma )=(17\eta /16)/(1\eta )^2`$ for a two-dimensional gas with packing fraction $`\eta =\pi \sigma ^2n/4`$), and $`\mu _2`$ is the second moment of the dimensionless collision integral, which may be written in the form

$$\mu _2=\frac{1}{4}𝑑\stackrel{}{c}_1𝑑\stackrel{}{c}_2𝑑\stackrel{}{e}\mathrm{\Theta }(\stackrel{}{c}_{12}\stackrel{}{e})\stackrel{~}{f}(\stackrel{}{c}_1)\stackrel{~}{f}(\stackrel{}{c}_2)(1\epsilon ^2)\left|\stackrel{}{c}_{12}\stackrel{}{e}\right|^3.$$
(8)
The integration is performed over the velocities $`\stackrel{}{c}_1`$, $`\stackrel{}{c}_2`$ of the colliding pair and over the unit vector $`\stackrel{}{e}=\stackrel{}{e}_{12}`$ as introduced in Eq. (1). The unit step function $`\mathrm{\Theta }(x)`$ guarantees that collisions occur only between approaching particles.

Approximating the velocity distribution function by the Maxwell distribution, $`\stackrel{~}{f}(\stackrel{}{c})=\pi ^{d/2}\mathrm{exp}(c^2)`$, the above moment may be calculated. For two-dimensional systems it reads

$$\mu _2=\frac{\sqrt{2\pi }}{2}\left(1\epsilon ^{\mathrm{\hspace{0.17em}2}}\right)\left(1+\frac{g_0^2}{2u}\right)\mathrm{exp}\left(\frac{g_0^2}{2u}\right),$$
(9)
where $`u=u(t)=T(t)/T_0`$ is the reduced temperature ($`T_0=T(0)`$) and $`g_0=g^{}/v_T(0)`$. Introducing the variable $`x=2u/g_0^2`$ we obtain from Eq. (7)

$$\dot{x}=\alpha \left(x^{3/2}+x^{1/2}\right)e^{1/x}$$
(10)
where

$$\alpha =g_0\left(1\epsilon ^{\mathrm{\hspace{0.17em}2}}\right)g_2(\sigma )\sigma n\sqrt{\frac{\pi T_0}{2m}}.$$
(11)
If $`v_T(0)g^{}`$, i.e. $`g_01`$, then $`x1`$ at the initial stage of the temperature evolution. Equation (10) reduces then to $`\dot{x}=\alpha x^{3/2}`$, yielding Haff’s law (5) with

$$\tau _H=\frac{\alpha }{g_0\sqrt{2}}=\frac{1}{2}\left(1\epsilon ^{\mathrm{\hspace{0.17em}2}}\right)g_2(\sigma )\sigma n\sqrt{\frac{\pi T_0}{m}}.$$
(12)
For the late stage of the temperature evolution, i.e. for sufficiently small $`x`$, the general solution of Eq. (10) can be found as a series

$$\frac{x^{3/2}e^{1/x}}{(1+x)}\left[1+\frac{1}{2}\frac{(3+x)}{(1+x)}x+\frac{1}{4}\frac{(15+10x+3x^2)}{(1+x)^2}x^2+\mathrm{}\right]=\alpha t+C,$$
(13)
where $`C`$ is an integration constant. For large times, i.e., $`x1`$ it simplifies to

$$x^{3/2}e^{1/x}=\alpha t$$
(14)
Applying the (implicit) definition of Lambert’s $`W`$ function

$$W(x)e^{W(x)}=x$$
(15)
we find the analytical solution of Eq. (14)

$$\frac{1}{x}=\frac{3}{2}W\left(\frac{2}{3(\alpha t)^{2/3}}\right);x=\frac{2}{g_0^2}\frac{T}{T_0}$$
(16)
For large values of $`t`$ we get the approximation (choosing the appropriate branch of $`W(x)`$)

$$\frac{1}{x}=\frac{3}{2}\left[\mathrm{ln}\left(\frac{3}{2}(\alpha t)^{2/3}\right)+\mathrm{ln}\mathrm{ln}\left(\frac{3}{2}(\alpha t)^{2/3}\right)+\mathrm{}\right]$$
(17)
which yields the asymptotics for $`t\mathrm{}`$

$$T(t)=\frac{g_0^2}{2}\frac{T_0}{\mathrm{log}\alpha t},$$
(18)
with the constants $`g_0`$ and $`\alpha `$ defined above.

At late times most of the particles collide elastically due to Eq. (3), therefore, only collisions of particles whose velocities belong to the high-energy tail contribute to the decay of temperature. For this reason the asymptotic temperature relaxation is extremely slow as compared with the power laws Eq. (5), see Fig. 2.

For Granular Gases of particles which collide with $`\epsilon =\text{const.}`$ it is known that the high energy tail of the scaled distribution function does not obey a Maxwell distribution ($`\mathrm{exp}(c^2)`$) but decays as $`\mathrm{exp}(\mathrm{const}.\times c)`$ . This overpopulation of the high energy tail arises due to a mechanism which is specific for gases of dissipative particles: the scaled velocity of a particle, $`c=v/v_T`$ can increase either by collisions with other particles or by keeping its velocity $`v`$ while the thermal velocity $`v_T`$, i.e. temperature, decreases. Both processes cause increase of the number of fast particles, i.e., they enhance the population of the high velocity tail. For very fast particles, $`c1`$, the former process may be neglected and the latter one leads to the exponential distribution function $`\mathrm{exp}(\mathrm{const}.\times c)`$ . We wish to stress that for $`\epsilon =\text{const.}`$ the Granular Gas has no inherent time scale and the shape of the scaled velocity distribution function $`\stackrel{~}{f}(c)`$ does not depend on time. Hence, the exponential overpopulation of the high-velocity tail is a time-independent effect.

For the coefficient of restitution according to Eq. (3) the overpopulation of the tail will not be observed at late times. The temperature decays only logarithmically slow as compared with the power law for $`\epsilon =\text{const}`$. Therefore, the process of depopulation of the tail due to collisions is not balanced anymore by the decrease of the temperature of the gas, which enhances the tail population. According to this argument in difference to the case $`\epsilon =\text{const.}`$ we expect an underpopulation of the high-energy tail. There is a non-negligible probability that particles collide in such a way that the after-collision velocity of one of the collision partners exceed $`g^{}`$ even if the pre-collision velocities of both partners are smaller than $`g^{}`$. The probability of such collisions, which are responsible for the population of the velocity range $`v^{}>g^{}`$, decays rapidly as a function of the after-collision velocity $`v^{}`$. To estimate the affected range of particle velocities we consider the collision of elastic particles: The maximal after-collision velocity which may emerge in a collision of elastic particles with pre-collision velocity $`g^{}`$ is $`\sqrt{2}g^{}`$. Hence we expect that only the interval of velocities $`(g^{},\sqrt{2}g^{})`$ will be significantly populated.

At present it remains unclear whether the deviation from the Maxwell distribution at high velocities does change the behavior of the gas at late times significantly. This intriguing problem is subject of current research.

## 3 Cluster formation

The spontaneous formation of clusters in a force free cooling Granular Gas can be understood by simple arguments : Consider fluctuations of the density in an otherwise homogeneous Granular Gas. In denser regions the particles collide more frequently than in more dilute regions, therefore, dense regions of the system cool faster than dilute regions and the local pressure decays in these colder regions. The resulting pressure gradient causes a flux of particles into the denser region, which leads to further increase of the density. Hence, small fluctuations of the density are enhanced which leads to the formation of clusters.

These arguments are certainly valid for the case $`\epsilon =\text{const.}`$ For viscoelastic material properties, however, with decaying thermal velocity the collisions occur less and less dissipatively since $`\epsilon _v`$ approaches one with propagating time. Hence the question arises whether cluster formation in Granular Gases of viscoelastic particles occurs as a transient phenomenon only. This question is difficult and there is no conclusive answer yet, although the linear stability analysis  seems to support the hypothesis that cluster formation is only a transient phenomenon in Granular Gases of viscoelastic particles.

The numerical simulation of clustering Granular Gases is problematic since clusters grow rapidly and at a certain time their extension is comparable with the system size. From this instant on the simulation becomes invalid since the cluster interacts with itself via the periodic boundary conditions. For the case of viscoelastic particle properties $`\epsilon =\epsilon (g)`$ we expect that clusters dissolve at late times when the typical coefficient of restitution $`\epsilon \left(v_T\right)`$ approaches one. To prove this by MD simulations we need a system which is large enough to assure that at no time the cluster size reaches the system dimensions. At the moment we are able to simulate systems up to about $`N=3\times 10^6`$ particles, which seems to be far below the necessary size.

We have simulated a system of $`n=10^5`$ particles which collide according to Eq. (3). This collision rule reveals the same asymptotic behavior as the rule (2) for the collision of viscoelastic spheres, $`lim_{g0}\epsilon =1`$, nevertheless, the law is drastically simplified. Snapshots of the simulation are shown in Fig. 3.

Starting with a uniform distribution after a short period of homogeneous cooling density inhomogeneities appear which grow up to a certain typical size. The collision frequency in dense regions is larger than in dilute regions which leads to decreasing thermal velocity. In the due of time more and more collisions inside the clusters occur with relative velocities $`g<g^{}`$, i.e., elastically. Therefore, the clusters dissolve. At late times the system returns to the homogeneous cooling state.

Why do we not observe system spanning clusters as in the case of $`\epsilon =\text{const.}`$? The reason is the dynamics of cluster growth. To grow a cluster of a given size the system needs a certain time which is determined by the density and by the value of the coefficient of restitution. We start our simulation at a temperature which corresponds to a thermal velocity which is well above $`g^{}`$, i.e., the gas has essentially the properties of a Granular Gas with $`\epsilon =\text{const.}`$ Therefore, clusters begin to form and grow due to the described pressure instability. At a certain time the typical velocity inside the clusters falls below $`g^{}`$, i.e., from this moment on the clusters start to dissolve. This time corresponds to a typical cluster size which is significantly smaller than the system size. Therefore, we do not observe clusters beyond a certain typical size.

According to a similar argumentation we expect clusters in Granular Gases of viscoelastic particles to dissolve too. Due to numerical limitations at present time it is not possible to perform the corresponding simulations without encountering the problem of boundary conditions.

In the recent study , where the restitution coefficient of the form (3) was used to avoid the inelastic collapse, the dissolution of clusters has not been detected. The reason for the persistent development of clusters in  was a very small value of $`g^{}`$, which was at least 100 times smaller than the thermal velocity of particles during the simulations. In our numerical experiments the clusters started to dissolve when the thermal velocity became comparable with $`g^{}`$. This regime has not been reached in simulations of .

## 4 Conclusion

We have investigated a model Granular Gas of particles which collide with a stepwise coefficient of restitution which mimics the coefficient of restitution of viscoelastic particles: The particles collide elastically ($`\epsilon =1`$) for small impact velocities and with $`\epsilon =\text{const.}<1`$ if the relative particle velocity exceeds a certain value ($`g>g^{}`$). This simplified dependence of the restitution coefficient reflects the main feature of viscoelastic particles that collisions tend to occur elastically for decreasing impact velocity. We analyze the temperature decay and the dynamics of cluster formation. We observe that clustering occurs in this model only as a transient process. The important property of the collision model is that the clusters dissolve before system spanning clusters appear. Hence, this model allows to investigate the total process without the drastic influence of the boundary conditions. The results of our study supports the hypothesis that clustering occurs as a transient process for realistic gases of viscoelastic particles, as was predicted before by means of the linear stability analysis.


text:
# Exchange interaction effects on the optical properties of LuMnO3

## Abstract

We have measured the optical conductivity of single crystal LuMnO3 from 10 to 45000 cm-1 at temperatures between 4 and 300 K. A symmetry allowed on-site Mn $`d`$-$`d`$ transition near 1.7 eV is observed to blue shift ($``$0.1 eV) in the antiferromagnetic state due to Mn–Mn superexchange interactions. Similar anomalies are observed in the temperature dependence of the TO phonon frequencies which arise from spin-phonon interaction. We find that the known anomaly in temperature dependence of the quasi-static dielectric constant $`ϵ_0`$ below the $`T_N90`$ K is overwhelmingly dominated by the phonon contributions.

The colossal magnetoresistance compounds based on doped pseudo-cubic LaMnO3 have excited much attention because of their interesting physical properties and potential applications. Another series of $`R`$MnO3 materials ($`R`$=Ho,Er,Tm,Yb,Lu, or Y,Sc,In) have smaller radius $`R^{3+}`$ ions and crystallize in the hexagonal lattice. The hexagonal manganites are interesting as examples of multiferroics (or ferroelectromagnets) Smolenskii and Chupis (1982 \[Sov. Phys. Usp. 25, 475 (1982) – they are both ferroelectric ($`T_c`$$``$900 K) and strongly frustrated antiferromagnets ($`T_N`$$``$90 K). The coupling between ferroelectric and magnetic order parameters provides the prospect of manipulating electrical properties through magnetic fields and vice versa which, in turn, gives these compounds potential for applications in electronics.

There have been several reports on the magnetic structure Bertaut and Mercier (1963); Munoz et al. (2000); Fiebig et al. (2000) and aspects of the electromagnetic response in hexa-manganites Fiebig et al. (2000); Kritayakirana et al. (1969); Penney et al. (1969); Iliev et al. (1997); Yi et al. (2000); Takahashi et al. (2002). Of particular interest is the observation of a temperature anomaly of the static dielectric constant, $`ϵ_0`$, below $`T_N`$. In this letter, we present a systematic study of the linear optical response of LuMnO3 which elucidates the origin of the $`ϵ_0`$ anomalies and the effects of the Mn–Mn exchange energy on the electrodynamics of this system. We show that the exchange interaction manifests itself in an antiferromagnetic resonance, spin-phonon coupling, and the temperature dependence of a 1.7 eV on-site Mn $`d`$-$`d`$ optical transition. The results give a comprehensive view of the magnetic and electronic structure of this interesting ferroelectric and strongly frustrated antiferromagnetic material.

The high-temperature paraelectric phase in the hexagonal manganites may be considered as layers of corner-sharing MnO5 triangular bipyramids connected by a layer of $`R`$-ions. In the ferroelectric phase, $`R`$-ions alternate their $`c`$-axis coordinates, producing a net electric moment of the unit cell, and the MnO5 pyramids tilt away from the $`c`$-axis. The $`d`$-orbitals of the Mn3+ ion in a triangular bipyramid of five O2- ions are split by the crystal field into three groups: $`d_{xz,yz},d_{xy,x^2y^2}`$, and $`d_{3z^2r^2}`$, in order of increasing energy. The latter orbital has the highest energy as the Mn-apical O bond lengths are shorter than in-plane distances. The in-plane $`d_{xy,x^2y^2}`$ orbitals are strongly hybridized with oxygen $`p`$ orbitals. In the ground state of the Mn3+ ion four electrons occupy four lowest orbitals giving $`<L>=0`$ and $`S=2`$. Below $`T_N`$, spins are ordered antiferromagnetically in-plane Munoz et al. (2000); Fiebig et al. (2000). The Neel temperature in hexa-manganites is 6–10 times smaller than the Weiss temperature which is due to the spin frustration in the triangular planar lattice and weak inter-plane exchange interaction.

Single crystals of LuMnO3 were grown using the travelling floating zone method and characterized by magnetization, resistivity, and x-ray powder diffraction. The lattice constants as well as the observed macroscopic properties agree well with measurements reported in literature. Platelet samples were cleaved perpendicular to the $`c`$-axis. In our experiments, we studied a 25 $`\mu `$m thick $`5\times 1.5`$ mm2 sample and a $`1.5\times 2\times 4`$ mm2 sample (for $`c`$-axis response). Transmittance and reflectance measurements were performed using a Fourier transform spectrometer in a frequency range from 10 to 45000 cm-1 (1.2 meV to 5.6 eV). Temperature dependence from 4 to 300 K is achieved using liquid He in a continuous flow cryostat (sample in vacuum) with optical access windows.

The optical spectrum of LuMnO3 is dominated by phonons in the 28–100 meV photon energy range and electronic transitions starting at $``$1.1 eV. The material is transparent below 28 meV and between 0.1–1.2 eV.

Figure 1 shows the electronic part of the optical conductivity spectrum of LuMnO3 calculated by Kramers-Kronig relations from the reflectance spectrum measured at 10 and 300 K. The lowest electronic excitation centered at $``$1.7 eV is seen to be strongly temperature dependent. A similar effect has been reported for hexa-YMnO3 Kritayakirana et al. (1969). Whereas, in ortho-YMnO3 this peak is not observed or is very weak Yi et al. (2000). From the Inset in Fig. 1 it is seen that the resonance energy decreases monotonically with temperature having an inflection point at $`T_N`$ 90 K while the integrated spectral intensity of this excitation is constant to within experimental accuracy of about 2 %. This feature is sitting near the foot of a large spectral feature which has a weak but sharp onset at 1.1 eV and is maximum near $``$5 eV. The 1.1 to 5 eV conductivity band is independent of temperature to within our measurement accuracy.

We understand the electronic conductivity spectrum in terms of a broad band of charge transfer transitions from the hybridized oxygen $`p`$-levels to the Mn $`d_{3z^2r^2}`$ levels centered at $``$5 eV and an on-site Mn $`d`$-$`d`$ transition centered at 1.7 eV. The fact that this feature is observed in the hexagonal phase but not ortho-YMnO3 is consistent with selection rules for the on-site Mn $`d`$-$`d`$ transitions in hexagonal and (near) cubic symmetry. Its interpretation as a charge transfer transition between Mn neighbors is ruled out because the spin dependence of the charge transfer matrix elements would lead to a strong temperature dependence of its oscillator strength which is not observed.

We have calculated the energies and spectral intensities of electronic transitions in the MnO5 complex in the framework of ligand-field theory. The electronic states are taken to be the Mn $`d`$-orbitals coupled to the oxygen $`p_\sigma `$ orbitals. Taking the energy of the $`d_{x^2y^2,xy}`$ pair as zero, the energy of O $`p`$-orbitals is $`\mathrm{\Delta }=3`$ eV, and the crystal field energy of the $`d_{3z^2r^2}`$ orbital is $`\mathrm{\Delta }_{CF}=0.7`$ eV. For the symmetric (paraelectric) structure, the hybridization is $`t_1`$$`=`$$`t_2`$$`=`$1.9 eV for apical and $`t_3`$$`=`$$`t_4`$$`=`$1.7 eV for in-plane oxygens. The allowed optical transitions for the $`Ec`$ polarization of light are: $`d_{(x^2y^2)^{},(xy)^{}}d_{(3z^2r^2)^{}}`$ at 1.58 eV and from the two in-plane bonding O $`pd_{(3z^2r^2)^{}}`$ at $``$6.3 eV, where the prime denotes the corresponding $`d`$-states as hybridized with the oxygen orbitals. The calculated spectral weights of both types of transitions are a factor of four smaller than observed but their ratio is correct. This discrepancy may be attributed to the neglect of the Mn $`p`$\- and $`s`$-orbitals which can also couple to the O $`p`$ states in the hexagonal symmetry. In the case of $`E||c`$ polarization of incident light the optical matrix elements for the $`d_{(x^2y^2)^{},(xy)^{}}d_{(3z^2r^2)^{}}`$ transitions are zero as we find experimentally.

What is the origin of the $``$0.15 eV shift of the 1.7 eV peak with temperature shown in Fig. 1? The anomaly at $`T_N`$ indicates that at least part of the shift is associated with the magnetic phase transition. A spectral shift due to thermal expansion is also expected. However, the observed shift is larger than is typically observed in interband features in solids ($``$ 0.1 eV). We attribute the magnetic part of the shift to the effects of the exchange interactions between the Mn ions. Level shifts of the Mn $`d`$-levels due to superexchange between Mn neighbors leads to a lowering of the $`d_{x^2y^2,xy}`$ levels in the antiferromagnetic state while the relatively isolated $`d_{3z^2r^2}`$ orbital is little affected. We believe that the shift in the resonance energy between 4 K and $`T_N`$ ($``$0.05 eV) underestimates the exchange energy. This is partly because short range antiferromanetic correlations in this frustrated magnetic system are expected and are observed to persist to higher temperatures Sat ; Katsufuji et al. (2001). Also, the absence of a shift in the 5 eV feature argues against such a large thermal shift in the 1.7 eV feature since they share the same final state. In addition, the extremely weak decrease in the oscillator strength of the 1.7 eV feature also argues against a strong thermal expansion effect since this optical transition, which is allowed only due to the hybridization with the O $`p`$ states, is more sensitive to the lattice constant than are the level shifts. At present, however, we cannot separate the exchange effects from the effects of the thermal expansion on the optical transition energy. This will require a detailed analysis of the thermal effects. For completeness we note that magnetostriction can also produce level shifts due to the changes in the lattice constant. However, these effects are much smaller than the effects of thermal expansion and can be safely neglected.

We can also estimate the exchange energy of the Mn spins from the Weiss temperature in the susceptibility. From molecular field theory: $`k_B\theta =zJ_{NN}S(S+1)/3`$ where $`\theta `$ is the Weiss temperature, $`z`$=6 is the number of nearest neighbors. Therefore the exchange energy for the manganese ion is $`E_{ex}3k_B\theta `$. In the literature $`\theta `$ data for LuMnO3 range between -519 K Tomuta et al. (2001) and -887 K Katsufuji et al. (2001), so that $`E_{ex}`$ is in the range 140 to 240 meV, which is somewhat larger than our estimate from the optical shifts. As noted, however, this estimate ignored the additional shift expected above $`T_N`$ due to the effects of frustration. Also, these two estimations represent different manifestations of the exchange interaction. The Weiss temperature represents the ground state exchange energy between the manganese moments. While the optical shift represents the change in exchange energy of the Mn ion between the ground and excited state of the Mn3+ ion. Therefore, the two estimates are in satisfactory agreement and our assignment of the magnetic shift of the 1.7 eV peak with temperature to the differences in exchange interaction in excited and ground optical states of a given Mn ion is reasonable.

The phonon part of the conductivity spectrum of LuMnO3 is shown in Fig. 2 (a) and (c). Mode assignments can be made by analogy with the assignment in YMnO3 Iliev et al. (1997). Panels (b) and (d) in Fig. 2 show the temperature dependence of the frequency shifts. The two low frequency vibrational modes in the $`Ec`$ polarization display the strongest absolute frequency shifts ($`\mathrm{\Delta }\omega _0`$), relative frequency shifts ($`\mathrm{\Delta }\omega _0/\omega _0`$), and inflection points at $`T_N`$. This observation suggests that these modes are coupled to the spin system. Measurements of the phonon spectrum in the $`E||c`$ polarization on a $`ac`$-plane sample show only a thermal shift of the phonons (Fig. 2d). The frequency shift of the spin-coupled phonon appears very similar to the temperature dependence of the 1.7 eV feature. This suggests that both are related to the same, nearest-neighbor spin correlation function: $`<𝐒_i𝐒_j>(T)`$.

The phonon shifts are understood in terms of the phonon induced modulation of the exchange energy produced by the ion modal displacements. The change in exchange energy produces a corresponding change in the effective restoring force for the phonon. The resulting shift will be $`\mathrm{\Delta }\omega ^2(T)=\mathrm{\Delta }k/M`$, where $`\mathrm{\Delta }k<𝐒_i𝐒_j>`$ is the exchange energy contribution to the force constant and $`M`$ is the reduced mass of the phonon mode. In addition to frequency shifts due to the exchange effects there will also be shifts due to the thermal expansion of the lattice. The exchange energy is dominated by superexchange between nearest neighbor Mn ions separated by oxygens. However, due to the complexity of the hybridized Mn $`d`$-states there are both ferromagnetic and antiferromagnetic contributions that differ for in- and out-of-plane neighbors.

The anomaly in temperature dependence of the static in-plane dielectric constant below the N$`\stackrel{`}{e}`$el temperature is considered to be one of the manifestations of coupling between magnetic and ferroelectric order parameters in the ferroelectromagnets Smolenskii and Chupis (1982 \[Sov. Phys. Usp. 25, 475 (1982). On the other hand, $`ϵ_0`$ is determined by all of the oscillators present in the optical response of the system. We have examined the contributions to the quasi-static dielectric constant of LuMnO3. Three groups of oscillators contribute to $`ϵ_0`$: ferroelectric domains, phonons, and electronic transitions. The contribution from the antiferromagnetic resonance, centered at 50 cm-1 at 10 K, is negligible. The contributions from ferroelectric domains falls off with frequency becoming negligible at the MHz frequencies of the quasi static measurements and, even more so, at the far-infrared frequencies of our measurements. Figure 3 shows the temperature dependence of the real part of the dielectric constant of LuMnO3. The top curve in panel (a) reproduces quasi-static data of Katsufuji et al. Katsufuji et al. (2001). The bottom curve in panel (a) was calculated from the measured frequency shift of the first interference maximum in the transmittance spectrum centered at 53 cm-1. The value of $`ϵ_1`$ at this frequency is determined by all optical phonons and all electronic oscillators. Panel (b) contains $`ϵ_1`$ determined from the interference fringes in the mid-infrared transparency region midway between the phonon and electronic absorption bands. The value of $`ϵ_1ϵ_{\mathrm{}}`$ at this frequency is due to the electronic transitions. Panel (c) shows an estimation of the input to the dielectric constant from the 1.7 eV electronic peak using $`\mathrm{\Delta }ϵ_1=2/\pi _{1eV}^{2.5eV}𝑑\omega ϵ_2(\omega )/\omega `$. From the data of Fig. 3, we conclude that practically all of $`\mathrm{\Delta }ϵ`$ below $`T_N`$ comes from phonon hardening. Comparing panels (a) and (b), only $``$5 % of this change is due to the shift of the 1.7 eV electronic peak. We can also conclude from panels (b) and (c) that the temperature dependence of $`ϵ_{\mathrm{}}`$ is almost entirely due to the 1.7 eV feature.

The absence of a magnetic anomaly in phonons in the $`E||c`$ polarization is consistent with the observed anisotropy in the static dielectric constant Katsufuji et al. (2001). Toward developing an understanding of this anisotropy of the spin lattice coupling we note that there are important distinctions between the $`E||c`$ and $`Ec`$ polarizations in the exchange modulation for this layered system. In the case of $`E||c`$ polarization the relative displacement of the in-plane Mn and O ions is transverse to their bond and bond stretching occurs only with the apical oxygen, primarily affecting the Mn $`z^2`$ orbital which is empty. Whereas for $`E`$ in the $`ab`$-plane both transverse and bond stretching Mn-O displacements are induced and the occupied in-plane orbitals are involved. This observation suggests that the in-plane bond stretching displacements dominate the spin-phonon interaction effects. The role of dynamic effective charge transfer for the bond stretching modes in transition metal compounds for ferroelectricity and the electron-phonon interaction has been recently discussed (see, for example, Resta (1994)). We suggest that the observed anisotropy of the spin-phonon coupling and $`ϵ_0`$ might be a consequence of dynamic charge redistribution for in-plane vibrations of Mn ions, which modulate the in-plane partially covalent bonds, and the absence of the effective in-plane charge redistribution for the out-of-plane modes.

More generally the question of a coupling between the two order parameters of this material is interesting. Within the Landau theory of phase transitions there are symmetry allowed terms in the Landau free energy describing the coupling between the magnetic and ferroelectric order Smolenskii and Chupis (1982 \[Sov. Phys. Usp. 25, 475 (1982). The experimental data on the quasi-static dielectric constant implies that this term is of the form $`\delta FL^2P^2(E_x^2+E_y^2)`$, where $`P`$ and $`L`$ are, respectively, the ferroelectric and antiferromagnetic order parameters. Establishing this coupling and relating it to the microscopic physics is a key issue in the study of this class of materials.

In conclusion, we have observed a strong coupling of the antiferromagnetism in LuMnO3 to a sharp low energy interband transition and to the infrared phonon spectrum. The optical feature has a large blue shift associated with the antiferromagnetism which is caused by the effects of the exchange interaction on the on-site Mn $`d`$-$`d`$ transition. A similar anomaly in the temperature dependence of the phonon frequencies is attributed to effects of spin-phonon coupling. These results demonstrate that optical spectroscopy is a powerful tool in the study of exchange interaction in the strongly frustrated magnetic system of the hexagonal manganites.

###### Acknowledgements.

This work supported in part by NSF-MRSEC DMR 0080008. Fruitful discussions with P. B. Allen and D. B. Romero are gratefully acknowledged. The work of J. S. Ahn is supported by the Korean Science and Engineering Foundation(KOSEF) through the Center for Strongly Correlated Materials Research (CSCMR) at Seoul National University.


text:
# Correlated two-electron transport: a principle for a novel charge pump

## Abstract

By considering a correlated two-electron transport process (TET) and using a diagrammatic analysis within the Keldysh nonequilibrium Green’s function formalism, we discuss a novel charge pump by which carriers are pumped from a contact with low chemical potential to another contact with a higher potential. The TET process involves two correlated incident electrons scattering and exchanging energy with each other. The process can significantly affect charge current density and it involves high empty states and/or low filled states of the Fermi liquid of the leads.

In the past two decades, quantum electron transport through meso- and nanoscopoic structures have received considerable attention. Charge transport theories for nanostructures typically consider the situation of an incident electron from a lead that is scattered in the device scattering region and transmitted to other leads or reflected to where it came from. The scattering region can be rather complicated and may involve electron-electron, electron-phonon, electron-impurity, and other interactions and scattering mechanisms. The scattering region itself can be a semiconductor nanostructure such as a quantum dot (QD), a carbon nanotube, or a single molecule. The scattering processes may involve exchange of energy or spin. While situations and physics can vary in a wide range of ways, charge transport in nanostructures has largely, so far, been considered as involving a single electron incident from a non-interacting lead of the device, it traverses through a complicated scattering region in which interactions occur, and finally exits the device through the non-interacting leads. We will refer them as single electron tunneling or single electron transport (SET) processes.

In this paper, we go one step further by investigating correlated two-electron tunneling and transport (TET) processes which, to the best of our knowledge, has not been studied before. In this process, two incident electrons with energies $`ϵ_1`$ and $`ϵ_2`$ from the leads tunnel into the scattering region of a device, they scatter with each other and exchange energy (they could also scatter and exchange energy with other particles in the scattering region). Following the scattering event, they transmit to the outside world through the contacts, but with energies $`ϵ_1^{^{}}`$ and $`ϵ_2^{^{}}`$ which are different from the initial values $`ϵ_1`$, $`ϵ_2`$ due to energy exchange. If $`ϵ_1+ϵ_2=ϵ_1^{^{}}+ϵ_2^{^{}}`$, i.e. the two incident electrons do not exchange energy with other particles, the transport process as a whole can be considered as elastic. In the TET process, the variation of electron number in the scattering region is always two: a situation that is different from SET. We show that such TET processes can indeed occur, and it affects measurable physical quantities such as charge current. In particular, TET induces a transport vertical flow—meaning carriers incident with energy $`ϵ_1`$ exit at a different energy $`ϵ_1^{}`$, resulting to a non-conserved current density, i.e. $`_nj_n(ϵ)0`$. Using the characteristics of the correlated two-electron transport process, we design an interesting charge pump so that incident carriers can be pumped from a contact with low chemical potential to another contact having a higher potential.

We consider a QD coupled to two or three leads, described by the following Hamiltonian: $`H=_\alpha ϵ_\alpha d_\alpha ^{}d_\alpha +Ud_{}^{}d_{}d_{}^{}d_{}+_{n,k,\alpha }ϵ_{nk}a_{nk\alpha }^{}a_{nk\alpha }+_{n,k,\alpha }[t_{nk}a_{nk\alpha }^{}d_\alpha +H.c.]`$, where $`a_{nk\alpha }^{}`$($`a_{nk\alpha }`$) and $`d_\alpha ^{}`$($`d_\alpha `$) are creation (annihilation) operators in lead $`n`$ and in QD, respectively. The QD includes two states and has an intradot Coulomb interaction $`U`$. Subscript $`\alpha `$ is the spin index, it may also indicate other quantum numbers. To account for a magnetic field, we let $`ϵ_{}ϵ_{}`$.

In the following we focus on investigating the elastic TET by analyzing the behavior of current density $`j_n(ϵ)`$ from lead $`n`$ to the QD. $`j_n(ϵ)`$ relates to the current $`I_n`$ through $`I_n=j_n(ϵ)𝑑ϵ`$. We also define the electron occupation number density operator $`\widehat{N}_n(ϵ,\tau )=_{k,\alpha }e^{iϵt}a_{nk\alpha }^{}(\tau )a_{nk\alpha }(\tau +t)\frac{dt}{2\pi }`$. $`\widehat{N}_n(ϵ,\tau )dϵ`$ describes the electron occupation number in lead $`n`$ in the energy range $`ϵ`$ to $`ϵ+dϵ`$ at time $`\tau `$. Current density $`j_n(ϵ)`$ can be calculated from the time evolution of $`\widehat{N}_n(ϵ,\tau )`$: $`j_n(ϵ)=e<\frac{d}{d\tau }\widehat{N}_n(ϵ,\tau )>`$. This leads to (in units of $`\mathrm{}=1`$):

$$j_n(ϵ)=eIm\underset{\alpha }{}\frac{\mathrm{\Gamma }_n(ϵ)}{2\pi }\left[2f_n(ϵ)G_\alpha ^r(ϵ)+G_\alpha ^<(ϵ)\right]$$
(1)
where linewidth function $`\mathrm{\Gamma }_n(ϵ)2\pi _k|t_{nk}|^2\delta (ϵϵ_{nk})`$; $`f_n(ϵ)`$ is the Fermi distribution function for lead $`n`$; $`G_\alpha ^{r,<}(ϵ)`$ are the retarded and Keldysh Green’s functions of the QD. Using the standard equation of motion technique, $`G_\alpha ^r(ϵ)`$ has already been solved in previous work:

$$G_\alpha ^r(ϵ)=\frac{1+UA_\alpha n_{\overline{\alpha }}}{ϵϵ_\alpha \mathrm{\Sigma }_\alpha ^{(0)}(ϵ)+UA_\alpha (\mathrm{\Sigma }_{\overline{\alpha }}^a+\mathrm{\Sigma }_{\overline{\alpha }}^b)}$$
(2)
where $`A_\alpha (ϵ)=[ϵϵ_\alpha U\mathrm{\Sigma }_\alpha ^{(0)}(ϵ)\mathrm{\Sigma }_{\overline{\alpha }}^{(1)}(ϵ)\mathrm{\Sigma }_{\overline{\alpha }}^{(2)}(ϵ)]^1`$; $`\mathrm{\Sigma }_\alpha ^{(0)}`$ is the lowest-order self-energy, $`\mathrm{\Sigma }_\alpha ^{(1),(2),a,b}`$ are the higher-order self-energies; and $`n_\alpha =Im\frac{dϵ}{2\pi }G_\alpha ^<(ϵ)`$ is the intradot electron occupation number of state $`\alpha `$. It is worth to mention that if temperature $`T`$ is lower than the Kondo temperature $`T_K`$, the solution of Eq.(2) has a Kondo resonance at the Fermi level which was the subject of many previous studies.

We investigate TET processes at a temperature higher than $`T_K`$. To this end we need to solve the the Keldysh Green’s function $`G_\alpha ^<(ϵ)`$. Note that if one applies the commonly used ansatz for interacting lesser and greater self-energies, or using the large-$`U`$ limit non-crossing approximation to solve $`G_\alpha ^<(ϵ)`$, the two-electron scattering will be lost in these approximations. Therefore a more precise analysis is needed in our problem and we proceed as follows. Introducing the intradot electron occupation number density operator $`\widehat{N}_\alpha (ϵ,\tau )=e^{iϵt}d_\alpha ^{}(\tau )d_\alpha (\tau +t)\frac{dt}{2\pi }`$ and due to the steady state condition $`<\frac{d}{d\tau }\widehat{N}_\alpha (ϵ,\tau )>=0`$, we have:

$$e\underset{n}{}\frac{i\mathrm{\Gamma }_n}{2\pi }\left[G_\alpha ^<+f_n(G_\alpha ^rG_\alpha ^a)\right]=i_\alpha (ϵ)$$
(3)
where $`i_\alpha (ϵ)`$ gives the intradot vertical flow: $`i_\alpha (ϵ)ieUe^{iϵt}\left\{<d_\alpha ^{}(0)d_{\overline{\alpha }}^{}(0)d_{\overline{\alpha }}(0)d_\alpha (t)><d_\alpha ^{}(0)d_{\overline{\alpha }}^{}(t)d_{\overline{\alpha }}(t)d_\alpha (t)>\right\}\frac{dt}{2\pi }`$. The quantity $`i_\alpha (ϵ)`$ can be viewed as an intradot “vertical” current density at energy $`ϵ`$ that is contributed by carriers with other energies due to the e-e scattering process. Summing over the index $`\alpha `$, Eq.(3) reduces to $`_nj_n(ϵ)+_\alpha i_\alpha (ϵ)=0`$, which is exactly the steady state current conservation equation so that the total current density—including the vertical flow, through the intradot energy level $`ϵ`$ is zero. If $`G_\alpha ^r(ϵ)`$ and $`i_\alpha (ϵ)`$ have been solved, from Eq.(3) $`G_\alpha ^<(ϵ)`$ can be obtained immediately.

Next, we solve the vertical flow quantity $`i_\alpha (ϵ)`$. We introduce two particle contour-ordered Green’s functions $`B_\alpha (t,0)`$ and $`B_c(t,0)`$, $`B_\alpha (t,0)<T_C[d_\alpha (t)d_{\overline{\alpha }}(0^+)d_\alpha ^{}(0)d_{\overline{\alpha }}^{}(0)]>`$ and $`B_c(t,0)<T_C[d_\alpha (t)d_{\overline{\alpha }}(t)d_\alpha ^{}(0)d_{\overline{\alpha }}^{}(0)]>`$. Although there are four operators in the definition $`B_{\alpha /c}(t,0)`$, only two time indices at $`(t,0)`$ appear. We can therefore write the contour ordered quantities as $`B_{\alpha /c}^{++}(t,0)`$, $`B_{\alpha /c}^+`$, $`B_{\alpha /c}^+`$, and $`B_{\alpha /c}^{}`$. Their Fourier transformations can be defined as $`B(ϵ)e^{iϵt}B(t,0)𝑑t`$. Using the two-particle Green’s function, the vertical flow $`i_\alpha (ϵ)`$ is reduced to $`i_\alpha (ϵ)=\frac{eU}{\pi }ImB_\alpha ^+(ϵ)`$.

The contour Green’s function $`B_{\alpha /c}`$ is solved by means of a Feynman diagram expansion using the Wick’s theorem. We take the interacting part of the Hamiltonian as $`H_I=Ud_{}^{}d_{}d_{}^{}d_{}`$, and the remaining part $`HH_I`$ as the non-interacting Hamiltonian $`H_0`$. The first-order irreducible self-energy which we consider is shown in Fig.1. It is clear that the first-order graph describes a two-particle propagation involving the exchange of an interacting energy $`U`$. Two important points should be mentioned here. (i) The vertex we calculate sums up all reducible diagrams constructed by the irreducible self-energy of Fig.1 (also see Eq.(4)). This level of approximation is equivalent to that of a typical random phase approximation, i.e. we compute the irreducible self-energy up to order $`U^1`$. It is reasonable to neglect other higher order irreducible diagrams which are of orders $`U^2`$ and higher, because the interaction energy $`U`$ is large. (ii) Single solid lines in Fig.1 stand for the intradot contour-ordered Green’s functions $`G_\alpha `$ of the Hamiltonian $`H`$ (not $`H_0`$). This means that we have summed over all terms in the Feynman diagram expansion of $`G_\alpha `$.

The corresponding equation for Fig.1 is

$`B_{c/\alpha }(t,0)=G_\alpha (t,0)G_{\overline{\alpha }}(t/0^+,0)`$ (4)
$`+iU{\displaystyle _C}𝑑t_1G_\alpha (t,t_1)G_{\overline{\alpha }}(t/0^+,t_1)B_c(t_1,0).`$ (5)
Finally, we get $`B_c^+(t,0)`$, $`B_c^{}(t,0)`$, and $`B_\alpha ^+(t,0)`$ from Eq.(4), and upon taking a Fourier transformation, $`B_\alpha ^+(ϵ)`$ can be expressed in terms of $`G_\alpha ^{++}(ϵ)`$, $`G_\alpha ^+(ϵ)`$, $`G_\alpha ^+(ϵ)`$, and $`G_\alpha ^{}(ϵ)`$—these four Green’s functions are directly related to $`G_\alpha ^r`$ and $`G_\alpha ^<`$. This completes the analytical derivations.

From Eqs. (2,3,4), the intradot occupation number $`n_\alpha `$ and the vertical flow $`i_\alpha (ϵ)`$ are determined self-consistently. Let’s consider a two-probe ($`n=L,R`$) device and the wideband approximation in which $`\mathrm{\Gamma }_n(ϵ)`$ is independent of energy $`ϵ`$. Fig.2 plots the vertical flow $`i(ϵ)=_\alpha i_\alpha (ϵ)`$ (solid line) as well as the current density $`j_{L/R}(ϵ)`$ at a high bias. Clearly, the vertical flow $`i(ϵ)`$ is non-zero due to the TET processes, and its value $`|i(ϵ)|`$ has four peaks at energies $`ϵ_{}`$, $`ϵ_{}`$, $`ϵ_{}+U`$, and $`ϵ_{}+U`$, respectively. This means that the incident electron indeed can vary its energy in QD by e-e scattering. In contrast, in a typical SET process the electron keeps its energy and does not induce any vertical flow $`i(ϵ)`$.

The physics of the TET process that induces the vertical flow is shown by inset (a) of Fig.2. To start, two incident electrons from left lead having energies $`ϵ_{}`$ and $`ϵ_{}+U`$ tunnel into the QD. They scatter with each other inside the QD and exchange energy to final states $`ϵ_{}`$ and $`ϵ_{}+U`$. Afterwards they tunnel out of the QD. In the vertical flow curve of Fig.2, two peaks at $`ϵ_{}`$ and $`ϵ_{}+U`$ are negative (dips), and the other two peaks at $`ϵ_{}`$ and $`ϵ_{}+U`$ are positive: precisely indicating the transfer of states from the initial ones at $`ϵ_{}`$ and $`ϵ_{}+U`$ to the final ones at $`ϵ_{}`$ and $`ϵ_{}+U`$. It is also worth mentioning that besides the new TET process, the usual SET processes also exist in charge transport through the QD in the present case.

So far we have demonstrated that TET processes can exist. In the following we analyze several important questions concerning TET. What is its consequence? (i) TET makes current density a non-conserved quantity, i.e. $`_nj_n(ϵ)0`$. Of course, the total current is still conserved, i.e. $`_nI_n=0`$. This can be easily proved from the definition of the vertical flow $`i(ϵ)`$, namely $`i(ϵ)`$ has the property $`i(ϵ)𝑑ϵ=0`$. (ii) TET can involve high energy empty states, namely states which are higher by about $`U`$ than the highest chemical potential $`max(\mu _L,\mu _R)`$ (see inset (a) of Fig.2); it may also involve electrons deep inside the Fermi sea, namely states which are lower by about $`U`$ than the lowest chemical potential $`min(\mu _L,\mu _R)`$ (see TET process shown in inset (b) of Fig.2). (iii) TET may induce a current density that is flowing out from the high voltage terminal, i.e. the left lead (indicated by the negative peak at $`ϵ_{}+U`$ in $`j_L(ϵ)`$ curve, and by the arrow A in inset (a) of Fig.2). Similarly, TET may also induce a current density that is flowing in from the low voltage terminal, i.e. the right lead (see arrow A in inset (b) of Fig.2). These characteristics are rather different from the typical elastic SET processes.

Under what conditions does TET or the vertical flow $`i(ϵ)`$ exist? (i) We found that an increase (decrease) of temperature $`T`$ or linewidth $`\mathrm{\Gamma }`$ will widen (narrow) the peaks of vertical flow, but does not affect peak positions and heights significantly. (ii) If $`U=0`$, the vertical flow $`i(ϵ)=0`$ identically: TET crucially depends on this parameter. (iii) If $`U\mathrm{}`$, $`i(ϵ)`$ tends to zero. This is because at large $`U`$, the intradot two-electron occupation is prohibited therefore TET is blockaded. In this case only SET processes occur. (iv) When bias potential $`eV=\mu _L\mu _R`$ is less than $`U`$, $`i(ϵ)`$ decreases drastically. In the limit of $`eV=0`$, $`i(ϵ)=j_n(ϵ)=0`$.

Are observable quantities of charge transport affected by the TET process? (i) Clearly the current density $`j_{L/R}(ϵ)`$ is affected significantly as already discussed above. (ii) In general, the current, conductance, and $`n_\alpha `$ will be affected significantly by TET (see below). The current noise, which reflects the e-e time correlation, will increase due to TET processes. However, if one uses the wideband approximation, the TET dependence in charge current will be lost.

In the rest of this paper, we apply the property of TET to design a device so that electrons can be pumped from a lead with lower chemical potential to another lead having a higher chemical potential. Consider a device with three leads ($`n=1,2,3`$) and consider the non-wideband case. We use a model of quasi-square bands where the coupling $`\mathrm{\Gamma }_n(ϵ)=\mathrm{\Gamma }/\{exp[(|ϵc_n|W)/0.05]+1\}`$, the width of the band is set by $`2W=1`$ and its center at $`c_n`$ which is dependent on the terminal voltage $`eV_n=\mu _n`$ but $`\mu _nc_n`$ is kept fixed. More specifically, let’s assume lead 1 to be a p-type semiconductor with $`\mu _1c_1=0.4`$; lead 2 an n-type semiconductor with $`\mu _2c_2=0.4`$; and lead 3 a metal with $`\mu _3c_3=0`$ (see Fig.4). The energy diagram of the device is set by external voltages as that shown in Fig.4 so that $`\mu _2>\mu _1>\mu _3`$. The current density $`j_n(ϵ)`$ in this case is shown in Fig.3. We note that $`j_1(ϵ)`$ (dotted line) has two positive peaks at $`ϵ_{}`$ and $`ϵ_{}+U`$; $`j_2(ϵ)`$ and $`j_3(ϵ)`$ each has one negative peak at $`ϵ_{}+U`$ and $`ϵ_{}`$, respectively; and $`i(ϵ)=_nj_n(ϵ)`$ has two negative peaks at $`ϵ_{}`$ and $`ϵ_{}+U`$, and two positive peaks at $`ϵ_{}+U`$ and $`ϵ_{}`$. The current $`I_n=j_n(ϵ)𝑑ϵ`$ is quite large. We emphasize two points for this pump. (i) In this device the SET process almost does not occur because bands of different leads do not overlap. Then, clearly, the large current $`I_n`$ originates from the TET process: $`I_n0`$ if the vertical flow $`i(ϵ)0`$. This demonstrates that TET can significantly affect charge current in the general case of non-wideband coupling. (ii) The charge current in the terminal with the highest bias voltage, e.g. lead 2, is negative (solid line in the inset of Fig.3), which demonstrates the pump effect. The pumps works because when an electron tunnels from lead 1 to 3, it emits energy $`U`$ to pump another electron from lead 1 to lead 2, through the TET process.

More clearly, the working principle of the TET pump is summarized in Fig.4. (a) We start from the situation where no charge is in the QD so that levels $`ϵ_{}`$ and $`ϵ_{}`$ are empty. In this situation, an electron in the Fermi sea of lead 1 having energy $`ϵ_{}`$ can easily tunnel into the QD (Fig.4a). (b) After this electron tunnels into the QD and occupies the QD level of $`ϵ_{}`$, the other intradot level $`ϵ_{}`$ is raised up to $`ϵ_{}+U`$, so that another electron with energy $`ϵ_{}+U`$ in lead 1 tunnels into the QD (Fig.4b). (c) When the second electron comes into QD, due to the e-e Coulomb interaction $`U`$, the level $`ϵ_{}`$ with its electron is raised up to $`ϵ_{}+U`$, leading to a negative peak at $`ϵ_{}`$ and a positive peak at $`ϵ_{}+U`$ in the vertical flow curve $`i(ϵ)`$. Now the intradot two-electron system has total energy $`ϵ_{}+ϵ_{}+U`$ (Fig.4c). Afterwards the first electron in state $`ϵ_{}+U`$ easily tunnels to lead 2 and takes away energy $`ϵ_{}+U`$. The net effect is that the two electrons exchanged energy $`U`$, which is the TET process discussed above. When the first electron leaves the QD, the other electron at $`ϵ_{}+U`$ falls down to $`ϵ_{}`$, leading to a negative peak at $`ϵ_{}+U`$ and a positive peak at $`ϵ_{}`$ in the curve of $`i(ϵ)`$. (d) Finally, the second electron tunnels to lead 3 (Fig.4d) and our device returns to its initial conformation of Fig.4a. This way an electron is pumped from lead 1 to lead 2, where $`\mu _1<\mu _2`$, via the TET process. We emphasize that each tunneling event from Fig.4a-d is a first-order normal tunneling event in which tunneling occurs at two aligning states (not like higher-order virtual co-tunneling process). We therefore conclude that the TET process should have large probability to occur so that $`i(ϵ)`$ can be near the unit value $`e/h`$ (see Fig.3). We have also investigated the terminal voltage (e.g. $`V_3`$) dependence of current $`I_n`$, shown in the inset of Fig.3. As $`V_3`$ is increased so that $`eV_3=\mu _3`$ passes the lowest resonance state $`ϵ_{}`$, the tunneling event in Fig.4d can not occur and the TET process is blockaded, leading to a significant reduction of all currents $`I_n`$ (including the pumping current $`I_2`$) (see inset of Fig.3).

In summary, we have investigated the two-electron correlated scattering process in mesoscopic system. TET induces a vertical flow in the scattering region so that electrons enter and exit the device with different energies. TET is found to affect current density significantly, and the process can involve high empty states or/and low filled states of the leads. The properties of TET suggests an interesting working principle of an electron pump which pumps charge carries from a lead with low chemical potential to another lead with a higher chemical potential. In fact, if the bands of lead 1 is full and bands of leads 2 and 3 are empty, these results are not affected. Our proposed pump is very different as compared to the electron-photon or parametric pumps. It should be experimentally feasible even for devices fabricated in two-dimensional electron gas. In that case, the bands in our theory can be replaced with Landau levels.

Acknowledgments: We gratefully acknowledge financial support from NSERC of Canada, FCAR of Quebec (Q.S., H.G), and a RGC grant from the SAR Government of Hong Kong under grant number HKU 7091/01P (J.W.).


text:
# Boundary conditions as constraints

## 1 Introduction

In recent years, the study of the boundary conditions in Quantum Theory has produced several important results. For example, an interesting idea is to consider the black hole event horizon as a physical boundary , this induce an extra term in the action, having as consequence the existence of a central charge in the algebra of generators of gauge transformations . Using this central charge it is possible to determine the asymptotic behavior of the density of states and in this way to get the entropy for a black hole in 2+1-dimensions .

In the context of string theory, the D-branes are natural boundaries for the open strings. These boundaries have very interesting effects in the theory among them is the non-commutativity of the D-brane coordinates . This non commutativity appears due to the presence of a constant B-field on the boundary that implies a change of the symplectic structure on the boundary. This change of the symplectic structure has been computed using different methods, the most simple is the direct solution of the Field equations subject to the Dirichlet boundary conditions on the 9-p transverse directions and Neumann boundary conditions on p+1 directions parallel to the Dp-brane. However, for some interesting systems is not possible to solve exactly this problem and the question appears if is there exist other method to compute the symplectic structure.

An alternative to solve this problem was to consider that the boundary conditions in a Quantum Field Theory can be interpreted as second class constraints in the sense of the Dirac’s Method . This procedure has some interesting characteristics, the non commutativity appears naturally due to a modification of the symplectic structure given by the Dirac’s brackets. These brackets are constructed with the second class constraints that arise from the boundary conditions and the time evolution of these conditions. However, the procedure have some problems, see for example . First the Lagrange multipliers are fixed by hand and not following the standard Dirac’s method where the Lagrange multipliers are fixed in the case of second class constraints by the time evolution of the constraints. This has as consequence that appears an infinite number of constraints and then a minus infinite number of degrees of freedom in the boundary. To solve these problems we propose an alternative procedure that follows in many aspects the previous proposal. However we have a very different interpretation that allows solve most of the shortcomings. The key point of our procedure is to perform a mapping from the original problem with a given Hamiltonian and boundary conditions to another problems with the same Hamiltonian, but now with second class constraints equal to the boundary conditions and no boundary. Our conjecture is that the results that we get for the Dirac’s brackets projected to the interior of the new problem, is the symplectic structure of the original problem. We show that this conjecture is valid in all examples that we know and inclusive we check that in the example of PP-waves  our result is fully consistent whereas the previous result it is not. In the literature exist several proposals  to solve the problems associated to the procedure of consider the boundary conditions as Dirac’s constraints. However, all the previous proposals are not useful in at least in one of the examples that we present, whereas our procedure is fully consistent and systematic in all these examples.

In section 2 we present an outline of our procedure, in the next section we analyze the case of the scalar field with Neumann and Dirichlet boundary conditions. In section 4 we present the case of the bosonic string with mixed boundary conditions and in section 5 we study the case of the bosonic string in a PP-wave background.

## 2 Boundary conditions as constraints

### 2.1 Boundary conditions and the action principle

Let $`M`$ be a $`(d+1)`$-dimensional manifold, with topology $`\mathrm{\Sigma }\times `$, where $`\mathrm{\Sigma }`$ is an oriented $`d`$-manifold with boundary $`\mathrm{\Sigma }`$. In $`M`$ we consider a Field Theory given by the action

$`S={\displaystyle _M}𝑑x^{d+1}L(\varphi _a(x),_\alpha \varphi _a(x)),`$ (1)
for the fields $`\varphi _a`$. The conditions that the integral (1) be stationary implies

$$\begin{array}{cc}\hfill \delta S=& _M𝑑x^{d+1}[\frac{L}{\varphi _a(x)}_\alpha (\frac{L}{(_\alpha \varphi _a(x))})]\delta \varphi _a(x)\hfill \\ & +_M𝑑x^{d+1}_\alpha [(\frac{L}{(_\alpha \varphi _a(x))})\delta \varphi _a(x)]=0.\hfill \end{array}$$
(2)
For a system without boundary the second term cancel automatically, but in our case we need to impose boundary conditions. There are three different ways to do that:

Dirichlet conditions:

$$\delta \varphi _a(x)|_{xϵ\mathrm{\Sigma }}=0$$
(3)
Neumann conditions:

$$\left(\frac{L}{(_\alpha \varphi _a(x))}\right)(x)|_{xϵ\mathrm{\Sigma }}=0,$$
(4)
or one combination of both types for the components of the field $`\varphi _a`$.

On the other hand, we have the canonical formalism defined by the Hamiltonian,

$`H_c={\displaystyle _\mathrm{\Sigma }}𝑑x^d_c(\varphi ,\varphi ,\mathrm{\Pi },\mathrm{\Pi }),`$ (5)
the symplectic structure,

$`\{\varphi _a(x,t),\varphi _b(x^{},t)\}=\{\mathrm{\Pi }_a(x,t),\mathrm{\Pi }_b(x^{},t)\}=0`$ (6)
$`\{\varphi ^a(x,t),\mathrm{\Pi }_b(x^{},t)\}=\delta _b^a\delta (xx^{}),`$ (7)
and the boundary conditions,

$$F_a(\varphi ,\varphi ,\mathrm{\Pi },\mathrm{\Pi })|_{xϵ\mathrm{\Sigma }}=0.$$
(8)
This set of boundary conditions is equivalent in the canonical formalism to the equations (3) or (4). In the definition of the canonical formalism we are assuming that the boundary conditions are not in contradiction with the local symmetries of the system. So from the beginning we consider that our system does not have gauge degrees of freedom or that a gauge choice consistent with the boundary conditions has been done. Also, we are assuming that the symplectic structure (6)-(7) is defined until a zero measure set, so this structure can be changed in the boundary without affect the rest of the theory. So, if the boundary conditions (8) are not consistent with the symplectic structure, it is possible to modify the structure in the boundary in such way that be consistent with the boundary conditions. For example suppose that we have the boundary conditions,

$$F_a=[\mathrm{\Pi }_a(x)\varphi _a(x)]|_{xϵ\mathrm{\Sigma }}=0,$$

these conditions are clearly inconsistent with (7), but we can modify the structure in the boundary, so we can introduce

$$\{\varphi ^a(x,t),\mathrm{\Pi }_b(x^{},t)\}|_{xϵ\mathrm{\Sigma }}=0.$$
(9)
This symplectic structure is consistent with the boundary conditions. In principle, there are many ways to modify the symplectic structure in such way that these be consistent with the boundary conditions. However, no all of these modified structures are consistent with the exact solution of the equations of motion of the system. The main result of our paper is to propose a new procedure to obtain a symplectic structure that agrees with the symplectic structure obtained by the exact solution of the problem.

### 2.2 Outline of the method to compute the symplectic structure in the boundary

Following the reference , we assume that the boundary conditions are primary constraints in the sense of the Dirac’s method. However, we want to construct a procedure fully consistent with the Dirac’s method, avoiding problems with regularization. To achieve this we consider that the primary constraints are not only valid in the boundary $`\mathrm{\Sigma }`$. Then, we suppose that these constraints are valid on all $`\mathrm{\Sigma }`$, i.e.,

$$F_a(\varphi ,\varphi ,\mathrm{\Pi },\mathrm{\Pi })|_{xϵ\mathrm{\Sigma }}0.$$
(10)
So, it make sense to write these constraints in smeared form

$$F[𝐍]=_\mathrm{\Sigma }𝑑x^dN^a(x)F_a(\varphi ,\varphi ,\mathrm{\Pi },\mathrm{\Pi })0,$$
(11)
for $`N^a(x)`$ smeared functions with compact support. In consequence, we are mapping our original problem with boundary conditions (8) and canonical Hamiltonian (5) to another problem where we have primary constraints (10) an a total Hamiltonian, given by

$$H_T=H_c+_\mathrm{\Sigma }𝑑x^d\lambda ^aF_a.$$
(12)
Using the total Hamiltonian (12) we evaluate the time evolution of the primary constraints,

$$\dot{F}[𝐍]=[F[𝐍],H_T].$$
(13)
From this procedure we can obtain new constraints. To this secondary constraints, we applied the same procedure, we evolve the constraints and see if we get new constraints. We finish the procedure when we don’t get new constraints and determine the Lagrange’s multipliers. This must be so, because we are considering that in the original problem (1) does not have gauge degrees of freedom and furthermore that the boundary condition (or primary constraints), do not generate gauge freedom in the boundary. We conclude that the system only have second class constraints, then all Lagrange’s multipliers associated to the primary constraints must be determined. The set of second class constraints that we obtain are denote by,

$$\chi _\alpha (\varphi ,\varphi ,\pi ,\pi )|_{x\mathrm{\Sigma }}0.$$
(14)
With these constraints we obtain the invertible matrix

$$C_{\alpha \beta }(x,x^{})=\{\chi _\alpha (x),\chi _\beta (x^{})\}$$
(15)
and using (15) we construct the Dirac bracket,

$$\{A(x),B(x^{})\}^{}=\{A(x),B(x^{})\}\{A(x),\chi _\alpha (y)\}C^{\alpha \beta }(y,z)\{\chi _\alpha (z),B(x^{})\}.$$
(16)
We notice that the Dirac bracket (16) is valuated on the surface $`\mathrm{\Sigma }`$, and it is not smeared. So, to obtain a Dirac bracket that is valid only in the interior of sigma $`\mathrm{\Sigma }\mathrm{\Sigma }`$, we need to eliminate from the computation of the bracket the part that is different of zero only in the boundary $`\mathrm{\Sigma }`$. Finally, we affirm that the correct symplectic structure in the boundary is obtained from the projection of this bracket to the boundary. Our procedure is in some sense equivalent to a Wick rotation, where one maps the original problem to the Euclidean space to make sense to the integrals, here we map the original problem defined with a boundary to another problem without boundary but with primary constraints that reflects the presence of the boundary. We do that in order to make sense of the Poisson brackets and in this way avoiding the problems with the regularization and obtaining a Dirac’s method fully consistent. Finally, when we obtain the symplectic structure in the new problem it is possible to obtain the symplectic structure of the original problem by projecting the structure resulting from the Dirac brackets to the boundary. We see, that our procedure is in some sense equivalent to an inverse of the Holographic map, since we obtain information about the boundary using computations in the bulk.

We don’t have a rigorous proof that the outlined procedure is always correct. However, we show in the following section that for several examples that our results coincide with the exact results obtained using the exact solution of the problem. In the example of the bosonic string in a background of PP-waves we disagree with one part of the result recently published. Nevertheless, we show that our symplectic structure is fully in the boundary whereas the result of  is inconsistent in this region.

## 3 The scalar field with Dirichlet and Neumann conditions

We introduce our procedure in a very simple example, the scalar field. This example is quite useful because we have an exact solution for boundary conditions of the type of Dirichlet and Neumann, and we show that our method reproduce correctly the symplectic structure in both cases. We start from the action

$`S={\displaystyle \frac{1}{2}}{\displaystyle _0^\pi }{\displaystyle _{t_1}^{t_2}}𝑑x𝑑t[(_t\varphi (x,t))^2(_x\varphi (x,t))^2].`$ (17)
From this action we get the equations of motion

$$(_t^2_x^2)\varphi (x,t)=0,$$
(18)
and we have two possible choices of boundary conditions

$`\varphi _i(x,t)|_{x=\pi ,0}`$ $`=0`$ $`\mathrm{Dirichlet},`$ (19)
$`(_x\varphi _i)(x,t)|_{x=\pi ,0}`$ $`=0`$ $`\mathrm{Neumann}.`$ (20)
and the canonical Hamiltonian is given by

$`H_c={\displaystyle \frac{1}{2}}{\displaystyle _0^l}𝑑x[(\mathrm{\Pi }(x,t))^2+(_x\varphi (x,t))^2].`$ (21)
### 3.1  Scalar Field with Dirichlet boundary conditions

For the Dirichlet boundary condition we have the primary constraint

$$\varphi (x,t)0.$$
(22)
That we can rewrite in densityzed form as,

$$\mathrm{\Phi }_D^{(1)}(N)=_0^\pi 𝑑xN(x)\varphi (x,t)0,$$
(23)
and the total Hamiltonian is,

$$H_T=H_c+_0^\pi 𝑑x\lambda (x)\varphi (x,t).$$
(24)
From, the time evolution of the primary constraint (22), we get

$$\dot{\mathrm{\Phi }}_D^{(1)}(N)=\{\mathrm{\Phi }^{(1)}(N),H_T\}=_0^\pi 𝑑xN(x)\pi (x)0,$$
(25)
which implies that we have the secondary constraint,

$$\mathrm{\Phi }_D^{(2)}(M)=_0^\pi 𝑑xM(x)\mathrm{\Pi }(x,t)0.$$
(26)
The time evolution of the constraint (26) produce no new constraints and we get the Lagrange multiplier,

$$\lambda =_x^2\varphi (x,t).$$
(27)
Computing the Poisson bracket of the constraints (22) and (26),

$$\{\mathrm{\Phi }_D^{(1)}(N),\mathrm{\Phi }_D^{(2)}(M)\}=_0^\pi 𝑑xM(x)N(x)0,$$
(28)
we see that the constraints are second class, we have only two constraints, and if we make a count of degrees of freedom in the boundary we obtain, zero degrees of freedom, that is the correct result. Furthermore, the Dirac brackets of the variables are,

$`\{\varphi (x,t),\varphi (x^{},t)\}^{}=\{\mathrm{\Pi }(x,t),\mathrm{\Pi }(x^{},t)\}^{}=0,`$ (29)
$`\{\varphi (x,t),\mathrm{\Pi }(x^{},t)\}^{}=0.`$ (30)
Then the projection of this brackets to the boundary, results in the symplectic structure

$`\{\varphi (x,t),\varphi (x^{},t)\}|_{x=0,\pi }=\{\mathrm{\Pi }(x,t),\mathrm{\Pi }(x^{},t)\}|_{x=0,\pi }=0,`$ (31)
$`\{\varphi (x,t),\mathrm{\Pi }(x^{},t)\}|_{x=0,\pi }=\delta (xx^{})|_{x=0,\pi }=0.`$ (32)
This structure agree with the results of the exact solution, see appendix.

### 3.2  Scalar Field with Neumann boundary conditions

For the Neumann boundary condition (20), we have the primary constraint

$$_x\varphi (x,t)0,$$
(33)
that we can rewrite in the form

$$\mathrm{\Phi }_N^{(1)}(N)=_0^\pi 𝑑xN(x)_x\varphi (x,t)=0.$$
(34)
Now the total Hamiltonian is given by

$$H_T=H_c+_0^\pi 𝑑x\lambda (x)\varphi (x,t).$$
(35)
From the evolution of the primary constraint $`\mathrm{\Phi }_N^{(1)}(N)`$ we obtain the secondary constraint,

$$\mathrm{\Phi }_N^{(2)}(M)=_0^\pi 𝑑xM(x)_x\mathrm{\Pi }(x,t)0.$$
(36)
From the stabilization of this constraint results that we don’t have more constraints and the Lagrange multiplier is,

$$\lambda (x)=0.$$

Computing the Poisson bracket between the constraints

$$\{\mathrm{\Phi }_N^{(1)}(N),\mathrm{\Phi }_N^{(2)}(M)\}=_0^\pi 𝑑xM(x)_x^2N(x)0,$$

we see that are second class constraints, then we construct the matrix (15) that in this case have the form,

$`C_{\alpha \beta }(x,x^{})=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)_x_x^{}\delta (xx^{})\mathrm{and}C^{\alpha \beta }=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)F(x,x^{}),`$ (41)
where $`F(x,x^{})`$ is a function with compact support that satisfies

$`_x^{}^2F(x,x^{})=\delta (xx^{}).`$ (42)
The solution to this equation is

$`F(x,x^{})={\displaystyle \underset{n1}{}}{\displaystyle \frac{1}{n^2\pi }}\mathrm{sin}(nx)\mathrm{sin}(nx^{}).`$ (43)
With the inverse matrix (41), we get the Dirac brackets,

$`\{\varphi (x,t),\varphi (x^{},t)\}^{}=\{\mathrm{\Pi }(x,t),\mathrm{\Pi }(x^{},t)\}^{}=0,`$ (44)
$`\{\varphi (x,t),\mathrm{\Pi }(x^{},t)\}^{}=\mathrm{\Delta }(xx^{})`$ (45)
$`\mathrm{with}\mathrm{\Delta }(xx^{})=\delta (xx^{})_x_x^{}F(x,x^{}).`$ (46)
To the function $`\mathrm{\Delta }(xx^{})`$ we call the Dirac’s delta transverse, because this have the property

$$_x\mathrm{\Delta }(xx^{})=_x^{}\mathrm{\Delta }(xx^{})=0.$$

This implies

$$\{_x\varphi (x,t),\mathrm{\Pi }(x^{},t)\}^{}=\{\varphi (x,t),_x^{}\mathrm{\Pi }(x^{},t)\}^{}=\{_x\varphi (x,t),_x^{}\mathrm{\Pi }(x^{},t)\}^{}=0.$$

So, the symplectic structure that we get for the scalar field with Neumann boundary conditions is given by,

$`\{\varphi (x,t),\varphi (x^{},t)\}|_{x=0,\pi }=\{\mathrm{\Pi }(x,t),\mathrm{\Pi }(x^{},t)\}|_{x=0,\pi }=0,`$ (47)
$`\{_x\varphi (x,t),\mathrm{\Pi }(x^{},t)\}|_{x=0,\pi }=_x\delta (xx^{})|_{x=0,\pi }=0.`$ (48)
Comparing (47-48) with the result obtained in the appendix from the exact solution (94\- 95). We see that both agree completely.

## 4 Bosonic string with mixed boundary conditions

Now we consider the case of the bosonic string in a constant magnetic field. The action for this string in the conformal gauge is,

$$S=\frac{1}{2}_\pi ^\pi _{t_1}^{t_2}𝑑x𝑑t[(_t\varphi _i(x,t))^2(_x\varphi _i(x,t))^2+F_{ij}_t\varphi _i(x,t)_x\varphi _j(x,t)],$$
(49)
with $`F_{ij}`$ an antisymmetric constant matrix, and $`i=1,..,N`$. The equations of motion for the system are,

$$(_t^2_x^2)\varphi _i(x,t)=0.$$
(50)
From the action (49) the canonical Hamiltonian is given by,

$$H_c=\frac{1}{2}_\pi ^\pi 𝑑x[\{\mathrm{\Pi }_i(x,t)T_{ij}_x\varphi _i(x,t)\}^2+(_x\varphi _i(x,t))^2]$$
(51)
and we choose Neumann boundary conditions,

$`(_x\varphi _i+T_{ij}_t\varphi _j)(x,t)|_{x=\pm \pi }=0\mathrm{with}T_{ij}={\displaystyle \frac{F_{ij}}{2}}.`$ (52)
Then, the primary constraints that we have are,

$$\mathrm{\Theta }_i^{(1)}(x,t)=M_{ij}_x\varphi _j(x,t)+T_{ij}\mathrm{\Pi }_j(x,t)0,\mathrm{with}M_{ij}=(IT^2)_{ji}.$$
(53)
In the smeared version we have

$$\mathrm{\Theta }_i^{(1)}[𝐍]=_\pi ^\pi 𝑑xN^i(x)\mathrm{\Theta }_i^{(1)}(x,t).$$
(54)
The Poisson brackets between these constraints are,

$$\{\mathrm{\Theta }^{(1)}[𝐍],\mathrm{\Theta }^{(1)}[𝐌]\}=0,$$
(55)
where we use that the smeared vectors $`N^i(x)`$ and $`M^i(x)`$ have compact support. From the evolution of these constraints with the total Hamiltonian, we obtain the secondary constraints,

$$\mathrm{\Theta }^{(2)}[𝐌]=_\pi ^\pi 𝑑xM^i(x)_x\mathrm{\Pi }_i(x,t)0.$$
(56)
The stabilization of these constraints, don’t imply new constraints. So, the complete set of constraints is,

$$\chi _\alpha (x,t)=(M_{ij}_x\varphi _j(x,t)+T_{ij}\mathrm{\Pi }_j(x,t),_x\mathrm{\Pi }_k(x,t))\mathrm{with}\alpha =1,\mathrm{},2N$$
(57)
The algebra algebra of these constraints is

$`\{\chi _\alpha (x,t),\chi _\beta (y,t)\}=C_{\alpha \beta }(x,y),`$ (58)
where the invertible matrix $`C_{\alpha \beta }(x,y)`$ is given by,

$`C_{\alpha \beta }(x,y)=\left(\begin{array}{cc}0& M\\ M& 0\end{array}\right)_x_y\delta (xy).`$ (61)
For the inverse matrix $`C^{\alpha \beta }(x,y)`$ we obtain,

$`C^{\alpha \beta }(x,y)=\left(\begin{array}{cc}0& M^1\\ M^1& 0\end{array}\right)F(x,y),`$ (64)
with $`F(x,y)`$ defined in (43). Using (64) we obtain the Dirac brackets,

$`\{\varphi _i(x,t),\varphi _j(y,t)\}^{}=(TM^1)_{ij}[_yF(x,y)+_xF(x,y)],`$ (65)
$`\{\mathrm{\Pi }_i(x,t),\mathrm{\Pi }_j(y,t)\}^{}=0,`$ (66)
$`\{\varphi _i(x,t),\mathrm{\Pi }_j(y,t)\}^{}=\delta _{ij}\mathrm{\Delta }(xy).`$ (67)
Notice that the Dirac bracket (65) vanishes in the boundary, but we have still an extra step in our procedure that implies to remove from the brackets the contribution that is different of zero only in the boundary. In order to obtain a Dirac bracket that is valid only in the bulk we take into account (43) and comparing with $`A(x,y)`$ defined in (102), we have

$`_yF(x,y)+_xF(x,y)=A(x,y){\displaystyle \frac{(x+y)}{2\pi }}.`$ (68)
To compare with the exact solution of the original problem we take the limit of our Dirac brackets to boundary, in this case $`xx^{}\pm \pi `$. So, we have

$`\{\varphi _i(x=\pm \pi ,t),\varphi _j(x^{}=\pm \pi ,t)\}^{}=\pm (TM^1)_{ij}.`$ (69)
In consequence using our procedure for the boundary we have the results

$`\{\varphi _i(x,t),\varphi _j(x^{},t)\}|_{x^{}=x=\pm \pi }=\pm (TM^1)_{ij},`$ (70)
$`\{\mathrm{\Pi }_i(x,t),\mathrm{\Pi }_j(x^{},t)\}|_{x^{}=x=\pm \pi }=0,`$ (71)
$`\{_x\varphi _i(x,t),\mathrm{\Pi }_j(x^{},t)\}|_{x^{}=x=\pm \pi }=\delta _{ij}_x\delta (xx^{})|_{x^{}=x=\pm \pi }=0.`$ (72)
with the first derivatives of the deltas null in the boundary. This result for the symplectic structure in the boundary agrees with  and the exact solution (100)-(101) projected in the boundary.

## 5 Klein-Gordon equation with mixed boundary conditions

In this section we analyze the case of the open string ending on a D-brane in the pp-wave background. This example has been recently studied in . Where the bosonic action for the open string in the light cone gauge is

$$S=\frac{1}{2}_\pi ^\pi _{t_1}^{t_2}𝑑x𝑑t[(_t\varphi _i(x,t))^2(_x\varphi _i(x,t))^2m^2(\varphi _i(x,t))^2+F_{ij}_t\varphi _i(x,t)_x\varphi _j(x,t)].$$
(73)
From this action the canonical Hamiltonian is

$$H_c=\frac{1}{2}_\pi ^\pi 𝑑x[(\mathrm{\Pi }_i(x,t)T_{ij}_x\varphi _i(x,t))^2+(_x\varphi _i(x,t))^2+m^2(\varphi _i(x,t))^2],$$
(74)
with $`T_{ij}`$ given in (52). The primary constraints are the Neuman boundary conditions,

$$\mathrm{\Theta }_i^{(1)}(x,t)==M_{ij}_x\varphi _j(x,t)+T_{ij}\mathrm{\Pi }_j(x,t)0,$$
(75)
with $`M_{ij}`$ given in (53). The variation in time of these constraints produce the secondary constraints

$`\mathrm{\Theta }^{(2)}[𝐌]={\displaystyle _\pi ^\pi }dxM(x)[_x\mathrm{\Pi }_i(x,t)m^2T_{ib}\varphi _b(x,t))]0.`$ (76)
These constraints in the same way that $`\mathrm{\Theta }^{(1)}[𝐍]`$ also satisfy

$$\{\mathrm{\Theta }^{(2)}[𝐍],\mathrm{\Theta }^{(2)}[𝐌]\}=0.$$

However, the matrix $`C_{\alpha \beta }(x,y)`$ given by the Poisson bracket of all constraints is invertible and have the form

$$C_{\alpha \beta }(x,y)=\left(\begin{array}{cc}0& (M_y_xm^2T^2)_{ij}\\ (M_y_xm^2T^2)_{ij}& 0\end{array}\right).$$
(77)
Then the inverse matrix $`C^{\alpha \beta }(x,y)`$ is given

$$C^{\alpha \beta }(x,y)=\left(\begin{array}{cc}0& R_{ij}(x,y)\\ R_{ij}(x,y)& 0\end{array}\right),$$
(78)
with $`R_{ij}(x,x^{})`$ a matrix with compact support, that satisfies

$`(M_x^{}^2+m^2T^2)_{ib}R_{bj}(x,x^{})=\delta _{ij}\delta (xx^{}).`$ (79)
The solution to this boundary problem is,

$`R(x,x^{})={\displaystyle \underset{n1}{}}{\displaystyle \frac{\mathrm{sin}(nx)\mathrm{sin}(nx^{})}{\pi [Mn^2m^2T^2]}}.`$ (80)
Then the Dirac brackets are,

$$\begin{array}{cc}\hfill \{\varphi _i(x,t),\varphi _j(x^{},t)\}^{}=& T_{ia}[_x^{}R_{a,j}(x,x^{})+_xR_{a,j}(x,x^{})]0,\hfill \\ \hfill \{\mathrm{\Pi }_i(x,t),\mathrm{\Pi }_j(x^{},t)\}^{}=& m^2(MT)_{ia}[_xR_{aj}(x,x^{})+_x^{}R_{aj}(x,x^{})],\hfill \\ \hfill \{\varphi _i(x,t),\mathrm{\Pi }_j(x^{},t)\}^{}=& \delta _{ij}(xx^{})+m^2T_{ib}^2R_{bj}(x,x^{})_x_x^{}R_{ib}(x,x^{})M_{bj}.\hfill \end{array}$$
(81)
Notice that, the massive term is included in all brackets. To compare with the exact solution, we see that in the term,

$$\begin{array}{cc}\hfill _xR(x,x^{})+_x^{}R_{aj}(x,x^{})=& \underset{n1}{}\frac{n\mathrm{sin}n(x+x^{})}{\pi [Mn^2m^2T^2]}\hfill \\ \hfill =& M^1\left[\underset{n1}{}\frac{1}{n\pi }\mathrm{sin}n(x+x^{})+\underset{n1}{}\frac{m^2M^1T^2\mathrm{sin}n(x+x^{})}{n\pi [n^2m^2M^1T^2]}\right].\hfill \end{array}$$
(82)
That, using $`A(x,x^{})`$ defined in (102)

$$_xR(x,x^{})+_x^{}R_{aj}(x,x^{})=M^1\left[\frac{(x+x^{})}{2\pi }+A(x,x^{})+\underset{n1}{}\frac{m^2M^1T^2\mathrm{sin}n(x+x^{})}{n\pi [n^2m^2M^1T^2]}\right].$$

Let us remember that before to take the projection to the boundary, we eliminate all terms that take values only in the boundary. Then to compare (81) with the exact solution we don’t take into account the term $`A(x,x^{})`$. Then taking the limit $`xx^{}\pm \pi `$ we get

$`[_xR(x,x^{})+_x^{}R_{aj}(x,x^{})]|_{\pm \pi }=\pm M^1.`$ (83)
In this way we arrive to the Dirac brackets projected in the boundary

$`\{\varphi _i(x=\pm \pi ,t),\varphi _j(x^{}=\pm \pi ,t)\}^{}`$ $`=`$ $`\pm (TM^1)_{ij},`$
$`\{\mathrm{\Pi }_i(x=\pm \pi ,t),\mathrm{\Pi }_j(x^{}=\pm \pi ,t)\}^{}`$ $`=`$ $`m^2T_{ij},`$ (84)
$`\{_x\varphi _i(x=\pm \pi ,t),\mathrm{\Pi }_j(x^{}=\pm \pi ,t)\}^{}`$ $`=`$ $`m^2(TM^1)_{ij}.`$
Now, we can compare with the results of , where the use the exact solution. They get for $`i=1,2`$

$`\{\varphi _i(x,t),\varphi _j(x^{},t)\}|_{x^{}=xϵ\mathrm{\Sigma }}`$ $``$ $`\pm (TM^1)_{ij},`$ (85)
$`\{\mathrm{\Pi }_i(x,t),\mathrm{\Pi }_j(x^{},t)\}|_{x^{}=xϵ\mathrm{\Sigma }}`$ $``$ $`\pm m^2T_{ij},`$ (86)
$`\{\varphi _i(x,t),\mathrm{\Pi }_j(x^{},t)\}`$ $`=`$ $`\delta _{ij}\delta (xx^{})`$ (87)
with $`_x\delta (xx^{})|_{xϵ\mathrm{\Sigma }}=0.`$ Here we see that our results agree in (85) and (86) but not for (87). However is easy to see that the result (87) is inconsistent with the boundary condition. Taking the partial derivative with respect to $`x`$ in (87), we get

$$\{_x\varphi _i(x,t),\mathrm{\Pi }_j(x^{},t)\}|_{xϵ\mathrm{\Sigma }}=\delta _{ij}_x\delta (xx^{})|_{xϵ\mathrm{\Sigma }}=0.$$

On the other hand using the boundary condition (75) in (86) we get

$$\{_x\varphi ^i(x,t),\mathrm{\Pi }^j(x^{},t)\}|_{xϵ\mathrm{\Sigma }}\pm m^2(M^1T^2)^{ij}.$$

Then, the results (86) and (87) are inconsistent, whereas (84) is fully consistent with the boundary conditions.

## 6 Conclusions

In the present paper was developed a procedure to compute the symplectic structure for a field theory with boundary. We consider that the boundary conditions can be interpreted as Dirac’s constraints and we construct a procedure that following the ordinary steps of the Dirac’s method, produce the symplectic structure on the boundary, for several examples. In our procedure, we have the problem that does not exist a definition for a smeared Dirac bracket so to compute this bracket we consider the prescription that the Dirac bracket in the interior of $`\mathrm{\Sigma }`$ is given by the ordinary Dirac bracket valued on $`\mathrm{\Sigma }`$ minus the contributions that are no null only in the boundary $`\mathrm{\Sigma }`$. Using this prescription we are capable to compute the symplectic structure in several examples always obtaining fully consistent results. In a future we will try to extend our procedure to the case of General Relativity and compare our results with the recently published papers

###### Acknowledgments.

The authors acknowledge partial support from DGAPA-UNAM grant IN117000.

## Appendix A Appendix

In this appendix, we show the results obtained by using the exact solutions.

### A.1 Scalar field

The general solution of the equation of motion for the scalar field (18) with Dirichlet boundary conditions (19) is

$`\varphi (x,t)={\displaystyle \underset{n1}{}}q_n(t)\sqrt{{\displaystyle \frac{2}{\pi }}}\mathrm{sin}(nx),`$ (88)
where $`q_n(t)`$ satisfies

$$\ddot{q}_n(t)=n^2q_n(t).$$

The inverse relation to (88) is given by

$`q_n(t)=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle _0^l}𝑑x\varphi (x,t)\mathrm{sin}(nx).`$ (89)
Using (89) in (88) we found

$$\varphi (x,t)=_0^\pi 𝑑x^{}\varphi (x^{},t)\frac{2}{\pi }\underset{n1}{}\mathrm{sin}(nx^{})\mathrm{sin}(nx)=_0^l𝑑x^{}\varphi (x^{},t)\delta (xx^{}),$$

that means that the Dirac delta is,

$$\delta (xx^{})=\frac{2}{\pi }\underset{n1}{}\mathrm{sin}(nx^{})\mathrm{sin}(nx).$$

On the other hand, the canonical momentum is,

$`\mathrm{\Pi }(x,t)={\displaystyle \frac{L}{\dot{\varphi }(x,t)}}=\dot{\varphi }(x,t)={\displaystyle \underset{n1}{}}\dot{q}_n(t)\sqrt{{\displaystyle \frac{2}{\pi }}}\mathrm{sin}(nx),`$ (90)
form this follows the notation $`\dot{q}_n(t)=p_n(t)`$, and then we have

$`\mathrm{\Pi }(x,t)=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \underset{n1}{}}p_n(t)\mathrm{sin}(nx)`$ (91)
with inverse given by

$`p_n(t)=\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle _0^\pi }𝑑x\mathrm{\Pi }(x,t)\mathrm{sin}(nx).`$ (92)
Now, considering that the Poisson brackets between the normal modes are,

$$\{q_n(t),q_m(t)\}=\{p_n(t),p_m(t)\}=0\mathrm{and}\{q_n(t),p_m(t)\}=\delta _{nm}.$$

In this way, the Poisson brackets that follows from the exact solution are,

$$\begin{array}{cc}\hfill \{\varphi (x,t),\mathrm{\Pi }(x^{},t)\}& =\frac{2}{\pi }\underset{n1}{}\mathrm{sin}(nx^{})\mathrm{sin}(nx)=\delta (xx^{})\hfill \\ \hfill \{\varphi (x,t),\varphi (x^{},t)\}& =\{\mathrm{\Pi }(x,t),\mathrm{\Pi }(x^{},t)\}=0.\hfill \end{array}$$
(93)
The projection of this result to the boundary, gives exactly the same result, that we get using our procedure, see (31) and (32).

For the Neumann boundary conditions, using the exact solution we get

$`\{\varphi (x,t),\mathrm{\Pi }(x^{},t)\}`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}+{\displaystyle \frac{2}{\pi }}{\displaystyle \underset{n1}{}}\mathrm{cos}(nx^{})\mathrm{cos}(nx)=\delta (xx^{})`$ (94)
$`\{_x\varphi (x,t),\mathrm{\Pi }(x^{},t)\}|_{x=\pi ,0}`$ $`=`$ $`0=_x\delta (xx^{})|_{x=\pi ,0}.`$ (95)
If we compare this result with our method of computation (47) and (48) we see that both result agree. So, we see that in the case of the scalar field the projection of the Dirac brackets obtained from our procedure agrees completely with the result obtained using the exact solution.

### A.2 Mixed Conditions

The general solution for the equations of motion in this case is of the form,

$$\varphi _i(x,t)=M_i+B_it+T_{ij}B_jx+\underset{n1}{}\frac{q_{ni}(t)}{\sqrt{\pi }}\mathrm{cos}(nx)\frac{T_{ij}p_{nj}(t)}{n\sqrt{\pi }}\mathrm{sin}(nx).$$

With $`\ddot{q}_i(t)=q_i(t)`$ and $`\dot{q}_i=p_i(t)`$. Furthermore, from the definition of the canonical momenta we have

$$\mathrm{\Pi }_i(x,t)=M_{ij}[B_j+\underset{n1}{}\frac{p_{nj}(t)}{\sqrt{\pi }}\mathrm{cos}(nx)].$$

Then, follows that the Fourier coefficients are,

$`q_{in}(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle _\pi ^\pi }𝑑x\varphi (x,t)\mathrm{cos}(nx),`$ (96)
$`p_{in}(t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle _\pi ^\pi }𝑑xM_{ij}^1\mathrm{\Pi }_j(x,t)\mathrm{cos}(nx),`$ (97)
$`B_i`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _\pi ^\pi }𝑑xM_{ij}^1\mathrm{\Pi }_j(x,t),`$ (98)
$`A_i`$ $`=`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _\pi ^\pi }𝑑x[\varphi _i(x,t)M_{ij}^1\mathrm{\Pi }_j(x,t)].`$ (99)
For the Poisson brackets we get,

$$\{q_{in}(t),p_{jm}(t)\}=M_{ij}^1\delta _{nm},\{A_i,B_j\}=\frac{M_{ij}^1}{2\pi },$$

and the brackets vanish for the other cases. Using, the exact solution we obtain

$$\{\varphi _i(x,t),\mathrm{\Pi }_j(x^{},t)\}=\delta _{ij}\left[\frac{1}{2\pi }+\frac{1}{\pi }\underset{n1}{}\mathrm{cos}(nx^{})\mathrm{cos}(nx)\right]=\delta _{ij}\delta (xx^{}),$$
(100)
and

$$\{\varphi _i(x,t),\varphi _j(x^{},t)\}=(M^1T)_{ij}\left[\frac{(x+x^{})}{2\pi }+\underset{n1}{}\frac{\mathrm{sin}n(x+x^{})}{n\pi }\right],$$
(101)
In particular on the boundary

$$\{\varphi _i(\pm \pi ,t),\varphi _j(\pm \pi ,t)\}=\pm (M^1T)_{ij}.$$

and we see that the fields do not commute on the boundary, whereas the canonical momenta are commutative.

On the other hand, from the consistency with the Poisson brackets we see that

$`A(x,x^{})=\left[{\displaystyle \frac{(x+x^{})}{2\pi }}+{\displaystyle \underset{n1}{}}{\displaystyle \frac{\mathrm{sin}n(x+x^{})}{n\pi }}\right],`$ (102)
must be vanish in all interval $`(\pi ,\pi )`$, except a zero measure set, in the boundary satisfies

$$A(\pm \pi ,\pm \pi )=\pm 1.$$

In fact we see that $`A`$ is a distribution given by

$$_xA(x,x^{})=_x^{}A(x,x^{})=\delta _a(xx^{})\delta _b(xx^{}).$$

With

$$\delta _a(xx^{})=\frac{1}{2\pi }+\frac{1}{\pi }\underset{n1}{}\mathrm{cos}(nx^{})\mathrm{cos}(nx),$$

and

$$\delta _b(xx^{})=\frac{1}{\pi }\underset{n1}{}\mathrm{sin}(nx^{})\mathrm{sin}(nx).$$