text:
# Coherent electron-phonon coupling and polaron-like transport in molecular wires

## I Introduction

Developments in nanofabrication, including both ‘bottom-up’ approaches and ‘top-down’ methods are focussing renewed attention on the properties of molecular-scale entities and their potential as electronic device components . One of the most basic theoretical questions that can be asked in this context is: what determines the conductance of a molecule if it is used as a current-carrying element bridging two reservoirs of differing electron chemical potential? This question now has an immediate relevance for experiments in which such conductances are measured, using either (i) scanning probe tips for individual molecules adsorbed on surfaces , for molecular wires adsorbed at step edges , embedded in self-assembled monolayers , or (ii) (macroscopic) electrodes obtained by nanolithography  or from a mechanically controllable break junction. .

Since the seminal work of Aviram and Ratner  concerning the electron transfer rate between acceptor and donor groups linked by a conjugated molecular bridge, numerous theoretical studies on electron transfer and transport through molecular systems have been performed. In the following, we briefly review some contributions on the electron transport through a single organic molecule (or a few molecules) whose ends are connected to electron reservoirs. Calculations of the electronic transmission through such systems have been done for purely one-dimensional models  and two-dimensional models . More realistic descriptions of the electrode/molecule system have also been developed. Combining elastic electron scattering theories with three-dimensional tight-binding-like Hamiltonians, models have been developed for molecular wires connected to two semi-infinite surfaces  or to two semi-infinite ‘rods’ , or to cluster-like leads where imaginary parts are introduced in the Hamiltonian to take into account the fact that electrons can leak into the metallic reservoirs . Within a framework equivalent to the latter model, calculations have been extended to the Hartree-Fock level for a molecule attached to gold clusters . More recently, density functional theory has been applied to a molecular wire (described by atomic pseudopotentials) connected to two jellium surfaces .

These theoretical studies have clarified the importance of three major points crucial for the transport properties of the molecular wire connected to the reservoirs. First, they have shown the importance of the electronic and chemical interaction between the ends of the molecular wire and the reservoirs. The larger the Hamiltonian matrix elements between the delocalized electronic states of the electron reservoirs and those molecular electronic states that extend along the wire, the better the conductance properties will be; furthermore, these matrix elements should be large compared with the characteristic intra-molecular Coulomb interaction between electrons in order to avoid the Coulomb blockade . Second, the theories show that in the limit of small applied voltage and away from the Coulomb blockade regime, the transport is dominated by charge carrier tunneling inside the HOMO-LUMO gap of the molecule. This gap is another crucial parameter for control of the conductance of the wire. The smaller the gap is, the larger the tunneling transmission will be. More generally, the gap of a molecule depends on the chemical nature and atomic structure of the system. This gap can also be modified by the electron-electron interactions or by a change of the structure of the molecule due to (i) external forces, (ii) (thermal) lattice fluctuations or (iii) electron-lattice interaction. Third, the calculations highlight the importance of the position of the molecular electronic levels with respect to the Fermi levels of the reservoirs in the presence of an applied voltage, and the related issue of where the potential drops occur inside the junction. This has been done both empirically  and by using approximate  or exact  self-consistent schemes.

Although the above theoretical work has shed light on several important physical processes for the transport in molecular wires and has matched to a certain level of accuracy with some experimental measurements, all the models previously cited are based on elastic electron scattering through the rigid lattice of the wire. However, in such highly confined electron systems, the coupling between electron and other excitations (phonons for instance) is strongly enhanced because of the size of the system and its quasi one-dimensionality. This makes the rigid lattice approximation questionable—particularly so since a one-dimensional metal is generically unstable to a Peierls transition at low temperature . In an infinite system, such a transition typically produces a semiconductor in which the states near the band extrema are very strongly coupled to distortions of the system; in a conjugated organic molecule, the corresponding phenomenon is a strong coupling of the $`\pi `$-electrons occupying the HOMO and LUMO states to the bond-alternation pattern. This coupling means that the low-lying states of a charged molecule (via which any net transport of charge through the molecule must proceed) involve an intimate coupling of electronic and lattice degrees of freedom, to produce excitations such as polarons or solitons . These coupled excitations can be thought of as conspiring to lower the energy gap locally around a charge carrier when it is introduced into the system. Such polaronic and solitonic phenomena have been studied in bulk or thin film samples of conducting polymers for decades .

The importance of this electron-lattice coupling means that the conventional manner of introducing lattice vibrations within a Landauer-type approach to conductance, as an extra broadening of the electronic levels (extra imaginary part in the corresponding Hamiltonian, see for example Refs. ), is not sufficient to describe the coherent lattice distortion due to charge injection. To our knowledge, the explicit nature of the distortion accompanying charge injection has only so far been partially addressed in two simplified limits . In the first case, a molecular wire was treated as a rigid lattice in which a static soliton-like defect is present . Although it possesses a mid-gap electronic state, this model does not permit the study of the dynamics of formation and transport of charge-induced lattice distortions. In the second case, the atoms of a conjugated molecule were assumed to respond classically to the injection of a electron wave-packet , via forces calculated from expectation values of the electron wave-packet and other electronic states. In this model, the lattice is able to respond to the injected charge, but not in the physically correct manner: within a wave-packet approach to tunneling, only a small part of the electronic charge enters the tunnel barrier. Therefore the lattice responds with probability unity to a small fraction of the charge of the injected particle, rather than responding with a small probability to the total charge of the injected particle.

The only way to overcome these limitations is to perform transport calculations in which the full dynamical correlation between charge carriers and quantum phonons is retained. We report the results of such calculations in this paper. In order to focus on this particular mechanism for charge transport, we use a simple tight-binding model of a conducting polymer (the Su-Schrieffer-Heeger (SSH) model for trans-polyacetylene ) that does not explicitly include any electron-electron interactions. However, our calculations cover a range of transport regimes that includes tunneling transport (virtual electrons) and resonant transport (where there is sufficient energy to inject real electrons into the system). In contrast to the more usual ‘phase-breaking’ approaches to the electron-phonon interaction in transport problems, we explicitly retain the phase coherence between elastic and inelastic processes. In this paper we concentrate on systems where the boundary conditions on the molecular wires force them to be semiconducting in the zero-bias limit. Within the SSH model, this correponds to molecular chains containing an even number of monomers. Calculations on odd-length chains, which incorporate mobile solitonic defects with associated mid-gap states that make an additional contribution to the transport, will be reported separately.

The paper is organized as follows. In Section II, we present the multichannel scattering technique used to calculate the transport properties of the molecular wires. This section also includes a detailed analysis of different approximations used for modelling the molecular wires (involving harmonic phonons) and for reducing the computational cost of the calculations (involving reduction of the parameter space). The results obtained, in the limit of low temperatures, for the response of the molecular wires to charge carrier injection are given in Section III. We show how the lattice is distorted by the injection of a tunneling electron and how the coherent coupling between the tunneling electron and the quantum phonons affects the transmission properties through the wires. We also compare our transport results with the effect of straightforward static fluctuations in the harmonic lattice, excluding the dynamical correlation of the electrons and the phonons. Finally we summarize the most important results and propose further developments of the present study (Section IV). Additionally, in an appendix we recall briefly the methods used to get the ground state and the harmonic phonon modes from the original SSH model. We also derive in the appendix the quantum electron-phonon Hamiltonian used for the molecular wires. A brief account of part of this work has already appeared .

## II Model

We are interested in modelling the coherent electron (or hole) transport through a finite size system (the molecular wire) connected to two leads which inject or collect the charge carriers. Within the wire, the charge carriers interact with the atomic motion which originally drives the Peierls transition in the molecule.

The interaction between the ends of the molecular wire and the leads is supposed to be strong enough to permit a good overlap between the electronic states of the wire and the surface electronic wave functions of the leads. In most of the practical applications, molecular wires end in “active” chemical groups, like thiol (S-H) for example, which are known to react easily in presence of a gold surface to form chemical Au-S bonds . We therefore assume that the electron transfer rate at the molecule/lead interface is such that we can consider the electron (hole) transport as being a coherent process throughout the nanojunction, rather than a sequential incoherent two-step process.

In the coherent transport regime, a stationary state wave function scattering technique can be used to calculate the electron transfer through the molecular wire. Since we assume that the basis sets used to describe the leads and the molecular wire form a complete set, there are basically two ways to solve the scattering problem for a single incident charge carrier. The technique is reminiscent of the Löwdin transformation . If one projects out the basis set associated to the molecular wire, the problem is reformulated as a “single impurity” with on-site energies and coupling matrix elements to the leads depending on the injection energy of the charge carrier . If one choose to project out the basis set associated to the leads, one can effectively remove the leads from the problem. This technique is identical to the embedding technique where a (finite size) effective Hamiltonian describing the region of interest is obtained by introducing complex embedding potentials . The embedding potentials characterise the matching of electronic spectrum of the wire to the continuum of states of the semi-infinite leads. These potentials also depend on the charge injection energy.

As we wish to obtain the response of the molecular wire to charge injection as well as the transport properties through the junction, we choose the embedding approach to solve the electron (hole) transport in the system. For this, we use a technique which permits us to map the many-body electron/phonon problem onto a single-particle problem with many channels .

In the remain of this section, we present the basis of the many-channel scattering technique and discuss the different approximations introduced to reduce the computational cost of the calculations. Details of the construction of the quantum Hamiltonian for the molecular wire from the model originally proposed by Su, Schrieffer and Heeger (SSH)  are given in Appendix A.

Starting from this model, we have derived a quantum electron-phonon Hamiltonian:

$$H_\mathrm{w}=\underset{n}{}ϵ_nc_n^{}c_n+\underset{q}{}\mathrm{}\omega _qa_q^{}a_q+\underset{q,n,m}{}\gamma _{qnm}(a_q^{}+a_q)c_n^{}c_m,$$
(1)
where $`c_n^{}`$ creates an electron in the $`n`$th one-electron state of the molecular wire with energy $`ϵ_n`$ and $`a_q^{}`$ creates an excitation in the $`q`$th eigenmode of vibration (phonon) of the molecule with energy $`\mathrm{}\omega _q`$. The Hamiltonian Eq. (1) goes beyond the Holstein and Fröhlich model for the electron-phonon ($`e`$-ph) interaction, in the sense that the electron couples to different non-local eigenmodes of vibration, each mode having a different frequency. The electron-phonon coupling is linear in the phonon field displacement and involves electronic transitions via a general form for the $`e`$-ph coupling matrix elements $`\gamma _{qnm}`$.

The electronic eigenstates and eigenvalues are determined self-consistently with the atomic configuration for the ground state of the neutral dimerized molecular chain taken to be the reference system (App. A). From the atomic and electronic structures of the reference system, the phonon modes and frequencies are calculated within the harmonic approximation (App. A). The $`e`$-ph matrix elements are derived from the SSH model by expanding the atomic displacements induced by adding a charge onto the phonon modes of the neutral molecule (App. A). We have checked the accurary and the validity of the harmonic approximation for the phonons (see below section II B).

### A Multichannel scattering technique

We are interested in solving the problem of electron transport through a ‘nanojunction’ within a two-terminal device. We are mostly interested in the coherent regime for the electron transport, that is the regime where no random phase breaking is arbitrarily introduced between different electron scattering states. Furthermore, we wish to use a formalism that can treat (i) different transport regimes (pure tunneling, resonant tunneling, eventually ballistic transport) on an equal footing or in a transparent way, and (ii) the coupling of an electron with other degrees of freedom within the ‘nanojunction’. In principle, to study the electronic transport in such open systems, one would have to deal with the non-equilibrium Green’s functions formalism .

In this paper, we use a model which maps a many-body problem (to be accurate a one electron/many bosons problem) onto a single-particle problem with many channels . In a such model, one deals directly with the multichannel scattering states, although the formalism can be reformulated in terms of Green’s functions. This multichannel scattering technique has already been used to study electron transport through one-dimensional models of (i) double-barrier resonant tunneling junctions with electron coupled to localised single phonon mode , (ii) Holstein phonon model in presence of an electric field , (iii) mesoscopic structures  and Aharonov-Bohm rings with on site phonon coupling , (iv) tunneling barriers with the electron coupled to surface plasmon modes . More recently such a technique has also been used to study inelastic electron tunneling through small molecules in a STM tunneling barrier .

We start with the following heterojunction: a molecular wire containing $`N_a`$ atomic sites, described by the Hamiltonian $`H_\mathrm{w}`$ in Eq. (1) is connected to ideal one-dimensional right (R) and left (L) metallic leads, whose Hamiltonians are

$$H_\mathrm{R}=\underset{l=N_a+1}{\overset{+\mathrm{}}{}}ϵ_\mathrm{R}d_l^{}d_l+\beta _\mathrm{R}(d_l^{}d_{l1}+d_{l1}^{}d_l)$$
(2)
and

$$H_\mathrm{L}=\underset{l=\mathrm{}}{\overset{0}{}}ϵ_\mathrm{L}d_l^{}d_l+\beta _\mathrm{L}(d_l^{}d_{l1}+d_{l1}^{}d_l),$$
(3)
via the coupling matrices,

$$T_\mathrm{R}=v_\mathrm{R}(d_{N_a+1}^{}c_{N_a}+c_{N_a}^{}d_{N_a+1})$$
(4)
and

$$T_\mathrm{L}=v_\mathrm{L}(d_0^{}c_1+c_1^{}d_0).$$
(5)
The operators $`d_l^{}`$ ($`d_l`$) create (annihilate) an electron on site $`l`$ inside the leads, with on-site energy $`ϵ_{\mathrm{L},\mathrm{R}}`$ and nearest neighbours hopping integrals $`\beta _{\mathrm{L},\mathrm{R}}`$, the operators $`c_i^{}`$ ($`c_i`$) are the corresponding operators for an electron on site $`i`$ within the molecule, and $`v_\mathrm{R}`$ and $`v_\mathrm{L}`$ are the hopping matrix elements between the ends of the molecule and the right and left leads respectively. Within the molecular wire, the transformation from site representation to eigenstate representation is easily performed knowing that $`c_n^{}=_{i=1}^{N_a}Z_n^{}(i)c_i^{}`$, where $`Z_n(i)`$ are the components of the $`n`$th electronic eigenstate of the wire (see Appendix A).

The procedure for mapping the problem onto a single-particle system with many channels is performed by writing the total scattering wave function $`|\mathrm{\Psi }(E)`$, for the total energy $`E`$ of the electron-phonon system as

$$|\mathrm{\Psi }(E)=\underset{l}{}\underset{\{n_q\}}{}\alpha _{l,\{n_q\}}(E)|l,\{n_q\},$$
(6)
where the basis set used to expand the scattering waves is defined (in the case of electron transport) as $`|l,\{n_q\}=c_l^{}_q\frac{(a_q^{})^{n_q}}{\sqrt{n_q!}}|0`$ (for $`1lN_a`$) and $`|l,\{n_q\}=d_l^{}_q\frac{(a_q^{})^{n_q}}{\sqrt{n_q!}}|0`$ (for other values of $`l`$), $`|0`$ being the vacuum state and $`\{n_q\}`$ being the phonon occupations. The vacuum state $`|0`$ is taken to be the neutral ground state of the system, with a definite number of electrons in each of the left lead, right lead and molecule. The electronic states we consider therefore involve adding a single electron to this neutral state; the added electron may be anywhere in the system (in the left lead, the right lead, or the molecule). For hole transport we use an identical basis, except that electron creation operators are replaced by annihilation operators. By writing the wave function coefficients as in Eq.(6), no explicit separation between the electronic and phonon degrees of freedom has been assumed.

As far as the electron (hole) propagation is concerned, each different channel is associated with a different set of phonon occupation numbers $`\{n_q\}`$. The total wave function $`|\mathrm{\Psi }(E)`$ is the eigenstate of the total Hamiltonian $`H=H_\mathrm{L}+T_\mathrm{L}+H_\mathrm{w}+T_\mathrm{R}+H_\mathrm{R}`$, $`H|\mathrm{\Psi }(E)=E|\mathrm{\Psi }(E)`$, with the full scattering boundary conditions applied (i.e. an incident electron or hole with the molecule in a given vibrational state).

In absence of dissipation, the total energy $`E`$ of the system is conserved during the scattering process, i.e

$$E=E_{\mathrm{in}}+\underset{q}{}n_q\mathrm{}\omega _q=E_{\mathrm{out}}+\underset{q}{}m_q\mathrm{}\omega _q,$$
(7)
where $`E_{\mathrm{in}}`$ is the energy of the incoming electron and $`\{n_q\}`$ is the initial set of phonon occupation numbers. $`E_{\mathrm{out}}`$ is the energy of the outgoing reflected or transmitted electron with the corresponding set of phonon occupancies $`\{m_q\}`$. For inelastic scattering processes $`\{n_q\}\{m_q\}`$; for elastic scattering the phonon distribution is conserved.

The wave function coefficients $`\alpha _{l,\{m_q\}}`$ take an asympotic form inside the leads. The form corresponds to propagating Bloch waves inside the left and right leads whose amplitudes are the reflection $`r_{\{m_q\}}`$ and transmission $`t_{\{m_q\}}`$ coefficients of the electron in the different channels. For an injected electron from the left lead (for example), we have

$`\alpha _{l,\{m_q\}}=e^{\mathrm{i}k_{\{n_q\}}^\mathrm{L}l}\delta _{\{n_q\},\{m_q\}}+r_{\{m_q\}}e^{\mathrm{i}k_{\{m_q\}}^\mathrm{L}l}`$ (8)
inside the left lead ($`l1`$) and

$`\alpha _{l,\{m_q\}}=t_{\{m_q\}}\mathrm{e}^{\mathrm{i}k_{\{m_q\}}^\mathrm{R}l}`$ (9)
inside the right lead ($`lN_a+2`$). The $`k_{\{n_q\}}^{\mathrm{L},\mathrm{R}}`$ are the dimensionless wave-vectors of the Bloch waves in the different channels. For a given total energy $`E`$, the wave-vectors depend on the phonon occupation numbers. Figure 1 shows a simplified sketch of the multichannel technique for different phonon excitations inside the molecular wire.

The solutions of the Schrörindger equation $`l|H|\mathrm{\Psi }(E)=El|\mathrm{\Psi }(E)`$ inside the leads give the dispersion relations for the electron wave-vectors inside the different channels, $`E_{\mathrm{in}}=ϵ_\mathrm{L}+2\beta _\mathrm{L}\mathrm{cos}k_{\{n_q\}}^\mathrm{L}`$ for the incoming wave and $`E_{\mathrm{out}}=ϵ_\mathrm{L}+2\beta _\mathrm{L}\mathrm{cos}k_{\{m_q\}}^\mathrm{L}`$ for the outgoing reflected wave and to $`E_{\mathrm{out}}=ϵ_\mathrm{R}+2\beta _\mathrm{R}\mathrm{cos}k_{\{m_q\}}^\mathrm{R}`$ for the transmitted wave to the right lead.

In this paper, we report calculations valid in the limit of low temperatures assuming that $`k_\mathrm{B}T\mathrm{}\omega _q`$. Therefore we take for the initial set of phonon occupation numbers, the set corresponding to all phonon modes in the ground state $`\{n_q\}=\{0\}`$. The dispersion relations become simply $`E=ϵ_{\mathrm{L},\mathrm{R}}+2\beta _{\mathrm{L},\mathrm{R}}\mathrm{cos}k_{\{0\}}^{\mathrm{L},\mathrm{R}}`$ for the elastic channels and $`E=ϵ_{\mathrm{L},\mathrm{R}}+2\beta _{\mathrm{L},\mathrm{R}}\mathrm{cos}k_{\{m_q\}}^{\mathrm{L},\mathrm{R}}+_qm_q\mathrm{}\omega _q`$ for the inelastic channels. Note that in this case, the total energy $`E`$ also represents the electron injection energy.

Solving $`l|H|\mathrm{\Psi }(E)=El|\mathrm{\Psi }(E)`$ at the interfaces of the heterojunction is the next step to perform in order to resolve the unknown wave function coefficients inside the molecular wire. For $`l=0`$ and $`l=N_a+1`$, the relations between the reflection and transmission coefficients with the wave function coefficients at the ends of the wire are obtained:

$$t_{\{m_q\}}=\frac{v_\mathrm{R}}{\beta _\mathrm{R}}\alpha _{N_a,\{m_q\}}e^{\mathrm{i}k_{\{m_q\}}^\mathrm{R}N_a},$$
(10)
and

$$r_{\{m_q\}}=\delta _{\{0\},\{m_q\}}+\frac{v_\mathrm{L}}{\beta _\mathrm{L}}\alpha _{1,\{m_q\}},$$
(11)
for the elastic $`\{0\}`$ and inelastic $`\{m_q\}\{0\}`$ channels.

Finally, solving the Schrödinger equation on sites $`l[1,N_a]`$ permits one to effectively remove the leads by introducing complex embedding potentials . Then the solution of the full scattering problem is obtained by solving the following complex linear system

$$\left[EH_\mathrm{w}\mathrm{\Sigma }^\mathrm{L}(E)\mathrm{\Sigma }^\mathrm{R}(E)\right]|\alpha (E)=|s(E),$$
(12)
where $`H_\mathrm{w}`$ is the molecular wire Hamiltonian Eq. (1), the components of $`|\alpha `$ are the scattering wave function coefficients inside the molecular wire expressed in the molecule eigenstate representation, i.e. $`\alpha _{n,\{n_q\}}=_iZ_n(i)\alpha _{i,\{n_q\}}`$, and $`\mathrm{\Sigma }^{\mathrm{L},\mathrm{R}}`$ are the embedding potentials due to the left and right lead respectively. The embedding potentials are diagonal matrices in the $`|n,\{n_q\}`$ basis set with components

$$\mathrm{\Sigma }_{n,\{n_q\}}^\mathrm{L}(E)=Z_n(1)v_\mathrm{L}g_{\{n_q\}}^\mathrm{L}(E)v_\mathrm{L}Z_n(1),$$
(13)
and

$$\mathrm{\Sigma }_{n,\{n_q\}}^\mathrm{R}(E)=Z_n(N_a)v_\mathrm{R}g_{\{n_q\}}^\mathrm{R}(E)v_\mathrm{R}Z_n(N_a),$$
(14)
where $`g_{\{n_q\}}^{\mathrm{L},\mathrm{R}}`$ are the surface Green’s function of the isolated left and right leads for the different channels given by

$$g_{\{n_q\}}^{\mathrm{L},\mathrm{R}}(E)=\mathrm{exp}(\mathrm{i}k_{\{n_q\}}^{\mathrm{L},\mathrm{R}}(E))/\beta _{\mathrm{L},\mathrm{R}}.$$
(15)
Finally $`|s(E)`$ represents the source term (i.e. the injected electron or hole at energy $`E`$) with components given by

$$s_{n,\{n_q\}}(E)=\delta _{\{0\},\{n_q\}}(2\mathrm{i}v_\mathrm{L}\mathrm{sin}k_{\{0\}}^\mathrm{L})Z_n(1).$$
(16)
In the present case, the boundary condition chosen for the source term corresponds to the injection of the charge, into the elastic channel $`\{0\}`$, from the left lead towards the right lead.

The solution of Eq. (12) can be obtained by several means using algorithms for sparse matrices. In the present work, we explicitly separate the real and imaginary parts of vectors and matrices as follows

$`\left[\begin{array}{cc}EH_\mathrm{w}\mathrm{}e\mathrm{\Sigma }(E)& \mathrm{}m\mathrm{\Sigma }(E)\\ \mathrm{}m\mathrm{\Sigma }(E)& E+H_\mathrm{w}+\mathrm{}e\mathrm{\Sigma }(E)\end{array}\right]\left[\begin{array}{c}\mathrm{}e\alpha (E)\\ \mathrm{}m\alpha (E)\end{array}\right]=\left[\begin{array}{c}\mathrm{}es(E)\\ \mathrm{}ms(E)\end{array}\right],`$ (23)
where $`\mathrm{}e`$ and $`\mathrm{}m`$ denote the real and imaginary parts of the different quantitites and $`\mathrm{\Sigma }=\mathrm{\Sigma }^\mathrm{L}+\mathrm{\Sigma }^\mathrm{R}`$. The matrix in the left-hand-side of Eq. (23) is therefore real and symmetric. We then use standard conjugate gradient (CG) technique to solve the linear system $`A|x=|b`$ where $`A`$ is a real and symmetric matrix . An efficient algorithm has been devised for computing the products $`H_\mathrm{w}|x_i`$ generated during the iterative CG steps. It is based on an optimal adressing of the vector components and uses the selection rules for the $`\gamma _{qnm}`$ matrix elements (see Appendix A 3).

### B Isolated molecular wire : harmonic phonons and reduced parameter space

We have shown that the solution of the scattering problem is obtained by solving Eq. (12) for the value of the wave functions inside the molecular wire. Assuming a truncated phonon space up to a finite number of excitations $`n_{\mathrm{occ}}^{\mathrm{max}}`$, the size of the basis set is given by $`N_{\mathrm{size}}=N_e\times (n_{\mathrm{occ}}^{\mathrm{max}}+1)^{N_{\mathrm{ph}}}`$ with $`N_e`$ ($`N_{\mathrm{ph}}`$) being the number of electronic states (phonon modes). Even for relatively short wires, the size $`N_{\mathrm{size}}`$ of the basis set quickly becomes too large for tractable numerical calculations and/or reasonable computing times. For example for $`N_a=20`$ atomic sites ($`N_{\mathrm{ph}}=19`$ acoustic and optic phonon modes) with only $`n_{\mathrm{occ}}^{\mathrm{max}}=2`$, we obtain $`N_{\mathrm{size}}10^6`$. In the following we show how to reduce the basis set size and prove the validity of the approximations introduced. The corresponding Hilbert space can be reduced by (i) considering only valence (conduction) band electronic states (this will correspond to hole (electron) transport respectively), by (ii) considering a limited but sufficient set of phonon modes and finally by (iii) using only a few excitations in each phonon modes (truncated harmonic oscillator approximation).

In order to determine which phonon modes are mainly contributing to the charge-induced deformation of the chain, we calculate the ground state atomic configurations for a neutral chain $`u_i^0`$ and for a chain charged with one additional electron $`u_i^\mathrm{c}`$. Note that because of the charge-conjugation symmetry (which holds exactly for the SSH model) adding an extra electron or removing an electron (i.e. adding a hole) produces exactly the same lattice distortion. Such a lattice distortion is known as an electron (or hole) polaron (c.f. Fig. 3). The lattice distortion is then projected onto the harmonic eigenmodes of vibration $`V_q`$ of the neutral chain and the corresponding Huang-Rhys factors $`S_q`$ are determined by

$$S_q=\frac{\frac{1}{2}M\omega _q^2\mathrm{\Delta }_q^2}{\mathrm{}\omega _q},$$
(24)
where $`\mathrm{\Delta }_q=_iV_q(i)(u_i^\mathrm{c}u_i^0)`$. The Huang-Rhys factors give the averaged number of quantum phonons that would be needed to achieve the corresponding elastic energy of the lattice distortion. We have checked that the largest values for those factors are obtained for the lowest energy (longest wave-length) optical phonon modes of the molecular wires . In particular, the optical phonon mode having the lowest energy also has the most important contribution (i.e. $`S_q>0.8`$ as can be seen for different chain lengths in Table I).

From these results, we can already reduce the number of phonon modes necessary by considering only the longest wave-length optical modes to describe the lattice deformation induced by adding a charge into the chain. As typical set of these modes is shown on Figure 2.

Furthermore, we can also check the validity of the harmonic approximation used to determine the phonon modes of the wire. Note that in section A 2, the elastic energy is expanded up to second order for small displacements around the equilibrium atomic positions in order to obtain the dynamical matrix from which the eigenmodes of vibration are determined. The deformation of the lattice due to charge addition is not purely harmonic because of the distance dependence in the electron hopping matrix elements, cf. Eq. (A1). However there are good reasons to believe that the potential surface can be fairly well described by the harmonic approximation to the neutral chain. To confirm this, we compare the deformation energy $`\mathrm{\Delta }E`$ with the harmonic distortion energy $`E_{\mathrm{harm}}=_q\frac{1}{2}M\omega _q^2\mathrm{\Delta }_q^2=_q\mathrm{}\omega _qS_q`$. The deformation energy is obtained from $`\mathrm{\Delta }E=E_0[u_i^\mathrm{c}]E_0[u_i^0]`$ where $`E_0[u_i^0]`$ is the selfconsistent total energy of the SSH Hamiltonian for a neutral chain with the corresponding atomic positions $`u_i^0`$. $`E_0[u_i^\mathrm{c}]`$ is the non-selfconsistent total energy of the neutral chain where the atomic positions $`u_i^\mathrm{c}`$ are taken to be those of the charged chain. In $`E_0[u_i^\mathrm{c}]`$ the effects of the imposed lattice distortions on the electronic spectrum are taken into account. Typical values of $`\mathrm{\Delta }E`$ and $`E_{\mathrm{harm}}`$ are given in Table I for different wire lengths. For short chains, the values of $`\mathrm{\Delta }E`$ and $`E_{\mathrm{harm}}`$ are almost indentical (less than $``$ 5% difference). For longer chains, the difference between $`\mathrm{\Delta }E`$ and $`E_{\mathrm{harm}}`$ increases but never exceed $``$ 15%. We also show below that the harmonic expansion of the elastic energy is sufficient to describe the formation of a static polaron by adding a permanent extra charge in the molecular chains.

Now we turn on the reduction of the Hilbert space in relation to the electron states. We want to check the validity of using only half the electronic spectum on the values of the atomic displacements induced by adding a charge (electron or hole) into the molecular chain. As mentioned in Ref. , we consider the action of the model Hamiltonian Eq. (1) on the ($`N_a\pm 1`$)-electron Hilbert space obtained by adding a charge into the neutral chain. It is thought that it is sufficient to consider only the ($`N_a\pm 1`$)-electron because for molecular wires strongly (electronically) coupled to the leads, the mean time between charge passages is $`10^7`$ s (for a corresponding current of 1 pA), a value orders of magnitude bigger than a typical residence time ($`10^{15}`$ s). Then we project out the addition of a charge into the electronic states by working with the electronic eigenstates representation, i.e. the sums over the eigenstates will only include the occupied valence band states when adding a hole (and later for hole transport) and only the empty conduction band states when adding an electron (for electron transport). In the following, we show that the static lattice distortions due to adding a charge are well reproduced by considering only one half of the electronic spectrum.

We start from the ground state of the reference system (neutral chain with atomic positions $`u_i^0`$). The lattice distortions due to charging are expanded onto the harmonic phonon modes $`V_q(i)`$ as $`u_i=u_i^0+_qV_q(i)\mathrm{\Delta }_q`$. Then the influence of the lattice distortions on the electronic Hamiltonian $`H^0ϵ_n\delta _{nm}`$ are taken into account by introducing the corresponding $`e`$-ph coupling off-diagonal elements $`\gamma _{qnm}`$ in $`H^0`$. The electronic Hamiltonian of the isolated molecule is

$$H_{nm}^{\mathrm{el}}=H_{nm}^0+H_{nm}^{e\mathrm{ph}}=ϵ_n\delta _{nm}+\underset{q}{}\stackrel{~}{\mathrm{\Delta }}_q\gamma _{qnm},$$
(25)
where the matrix elements $`\gamma _{qnm}`$ are given by Eq. (A9) in Appendix A, and $`\stackrel{~}{\mathrm{\Delta }}_q`$ is the dimensionless displacement $`\mathrm{\Delta }_q=\stackrel{~}{\mathrm{\Delta }}_q\sqrt{\mathrm{}/(2M\omega _q)}`$. The total energy of the distorted lattice is

$$E_{\mathrm{half}}(\{\mathrm{\Delta }_q\})=\underset{q}{}\frac{1}{2}M\omega _q^2\mathrm{\Delta }_q^2+\mathrm{Tr}[\rho ^{\mathrm{el}}H^{\mathrm{el}}(\{\mathrm{\Delta }_q\})],$$
(26)
where $`\rho ^{\mathrm{el}}`$ is the electronic density operator. To find the corresponding ground state, $`E_{\mathrm{half}}(\{\mathrm{\Delta }_q\})`$ has to be minimised versus the classical lattice phonon displacements $`\{\mathrm{\Delta }_q\}`$.

Within the half spectrum approximation, the trace in Eq. (26) runs only over the conduction (or valence) band eigenstates when we consider the wire charged by an extra electron (or hole). Adding a charge to the system involves only the LUMO or HOMO electronic state for an electron or hole respectively, therefore the functional to be minimized is actually $`E_{\mathrm{half}}(\{\mathrm{\Delta }_q\})=_q\frac{1}{2}M\omega _q^2\mathrm{\Delta }_q^2+\lambda _e(\{\mathrm{\Delta }_q\})`$ where $`\lambda _e=\phi _e|H^{\mathrm{el}}(\{\mathrm{\Delta }_q\})|\phi _e`$ is the lowest (or highest) energy eigenvalue of the electronic half-space Hamiltonian $`H^{\mathrm{el}}(\{\mathrm{\Delta }_q\})`$ and corresponds to the LUMO (or HOMO) as modified by the atomic distortion. The minimisation of $`E_{\mathrm{half}}(\{\mathrm{\Delta }_q\})`$ is obtained when the forces $`F_q=M\omega _q^2\mathrm{\Delta }_q+\phi _e|\widehat{\gamma }_q|\phi _e\sqrt{(2M\omega _q)/\mathrm{}}`$ are zero for all the phonon modes $`q`$, $`\widehat{\gamma }_q`$ being the $`e`$-ph coupling matrix with components $`\gamma _{qnm}`$. With this procedure, we can also study the contribution of the different phonon modes $`q`$ and check the validity of using only a limited number of optical phonon modes to create the distortions as we have already suggested from the values of the Huang-Rhys factors.

Figure 3 shows the lattice distortions calculated within different approximation from the ground state of a charged chain. The lattice distortions are best represented by the staggered difference $`d_i`$ between adjacent bond lengths. The quantity $`d_i=(1)^i(u_{i+1}2u_i+u_{i1})`$ is known as the dimerization. A constant dimerization pattern indicates a perfect bond length alternation in the chain, while a decrease of the dimerization indicates a deformation of the bond lengths (i.e. an increase of the short bonds and a decrease of the long bonds). These lattice distortions are localized around the charge added and are characteristic of the polaron defect in the molecular chain. The general shape of the dimerization pattern in Fig. 3 indicates that a static polaron has formed in the chain and it extents over several atomic sites.

The dimerization patterns obtained from the original SSH model are shown on Figure 3 with the dimerizations patterns calculated by considering half the electronic spectrum and classical phonon modes (all the modes and the reduced set of long wave-length optical modes). We can see that the lattice distortions representing a static polaron defect in the chain are very well reproduced by considering only the long wave-length optical phonons (for example those shown on Fig. 2 for the $`N_a=100`$). The discrepancies between the dimerization amplitudes obtained for the full and half electronic spectrum are more important for long chains than for short chains. For the short chains, the difference in dimerization amplitude does not exceed $``$ 10%. In the extreme case of the shortest two-atom chain, working with half the electronic spectrum gives the exact results. For the long chains, the difference in dimerization amplitude increases and is $``$ 17% for the $`N_a=100`$ chain length. We attribute the origin of these differences to the fact that the gap of the molecular wires decreases with the chain length. However, as we show in the next section, an asymptotic regime is reached for chain length $`N_a100`$, where the gap becomes independent of the chain length. We therefore can assume that for longer chains, the difference in dimerization amplitude should not increase further.

We also calculated the corresponding lattice distortions using quantum phonons. The calculations were done by determining the ground state of the fully quantum electron/phonon Hamiltonian Eq. (1) where we introduced a cut-off for the number of possible phonon excitations. Each harmonic quantum phonon mode can contain only up to $`n_{\mathrm{occ}}^{\mathrm{max}}`$ quanta. Note again, that the calculations were performed with the $`n,m`$ sums running only over half of the electronic spectrum. Then, for the sums running over the originally empty conductance band states, the ground state $`|\mathrm{\Psi }_0`$ corresponds to the situation where an extra electron has been added to the chain. The equivalent situation corresponding to removing an electron is obtained by summing the electronic eigenstates only over the originally occupied valence band states and considering among these eigenstates of Eq. (1) that with the highest eigenvalue.

In practice, here we choose to calculate the situation corresponding to a chain charged with one extra electron. Due to charge-conjugation symmetry, the results for hole injection will be identical. Once the ground state $`|\mathrm{\Psi }_0`$ of Eq. (1) is obtained, we can calculate the quantum average $`\delta _q=\mathrm{\Psi }_0|\delta _q|\mathrm{\Psi }_0`$ for the mean displacement of the phonon mode $`q`$, where $`\delta _q`$ is given by Eq. (A7). The atomic displacements $`u_i`$ induced by charging are obtained from $`u_i=_qV_q(i)\delta _q`$. The resulting dimerization patterns are shown on Figure 3. Convergence of the results is obtained for a small and finite number of quanta in each mode, roughly $`n_{\mathrm{occ}}^{\mathrm{max}}4,5,6`$. We will show in the next section that the results for the transport properties of the molecular wire coupled to the electrodes converge faster with respect to the values of $`n_{\mathrm{occ}}^{\mathrm{max}}`$, especially when one considers charge transport in the tunneling regime. As expected, the dimerization patterns are very close to those obtained from the classical phonon model with half electronic spectrum. The slight differences may be due to quantum delocalization of the eigenstate of Eq. (1).

Finally, it can be noticed that the general shape of the dimerization induced by charging and the spatial extent of the corresponding polaron are well reproduced by the different approximations (half electronic spectrum, limited set of optical phonon modes, finite number of phonon excitations) introduced to reduce the Hilbert space. The reduction of the parameter space is therefore entirely justified. This reduction permits us to treat a great range of molecular wire lengths and to determine the transport properties of the wires presented in the next section, with reasonable computing times.

## III Results

In this section, we present results for the electron transport through typical heterojunctions made of a molecular wire connected to two electron reservoirs. In our calculations, the dynamical (i.e. energy-dependent) correlation between the electron and the phonon degrees of freedom is kept. Such quantum coherence between the electron and the lattice is essential to treat the propagation of polarons through the wire in the tunneling regime (i.e. for an injection energy $`E`$ inside the gap of the molecular wire). We use the term coherence because there is no random dephasing of the wave functions introduced by the electron-phonon interaction. Instead each wave vector $`|j,\{n_q\}`$ has both a definite amplitude and phase, which are both dependent on the energy $`E`$ .

### A Lattice distortions induced by a tunneling electron

In order to analyse the response of the molecular wire lattice to a tunneling electron in the stationary state, we calculate the expectation value of a correlation function between the phonon field displacements and the electron density. Such a correlation function is defined as the following quantum average

$$\delta _q^{[i]}=\frac{P_i\sqrt{\frac{\mathrm{}}{2M\omega _q}}(a_q+a_q^{})P_i}{P_i},$$
(27)
where $`P_i=c_i^{}c_i`$ is the electron wave function projector onto atomic site $`i`$ (the electron density operator on site $`i`$). The quantum average of any electron and/or phonon operator $`𝒪`$ is given by $`𝒪=𝒪(E)=\mathrm{\Psi }(E)|𝒪|\mathrm{\Psi }(E)`$. The correlation function $`\delta _q^{[i]}`$ represents the mean displacement of the phonon mode $`q`$ when the electron is on site $`i`$. We can then define the conditionally averaged atomic displacement $`x_j^{[i]}=_qV_q(j)\delta _q^{[i]}`$ representing the atomic displacement on site $`j`$ on the condition that the electron is on site $`i`$. From these atomic displacements, we can obtain the dimerization pattern $`d_j^{[i]}=(1)^j(x_{j+1}^{[i]}2x_j^{[i]}+x_{j1}^{[i]})`$.

We have calculated the response of the molecular lattice to the tunneling electron for different wire lengths and different injection energies $`E`$. For example, Figure 4 shows a two-dimensional map (atomic position $`j`$/electron position $`i`$) of the dimerization $`d_j^{[i]}`$ for a $`N_a=100`$ chain length. The tunneling electron is injected at mid-gap ($`E=0.0`$) from the left ($`j=1`$) and propagates to the right ($`j=100`$). The bright part around the first diagonal in Fig. 4 represents a dip in the dimerization pattern along the atomic positions. As shown in the previous section, such a lattice distortion is the signature of the formation of a polaron inside the chain. We name it a virtual polaron because of the transient nature of the tunneling electron injected inside the molecular wire gap. The induced atomic distortions are located around the electron position $`i`$ (vertical axis). As can be seen in Figs. 4 and 5, the virtual polaron has its own intrinsinc width (estimated around $``$ 15 atomic sites for the present model of molecular chains which is consistent with the correlation length $`\xi `$ of the continuum model ). Its formation is therefore not possible for the electron positions at the ends of the molecular wire. For short chains ($`N_a<15`$), the lattice distortion associated with the polaron cannot be accommodated in the wire. This has an important consequence on the electronic transport properties, as we shall see in the next section.

The atomic displacements are, for most electron positions, slightly less than those for an isolated chain (Fig. 5) because the lattice does not respond fully to the tunneling electron as in the case of the static charge added into the chain. The width of the lattice distortion around the electron is also smaller for the tunneling electron than for the static charge.

The amplitude of the virtual polaron slightly increases when the injection energy $`E`$ increases above mid-gap. Such a behaviour persists until $`E`$ reaches the first resonance peak in the electron transmission (see next section). This may be understood from the fact that far from mid-gap, the electron wave function gets more weight (the amplitude of the electron wave function gets larger) leading to a stronger coupling to the phonons and therefore to more important lattice distortions. Above the first resonance peak, the electron is not, strictly speaking, tunneling anymore through the molecular gap. For these energies, the transport regime is better described as tunneling resonantly through the coupled $`e`$-ph states of the system. The scattering wave function acquires more weight from these quasi-standing waves which in turn modify qualitatively the polaron-like nature of the corresponding atomic distortions. A detailed analysis of such distortions is out of the scope of this paper and will be presented elsewhere .

The coherent electron-lattice distortion leads to a modification of the electronic spectrum of the molecular wire compared to the spectrum of the undistorted chain. In all the cases studied, the polaron formation is associated with a reduction of the gap of the originally undistorted chain. Because the atomic distortions induced in the case of a static electron added in the chain are different from those obtained for a tunneling electron, we expect the electron transport properties through the corresponding spectra to be different. The behaviour of the transport properties in the molecular wires is shown and analysed in terms of electron transmission probabilities in the next section.

### B Transport properties

The current flowing through the different channels (for example in the right outgoing channels) is given by

$$j_{\{m_q\}}^\mathrm{R}(E)=\frac{2e}{h}\mathrm{}m\left(\alpha _{l,\{m_q\}}^{}\beta _\mathrm{R}\alpha _{l+1,\{m_q\}}\right),$$
(28)
for sites $`l>N_a+1`$ located inside the right lead. From the asymptotic form of the wave functions inside the leads Eqs. (8) and (9), it is found that $`j_{\{m_q\}}^\mathrm{R}`$ is related to the transmission probability by $`j_{\{m_q\}}^\mathrm{R}=2\frac{e}{h}\beta _\mathrm{R}\mathrm{sin}k_{\{m_q\}}^\mathrm{R}|t_{\{m_q\}}(E)|^2`$. Similarly, the current inside the left lead channels $`j_{\{m_q\}}^\mathrm{L}`$ is related to the reflection probability $`|r_{\{m_q\}}(E)|^2`$.

We can define an effective total transmission probability as

$$T(E)=\underset{\{m_q\}}{}|t_{\{m_q\}}(E)|^2\frac{\beta _\mathrm{R}\mathrm{sin}k_{\{m_q\}}^\mathrm{R}}{\beta _\mathrm{L}\mathrm{sin}k_{\{0\}}^\mathrm{L}}.$$
(29)
This takes the form of a sum of contributions from the different outgoing channels .

In the following, we present results for the injection of an electron and transport by tunneling effect inside the gap and by resonant tunneling through the levels of the system above the gap. The electron injection energy is defined as positive with respect to the reference energy $`E=0`$ for all the molecular wires. The Fermi energies of the leads are assumed to be pinned at mid-gap in absence of any applied bias.

Figure 6 shows typical results for the effective total transmission $`T(E)`$ and the contribution from the elastic channel and from one inelastic channel. The different transmission probabilities have, as expected, an exponential behaviour in the tunneling regime (below the first resonance peak) and present resonances through the spectrum of the system for higher energies.

It is interesting to note that in the tunneling regime when the injection energies are below the first resonance peak ($`E<0.51`$ eV for $`N_a=40`$ and $`E<0.39`$ eV for $`N_a=100`$), the most important contribution to the total transmission $`T(E)`$ comes from the elastic channel. Above the first resonance, the contribution of the other inelastic channels may become as important as the contribution of the elastic channel. Although the contribution of the elastic process is dominant in the tunneling regime, it should be noted that this contribution is quite different to that obtained from an elastic scattering treatment of the transport through the undistorted molecule. This is because the ‘elastic channel’ corresponds to processes in which, once the electron has crossed the molecule, there is no overall absorption or emission of phonons. Nevertheless, phonons are emitted and absorbed during the intermediate stages of the transport, and therefore this ‘elastic channel’ is quite different from a rigid-molecule calculation where no coupling to the phonons is present. We have already shown in Ref.  that the transmission is actually enhanced (in the tunneling regime and in the limit of low temperatures) because of the $`e`$-ph coupling. The transmission enhancement comes from an effective gap reduction associated with the virtual polaron formation. This gap reduction is also different from the one obtained by charging a (classical lattice) chain with a static electron.

Figure 7 represents half the HOMO-LUMO gap of isolated neutral and charged chains obtained from the SSH model and the effective gap of our coupled quantum electron-phonon model for different chain lengths. Half the gap corresponds to the value of the lowest unoccupied molecular states (LUMO) of the isolated (neutral or charged) chain. The effective gap of the molecular wire (including $`e`$-ph coupling) connected to the leads is obtained from the position of the first resonance peak in the transmission $`T(E)`$. We also plotted the position of the first resonance peak obtained from purely elastic calculations (ignoring the $`e`$-ph coupling). As expected, the positions of these resonances reproduce the values of half the gap of the corresponding neutral chains. An energy shift appears in the first resonance position because of the real part of the embedding potentials $`\mathrm{\Sigma }^\mathrm{L}`$ and $`\mathrm{\Sigma }^\mathrm{R}`$ that appear in the solution of the scattering problem Eq. (12).

For the fully inelastic scattering calculations, the first resonance in $`T(E)`$ occurs for injection energy $`E`$ smaller than those obtained from the (solely) elastic calculations. The behaviour characterises the effective gap reduction of the wire due to the $`e`$-ph interaction. The energy position of the first resonance in $`T(E)`$ is close to the charging energy $`E^{\mathrm{charg}}`$ of the isolated molecular wire. The charging energy $`E^{\mathrm{charg}}`$ is defined as the difference between the (selfconsistent) ground state total energy of the charged chain $`E_{+1}[u_i^\mathrm{c}]`$ (with the distorted lattice positions $`u_i^\mathrm{c}`$ corresponding to a static polaron) and the ground state total energy of the neutral chain $`E_0[u_i^0]`$ (with the undistorted, perfectly dimerized, lattice positions $`u_i^0`$). Some values of $`E^{\mathrm{charg}}`$ for different molecular wire lengths are given in Table I. Differences between the values of $`E^{\mathrm{charg}}`$ and the first resonance in $`T(E)`$ do occur, they are due to (i) the systematic energy shift introduced by the embedding potentials and also to (ii) the differences arising from treating the lattice distortions classically or with quantum phonons as pointed out in section III A.

To summarise, there are two distinct physical processes that affect the value of the band-gap of the molecular wire and therefore the transport properties through this wire. Firstly, the intrinsic band-gap dimishes with increasing wire lengths. The values reach a asymptotic regime for a length $`N_a\stackrel{>}{}100`$. The asymptotic band-gap value and the length above which the asympotic regime is obtained depend of the chemical nature of the molecule, i.e. on the SSH parameters used to model the molecular chain (cf. Appendix A). Secondly, an effective band-gap reduction occurs upon charging the molecular wire with a static charge or with a transient tunneling charge. The reduction is due to the the coupling between the charge and the lattice leading to the formation of static or virtual polaron respectively. Therefore, the electron transmission through the molecular wire increases because of the $`e`$-ph coupling for injection energies inside the gap. Although, the trends for the gap reduction are similar for the static and virtual polarons, they differ quantitatively over the whole range of molecular lengths studied here.

Finally, it should be noticed that for short wires ($`N_a10`$) the effective gap obtained from the fully inelastic calculations converges towards the gap obtained from elastic calculations. Therefore the values of the electron transmission are almost identical for both inelastic and purely elastic calculations. As mentioned in section III A, although the polaron cannot be accomodated in very short wires, in these conditions (short tunneling length and large gap), the tunneling process itself is too fast to get significant lattice distortions associated with the charge injection. In order to illustrate this point, we give in Table II the characteristic times associated with the electron and the phonons. We calculate the time domain $`\tau _{\mathrm{ph}}`$ of the phonons from the extremal phonon frequencies used in our calculations. We estimate a residence time $`\tau _{\mathrm{res}}`$ for the electron from the full width at half maximum $`\eta `$ of the first resonance peak using $`\tau _{\mathrm{res}}=\mathrm{}/2\eta `$. Following Büttiker and Landauer , a traversal time $`\tau _{\mathrm{BL}}`$ for the tunneling electron can be obtained; it is calculated from the elastic channel transmission coefficient as done in Ref. . The traversal time $`\tau _{\mathrm{BL}}`$ is by definition dependent on the electron injection energy $`E`$. We give in Table II the range of $`\tau _{\mathrm{BL}}`$ corresponding to energies inside the tunneling gap. We can see that for short wires ($`N_a2\xi `$), the tunneling time (and also $`\tau _{\mathrm{res}}`$) is smaller than the characteristic times associated with the phonons. The lattice dynamics is therefore not fast enough to respond to the tunneling electron. In this regime, the transport properties obtained from purely elastic scattering (rigid lattice) are not strongly different from those obtained by inelastic scattering (coherently distorted lattice by the $`e`$-ph coupling). For longer wires ($`N_a402\xi `$), the tunneling time (and $`\tau _{\mathrm{res}}`$) is comparable to or larger than the characteristic phonon times; the polaron can be formed inside the molecular wire. It is in this regime that we observe the most important differences between the electron transmission for a rigid lattice, and for a lattice that can be deformed by the tunneling electron.

### C Lattice fluctuations

The importance of lattice fluctuations on the electronic structure of conjugated molecules has already been considered . As the lattice fluctuations (even in the limit of zero temperature, i.e. the zero-point motion) are of the same order of magnitude than the distortions induced by charge injection, it is important to know if such species (polarons, solitons) survive the lattice fluctuations.

We present here results for the electron transmission through a molecular chain where disorder is introduced due to lattice fluctuations. In order to compare these results with the full quantum inelastic transmission, we consider the limit of low temperatures where the lattice fluctuations are due to the zero-point motion of each phonon mode. The averaged transmission is obtained by summing up the different elastic transmission probabilities $`T(E;\{\mathrm{\Delta }_q\})`$ associated with the displaced phonon configurations $`\{\mathrm{\Delta }_q\}`$, weighted by the Gaussian distribution probability of the ground state harmonic oscillator wave function. The average of a function $`f(\{\mathrm{\Delta }_q\})`$ depending on the phonon displacements $`\{\mathrm{\Delta }_q\}`$ is calculated as

$`f_{\mathrm{dis}}={\displaystyle 𝑑\mathrm{\Delta }_1\mathrm{}𝑑\mathrm{\Delta }_q\mathrm{}𝑑\mathrm{\Delta }_{N_{\mathrm{ph}}}f(\{\mathrm{\Delta }_q\})\underset{q}{}\frac{\mathrm{exp}(\frac{1}{\mathrm{}}M\omega _q\mathrm{\Delta }_q^2)}{\sqrt{2\pi }\sigma _q}},`$ (30)
where the width of the Gaussian distribution is given by the virial theorem for the zero-point fluctuation of each mode $`\frac{1}{2}M\omega _q^2X_q^2=\frac{1}{2}E_{n=0}=\frac{1}{2}(n+\frac{1}{2})\mathrm{}\omega _q=\frac{1}{4}\mathrm{}\omega _q`$ and $`\sigma _q=X_q^2^{1/2}`$.

When $`f(\{\mathrm{\Delta }_q\})=T(E;\{\mathrm{\Delta }_q\})`$, Eq. (30) is equivalent to the static lattice approximation expression obtained by Pazy and Laikhtman  using a path-integral formalism. In our calculations, the transmission $`T(E;\{\mathrm{\Delta }_q\})`$ is calculated from elastic scattering through the one-electron spectrum of the molecular chain distorted by the zero-point lattice fluctuations. The distorted molecule eigenstates are obtained from the Hamiltonian Eq. (A1) by replacing the equilibrium atomic positions $`u_i^0`$ by $`u_i=u_i^0+\delta u_i`$ where $`\delta u_i=_qV_q(i)\mathrm{\Delta }_q`$ for the different configurations $`\{\mathrm{\Delta }_q\}`$ of the distorted lattice. In practice, $`T(E;\{\mathrm{\Delta }_q\})`$ is calculated in a similar way as described in section II A but using the approximation of elastic scattering (i.e. no energy exchange between the electron and the phonons is allowed, setting $`n_{\mathrm{occ}}^{\mathrm{max}}=0`$ for all the phonon modes considered). Futhermore, the integrals over the phonon displacements are performed using an algorithm to generate random deviates with a normal Gaussian distribution with a given zero mean value and a variance $`\sigma _q`$.

The transmission Eq. (30) corresponds to an average over the phonon modes, treated classically (as in Appendix A 2) but with a mean-square displacement equal to the quantum zero-point motion, while the electronic transmission $`T(E;\{\mathrm{\Delta }_q\})`$ is obtained from quantum mechanics. In this average the influence of the phonon displacements on the electronic spectrum is taken into account but not vice versa. This is the fundamental difference with the transport calculations presented in section III B where both the electron and phonon degrees of freedom are treated at the quantum level and where the lattice distortion is induced by the injected electron. We therefore expect the results for zero-point motion to be different from the one due to the quantum coherent electron-phonon coupling.

Figure 8 shows the elastic transmission probability through the rigid undistorted lattice and the transmission averaged (in the manner discussed below) over the zero-point fluctuations. The total effective transmission obtained from the inelastic scattering calculations is also shown.

As expected, the lattice fluctuations substantially modify the electron transmission through the molecular wire. The transmission is reduced around the resonance peaks of the one-electron spectrum of the undistorted chain, and the peaks are smeared out. However, there is no shift in the position of the resonances corresponding to the formation of virtual polarons.

As in the fully quantum calculations, the transmission is enhanced for injection energies inside the gap of the undistorted chain. However, at this point it becomes important to know exactly how the average over the lattice fluctuations in Eq. (30) is calculated. It is known  that the ensemble average of the transmission probability in a disordered one-dimensional conductor is a statistically ill-defined quantity, in the sense that it is dominated by a very small number of exceptional configurations. In other words, the mean is in no way representative of the typical transmission probability to be expected from a randomly sampled member of the ensemble. This difficulty does not arise if $`\mathrm{log}T`$, rather then $`T`$, is averaged . In Figure 8 we therefore compare $`\mathrm{log}T`$ from our fully quantum results with $`\mathrm{log}T_{\mathrm{dis}}`$. In this case we see that fully including the electron-phonon coupling (allowing the possibility of polaron formation) gives a greater increase in the tunneling transmission than averaging over the disorder introduced by the zero-point motion. This shows that only the fully quantum coherent electron-phonon coupling calculations (exact in the limit of zero temperature) give the correct physics.

## IV Conclusion

In this paper, we have presented a method to calculate the inelastic effects on electron transport through one-dimensional molecular wires, including realistic electron-phonon coupling. The model for the wires is inspired by the Su-Schrieffer-Heeger model for trans-polyacetylene. The transport through the quantum electron-phonon system is solved by means of a multichannel scattering technique where each channel is associated with probabilities for the electron to be reflected or transmitted, given the phonon occupation number configuration.

The results show that the transport in the one-dimensional molecular wires does not occur as in traditional (three-dimensional) semiconducting molecular devices. It does not involve the propagation of free electron-like particles, but instead is due to the coherent propagation of “quasi-particles” : an electron (or hole) surrounded by a lattice distortion, characteristic of the formation of an electron (or hole) polaron. This object is called here a virtual polaron because of the transient nature of the tunneling electron (or hole) injected inside the HOMO-LUMO gap of the molecular wire. In this regime and in the limit of low temperatures, the tunneling transmission probability through the device increases, due to the electron-phonon coupling, in comparison with the transmission obtained from elastic scattering through the undistorted molecular wire. Lattice fluctuations also modify the electron transmission through the wire. However, the corresponding enhancement of the transmission in the tunneling regime is less than that produced by the virtual polarons.

The influence of other defects (such as a soliton-like defect existing in the chain) on the transmission has already been considered, and will be presented elsewhere. The transport for finite temperatures (different initial phonon occupations and proper statistical averages) is under study and results will be presented in the near future.

The results presented in this paper are general and can be applied to other types of one-dimensional atomic-scale wires subject to a Peierls transition. For example, it has been shown both experimentally and theoretically that a Peierls-like transition also occurs in dangling-bond (DB) lines fabricated on the H-passivated Si(001) surface . More recently, it has been shown theoretically that the injection of a static charge in the bands around the band-gap leads to a distortion of the atomic positions along the line . This distortion corresponds to the formation of a small polaron in the DB line. We therefore expect that carrier injection in the band-gap of the DB lines will lead to similar physical results to those presented here for molecular wires. Many other quasi-one-dimensional systems may be expected to show similar characteristics.

###### Acknowledgements.

The authors greatly appreciate enlightening discussions with L. Kantorovich, J. Gavartin, A.L. Shluger and A.M. Stoneham. HN thanks C. Joachim for stimulating discussions on inelastic tunneling in molecular systems. We acknowledge support from the U.K. Engineering and Physical Sciences Research Council in the form of an Advanced Fellowship (AJF), a Postdoctoral Fellowship (HN), a Research Studentship (SAS) and under Grant No. GR/M09193.

## A Isolated molecular wires

### 1 The SSH model

Su, Schrieffer and Heeger (SSH) modelled a chain of trans-polyacetylene ($`t`$-PA) as a purely one-dimensional atomic (CH)<sub>x</sub> chain . The model combines a classical ball-and-spring term for the distortion of the $`\sigma `$-bond backbone with a tight-binding representation for the delocalised $`\pi `$-electron orbitals along the chain. Furthermore the electron hopping integrals between adjacent (CH) groups are expanded linearly about some reference values. The corresponding Hamiltonian is

$$H_{\mathrm{SSH}}=\underset{i,s}{}\left(t_0\alpha \left(u_{i+1}u_i\right)\right)\left[c_{i,s}^{}c_{i+1,s}+c_{i+1,s}^{}c_{i,s}\right]+\frac{1}{2}K\underset{i}{}\left(u_{i+1}u_i\right)^2,$$
(A1)
where $`c_{i,s}^{}`$ ($`c_{i,s}`$) creates (annihilates) a $`\pi `$-electron of spin $`s`$ at site $`i`$. $`u_i`$ is the displacement of the $`i`$-th (CH) group from its place in the reference (undimerised) system with hopping integral $`t_0`$. $`K`$ is the spring constant corresponding to the $`\sigma `$-bond and $`\alpha `$ is the electron-lattice coupling constant. As explained in the introduction, the undimerised metallic chain is unstable with respect to a Peierls distortion, and the ground state of an infinite neutral chain has displacements given by $`u_i=(1)^iu_0`$, where $`u_0`$ is a constant depending on $`t_0`$, $`\alpha `$ and $`K`$. The values chosen for the different parameters ($`t_0`$=2.5 eV, $`\alpha `$=6.1 eV/Å and $`K`$=42.0 eV/Å<sup>2</sup> give for an infinite perfectly dimerised chain a total band width of 10 eV, a band gap of 1.4 eV, and a difference between long/short bond lengths of 0.1 Å in agreement with experiments.

For finite size chains containing $`N_a`$ atomic sites (i.e. (CH) groups), the ground state GS($`N_a`$,$`N_e`$) of Hamiltonian (A1) can be obtained for neutral ($`N_e=N_a`$) or charged ($`N_eN_a`$) chains, $`N_e`$ being the total number of $`\pi `$-electrons inside the chain. The ground state is obtained when the restoring forces of the springs balance the electronic forces by solving

$$2\alpha \underset{n}{\overset{\mathrm{occ}}{}}Z_n(i)\left(Z_n(i1)Z_n(i+1)\right)+K\left(2u_iu_{i1}u_{i+1}\right)=0,$$
(A2)
where $`Z_n(i)`$ are the components of the one-electron eigenstates (with energy $`ϵ_n`$) of $`H_{\mathrm{SSH}}`$ for a given atomic configuration $`\{u_i\}`$. The sum in Eq.(A2) runs only over the occupied states and the spin index $`s`$ is implicitly taken into account in the $`n`$ summations. For the finite size molecular chains we study here, the boundary conditions are $`Z_n(i)=0`$ and $`u_i=0`$ when $`i<1`$ or $`i>N_a`$. Furthermore in practice, in order to avoid the uniform translation of the chain through space, one atomic site is kept fixed, for instance one end of the molecular chain is kept fixed (i.e. $`u_{N_a}=0`$).

### 2 Classical phonons

In the limit of small displacements $`\{\delta u_i\}`$ around the equilibrium atomic positions $`\{u_i^0\}`$ of the ground state GS($`N_a`$,$`N_e`$), it is possible to derive an effective harmonic phonon Hamiltonian by partitioning Eq.(A1) as follows

$$H_{\mathrm{SSH}}[\{c_n\},\{u_i^0+\delta u_i\}]=H_{\mathrm{st}}[\{c_n\},\{u_i^0\}]+H_{\mathrm{ph}}[\{u_i^0\},\{\delta u_i\}]+H_{e\mathrm{ph}}[\{c_n\},\{\delta u_i\}],$$
(A3)
where $`H_{\mathrm{st}}`$, $`H_{e\mathrm{ph}}`$, $`H_{\mathrm{ph}}`$ are respectively the static, phonon and electron-phonon coupling parts of the total SSH Hamiltonian. $`H_{e\mathrm{ph}}`$ is usually treated as a perturbation up to second order to give a quadratic term in $`\delta u_i`$ in the effective phonon Hamiltonian $`_{i,j}K_{ij}\delta u_i\delta u_j`$. The dynamical matrix $`K_{ij}`$ is given by

$$K_{ij}=K_{ji}=2\alpha ^2\underset{n}{\overset{\mathrm{unocc}}{}}\underset{m}{\overset{\mathrm{occ}}{}}\left(F(i,n,m)F(i+1,n,m)\right)\left(F(j,n,m)F(j+1,n,m)\right)/\left(ϵ_mϵ_n\right),$$
(A4)
where $`F(j,n,m)=Z_n(j)Z_m(j1)+Z_m(j)Z_n(j1)`$. The $`n`$ ($`m`$) sum runs over the empty (occupied) electronic states. The eigenstates of dynamical matrix $`K_{ij}`$ give the eigenmodes $`V_q(i)`$ of vibration (phonons) of the finite size molecular chain, while the eigenvalues, $`M\omega _q^2`$, of $`K_{ij}`$ are related to the phonon frequencies $`\omega _q`$ ($`M`$ being the mass of the CH group).

In the present paper, we use the same boundary conditions as above to determine the phonon modes. However, different boundary conditions could be used : fixed ends (i.e. constant molecular chain length), free ends eventually coupled to different spring constants to simulate the effective coupling to the electrodes. Although, these different boundary conditions would in principle affect the electronic and the vibrational properties of the chain , we infer that the main physical results obtained in the present study will not be drastically modified. For instance, it appears that these different conditions will mostly affect the acoustic modes of the chain. Those modes have been however neglected in the present work because their contribution to the (virtual) polaron formation is negligible.

### 3 Quantum phonons

At this stage, we already have all the ingredients to derive a quantum version of the SSH Hamiltonian. From the reference system, chosen to be the neutral molecular chain of length $`N_a`$, we can write the Hamiltonians for the non-interacting electron and phonon degrees of freedom as

$$H_e=\underset{n}{}ϵ_nc_n^{}c_n,$$
(A5)
where $`c_n^{}=_iZ_n(i)c_i^{}`$ creates ($`c_n`$ annihilates) an electron in the $`n`$th electronic state of the reference system with energy $`ϵ_n`$. The harmonic phonon Hamiltonian (neglecting the zero-point energy) is

$$H_{\mathrm{ph}}=\underset{q}{}\mathrm{}\omega _qa_q^{}a_q,$$
(A6)
where $`a_q^{}`$ ($`a_q`$) creates (annihilates) a phonon mode $`q`$ with frequency $`\omega _q`$. The basis set associated to $`H_e+H_{\mathrm{ph}}`$ is formed by the eigenstates $`|n,\{n_q\}=c_n^{}_q\frac{(a_q^{})^{n_q}}{\sqrt{n_q!}}|0`$ with eigenvalues $`ϵ_{n,\{n_q\}}=ϵ_n+_qn_q\mathrm{}\omega _q`$, where $`|0`$ is the vacuum state and $`\{n_q\}`$ the set of phonon occupation numbers.

We expand the lattice deformations $`\delta x_i`$ induced by an additional charge introduced in the chain onto the phonon modes of the neutral chain: $`\delta u_i=_qV_q(i)\delta _q`$. The new lattice positions, displaced from the equilibrium position $`u_i^0`$, are $`u_i=u_i^0+\delta u_i`$. Then, the linear electron-phonon coupling term of the original SSH Hamiltonian is written in a quantum form by quantizing the phonon field displacements

$$\delta _q=\sqrt{\frac{\mathrm{}}{2M\omega _q}}\left(a_q+a_q^{}\right).$$
(A7)
Therefore the $`e`$-ph coupling Hamiltonian is

$$H_{e\mathrm{ph}}=\underset{q,n,m}{}\gamma _{qnm}\left(a_q^{}+a_q\right)c_n^{}c_m,$$
(A8)
where

$$\gamma _{qnm}=\underset{i=2}{\overset{N_a}{}}\lambda _q(i)\left[Z_n(i)Z_m(i1)+Z_n(i1)Z_m(i)\right],$$
(A9)
and

$$\lambda _q(i)=\alpha [V_q(i)V_q(i1)]\times \sqrt{\frac{\mathrm{}}{2M\omega _q}}.$$
(A10)
The total Hamiltonian $`H_\mathrm{w}`$ for the molecular wire with quantum phonons and linear electron-phonon coupling inspired by the SSH model is given by the sum $`H_\mathrm{w}=H_e+H_{\mathrm{ph}}+H_{e\mathrm{ph}}`$ as in Eq. (1).

Finally it should be noted that the $`e`$-ph coupling matrix elements $`\gamma _{qnm}`$ obey some selection rules. Generally, $`\gamma _{qnm}=0`$ unless the direct product $`\mathrm{\Gamma }_q\mathrm{\Gamma }_n\mathrm{\Gamma }_m`$ contains the identity representation ($`\mathrm{\Gamma }_n`$, $`\mathrm{\Gamma }_m`$ and $`\mathrm{\Gamma }_q`$ being the irreducible representation of the eigenstate $`n`$, $`m`$ and of $`\lambda _q`$ respectively). In practice, $`Z_n`$ and $`V_q`$ are even/odd functions with respect to the center of the molecule. The quantity $`Z_n(i)Z_m(i1)+Z_n(i1)Z_m(i)`$ is even (odd) when $`n+m`$ is an even (odd) integer (indexing the eigenvectors $`Z_n`$ by increasing eigenvalues and $`Z_{n=1}`$ being even). Whenever the quantity under the site $`i`$ summation in Eq. (A9) is odd, $`\gamma _{qnm}=0`$. Although it is not surprising to obtain selection rules for the $`e`$-ph coupling, their existence is very important in order to reduce the computing time of the product $`H_\mathrm{w}|\varphi `$ needed to solve Eq. (12).


format: text

text:
# An analogue of the Sommerfeld radiation condition for the Dirac operator.

## 1 Introduction

There are two principal ways to guarantee the unicity of solution of a not homogeneous differential equation of mathematical physics in the whole space. The first is the famous limiting absorption principle which in the case of some operators, for example the Helmholtz operator, has a very natural physical meaning reflected in its name (see, e.g., ). The medium is assumed to be slightly absorbent with a corresponding absorption parameter $`\epsilon .`$ For $`\epsilon >0`$ the inhomogeneous equation (under some additional natural conditions) has a unique solution $`u_\epsilon `$ decreasing at infinity, and the solution for $`\epsilon =0`$ is obtained as the limit of $`u_\epsilon `$ when $`\epsilon 0.`$ This approach, based on a deep understanding of the physical nature of the differential operator, is not applicable in all cases. Another approach is to impose a radiation condition at infinity which “prohibits” the existence of some unwanted solutions. Such a prohibition is related to the unicity of an integral representation for null-solutions of the differential operator. In the case of the Helmholtz operator the corresponding radiation condition was introduced by A. Sommerfeld .

For the Maxwell equations such condition was proposed by C. Müller  and S. Silver . Recently, A. McIntosh and M. Mitrea in  introduced a radiation condition for a special perturbed Clifford-Dirac operator. We use their approach to obtain a radiation condition for the operator $`D+M^\stackrel{}{\alpha }`$, where $`D`$ is the Moisil-Theodoresco operator and $`M^\stackrel{}{\alpha }`$ is the operator of multiplication by a complex quaternion $`\stackrel{}{\alpha }`$ from the right-hand side (for exact definitions see Section 2). Then using a simple matrix transformation proposed in , which converts the operator $`D+M^\stackrel{}{\alpha }`$ into the Dirac operator in its traditional form through the $`\gamma `$-matrices, we rewrite the radiation condition for the Dirac spinors in a form suitable for applications and obtain an integral representation for null-solutions of the classical Dirac operator in unbounded domains, based on the proposed radiation condition.

## 2 Preliminaries

### 2.1
Quaternions

A quaternion $`a`$ can be represented as a linear combination of its components ($`\{a_k\}`$ for real quaternions and $`\{a_k\}`$ for complex quaternions) and the elements of the orthonormal basis $`i_k`$ as follows

$$a=\underset{k=0}{\overset{3}{}}a_ki_k$$

where $`i_0:=1`$ and $`\{i_kk=1,2,3\}`$ are the quaternionic imaginary units, which have the following properties:

$$i_k^2=1,k=1,2,3;$$

$$i_1i_2=i_2i_1=i_3;i_2i_3=i_3i_2=i_1;i_3i_1=i_1i_3=i_2.$$
(1)
The imaginary unit in $``$ is denoted as usual by $`i`$ and by definition commutes with quaternionic imaginary units.

$$ii_k=i_ki,k=\overline{0,3}.$$

In some cases it is useful to represent a quaternion $`a`$ as the sum of a scalar part $`\mathrm{Sc}(a)=a_0`$ and a vectorial part $`Vec(a)=\stackrel{}{a}=a_1i_1+a_2i_2+a_3i_3`$. Using this representation we can define the product of two arbitrary quaternions $`a`$ and $`b`$ in the form

$$ab=a_0b_0\stackrel{}{a},\stackrel{}{b}+a_0\stackrel{}{b}+\stackrel{}{a}b_0+\left[\stackrel{}{a}\times \stackrel{}{b}\right],$$
(2)
where $`,`$ is the scalar product and $`[\times ]`$ is the vector product.

The quaternionic conjugation of a quaternion $`a=a_0+\stackrel{}{a}`$ is defined as $`\overline{a}=a_0\stackrel{}{a}.`$ From (2) it is easy to observe that

$$a\overline{a}=a_0^2+a_1^2+a_2^2+a_3^2=:|a|^2,$$
(3)
and to obtain an important property of quaternionic conjugation:

$$\overline{ab}=\overline{b}\overline{a}$$

Note that $`\left|ab\right|=\left|a\right|\left|b\right|.`$

One important conclusion which can be made using (3) is that each non-zero real quaternion $`a`$ is invertible and its inverse is given by

$$a^1=\frac{\overline{a}}{\left|a\right|^2}.$$

We denote the set of zero divisors in $`()`$ by $`𝔖`$:

$$𝔖:=\{a()a0;b0:ab=0\}.$$

###### Lemma 1 ()

(Structure of the set of zero divisors) Let $`a()`$ and $`a0.`$ The following statements are equivalent:

1. $`a𝔖`$.

2. $`a\overline{a}=0.`$

3. $`a_0^2=\stackrel{}{a}^2.`$

4. $`a^2=2a_0a=2\stackrel{}{a}a.`$

As the modulus introduced by (3) does not in general give information about the absolute values of the components of a complex quaternion when it is a zero divisor, another kind of modulus is frequently used

$$\left|a\right|_c:=\sqrt{\left|a_0\right|^2+\left|a_1\right|^2+\left|a_2\right|^2+\left|a_3\right|^2},$$
(4)
where $`\left|a_k\right|^2=a_ka_k^{}`$ and “\*” stands for the usual complex conjugation. (4) can be expressed also in the following ways

$$\left|a\right|_c^2=\left|\mathrm{Re}a\right|^2+\left|\mathrm{Im}a\right|^2=\mathrm{Sc}(a\overline{a}^{})=\mathrm{Sc}(\overline{a}^{}a).$$

### 2.2 Quaternionic differential operators

Let $`\mathrm{\Omega }^3`$ and $`fC^1(\mathrm{\Omega };())`$. The Moisil-Theodoresco operator is defined as

$$Df:=\underset{k=1}{\overset{3}{}}i_k_kf$$

where $`_k:=\frac{}{x_k}.`$ An important property of $`D`$ is obtained from (3)

$$D^2=\mathrm{\Delta },$$
(5)
where $`\mathrm{\Delta }=_1^2+_2^2+_3^2`$ is the usual Laplace operator.

We will consider the following operator

$$D_{\pm \nu }:=D\pm \nu I$$

where $`\nu `$ is a complex number, which was studied for the first time in  (for $`\nu `$ real). With the aid of the operator $`D_\nu `$ it is possible to factorize the Helmholtz operator

$$\mathrm{\Delta }+\nu ^2=(D+\nu )(D\nu )=D_\nu D_\nu .$$
(6)
Suppose that $`\theta _\nu `$ is a fundamental solution of the Helmholtz operator

$$(\mathrm{\Delta }+\nu ^2)\theta _\nu =\delta .$$

Then, using (6) we obtain the functions

$$𝒦_{\pm \nu }:=(D\nu )\theta _\nu $$
(7)
which are the fundamental solutions of $`D_{\pm \nu }`$, that is

$$D_{\pm \nu }𝒦_{\pm \nu }=\delta .$$

Let us assume that $`\mathrm{Im}\nu 0.`$ Under this condition, the fundamental solution of the Helmholtz operator is chosen as

$$\theta _\nu (x)=\frac{e^{i\nu \left|x\right|}}{4\pi \left|x\right|},$$
(8)
in which case it fulfills the Sommerfeld radiation condition at infinity and has a clear physical sense of an outgoing wave generated by a point source at the origin.

Substituting (8) into equality (7), we obtain

$$𝒦_{\pm \nu }(x)=D\theta _\nu (x)\pm \nu \theta _\nu (x)=\left(\pm \nu +\frac{x}{\left|x\right|^2}i\nu \frac{x}{\left|x\right|}\right)\theta _\nu (x),$$

where $`x:=_{k=1}^3x_ki_k.`$ These functions were obtained in , (see also \[6, Section 3\]).

Let us assume that $`\mathrm{\Omega }`$ is a bounded domain in $`^3`$ with a piecewise smooth boundary $`\mathrm{\Gamma }:=\mathrm{\Omega }`$.

###### Theorem 2

(Quaternionic Stokes’ formula) Let $`f`$ and $`g`$ belong to $`C^1(\mathrm{\Omega };())C(\overline{\mathrm{\Omega }};())`$. Then

$$_\mathrm{\Omega }\left((D_rf(y))g(y)+f(y)(Dg(y))\right)𝑑y=_\mathrm{\Gamma }f(y)\stackrel{}{n}(y)g(y)𝑑\mathrm{\Gamma }_y,$$
(9)
where $`\stackrel{}{n}(y)`$ is the quaternionic representation of the outward unit normal to the surface $`\mathrm{\Gamma }:\stackrel{}{n}(y)=_{k=1}^3i_kn_k(y)`$ and $`D_rf:=_{k=1}^3_kfi_k`$ is the right Moisil-Theodoresco operator.

The proof of this fact can be found in \[6, Chapter 4\] or \[3, p. 86\].

Let us consider the following integral operator

$$K_\nu [f](x):=_\mathrm{\Gamma }𝒦_\nu (xy)\stackrel{}{n}(y)f(y)𝑑\mathrm{\Gamma }_y,x^3\mathrm{\Gamma }$$
(10)
which can be considered as an analogue of the Cauchy integral operator from complex analysis due to the following fact.

###### Theorem 3

(Quaternionic Cauchy integral formula) Let $`fC^1(\mathrm{\Omega };())C(\overline{\mathrm{\Omega }};())`$ and $`f\mathrm{ker}D_\nu (\mathrm{\Omega })`$. Then

$$f(x)=K_\nu [f](x),x\mathrm{\Omega }.$$
(11)
### 2.3 Dirac operator in quaternionic form

Denote $`\stackrel{~}{\mathrm{\Phi }}:=\mathrm{\Phi }(t,x_1,x_2,x_3).`$ The domain $`\stackrel{~}{𝒢}`$ is assumed to be obtained from the domain $`𝒢^4`$ by the reflection $`x_3x_3.`$ Let us consider the transformation denoted by $`A`$ and introduced in the following way. A function $`\mathrm{\Phi }:𝒢^4^4`$ is transformed into a function $`\rho `$ $`:\stackrel{~}{𝒢}^4()`$ by the rule

$$\rho =A[\mathrm{\Phi }]:=\frac{1}{2}\left((\stackrel{~}{\mathrm{\Phi }_1}\stackrel{~}{\mathrm{\Phi }_2})i_0+i(\stackrel{~}{\mathrm{\Phi }_0}\stackrel{~}{\mathrm{\Phi }_3})i_1(\stackrel{~}{\mathrm{\Phi }_0}+\stackrel{~}{\mathrm{\Phi }_3})i_2+i(\stackrel{~}{\mathrm{\Phi }_1}+\stackrel{~}{\mathrm{\Phi }_2})i_3\right).$$

Then the definition of the inverse transformation $`A^1`$ is

$$A^1[\rho ]=(i\stackrel{~}{\rho }_1\stackrel{~}{\rho }_2,\stackrel{~}{\rho }_0i\stackrel{~}{\rho }_3,\stackrel{~}{\rho }_0i\stackrel{~}{\rho }_3,i\stackrel{~}{\rho }_1\stackrel{~}{\rho }_2).$$

These transformations can be represented in a more explicit matrix form:

$$\rho =A[\mathrm{\Phi }]=\frac{1}{2}\left(\begin{array}{cccc}\hfill 0& \hfill 1& \hfill 1& \hfill 0\\ \hfill i& \hfill 0& \hfill 0& \hfill i\\ \hfill 1& \hfill 0& \hfill 0& \hfill 1\\ \hfill 0& \hfill i& \hfill i& \hfill 0\end{array}\right)\left(\begin{array}{c}\stackrel{~}{\mathrm{\Phi }_0}\\ \stackrel{~}{\mathrm{\Phi }_1}\\ \stackrel{~}{\mathrm{\Phi }_2}\\ \stackrel{~}{\mathrm{\Phi }_3}\end{array}\right)$$

and

$$\mathrm{\Phi }=A^1[\rho ]=\left(\begin{array}{cccc}\hfill 0& \hfill i& \hfill 1& \hfill 0\\ \hfill 1& \hfill 0& \hfill 0& \hfill i\\ \hfill 1& \hfill 0& \hfill 0& \hfill i\\ \hfill 0& \hfill i& \hfill 1& \hfill 0\end{array}\right)\left(\begin{array}{c}\stackrel{~}{\rho }_0\\ \stackrel{~}{\rho }_1\\ \stackrel{~}{\rho }_2\\ \stackrel{~}{\rho }_3\end{array}\right).$$

The Dirac $`\gamma `$-matrices in the standard (, ) Dirac-Pauli form are

$$\begin{array}{cc}\gamma _0:=\left(\begin{array}{cccc}\hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\end{array}\right),& \gamma _1:=\left(\begin{array}{cccc}\hfill 0& \hfill 0& \hfill 0& \hfill 1\\ \hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\\ \hfill 1& \hfill 0& \hfill 0& \hfill 0\end{array}\right),\end{array}$$

$$\begin{array}{cc}\gamma _2:=\left(\begin{array}{cccc}\hfill 0& \hfill 0& \hfill 0& \hfill i\\ \hfill 0& \hfill 0& \hfill i& \hfill 0\\ \hfill 0& \hfill i& \hfill 0& \hfill 0\\ \hfill i& \hfill 0& \hfill 0& \hfill 0\end{array}\right),& \gamma _3:=\left(\begin{array}{cccc}\hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\\ \hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\end{array}\right);\end{array}$$

$$\gamma _5:=i\gamma _0\gamma _1\gamma _2\gamma _3=\left(\begin{array}{cccc}\hfill 0& \hfill 0& \hfill 1& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0& \hfill 1\\ \hfill 1& \hfill 0& \hfill 0& \hfill 0\\ \hfill 0& \hfill 1& \hfill 0& \hfill 0\end{array}\right).$$

Recall that:

1. $`\gamma _0^2=\gamma _5^2=E_4,`$ the identity matrix;

2. $`\gamma _k^2=E_4,`$ $`k=1,2,3`$;

3. $`\gamma _j\gamma _k=\gamma _k\gamma _j`$, $`j,k=0,1,2,3,5,`$ $`jk`$.

The transforms $`A`$ and $`A^1`$ have the following algebraic properties :

1. $`A\gamma _1\gamma _2\gamma _3\gamma _1A^1[\rho ]=i_1\rho ;`$

2. $`A\gamma _1\gamma _2\gamma _3\gamma _2A^1[\rho ]=i_2\rho ;`$

3. $`A\gamma _1\gamma _2\gamma _3\gamma _3A^1[\rho ]=i_3\rho ;`$

4. $`A\gamma _1\gamma _2\gamma _3\gamma _0A^1[\rho ]=\rho i_1;`$

5. $`A\gamma _1\gamma _2\gamma _3A^1[\rho ]=i\rho i_2.`$

The Dirac operator for a free massive particle of spin 1/2 is

$$𝔻[\mathrm{\Phi }]:=\left(\gamma _0_t\underset{k=1}{\overset{3}{}}\gamma _k_k+im\right)[\mathrm{\Phi }],$$

where $`m`$ and $`\mathrm{\Phi }`$ is a $`^4`$-valued function. Consider the quaternionic differential operator

$$R:=P_1^+(i_t+D)+P_1^{}(i_t+D)mM^{i_2},$$

where $`P_k^\pm :=\frac{1}{2}M^{(1\pm ii_k)}`$ and $`M^\lambda [f]:=f\lambda ,`$ $`\lambda ()`$. Using the algebraic properties of $`A`$ and $`A^1`$ we obtain the following equality :

$$R=A\gamma _1\gamma _2\gamma _3𝔻A^1.$$
(12)
We consider a time-harmonic null-solution $`\mathrm{\Phi }`$ of the operator $`𝔻`$ which has the form $`\mathrm{\Phi }(t,x)=q(x)e^{i\omega t}`$, where $`\omega `$ and $`q(x)`$ is a $`^4`$-valued function depending only on $`x`$. We have the following equation for $`q`$

$$𝔻_{\omega ,m}q(x)=0,$$

where $`𝔻_{\omega ,m}:=i\omega \gamma _0_{k=1}^3\gamma _k_k+imI`$. The corresponding quaternionic reformulation of $`𝔻_{\omega ,m}`$ is the operator

$$D_\alpha :=D+M^\alpha ,$$

where $`\alpha :=(i\omega i_1+mi_2).`$ That is

$$D_\alpha =A\gamma _1\gamma _2\gamma _3𝔻_{\omega ,m}A^1.$$
(13)
## 3 Radiation condition for the Dirac operator

In order to prove the Cauchy integral formula for null-solutions of the operator $`D_\nu `$ for the exterior domain $`^3\backslash \overline{\mathrm{\Omega }}`$ we need some appropriate radiation condition at infinity.

Let us consider the equation

$$D_\nu 𝒦(x)=\delta (x),\nu 0,x^3.$$

When $`\mathrm{Im}\nu =0`$ it admits two solutions decreasing at infinity obtained (according to (7)) by the application of the operator $`D_\nu `$ to the fundamental solutions of the Helmholtz operator $`u^\pm (x):=\frac{e^{\pm i\nu \left|x\right|}}{4\pi \left|x\right|}`$ respectively:

$$𝒦^\pm (x)=\left(\nu +\frac{x}{\left|x\right|^2}i\nu \frac{x}{\left|x\right|}\right)u^\pm (x).$$

In order to omit one of these possibilities we impose the following radiation condition

$$\left(\nu \frac{x}{\left|x\right|^2}+i\nu \frac{x}{\left|x\right|}\right)𝒦(x)=o\left(\frac{1}{\left|x\right|}\right),\text{when }\left|x\right|\mathrm{}.$$
(14)
It is easy to see that

$$\left(\nu \frac{x}{\left|x\right|^2}+i\nu \frac{x}{\left|x\right|}\right)𝒦^+(x)=O\left(\frac{1}{\left|x\right|^2}\right),\text{when }\left|x\right|\mathrm{}$$

and hence $`𝒦^+`$ fulfills (14). For the function $`𝒦^{}`$ this is not true. Note that $`𝒦^+`$ is precisely the fundamental solution $`𝒦_\nu `$.

We are now ready to prove the following

###### Theorem 4

(Quaternionic Cauchy integral formula for the exterior domain) Let

$$fC^1(^3\overline{\mathrm{\Omega }};())C(^3\mathrm{\Omega };()),f\mathrm{ker}D_\nu (^3\overline{\mathrm{\Omega }}),Im\nu 0$$

and let $`f`$ satisfy the radiation condition

$$\left(\nu \frac{x}{\left|x\right|^2}+i\nu \frac{x}{\left|x\right|}\right)f(x)=o\left(\frac{1}{\left|x\right|}\right),\text{when}\left|x\right|\mathrm{}.$$

Then

$$f(x)=K_\nu [f](x),x^3\overline{\mathrm{\Omega }}.$$
(15)
Proof. Let $`\mathrm{\Gamma }^R`$ be a sphere with center at the origin and radius $`R`$ sufficiently large so that $`\mathrm{\Omega }`$ is contained in the ball $`B^R`$ with boundary $`B^R=\mathrm{\Gamma }^R`$ (Fig. 3).

According to Theorem 3 in each point $`x`$ of the domain $`\mathrm{\Omega }^R:=B^R\overline{\mathrm{\Omega }}`$ we have the equality

$$f(x)=_\mathrm{\Gamma }𝒦_\nu (xy)\stackrel{}{n}(y)f(y)𝑑\mathrm{\Gamma }_y_{\mathrm{\Gamma }^R}𝒦_\nu (xy)\frac{y}{\left|y\right|}f(y)𝑑\mathrm{\Gamma }_y^R.$$

We now consider the limit of this equality when $`R\mathrm{}`$. We have the following asymptotic relation

$$_{\mathrm{\Gamma }^R}𝒦_\nu (xy)\frac{y}{\left|y\right|}f(y)𝑑\mathrm{\Gamma }_y^R_{\mathrm{\Gamma }^R}\theta (y)\left(\nu \frac{y}{\left|y\right|^2}+i\nu \frac{y}{\left|y\right|}\right)f(y)𝑑\mathrm{\Gamma }_y^R,R\mathrm{}.$$

Using the radiation condition we obtain that this integral tends to zero when $`R\mathrm{}`$. Thus,

$$f(x)=_\mathrm{\Gamma }𝒦_\nu (xy)\stackrel{}{n}(y)f(y)𝑑\mathrm{\Gamma }_y.$$

The introduction of the radiation condition in the form (14) allowed us to obtain a very simple proof of the Cauchy integral formula for the exterior domain. Nevertheless, as was shown in  this condition can be represented in a more elegant form. We note that the expression which appears in (14) can be rewritten as follows

$$\left(\nu \frac{x}{\left|x\right|^2}+i\nu \frac{x}{\left|x\right|}\right)=\nu \left(1+\frac{ix}{\left|x\right|}\right)+O\left(\frac{1}{\left|x\right|}\right),\left|x\right|\mathrm{}.$$

Thus, a natural idea is to introduce the radiation condition as

$$\left(1+\frac{ix}{\left|x\right|}\right)f(x)=o\left(\frac{1}{\left|x\right|}\right),$$
(16)
because the term $`x/\left|x\right|^2`$ apparently gives a faster decay. The problem here is that (16) does not imply the decay of the product $`\frac{x}{\left|x\right|^2}f(x)`$ when $`\left|x\right|\mathrm{}`$ due to the fact that $`(1+\frac{ix}{\left|x\right|})`$ is a zero divisor. The following example makes this clear. Consider $`f(x):=(1\frac{ix}{\left|x\right|})\left|x\right|^2`$. This function obviously satisfies (16) but $`\frac{x}{\left|x\right|^2}f(x)=O(\left|x\right|)`$, $`\left|x\right|\mathrm{}`$. Thus, in order to be able to apply (16) instead of (14), we should prove that if $`f`$ belongs to $`\mathrm{ker}D_\nu (^3\overline{\mathrm{\Omega }})`$ and satisfies (16) then it decreases at infinity. Following arguments from  we prove this.

###### Theorem 5

Let $`fC^1(^3\overline{\mathrm{\Omega }};())C(^3\mathrm{\Omega };()),f\mathrm{ker}D_\nu (^3\overline{\mathrm{\Omega }}),Im\nu 0`$, such that satisfies (16). Then

$$_{\left|x\right|=R}\left|f(x)\right|_c^2𝑑\mathrm{\Gamma }_x^R=O(1)\text{as }R\mathrm{}.$$
(17)
Proof. Consider the expression

$`\left|\left(1+{\displaystyle \frac{ix}{\left|x\right|}}\right)f\right|_c^2`$ $`=\mathrm{Sc}\left(\left(1+{\displaystyle \frac{ix}{\left|x\right|}}\right)f\overline{f}^{}\left(1+{\displaystyle \frac{ix}{\left|x\right|}}\right)\right)=2Sc\left(\left(1+{\displaystyle \frac{ix}{\left|x\right|}}\right)f\overline{f}^{}\right)=`$
$`=2\left(\left|f\right|_c^2ImSc(\overline{f}^{}{\displaystyle \frac{x}{\left|x\right|}}f)\right).`$
Consequently,

$`{\displaystyle _{\left|x\right|=R}}\left|\left(1+{\displaystyle \frac{ix}{\left|x\right|}}\right)f(x)\right|_c^2𝑑\mathrm{\Gamma }_x^R`$ $`=2({\displaystyle _{\left|x\right|=R}}|f(x)|_c^2d\mathrm{\Gamma }_x^R`$ (18)
$`ImSc({\displaystyle _{\left|x\right|=R}}\overline{f}^{}(x){\displaystyle \frac{x}{\left|x\right|}}f(x)d\mathrm{\Gamma }_x^R)).`$
Due to Theorem 2 we have (using the notation of Theorem 4)

$`{\displaystyle _{\left|x\right|=R}}\overline{f}^{}(x){\displaystyle \frac{x}{\left|x\right|}}f(x)𝑑\mathrm{\Gamma }_x^R`$ $`={\displaystyle _\mathrm{\Gamma }}\overline{f}^{}(x)\stackrel{}{n}(x)f(x)𝑑\mathrm{\Gamma }_x+`$
$`{\displaystyle _{\mathrm{\Omega }^R}}(D_r(\overline{f}^{}(x))f(x)+\overline{f}^{}(x)Df(x))𝑑x.`$
The first integral on the right-hand side is some constant $`C`$. In order to simplify the second we use $`Df=\nu f`$ and observe that

$$D_r\overline{f}^{}=\overline{Df^{}}=\overline{\nu ^{}f^{}}=\nu ^{}\overline{f}^{}.$$

Then

$`{\displaystyle _{\mathrm{\Omega }^R}}(D_r(\overline{f}^{}(x))f(x)+\overline{f}^{}(x)Df(x))𝑑x`$ $`={\displaystyle _{\mathrm{\Omega }^R}}(\nu ^{}\overline{f}^{}(x)f(x)\nu \overline{f}^{}(x)f(x))𝑑x=`$
$`=2iIm\nu {\displaystyle _{\mathrm{\Omega }^R}}\overline{f}^{}(x)f(x)𝑑x,`$
we obtain that

$$Sc\left(_{\left|x\right|=R}\overline{f}^{}(x)\frac{x}{\left|x\right|}f(x)𝑑\mathrm{\Gamma }_x^R\right)=C_02iIm\nu _{\mathrm{\Omega }^R}\left|f(x)\right|_c^2𝑑x.$$

Substituting this expression in (18) and using the radiation condition (16) we have

$$_{\left|x\right|=R}\left|f(x)\right|_c^2𝑑\mathrm{\Gamma }_x^R=ImC_02Im\nu _{\mathrm{\Omega }^R}\left|f(x)\right|_c^2𝑑x,R\mathrm{}.$$

The last term above vanishes when $`Im\nu =0`$ and is negative when $`Im\nu >0`$. In both cases (17) is proved.

This theorem establishes the equivalence of the radiation conditions (14) and (16) for the functions from $`\mathrm{ker}D_\nu `$.

We now obtain the radiation condition for the Dirac operator in quaternionic form, that is for $`D_\alpha `$, and prove the Cauchy integral formula for the exterior domain for this operator. After that, using the transform $`A`$ we will rewrite these results for the operator $`𝔻_{\omega ,m}.`$

Let $`\alpha `$ be a purely vectorial complex quaternion $`\alpha =\stackrel{}{\alpha }`$, $`\alpha 𝔖`$, $`\stackrel{}{\alpha }^20`$ and consider $`\nu :=\sqrt{\stackrel{}{\alpha }^2}`$ where the sign of the root is chosen in a way such that $`Im`$ $`\nu 0.`$

Let

$$K_\stackrel{}{\alpha }[f](x):=P^+K_\nu [f](x)+P^{}K_\nu [f](x)$$
(19)
where the operators $`K_{\pm \nu }`$ are defined by (10) and $`P^\pm :=\frac{1}{2\nu }M^{(\nu \pm \stackrel{}{\alpha })}`$ are two mutually complementary projection operators in the set of $`()`$-valued functions.

In the case when $`\alpha :=\stackrel{}{\alpha }𝔖`$, $`\nu =0,`$ we introduce the operator $`K_\stackrel{}{\alpha }`$ as follows

$$K_\stackrel{}{\alpha }f:=K_0fM^\stackrel{}{\alpha }V_0f$$
(20)
where

$$V_0[f](x):=\underset{\mathrm{\Gamma }}{}\theta _0(xy)\stackrel{}{n}(y)f(y)𝑑\mathrm{\Gamma }_y,x^3\mathrm{\Gamma }.$$
(21)
Here $`\theta _0(xy)=\frac{1}{4\pi \left|xy\right|}.`$

For all $`\alpha ()`$ (that is $`\alpha 𝔖`$ or $`\alpha 𝔖`$), $`K_\stackrel{}{\alpha }f`$ admits the following representation

$$\begin{array}{c}K_\stackrel{}{\alpha }[f](x)=\underset{\mathrm{\Gamma }}{}\{\theta _\nu (xy)(\frac{xy}{\left|xy\right|^2}\frac{i\nu (xy)}{\left|xy\right|})\stackrel{}{n}(y)f(y)+\hfill \\ \\ \hfill +\theta _\nu (xy)\stackrel{}{n}(y)f(y)\stackrel{}{\alpha }\}d\mathrm{\Gamma }_y,x\mathrm{\Gamma }.\end{array}$$
(22)
###### Theorem 6

Let $`fC^1(^3\overline{\mathrm{\Omega }};())C(^3\mathrm{\Omega };())`$, $`f\mathrm{ker}D_\alpha (^3\overline{\mathrm{\Omega }})`$, $`\alpha ()`$, $`\mathrm{Sc}\alpha =0`$ such that $`f`$ satisfies the radiation condition

$$\nu f(x)+\frac{ix}{\left|x\right|}f(x)\stackrel{}{\alpha }=o\left(\frac{1}{\left|x\right|}\right),\text{when }\left|x\right|\mathrm{},$$
(23)
where $`\nu :=\sqrt{\stackrel{}{\alpha }^2}`$ and $`\mathrm{Im}\nu 0.`$ If $`\alpha 𝔖`$ then we suppose additionally that $`f(x)=o(1).`$ Then

$$f(x)=K_\stackrel{}{\alpha }[f](x),x^3\overline{\mathrm{\Omega }}.$$

Proof. First we suppose that $`\stackrel{}{\alpha }𝔖`$. Multiplying (23) by $`\stackrel{}{\alpha }/\nu `$ from the right-hand side we obtain

$$f(x)\stackrel{}{\alpha }+\nu \frac{ix}{\left|x\right|}f(x)=o\left(\frac{1}{\left|x\right|}\right)$$
(24)
adding and subtracting (23) and (24) we obtain

$$\left(\nu f(x)+\frac{ix}{\left|x\right|}f(x)\stackrel{}{\alpha }\right)\pm \left(f(x)\stackrel{}{\alpha }+\nu \frac{ix}{\left|x\right|}f(x)\right)=$$

$$=\left(1\pm \frac{ix}{\left|x\right|}\right)f(x)\nu \pm \left(1\pm \frac{ix}{\left|x\right|}\right)f(x)\stackrel{}{\alpha }=$$

$$=\left(1\pm \frac{ix}{\left|x\right|}\right)f(x)\left(\nu \pm \stackrel{}{\alpha }\right)=o\left(\frac{1}{\left|x\right|}\right).$$
(25)
Which can be writen as follows

$$P^+\left(\left(1+\frac{ix}{\left|x\right|}\right)f(x)\right)+P^{}\left(\left(1\frac{ix}{\left|x\right|}\right)f(x)\right)=o\left(\frac{1}{\left|x\right|}\right)\text{when }\left|x\right|\mathrm{}.$$
(26)
Thus, (23) is equivalent to (26), from which it can be seen that $`P^+f`$ fulfills the radiation condition (16) and $`P^{}f`$ fulfills its conjugate which corresponds to the operator $`D_\nu .`$ Consequently the integrals $`K_{\pm \nu }P^\pm f`$ taken over the sphere $`\mathrm{\Gamma }^R`$ (see the proof of Theorem 4) decrease when $`R\mathrm{}.`$

Since $`K_{\pm \nu }`$ commute with $`P^\pm `$ we obtain that the integral $`K_\stackrel{}{\alpha }f(x)`$ taken over the sphere $`\mathrm{\Gamma }^R`$ also decreases when $`R\mathrm{}.`$

In the case when $`\alpha `$ is a zero divisor the radiation condition (23) becomes

$$\frac{ix}{\left|x\right|}f\stackrel{}{\alpha }=o\left(\frac{1}{\left|x\right|}\right),\text{when }\left|x\right|\mathrm{}.$$
(27)
Since the behavior of $`\theta _0`$ in (21) is of the type $`O\left(\frac{1}{\left|x\right|}\right)`$ when $`\left|x\right|\mathrm{}`$ and since the expression $`M^\stackrel{}{\alpha }V_0`$ contains the multiplication from the right-hand side by $`\stackrel{}{\alpha }`$, and in this case $`f\stackrel{}{\alpha }=o\left(\frac{1}{\left|x\right|}\right)`$, it can be seen that the integral (20) taken over the sphere $`\mathrm{\Gamma }^R`$ decreases at infinity also.

We now consider the Dirac operator $`D_\stackrel{}{\alpha }`$, $`\stackrel{}{\alpha }:=(i\omega i_1+mi_2)`$. Note that $`\stackrel{}{\alpha }𝔖`$ iff $`\omega ^2=m^2.`$ In the following lines we will obtain the radiation condition for the Dirac operator in its classical representation using $`\gamma `$-matrices.

Let us denote

$$q:=A^1f.$$
(28)
According to (13), if $`f\mathrm{ker}D_\stackrel{}{\alpha }`$ then $`q\mathrm{ker}𝔻_{\omega ,m}`$. Suppose that $`f`$ fulfills the radiation condition (23) then we have

$$\nu A\gamma _1\gamma _2\gamma _3A^1f+A\gamma _1\gamma _2\gamma _3A^1\frac{ix}{\left|x\right|}f\stackrel{}{\alpha }=o\left(\frac{1}{\left|x\right|}\right).$$

Let us denote

$$f\stackrel{}{\alpha }=:F$$

and consider the product

$$xF=(x_1i_1+x_2i_2+x_3i_3)F\text{.}$$

Using properties 1.-3. of the transforms $`A`$ and $`A^1`$ (Section 2.3), we obtain

$$xF=A\gamma _1\gamma _2\gamma _3\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)A^1[F].$$

Consequently, using the properties 4. and 5., we have

$$xF(x)=A\gamma _1\gamma _2\gamma _3\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)A^1[f(x)(i\omega i_1mi_2)]=$$

$$=A\gamma _1\gamma _2\gamma _3\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)A^1A\gamma _1\gamma _2\gamma _3\gamma _0A^1f(x)(i\omega )+$$

$$+A\gamma _1\gamma _2\gamma _3\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)A^1A\gamma _1\gamma _2\gamma _3A^1f(x)(im)=$$

$$=A\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)\gamma _0A^1f(x)(i\omega )+A\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)A^1f(x)(im).$$

Using (28) we can rewrite the previous equation as

$$xF\text{ }=A\left(\underset{k=1}{\overset{3}{}}x_k\gamma _k\right)(i\omega \gamma _0im)q.$$

Thus (23) can be rewritten as follows

$$\nu f+\frac{ix}{\left|x\right|}f\stackrel{}{\alpha }=\nu Aq+A\frac{1}{\left|x\right|}\underset{k=1}{\overset{3}{}}x_k\gamma _k(\omega \gamma _0+m)q=o\left(\frac{1}{\left|x\right|}\right),$$

and hence the radiation condition for the Dirac operator in its traditional form is obtained in the following form

$$\nu q(x)(\omega \gamma _0m)\frac{\stackrel{}{x_\gamma }}{\left|x\right|}q(x)=o\left(\frac{1}{\left|x\right|}\right),\left|x\right|\mathrm{}$$
(29)
where $`\stackrel{}{x_\gamma }:=_{k=1}^3x_k\gamma _k`$.

We introduce the following operator

$$𝕂_{\omega ,m}:=A^1K_\stackrel{}{\alpha }A,$$

where $`\stackrel{}{\alpha }=(i\omega i_1+mi_2)`$. Using this operator, in  an analogue of the Cauchy integral formula for the Dirac operator $`𝔻_{\omega ,m}`$ in bounded domains was proved. Here with the aid of the obtained radiation condition we obtain this fact for the unbounded domains.

###### Theorem 7

(Cauchy integral formula for $`𝔻_{\omega ,m}`$ in the exterior domain $`^3\backslash \overline{\mathrm{\Omega }}`$) Let $`qC^1(^3\backslash \overline{\mathrm{\Omega }};())C(^3\backslash \mathrm{\Omega };())`$, $`q\mathrm{ker}𝔻_{\omega ,m}(^3\backslash \overline{\mathrm{\Omega }})`$, and $`q`$ satisfies the radiation condition (29). Then

$$q(x)=𝕂_{\omega ,m}[q](x),x(^3\backslash \overline{\mathrm{\Omega }}).$$

The proof consists in the application of (13) and Theorem 4.

Acknowledgment: This work was partially supported by CONACYT, project 32424-E.