AC quantum transport: Non-equilibrium in mesoscopic wires due to time-dependent fields

TL;DR

This paper models the non-equilibrium energy distribution of quasi-particles in a mesoscopic normal metal wire under time-dependent ac bias, revealing how high-frequency irradiation and interactions influence charge transport and energy states.

Contribution

It introduces a Green function-based quantum diffusion model for ac-driven mesoscopic wires, analyzing photon absorption effects and the impact of electron-electron and electron-phonon interactions on energy distributions.

Findings

Photon steps depend on field amplitude and photon energy.

In slow field regime, photon absorption is time-dependent.

Strong interactions lead to Fermi distributions with bath or effective temperature.

Abstract

A model is developed describing the energy distribution of quasi-particles in a quasi-one dimensional, normal metal wire, where the transport is diffusive, connected between equilibrium reservoirs. When an ac bias is applied to the wire by means of the reservoirs, the statistics of the charge carriers is influence by the formed non-equilibrium. The proposed model is derived from Green function formalism. The quasi-particle energy distribution is calculated with a quantum diffusion equation including a collision term accounting for inelastic scattering. The ac bias, due to high frequency irradiation, drives the wire out of equilibrium. For coherent transport the photon absorption processes create multiple photon steps in the energy distribution, where the number of steps is dependent on the relation between the amplitude of the field eV and the photon energy \omega. Furthermore we observe that for the slow field regime, \omega \tau_D < 1, the photon absorption is highly time-dependent. In the fast field regime \omega \tau_D > 1 this time-dependency disappears and the photon steps in the distribution have a fixed value. When the wire is extended, the transport becomes incoherent due to interaction processes, like electron-electron interaction and electron-phonon interaction. These interactions give rise to a redistribution of the quasi-particles with respect to the energy. We focused on the fast field regime and concluded that the strong interaction limit for both mechanisms gives the expected result. Strong electron-phonon interaction forces the distribution function on every position in the wire to become a Fermi function with the bath temperature, while strong electron-electron interaction causes an effective temperature profile across the wire and the distribution function on every position in the wire is a Fermi function with an effective temperature.

Figures

Equations

$$\frac{1}{{\tau }_D}\frac{d^{2}f(x,E)}{d{x}^2}+I_{coll}(x,E,f)=0$$

$$i\hslash \frac{\partial \psi (\text{r},t)}{\partial t}=\hat{H}\psi (\text{r},t);$$

$$\Psi (\text{r},t)=e^{-iEt/\hslash }\psi (\text{r}).$$

$$I=\frac{2_{s}e}{2\pi \hslash }{\int }_0^{\infty }Tr[{\hat{t}}^{+}\hat{t}][f_L(E)-f_R(E)]$$

$$I=G_{Q}V\sum _n{T}_n(\mu ).$$

$${\rho }_D(T)=\frac{⟨G⟩}{2G_Q}\frac{1}{T\sqrt{1-T}}$$

$$⟨\delta n_T\delta n_T⟩$$

$$⟨(\delta n_T{)}^2⟩$$

$$⟨(\delta n_R{)}^2⟩$$

$$S=\frac{e^2}{\pi \hslash }[2k_{b}T\sum _n{T}_n^2+eVcoth\left(\frac{eV}{2k_{b}T}\right)\sum _n{T}_n(1-T_n)].$$

$$\psi (x,y,z,t)={\psi }_0(x,y,z)e^{-iEt/\hslash }$$

$$H=H_0+eVcos(\omega t).$$

$$\psi (x,y,z,t)$$

$$\tilde{\rho }(E)=\sum _{n=-\infty }^{\infty }\rho (E+n\hslash \omega )J_n^2\left(\frac{eV}{\hslash \omega }\right).$$

$$I_{AB}=C{\int }_{-\infty }^{\infty }\left[f(E-e{V}_0)-f(E)\right]{\rho }_A(E-e{V}_0){\rho }_B(E)dE.$$

$${\tilde{I}}_{AB}=c\sum _{n=-\infty }^{\infty }{J}_n^2\left(\frac{eV}{\hslash \omega }\right){\int }_{-\infty }^{\infty }\left[f(E-e{V}_0)-f(E+n\hslash \omega )\right]{\rho }_A(E-e{V}_0){\rho }_B(E+n\hslash \omega )dE.$$

$$\tilde{\Gamma }=\sum _{n=-\infty }^{\infty }{J}_n^2(z)\Gamma (E+n\hslash \omega )$$

$$S_I=4G_Q{k}_{b}T\sum _n{T}_n^2+2G\sum _n{T}_n(1-T_n)\sum _{\pm }\sum _{l=0}^{\infty }{J}_l^2(\alpha )(eV\pm l\hslash \omega )coth\left(\frac{eV\pm l\hslash \omega }{2k_{b}T}\right)$$

$$\text{j}=\frac{1}{A}\frac{dQ}{dt}=-ne{\text{v}}_{ave}.$$

$$\text{j}=\frac{n{e}^2\tau }{m}\text{E}.$$

$$\sigma =\frac{n{e}^2\tau }{m}$$

$$\frac{\kappa }{\sigma }=\frac{{\pi }^2}{3}{\left(\frac{k_b}{e}\right)}^{2}T.$$

$$i\hslash \frac{\partial \psi (x,t)}{\partial t}=H(x,t)\psi (x,t)=\left(\frac{p^2}{2m}+V(x,t)\right)\psi (x,t).$$

$${\psi }_{adiabatic}(x,t)\approx {\psi }_{V_0(t)}(x)e^{-i{E}_{V_0(t)}t/\hslash }.$$

$$i\hslash \frac{\partial {\psi }_{adiabatic}(x,t)}{\partial t}$$

$$E={\hslash }^2\sum _{\text{k}}\frac{{\text{k}}^2}{2m}n(\text{k})$$

$$\delta E={\hslash }^2\sum _{\text{k}}\frac{{\text{k}}^2}{2m}\delta n(\text{k})$$

$$\delta E=\sum _{\text{k}}{ϵ}_{\text{k}}^0\delta n(\text{k})+\frac{1}{2V}\sum _{\text{k},{\text{k}}^{\prime }}g(\text{k},{\text{k}}^{\prime })\delta n(\text{k})\delta n({\text{k}}^{\prime }).$$

$${ϵ}_{\text{k}}=\frac{\delta E}{\delta n(\text{k})}={ϵ}_{\text{k}}^0+\frac{1}{V}\sum _{\text{k},{\text{k}}^{\prime }}g(\text{k},{\text{k}}^{\prime })\delta n({\text{k}}^{\prime }).$$

$$v_F=\frac{p_F}{m^{\ast }},$$

$$\frac{\partial f(x,E)}{\partial t}+D\frac{{\partial }^{2}f(x,E)}{\partial x^2}=0$$

Full text

The proposed model is derived from Green function formalism. The quasi-particle energy distribution is calculated with a quantum diffusion equation including a collision term accounting for inelastic scattering. The ac bias, due to high frequency irradiation, drives the wire out of equilibrium. For coherent transport the photon absorption processes create multiple photon steps in the energy distribution, where the number of steps is dependent on the relation between the amplitude of the field $eV$ and the photon energy $\hbar\omega$. Furthermore we observe that for the slow field regime, $\omega\tau_{D}<1$, the photon absorption is highly time-dependent. In the fast field regime $\omega\tau_{D}>1$ this time-dependency disappears and the photon steps in the distribution have a fixed value.

When the wire is extended, the transport becomes incoherent due to interaction processes, like electron-electron interaction and electron-phonon interaction. These interactions give rise to a redistribution of the quasi-particles with respect to the energy. We focused on the fast field regime and concluded that the strong interaction limit for both mechanisms gives the expected result. Strong electron-phonon interaction forces the distribution function on every position in the wire to become a Fermi function with the bath temperature, while strong electron-electron interaction causes an effective temperature profile across the wire and the distribution function on every position in the wire is a Fermi function with an effective temperature.

So the complicated interplay between the effect of photon absorption, diffusive transport and inelastic scattering on the quasi-particle energy distribution seems to be accurately described by our model.

[TABLE]

**AC quantum transport:

Non-equilibrium in mesoscopic wires due to time-dependent fields**

1. 1 Introduction

2. 1.1 Non-equilibrium and mesoscopic systems

3. 1.2 AC quantum transport

4. 2 Phase coherent quantum transport

5. 2.1 Scattering theory

6. 2.1.1 Transport

7. 2.1.2 Shot noise

8. 2.2 Tien-Gordon theory

9. 2.3 Photon-assisted shot noise

10. 3 Diffusive transport

11. 3.1 Drude-Sommerfeld model

12. 3.2 Landau theory of Fermi liquids

13. 3.3 Transport in quasi-one dimensional metallic systems

14. 3.4 Quantum corrections to the conductivity

15. 4 Photon absorption and other energy exchange processes in diffusive wires

16. 4.1 Introduction

17. 4.2 Green function formalism

18. 4.3 Quantum diffusion equation for the distribution function

19. 4.4 Limit situations for the simple quantum diffusion equation

20. 4.4.1 Slow field limit

21. 4.4.2 Fast field limit

22. 4.5 Incorporating inelastic interactions

23. 4.6 Inelastic scatterering

24. 4.6.1 Electron-phonon interaction

25. 4.6.2 Electron-electron interaction

26. 4.7 Summary

27. 5 Numerical results

28. 5.1 Introduction

29. 5.2 Calculation method

30. 5.3 Simulation of realistic coherent and incoherent transport situations

31. 5.3.1 Coherent transport in diffusive wires

32. 5.3.2 Incoherent transport in diffusive wires

33. 6 Probing the quasi-particles energy distribution

34. 6.1 System description

35. 6.2 Fabrication

36. 6.3 Measurements on NIS junctions

37. 6.3.1 Differential conductance of a NIS junction

38. 6.3.2 NIS differential conductance for ac biased wire

39. 7 Conclusion and discussion

40. 7.1 Conclusion

41. 7.2 Discussion

42. A Shot noise

43. B MATLAB code of the simulation program

44. B.1 Script for the simulation of coherent transport

45. B.2 Functions for the used operators

46. B.3 Script for the simulation of incoherent transport

47. C Space dependency in the distribution function

48. D Differential conductance of a NIS junction

49. E MATLAB script for deconvolution

Chapter 1 Introduction

1.1 Non-equilibrium and mesoscopic systems

The last decades the non-equilibrium in mesoscopic systems is intensively studied by a part of the nano-scientific community. Despite all the efforts the physics of this is still not fully understood due to the complexity of these systems. The systems have length scales between microscopic and macroscopic. On one hand the system contains many particles, but on the other hand it can still exhibit quantum features. Because of the intermediate dimensions a specific approach is needed for calculating the physical properties. Pure quantum mechanics can not be used because the many particles complicate the quantum mechanical description in a horrible way and thermodynamics can not be used because of the significance of the quantum features in the system. Therefore often a quantum statistical approach is used which reveals the intriguing world of mesoscopic physics.

Before looking at mesoscopic systems, let’s look at macroscopic and microscopic systems and the meaning of equilibrium and non-equilibrium in this context. Consider a macroscopic resistor $R$ placed between electron reservoirs at equilibrium, which means that the electrons in the reservoirs obey Fermi statistics and the electrons with energy $E$ are distributed according to a Fermi function, $f(E)=(e^{(E-\mu)/(k_{b}T)}+1)^{-1}$, where $k_{b}$ is the Boltzmann constant and $T$ the temperature. At zero temperature this Fermi function is just a step function at the chemical potential $\mu$ of the material. For energies lower than the chemical potential all energy levels are occupied and for higher energies all levels are empty. When the temperature is increased the electrons become thermally excited, creating holes for energies below chemical potential and electrons for higher energies. This can be seen as a quasi-equilibrium situation. When we look at the unexcited resistor between the reservoirs, we see that the electrons are at the same equilibrium, or quasi-equilibrium, as the reservoirs. Now when a dc voltage $V$ is applied on the reservoirs a current will flow from one reservoir through the resistor to the other reservoir by the relation $I=V/R$. The resistance on the flowing electrons due to impurities causes dissipation, heating the resistor. The statistics of the electrons in the resistor are no longer the same statistics as that of the reservoirs and becomes spatial dependent. The heating of the resistor causes a local equilibrium in the resistor and the electron energy distribution is described by an effective electron temperature [1]. The applied power $P=V^{2}/R$ causes a temperature profile along the wire which is bounded by the temperature of the reservoirs. Such an effective temperature profile is shown in figure 1.1. The temperature of the reservoirs is held at 4.2 K and the effective temperature is at maximum in the middle of the wire. Figure 1.1 also shows the local equilibrium distribution function at the boundary of the wire and in the middle of the wire. The effect of the dissipated energy is a thermal smearing around the Fermi energy.

For the opposite case, a microscopic system, the situation is completely different. A scatterer is placed between two reservoirs and a voltage is applied. In this microscopic situation it becomes more convenient to evaluate the transport using scattering theory [3], so we do not speak anymore of distribution functions inside the transport region. The electron approaches the scatterer as a Fermi particle with a wave function. Because of the wave-particle duality the electron can be transmitted or reflected with a certain probability by the scatterer, whereafter the electron leaves the scatterer as a Fermi particle with a certain wave function. The transport of the electrons through the scatterer depends on the properties of the scatterer. These properties are described by the transmission distribution, which gives the probability of finding a transport channel in the scatterer with a certain transmission probability. A bias voltage applied to the reservoirs will only create a potential difference across the structure and the transport depends on this potential difference and the transmission distribution of the scatterer. The number of electrons involved in the transport is a measure of the non-equilibrium.

So the non-equilibrium of macroscopic systems is described by the temperature and resistance of the object and the non-equilibrium of microscopic systems is revealed by scattering theory. Now the intermediate regime between macroscopic and microscopic: mesoscopic. In this research we will focus on a diffusive wire, which shows the most resemblance with the macroscopic situation where a resistor was evaluated. However, the general idea that the electron energy distribution inside the wire can be described by an effective temperature appears to breaks down. Pothier et al. studied the effect of a dc voltage on a diffusive wire between electron reservoirs [4]. From this research it was concluded that the electron energy distribution obeys the time-independent Boltzmann equation when the driving term, i.e. the potential difference across the wire, is absorbed in the boundary conditions.

[TABLE]

In absence of inelastic interactions the collision integral vanishes and the solution is on every position in the wire a superposition of the boundary conditions which are the Fermi function of right reservoir and that of the left reservoir. When one reservoir is held at zero potential, the other reservoir is at maximum potential which shifts this Fermi function with $eU$. The superposition of these two distribution functions creates a two step function dependent on the position on the wire as shown in figure 1.2.

When inelastic interactions are involved the situation becomes a bit more complicated. The collision integral in the Boltzmann equation has to be evaluated. The energy that electrons gained from the electric field is redistributed during collisions on inelastic scatterers. These inelastic scattering processes are electron-electron and electron-phonon interactions. Depending on the characteristics of the diffusive wire and the dominant scattering processes, a relation for the energy relaxation time can be obtained which is self-consistently used in calculating the distribution function.

1.2 AC quantum transport

So far only non-equilibrium of dc quantum transport is considered. The study of non-equilibrium of ac quantum transport in mesoscopic systems is interesting for better understanding of the physics of many-body systems and how small electronic devices respond to high frequency irradiation. Previous studies on ac quantum transport focused mainly on coherent structures, where the phase of electrons is preserved. Different examples of study objects of ac quantum transport are SIS junctions, quantum point contacts (QPC), quantum dots (QD) and resonant tunneling diodes (RTD). Tien and Gordon successfully constructed a theory describing the tunneling current between two superconducting films separated by an insulating layer biased with an ac voltage [5]. The electrons involved in the transport can gain energy in discrete values from the ac field creating steps in the $I-V$ characteristics. The success of their theory reached further than the SIS and was also successfully applied to the QPC, QD and RTD.

Stimulated by the success of this theory for different structures Remco Schrijvers [6] tried to apply this theory to the reservoirs and use the Boltzmann equation to calculate the electron energy distribution in a diffusive wire excited by an ac voltage. The validity of this approach was a bit disappointing. The model was only valid for low frequencies in a wire without inelastic scattering. This was caused by the fact that Tien-Gordon theory assumes averaging over time and therefore the collision term of the Boltzmann equation can not be evaluated in a correct manner.

From the previous research on non-equilibrium due to time-dependent fields in diffusive wires it was concluded that the situation is still not completely understood. To avoid the deducted problems put forward by Remco Schrijvers, we derived from the Green function formalism a quantum diffusion equation for the electron energy distribution in a quasi-one dimensional diffusive wire subject to an oscillating electric field. The model is first derived for a coherent structure with elastic impurity scattering, whereafter this model is extended to account for inelastic scattering processes such as electron-electron and electron-phonon interactions.

Chapter 2 Phase coherent quantum transport

2.1 Scattering theory

2.1.1 Transport

Phase coherent quantum transport involves the transport of charge carriers where the phase of these charge carriers is preserved. Generally this means that scattering inside the structure is elastic, so that the energy of the charge carriers is not redistributed. Phase coherent transport of electrons in nanostructures is usually described with scattering theory. The nanostructure is defined as a scattering region between reservoirs and the wave function of the electrons subject to Hamiltonian $\hat{H}$ with potential $U(\textbf{r},t)$ obeys the Schrodinger equation

[TABLE]

The solution of the Schrodinger equation is a stationary space-dependent function multiplied by a time-dependent function dependent on the eigen energy $E$ of the Hamiltonian:

[TABLE]

The wave function $\psi(\textbf{r})$ obeys the time-independent Schrodinger equation $\hat{H}\psi(\textbf{r})=E\psi(\textbf{r})$. Due to the wave character of a charge carrier, an electron can contribute to the current through the scatterer between the reservoirs by either being reflected or being transmitted. The probability of being reflected or transmitted is dependent on the thickness and height of the barrier whereon the electron scatters. The potential difference across the structure is determined by the difference of the energy distribution of the two electron reservoirs. Landauer’s result for the current through the scatterer between reservoirs is proportional to the integral over energy of the trace of the product of the transmission matrix $\hat{t}$ and its conjugate transpose $\hat{t}^{+}$ and the difference between the energy distribution of the left and right reservoir [7]. An insightful derivation of the Landauer formula can be found in Ref [3].

[TABLE]

The factor $2_{s}$ accounts for the degeneracy of electrons with charge $e$. When a bias is applied to the reservoirs, creating across the structure a potential difference $V$, much smaller than the scale of energy dependence in the transmission eigenvalues $T_{n}$, equation 2.3 can be evaluated at the Fermi energy $\mu$. Introducing the conductance quantum $G_{Q}=2e^{2}/h$ gives for the current

[TABLE]

This expression for the current through a scattering structure clearly shows that the structure exists of different channels in which the electrons are transported with a certain probability from one reservoir to the other. The type of transport structure is characterized by the distribution of the transmission probabilities. This distribution is constructed by taking one specific nanostructure from an ensemble of identical design and counting the number of transmission eigenvalues of the transmission matrix in the interval of $T$ to $dT$. This is divided by the total number of nanostructures in the ensemble. For large enough ensembles, the result converges to $P(t)dT$, so that the transmission distribution is defined as $P(t)=\left\langle\sum_{p}\delta(T-T_{p}(E))\right\rangle$. For very short structures, where the wavelength of the electron exceeds the length of the structure, the conductance quantization is prominent present and the distribution of the transmission probabilities is sharply peaked on certain values. When the length of the structure increases, the resistance due to defects in the system becomes dominant. The diffusive behavior of the electrons in the scatterer is random and for a diffusive scatterer the distribution of transmission probabilities is universal, i.e. independent on the details of the scatterer [3].

[TABLE]

Here $\left\langle G\right\rangle$ is the average conductance due to many scattering events. Now with the increasing dimensions of the structure the describing picture becomes more and more complicated due to the fact that more charge carriers are involved and inelastic scattering processes affect the energy of the charge carriers. Therefore one has to let go the idea that the energy of electrons is unchanged by the scattering events. Pure scattering theory can no longer describe in an effective way the transport. Quantum statistical mechanics provides a way out as we will see later on. First we look at the statistical information of charge carriers that the noise due to finite transmission probabilities in scattering processes provides.

2.1.2 Shot noise

A physical phenomenon that contains information about statistics of charge carriers in a mesoscopic conductor is shot noise. Shot noise is caused by the quantization of charge [8]. When a single incident charge in a state with occupation 1 scatters on some potential barrier it has a probability $R$ of being reflected and a probability $T=1-R$ of being transmitted. Figure 2.2 shows how the incoming wave packet of an electron scattering on a barrier with transmission probability $T$ is splitted and only a part of the initial wave packet is transmitted with a probability $T$, causing fluctuations in the current.

When the initial state is occupied by the distribution function $f$, an incident particle is reflected with probability $fR$ and transmitted with probability $fT$, so the averaged occupation of the reflected state is $\left\langle n_{R}\right\rangle=fR$ and the averaged occupation of the transmitted state is $\left\langle n_{T}\right\rangle=fT$. By looking at many scattering processes the fluctuations from the average occupation can be determined. For the incident state the average occupation is just the Fermi distribution $\left\langle n_{in}\right\rangle=f$, so that the mean squared fluctuations in the incident state vanishes: $\left\langle(f-\left\langle n_{in}\right\rangle)^{2}\right\rangle=0$. The fluctuations in the reflected and transmitted state have a finite value. The fluctuations are expressed as a deviation from the average so $\delta n_{T}=n_{T}-\left\langle n_{T}\right\rangle$ and $\delta n_{R}=n_{R}-\left\langle n_{R}\right\rangle$. When we use these identities to calculate the mean squares of the correlations between reflected and transmitted state and of the reflected and transmitted state itself we find:

[TABLE]

From these expressions we can distinguish two limits. One limit is given by full transparency and the other limit is given by full reflectance. Both limits have the same outcome in the fluctuations. In a situation where the occupation of the initial state is given by a Fermi distribution at zero temperature, the mean square fluctuations vanish. However, for finite temperature this is not the case. The mean square fluctuations does not vanish, but fluctuates like the incident state with occupation $f$.

These mean square fluctuations contribute in the current and from the current expressions derived in appendix A the noise power can be obtained. When a multi-channel scatterer between two reservoirs is considered, the noise power can be evaluated at Fermi energy when the scale of energy dependence of the transmission coefficients is much larger than the thermal energy and the energy associated with the applied bias voltage on the reservoirs. The shot noise power is then [8]:

[TABLE]

As we will see later on, the shot noise for an ac bias has a bit different form than equation 2.9 and therefore the non-equilibrium due to the ac transport can be seen in the shot noise.

2.2 Tien-Gordon theory

When we make the switch from dc quantum transport to ac quantum transport, the needed describing theoretical frameworks become a bit more sophisticated. Approximately five decades ago Dayem and Martin observed interactions of electrons with photons in the tunneling current between the superconducting films A and B separated by an insulating layer, when the structure was illuminated with microwave radiation, causing an ac bias across the junction [10]. Figure 2.3 shows the clear difference between the $I-V$ characteristic with and without this oscillating electric field.

In order to explain these quantum interactions Tien and Gordon developed a describing theory for electric fields normal and parallel to the surface of the superconductor [5]. Here we will only consider the case where the field is normal to the surface of the superconductor.

The potential difference between the superconductors A and B due to the electric field is given by $Vcos(\omega t)$, where the bias is applied to one reservoir and the other reservoir is held at zero potential. When no field is present the wave functions of the charge carriers of energy $E$ satisfy the unperturbed Hamiltonian $H_{0}$.

[TABLE]

The perturbed Hamiltonian due to the oscillating electric field is given by

[TABLE]

This interaction Hamiltonian only effects the time-dependent part of the wave function given by equation 2.10. The new wave function under influence of the oscillating electric field becomes

[TABLE]

To come to the last line the identity $e^{zsin(\theta)}=\sum^{\infty}_{n=-\infty}J_{n}(z)e^{in\theta}$ is used, where $J_{n}(z)$ is the Bessel function giving the probability of the absorption of $n$ field quanta. The wave function in equation 2.12 is normalized, since $\left[\sum^{\infty}_{n=-\infty}J_{n}(z)\right]^{2}=1$. It appears that the wave function no longer has one energy variable. The energy variable is extended in a sum of multiples of the photon energy. This means that where a charge carrier in the situation without the oscillating field could only tunnel to a state with the same energy, now also could tunnel to states with energy $E\pm n\hbar\omega$. Basically the density of states of the superconductor is modulated by the electric field. The unperturbed density of states of the superconductor is $\rho(E)$. In the presence of the oscillating field the density of states $\widetilde{\rho}(E)$ becomes

[TABLE]

The tunnel current is calculated from the density of states. For an SIS junction biased with a dc voltage $V_{0}$ the tunnel current is

[TABLE]

Here $C$ is a proportionality constant depending on the junction resistance. When an additional ac voltage is applied to the SIS junction the tunnel current shows the multiple photon steps.

[TABLE]

When the tunnel current is explicitly calculated it shows indeed the photon steps as measured by Dayem and Martin. Figure 2.4 shows the difference between the measured tunnel current without oscillating electric field given by the solid lines and the calculated tunnel current with oscillating electric field between two superconducting films for two different ratios of $\frac{eV}{\hbar\omega}$ given by the dashed lines.

The energy diagram of an ac biased SIS junction in figure 2.5 shows explicitly how a photon assist the transport of an electron from the first superconductor through the insulating layer to the second superconductor. The gap in the density of states of the superconductor makes it impossible for an electron unaffected by the electric field to tunnel through the barrier to an unoccupied level in the second superconductor. The absorption of a photon can provide the required energy to make this possible.

As said in chapter 1 photon-assisted transport is, besides in SIS junctions, also observed in other nano-electronic systems. We won’t discuss all examples. Here we only have an additional look at the transport in a quantum dot illuminated with radiation, where the driving frequency exceeds the normal tunneling rate of electrons through the dot, since it provides great insight in the mechanism of photon-assisted transport.

A quantum dot is usually some island coupled by tunnel barriers to leads, the source and drain. The electronic properties of the island and the tunnel barriers can be controlled by gates. Figure 2.6 shows this schematically.

The energy levels on the island are assumed to be discrete with a spacing $\Delta E$ while the energy spectrum of the leads is assumed to be a continuum. The radiation is coupled to the island by the gate [13]. We will not go into detail about this, since we mainly want to focus on the transport from drain to source. The normal tunneling rates are modified by the radiation due to the modification of the wave function of the electrons given by equation 2.12.

[TABLE]

Here $z=e\widetilde{V}/\hbar\omega$ and $\widetilde{V}$ is the amplitude of the oscillation. The tunneling is assisted by the absorption of photon with energy $E+n\hbar\omega$ and emission of photons with energy $E-n\hbar\omega$. The possible tunneling processes in the dot with and without radiation are shown in figure 2.7. Only the upper energy diagram in the middle can contribute to a current through the dot without help of radiation. The remaining diagrams show the photon-assisted tunneling through the ground state $\epsilon_{0}$ and the first excited state $\epsilon_{1}$ of the dot. Electrons which normally do not have the right energy to tunnel to an unoccupied state can now absorb or emit a photon. This modifies their energy in such a way that tunneling becomes possible.

For both the SIS junction and the quantum dot the ac bias, due to radiation coupled on the structure, modulates the electronic properties making transport possible to energy states which are not accessible without the energy gain from the field. The photons from the field assist in the transport of charge carriers through the structure.

2.3 Photon-assisted shot noise

In section 2.1.2 the basic idea of shot noise in mesoscopic conductors for dc quantum transport is evaluated and we stated that the expression for the shot noise differs a bit for ac quantum transport. Here we will look how it differs and how this difference arises.

A general scatterer is placed between two reservoirs and an ac voltage is applied to the scatterer by the left reservoirs while the other reservoir is grounded. The transport of electrons can be divided into two regimes: transport of affected and unaffected electrons by the ac bias [9]. The unaffected electrons do not contribute to the shot noise, because the number of emitted, unaffected electrons from the right reservoir is the same as that of the left reservoir. Since according to the Pauli exclusion principle both left and right outgoing states can only be occupied by one electron, the current cancels and so does the fluctuation in current.

The affected electrons from the left reservoir can contribute to the shot noise. An electron with energy $\epsilon\leq\hbar\omega$ below the Fermi energy can get excited to an energy $\hbar\omega-\epsilon$. At energy $-\epsilon$ a hole is created. Since only the left reservoir can excite electrons in this way (the other reservoir is grounded), there is no counter current, so that this becomes the source of the fluctuations in the current. Now when also a dc voltage is applied to the scatterer, the shot noise expression becomes an extended version of equation 2.9 [15], where the photon-assisted features are presented by the Bessel functions like in the tunnel current calculated by Tien and Gordon.

[TABLE]

Here $\alpha=eV_{ac}/\hbar\omega$. For $V_{ac}=0$ the normal expression for shot noise is obtained. Now when we make the transition to a diffusive wire it appears that this description still holds. Schoelkopf et al. [16] investigated the photon-assisted shot noise experimentally for phase-coherent diffusive conductors and compared their results to the theoretical predictions for photon-assisted shot noise stated by Lesovik and Levitov [15]. The ac bias is applied on the conductor by bending the conductor between the reservoirs in a loop. A time-dependent magnetic field enters the loop, which induces a time-dependent electric field in the conductor. The situation is depicted in figure 2.8.

Lesovik and Levitov predicted theoretically the photon steps in the noise power for an ac biased diffusive conductor where the phase of electrons is preserved. The experiment of Schoelkopf verifies this model. Figure 2.9 shows the experimental results and the expected results from equation 2.17 of the differential noise power. The photon steps are not that clear in the first derivative of the noise power. The second derivative of the noise power however clearly shows at the expected energies the steps, indicating the photon-assisted mechanism in the shot noise.

The discrete steps in the shot noise shows the absorption of field quanta and give information about the statistics of the charge carriers in the diffusive wire. It reveals that the energy distribution of the charge carriers inside the wire is affected by the ac bias. This is a completely different point of view in comparison to the transport in the SIS junction and the quantum dot where the electronic properties of the reservoirs are affected by the ac bias. So apparently there arises some interesting physics in the diffusive wire. This is still a relatively simple model, where the electron transport is coherent, so that scattering theory still can be used to describe the transport. However, when the length of the diffusive wire is increased and not only diffusivity and photon absorption causes a change in statistics in the wire, but also inelastic scattering processes induce energy redistribution, scattering theory is no longer the most convenient describing theory. As said in section 2.1.1 we can proceed with quantum statistical theory to determine the statistics of the charge carriers described by the energy distribution function. In the next chapter we will evaluate the conditions for such an approach.

Chapter 3 Diffusive transport

3.1 Drude-Sommerfeld model

The model that was proposed in the 1900s by Drude describes the transport properties of electrons in metals on a microscopic level from a classical point of view. The electronic properties of a metal are then described by a gas of electrons bouncing on heavier positive charged ions. Because of the higher mass of the ions, they are seen as static potentials and the collisions of the electrons on these ions are purely elastic. The electrons involved in the transport are assumed to be free. Between two scattering events no forces act on the electron. In a situation where no electric field is applied on the metal conductor, the average velocity due to different electrons cancels, as the electrons move in a variety of directions. When an electric field is applied the average velocity and thus the net current becomes finite. If $n$ electrons per unit volume with charge $-e$ move with the average velocity $\textbf{v}_{ave}$ and move in a time $dt$ a distance $\textbf{v}dt$, then the net charge passing through a cross-section $A$ is $-ne\textbf{v}_{ave}Adt$ [17] [18]. The current density becomes

[TABLE]

Now when an electron is considered at time zero with velocity $\textbf{v}_{0}$, the velocity that this electron can gain from the electric field in time $t$ is $-e\textbf{E}t/m$ following from Newton’s laws of motion. The initial velocity $\textbf{v}_{0}$ of every electron does not contribute to the average velocity, due to the random collisions from which the electron emerges on time zero. From this it is also directly clear that the average time $t$ is the average time between collision $\tau$, so that the average velocity is $\textbf{v}_{ave}=-e\textbf{E}\tau/m$. Substituting this in the current density gives

[TABLE]

Ohm’s law is given by $\textbf{j}=\sigma\textbf{E}$, where $\sigma$ is the conductivity. Equating the current density of equation 3.2 and from Ohm’s law gives the final expression of the conductivity.

[TABLE]

Based on the observation that metals conduct heat better than insulators the assumption was made that the electrons involved in the electric conduction also carry the thermal current. The original Drude model used the Maxwell-Boltzmann distribution to account for the probability of finding an electron with a certain energy and thus a certain velocity. However, the ratio between thermal and electric conductivity observed in experiments was not explained in this way. Then the Pauli exclusion principle was put forward, which stated that two fermions can never occupy the same state. From this the conclusion was drawn that the Maxwell-Boltzmann distribution had to be replaced by the Fermi-Dirac distribution. Sommerfeld exchanged the Maxwell-Boltzmann distribution by the Fermi-Dirac distribution in the classical electron gas of Drude. This modified the expression for the electronic velocity and gave the correct expression of the ratio between thermal and electric conductivity, the Wiedemann-Franz law [17]:

[TABLE]

The idea that electrons form a gas in a metal is sufficient for cases where no energy exchange is present in all processes involving the electrons. However, this is not always the situation. When the collisions of the electrons are no longer purely elastic and they cause energy exchange, the Drude-Sommerfeld model breaks down. Fortunately Landau’s theory of Fermi liquids provides a strong replacement.

3.2 Landau theory of Fermi liquids

As said in the previous section, at a certain stage the transport of electrons can no longer be explained in a electron gas model where the interactions are purely elastic. The effect of inelastic interactions becomes significant and the energy exchange processes initiate the break down of the electron gas concept. Instead one considers the transport of electrons in a liquid model. This Fermi liquid model is developed by Lev Landau in 1956. The transport of one electron is affected by the surrounding electrons and its wave function is extremely complicated due to screening effects. It behaves however still very like a particle with a charge $e$. The screening can simply be seen as the modification of the relation between energy and wave vector, so $E(\textbf{k})=\hbar^{2}k^{2}/2m^{*}$, where $m^{*}$ deviates from the free electron mass $m$. The electrons are defined as quasi-particles which are stable near the Fermi level, but lose their stability far from the Fermi level [19].

The domain of validity for excitations near the Fermi surface in the Landau theory of Fermi liquids has its origin in the assumed one-to-one correspondence between states of a non-interacting system and states of an interacting system when the interaction is adiabatically turned on. Since the lifetime of a quasi-particle is proportional to $(\epsilon-\epsilon_{F})^{-2}$, the high energy quasi-particles are decayed before the interaction process is fully complete [20]. The adiabatic continuation leads to the assumption that the excited states of the interacting system are labeled with the same quantum numbers as the excited states of the non-interacting system. The validity of adiabatic continuation from a non-interacting system to a interacting system can be shown by looking at the wave function. An example is given by a particle trapped in an one-dimensional potential $V(x,t)=V_{0}(t)h(x)$ [21]. The wave function obeys the Schrodinger equation

[TABLE]

Now the potential changes slowly from initial value $V_{01}$ to a final value $V_{02}$. Because the potential varies slowly, the solution of the Schrodinger equation can be approximated by the solution of the static Schrodinger equation $H(x,t)\psi_{V_{0}(t)}(x)=E_{V_{0}(t)}\psi_{V_{0}(t)}(x)$ [22]. The adiabatic solution becomes

[TABLE]

By inserting equation 3.6 in equation 3.5 the accuracy of equation 3.6 is obtained.

[TABLE]

The adiabatic solution is a good approximation for the wave function in an one-dimensional potential $V(x,t)$ if the first term of equation 3.7 dominates the second term, which is true if the rate of change of $V_{0}(t)$ is small enough. Then the solution for the new potential $V_{0}=V_{02}$ is found from the old value of the potential $V_{0}=V_{01}$ from which it adiabatically rises. This implies that when the excited state of the initial potential is a bound state, the excited state of the final potential is also a bound state. A transition from a bound state to an un-bound state will never occur from an adiabatic continuation, no matter how small the rate of change in $V(x,t)$, because one is a decaying function while the other is an oscillatory function.

As said the interactions cause a modification of the relation between energy and momentum of a particle. The total energy of an unperturbed electron system is given by the kinetic energy of the electrons [23].

[TABLE]

Here $n(\textbf{k})$ is the occupation number of the state with momentum k. When a weak external field is coupled on the system, there will occur a change in occupation number and thus a change in total energy.

[TABLE]

If the system now is perturbed by a adiabatically turned on interaction, with interaction energy $g(\textbf{k},\textbf{k}^{\prime})$ between states of wave vector k and $\textbf{k}^{\prime}$, the system is taken away from its ground state energy and a change of occupation numbers is induced. Therefore the change of energy is

[TABLE]

Due to the interaction the electron is no longer a pure particle, but it is a quasi-particle. It behaves still like a particle, but it arises from the interactions with its local environment. A quasi-particle with wave vector k has an energy of

[TABLE]

In the above we have suppressed magnetic fields, so that spin dependency can be neglected, since $\epsilon(\textbf{k},\sigma)=\epsilon(\textbf{k})$ in absence of magnetic fields.

A fundamental parameter in the Landau theory of Fermi liquids is the effective mass. The interaction experienced by a quasi-particle changes its mass with respect to the mass in an environment free of interactions. The velocity and density of states at the Fermi surface can be calculated using this effective mass.

[TABLE]

The expressions for these quantities are similar to that of a non-interacting system which confirms the one-to-one correspondence between the states of a non-interacting system and an interacting system. So concluding this section, we can take interactions into account in calculating the electronic properties in quantum transport by considering the charge carriers being quasi-particles for low excited states. Therefore the total energy of the system is not the sum of the energy of the individual particles, but is function of the energy distribution among the quasi-particles. Also due to the one-to-one correspondence between the states of a non-interacting system and an interacting system, the energy distribution of the quasi-particles can be calculated from a diffusion equation, like the semi-classical Boltzmann equation.

3.3 Transport in quasi-one dimensional metallic systems

The Landau theory of Fermi liquids, discussed in the previous section, provides the justification of using a semi-classical Boltzmann equation to calculate the energy distribution of the quasi-particles in a diffusive wire. In this work we focus on a quasi-one dimensional metallic wire of mesoscopic dimensions where the transport of the quasi-particles is diffusive. We will first explain what we exactly understand when we talk about quasi-one dimensional, mesoscopic and diffusive. Then we discuss the non-equilibrium in such a system biased with a dc voltage by looking at the energy distribution of the quasi-particle involved in the transport.

Mesoscopic structures are defined by the relation between length scales defining the geometrics of the structure and defining microscopic processes in the structure.

The length scales defining the microscopic processes involving an quasi-particle are:

• •

The Fermi wavelength $\lambda_{F}=2\pi/k_{F}$, where $k_{F}$ is the Fermi wave vector,

• •

The elastic mean free path $l_{e}$, which is the average distance between elastic collisions on impurities for instance,

• •

The phase coherence length $l_{\phi}$, which is the distance that the phase of a quasi-particle is preserved,

• •

The energy relaxation length $l_{E}$, which is the distance the energy of a quasi-particle is preserved.

The length scales defining the geometrics of the structure are given by:

• •

The length of the structure $L$,

• •

The cross-section of the structure $S$.

When the length of the structure is significantly larger than the cross-section it is more natural to talk about the structure as being a wire. The wire is said to be diffusive if the length $L$ of the wire is significantly larger than the elastic mean free path $l_{e}$ of an quasi-particle in the wire. The wire is quasi-one dimensional for $\lambda^{2}_{F}<<S<<l^{2}_{e}$ provided that the width and the thickness are of the same order of magnitude. When we want to be able to apply the Landau theory of Fermi liquids we are bound to at least quasi-one dimensional systems. For purely one dimensional systems the Landau theory of Fermi liquids is no longer valid. This has its origin in the nesting property of the Fermi surface, which means that a part of the Fermi surface can be matched onto an other part by a translation of $2k_{F}$. Therefore there arises a divergence in calculating physical properties. A more detailed explanation can be found in Ref. [24].

Pothier et al. studied the quantum transport in dc biased diffusive wires by looking at the effect of the induced non-equilibrium on the quasi-particle energy distribution [4]. The diffusive wire is placed between large electron reservoirs where the electron energy distribution is described by an equilibrium Fermi function. A dc voltage $U$ is applied on one reservoir while the other reservoir is held at zero potential, creating a potential difference $U$ over the wire.

The distribution function can be calculated by using semi-classical kinetic theory, which is the Boltzmann equation extended with an interaction term. For wires where the diffusion time $\tau_{D}$ is shorter than the relaxation time $\tau_{E}$ the transport is coherent and the distribution is described by a Boltzmann equation without interaction term. When the driving term $eU$ due to the potential difference $U$ across the wire is absorbed in the boundary conditions at the reservoirs, one Fermi function is unchanged, while the other is shifted by $eU$. This leads to the equation:

[TABLE]

The explicit boundary conditions for this equation are given by $f(0,E)=f_{F}(E)$ and $f(L,E)=f_{F}(E+eU)$, where $f_{F}(E)=(1+e^{E/k_{b}T})^{-1}$ is just the Fermi distribution. The stationary solution on every position in the wire is a superposition of the two boundary conditions.

[TABLE]

Figure 3.2 shows the two step function of the electron energy distribution on every position in a dc biased wire with no interactions present.

If inelastic scattering is introduced the situation becomes a bit more sophisticated. Two main phase breaking mechanisms can be distinguished: electron-electron interactions and electron-phonon interactions, where the electrons are considered to be quasi-particles. We first consider electron-electron interactions and neglect electron-phonon interactions. Strong scattering induces a local equilibrium with temperature $T_{e}(x)$ and the distribution is described by

[TABLE]

where $\mu(x)=-eU\frac{x}{L}$ [25]. The effective temperature $T_{e}(x)$ in a wire with cross-section $S$ and resistance $R$ is calculated from the heat equation [25].

[TABLE]

The boundary conditions of this equation are $T_{e}(0)=T_{e}(L)=T$ and using the Wiedemann-Franz law (equation 3.4) for the heat conductivity $\kappa$ the effective temperature is [25]

[TABLE]

Now the electron-electron interactions are negligible and the electron-phonon scattering is the dominant phase breaking mechanism. For strong scattering the electrons thermalize with the temperature of the phonons. The distribution function is given by $f(x,E)=f_{F}(E-\mu(x),T)$ where $\mu=-eU\frac{x}{L}$ and $T$ is the phonon bath temperature [25]. The space dependence of the distribution functions is shown for both situations in figure 3.3.

For intermediate regimes where neither electron-electron scattering nor electron-phonon scattering is strong, but still present, the interaction term in the Boltzmann equation has to be evaluated. The interaction term can be calculated from the Fermi golden rule and the belonging kernel follows from a microscopic derivation [26]. We will come back to this later in chapter 4 where we calculate the interactions in a diffusive wire due to electron-electron scattering and electron-phonon scattering.

3.4 Quantum corrections to the conductivity

On quantum scale the conductance of a diffusive wire is not simply given by the Drude result of the conductivity in equation 3.3. Because an electron has a wave-character, the electron is not localized. Therefore, when no phase-breaking processes are present, an electron can interfere with itself when it returns to a certain initial position after multiple elastic scattering events. This modification of the conductance is called localization and is depicted in figure 3.4.

The probability for an electron of passing between A and B is given by a classical probability and additionally an interference term

[TABLE]

The phase gained by an electron while traveling through the diffusive medium is $\Delta\phi=\hbar^{-1}\int^{B}_{A}\textbf{p}d\textbf{l}$. For most of the trajectories this phase gain will be much larger than one and therefore vanish in the interference term. The self-crossings have the same phase gain when the direction of the traveled trajectory is reversed, i.e. $\textbf{p}\rightarrow-\textbf{p}$ and $d\textbf{l}\rightarrow-d\textbf{l}$. This results in two paths and the probability of self-crossing is

[TABLE]

The quantum interference doubles the result. So the probability of scattering is increased, which results in a decrease of conductance. To determine qualitatively the effect of weak localization on the conductance we shall follow a heuristic derivation which can be found in Ref. [27]. The de Broglie wavelength $\lambda_{F}=2\pi/k_{F}$ of the electron determines the scattering cross-section on site $O$. In time $t$ it travels diffusively a distance $\sqrt{Dt}$, where $D$ is the diffusion coefficient. The interference volume in $d$ dimensions becomes $(Dt)^{d/2}r^{3-d}$, where $r$ is the thickness of the system. The electron has to enter the interference volume to experience interference, which occurs with a probability of $\frac{v_{F}\lambda_{F}^{2}dt}{(Dt)^{d/2}r^{3-d}}$. This leads to a relative correction to the conductivity of

[TABLE]

The phase coherence time in the upper limit of the integral shows the condition for phase preservation. Now when we focus on the one dimensional situation for our quasi-one dimensional wire, the evaluation of the integral gives

[TABLE]

For the last expression we used

[TABLE]

If the elastic mean free path is much smaller than the phase coherence length we can neglect this term in the conductivity correction.

[TABLE]

The Drude conductivity can be expressed in terms of the elastic mean free path and the Fermi momentum.

[TABLE]

Substituting this in the relative correction expression, where we use the identities 3.22 and $\lambda_{F}\propto\hbar/p_{F}$ leads to

[TABLE]

To get the correction to the conductance we introduce $\Delta G=\Delta\sigma r^{2}/L$ and arrive at

[TABLE]

This quantum correction to the conductance is known as weak localization and arises due to a self-crossing in the diffusive transport of an electron. However, when $l_{\phi}<<L$ the correction becomes negligible.

At zero temperature the phase of an electron is not broken $(l_{\phi}\rightarrow L)$, so that the correction to the conductance is no longer negligible [27]. This is known as strong, or Anderson, localization. When again the conductance is implemented, the expression for this correction is obtained [6].

[TABLE]

The number of transverse channels available for conduction is determined by the ratio of Fermi wavelength and cross-section. Now the correction is negligible if $\Delta G/G<<1$, which is true for a large number of open conduction channels. Since we consider a diffusive wire, we can look at the distribution of transmission probabilities in equation 2.5, and see that if the average conductance increases the number of open channels increases.

We can conclude that for our quasi-one dimensional diffusive wire, we can neglect the quantum correction to the conductance due to interference effects when we consider wires with length much larger than the phase coherence length and a conductance significantly larger than the conductance quantum.

Chapter 4 Photon absorption and other energy exchange processes in diffusive wires

4.1 Introduction

The model proposed by Remco Schrijvers had the aim to describe the electron energy distribution in a diffusive wire subject to high frequency irradiation with energy relaxation present inside the wire [6]. Unfortunately this aim was not fully achieved. The assumption was made that the path traveled by the electron inside the wire is of no influence to the energy distribution, so that Tien-Gordon theory could be applied to the reservoirs and the distribution inside the wire was described by the Boltzmann equation. However, this turned out to be incorrect since Tien-Gordon theory assumes averaging over time and therefore the collision integral can not be evaluated in the correct manner. Therefore a different approach is required.

A.V. Shytov developed a theoretical framework to calculate the electron energy distribution for wires where the phase coherence time and energy relaxation time exceeds the diffusion time, so that the transport is fully coherent. We derive from Green function formalism an equivalent model. The insight we gain from this derivation is helpfull in the extension of the theoretical framework of Shytov with a term accounting for inelastic scattering, breaking the phase of the electrons. Since Green function formalism is not basic knowledge, the most important parts for our derivation are first shortly explained.

4.2 Green function formalism

The Green function formalism provides a strong calculation method which can be used to calculate a variety of properties of many-particle systems. In mathematics Green functions obey a inhomogeneous differential equation, where the inhomogeneity is singular. As we have seen in the previous chapters, the Schrodinger equation is the central equation in quantum mechanics. Since this is a differential equation the Green functions apply in describing many-body physics in both equilibrium and non-equilibrium situations. The basis of the formalism is the definition of the single-particle Green function by the wave function [28].

[TABLE]

So the Green function is based on the wave function $\psi_{0}$ of the ground state of the system with Hamiltonian $H$ and the time-evolving wave function $\psi_{H}$ of the system which evolves like $e^{iHt/\hbar}\psi(t=0)e^{-iHt/\hbar}$. The time-ordening operator $T$ is defined in such a way that it always moves the operator with the earlier time-argument to the right.

[TABLE]

The sign in the time-ordening is dependent on the nature of the considered particle. For fermions the sign is negative, so that the Pauli exclusion principle is not violated, and for bosons the sign is positive. In the following we shall only consider fermions. The equation of motion is now derived by differentiating the equation for the single particle Green function with respect to $t$.

[TABLE]

From second quantization it is know that the anticommutation of a wave function in the Heisenberg picture with its conjugate gives a delta-function, so that the first term on the right side of the equation of motion is a multiplication of a spatial and a temporal delta-function. For the second term we use the Heisenberg equation of motion $i\hbar\frac{\partial\psi_{H}}{\partial t}=[\psi_{H},H]$. When we consider a particle free of interactions subject to a Hamiltonian $H=-\frac{\hbar^{2}}{2m}(-i\nabla-\frac{e}{\hbar}A(t))^{2}$, where the vector potential $A(t)$ representing an electric field is integrated in the momentum operator by principle of minimal substitution, the equation of motion for the Green function $G_{0}$ of a free particle becomes

[TABLE]

Because the Hamiltonian is time-dependent in the vector potential we are already considering non-equilibrium. When now also a many-particle system is considered where the particles interact with eachother, the picture becomes a bit complicated. The wave functions, and thus the Green functions, are subject to both an external potential and an internal potential. To ease the calculations the operations are contour-ordered. This replaces the time-ordening operator $T$ in equation 4.1 with the contour-ordening operator $T_{C}$ which has the same properties, only not in time, but on the defined contour. Because in non-equilibrium the final state does not have to return to the initial state the contour, on which the particle is defined, lies in the complex plane depicted in figure 4.1. We won’t go into detail on this, but a insightful derivation can be found in Ref.[28] and Ref.[29].

The derivation in Ref.[28] and Ref.[29] is an approach from non-equilibrium statistical mechanics and leads to the Dyson equation for the Green function which consists of the free particle Green function $G_{0}$ and a self energy term responsible for the interactions.

[TABLE]

The complex contour integral in equation 4.5 is rather impractical in calculations. Fortunately analytic continuation provides a method to replace the contour integrals by real time integrals. The Green function is defined by different Green functions on the contour, the lesser and greater Green function, the time-ordered and anti-time-ordered Green function and the advanced and retarded Green function, dependent on the position of the time coordinates of the Green function on the contour. When the initial time $t_{0}$ is set to infinity and the interactions are coupled adiabatically, the complex part of the contour depicted in figure 4.1 vanishes. By doing this one neglect initial correlations, but in many situations the interactions in the process of reaching a steady state will wash out these initial correlations. In highly transient situations it can however cause problems.

When we consider the lesser Green function, which contains the information on the energy distribution, the first time coordinate is on the first half of the contour and the second time coordinate on the second half. The contour can be deformed to form two contours in the limit of initial time going to infinity as indicated in figure 4.2.

When we look at the product $C(t_{1},t_{1^{\prime}})=\int_{C}d\tau A(t_{1},\tau)B(\tau,t_{1^{\prime}})$, the lesser function becomes on the new deformed contour $C^{<}(t_{1},t_{1^{\prime}})=\int_{C_{1}}d\tau A(t_{1},\tau)B^{<}(\tau,t_{1^{\prime}})+\int_{C_{2}}d\tau A^{<}(t_{1},\tau)B(\tau,t_{1^{\prime}})$. The integration on the first contour can run from $-\infty$ to $t_{1}$ and from $t_{1}$ to $+\infty$ and on the second contour from $-\infty$ to $t_{1^{\prime}}$ and from $t_{1^{\prime}}$ to $+\infty$. By doing this all functions can be expressed in lesser functions (for $t_{1}<t_{1^{\prime}}$) and greater functions (for $t_{1}>t_{1^{\prime}}$) and when the relations $G^{a}(1,1^{\prime})=\theta(t_{1^{\prime}}-t_{1})[G^{<}(1,1^{\prime})-G^{>}(1,1^{\prime})]$ and $G^{r}(1,1^{\prime})=\theta(t_{1}-t_{1^{\prime}})[G^{>}(1,1^{\prime})-G^{<}(1,1^{\prime})]$ are used Langreth’s result for analytic continuation is obtained [30].

[TABLE]

In the next section we shall derive from a simplified Dyson equation a quantum diffusion equation. In the subsequent section equation 4.6 is used to derive from the complete Dyson equation a quantum diffusion equation with an interaction term accounting for inelactic scattering.

4.3 Quantum diffusion equation for the distribution function

To calculate the energy distribution function of electrons in a mesoscopic wire biased with an ac voltage induced by THz radiation on the reservoirs whereon the wire is coupled, we derive a quantum diffusion equation from the Dyson equation. First we derive an equation for a situation where inelastic interactions are neglected by neglecting the self energy term in the Dyson equation and introduce instead an elastic interaction term which will lead to a relaxation time approximation to account for the diffusivity of the system [31].

[TABLE]

Here $G(1,1^{\prime})$ is the non-equilibrium Green function of a particle at coordinates $x_{1}$ and $t_{1}$ provided that the particle arises from the coordinates $x_{1^{\prime}}$ and $t_{1^{\prime}}$ defined by equation 4.1. $G_{0}(1,1^{\prime})$ is the Green function of a free particle given by equation 4.4 and $I[G(1,1^{\prime})]$ is the collision term for elastic impurity scattering. By substituting the equation 4.4 in the Dyson equation we can obtain the differential form consisting of the two conjugate parts.

[TABLE]

[TABLE]

These two conjugate parts are subtracted from each other where the two collision terms are re-defined in a single collision term which will later provide the relaxation time approximation for elastic impurity scattering.

[TABLE]

Now the quadratic terms are expanded and we can use the fact that the vector potential is taken only time-dependent, so that according to commutation rules the operation $\nabla A$ is equivalent to $A\nabla$.

[TABLE]

For reasons of convenience we will proceed with this equation expressed in Wigner coordinates defined like:

[TABLE]

To introduce the Wigner coordinates the quadratic parts of equation 4.3 has to be expanded. The summation and difference of the vector potential can be replaced by an representive symbols: $A_{+}(t)=A(t_{1})+A(t_{1^{\prime}})$ and $A_{-}(t)=A(t_{1})-A(t_{1^{\prime}})$. Because we are interested in the distribution function we proceed with the lesser Green function in the equations. The interaction term is now just dependent on the lesser Green function. No analytic continuation procedures have to be followed, because in the end the interaction is given by a relaxation time approximation.

[TABLE]

Here we make the transition to proceed with the distribution function in a momentum representation of equation 4.3.

[TABLE]

The terms containing $\nabla_{r}$ operate first on the integral, so that the operator is replaced by $ip^{\prime}/\hbar$, and the terms are rearranged.

[TABLE]

Then the equation is multiplied by $e^{-ipr/\hbar}$ and a Fourier transform is performed by integrating over all $r$.

[TABLE]

The Fourier transform in $r$ of the exponent creates the delta function $\delta(p^{\prime}-p)$ and the integral over $p^{\prime}$ forces by means of the delta function all $p^{\prime}$ to $p$.

[TABLE]

The sum of the vector potential on time $t_{1}$ and $t_{1^{\prime}}$ modulates the momentum of the charge carrier. This is a second order effect so that the term in front of the momentum part of the equation above can be replaced by the velocity of the charge carrier. The vector potential is defined as $A(t)=U/(L\omega)cos(\omega t)$. The difference term in the vector potential is then expressed in the Wigner coordinates.

[TABLE]

This vector potential is substituted in equation 4.20 and the same procedure is followed for an energy representation as previous done for the momentum representation. A Fourier transform in $t$ is performed and this is integrated over $E^{\prime}$. For the terms without the vector potential this operation is trivial since it just replaces the variable $t$ in the distribution function by $E$. For the part containing the vector potential the situation is a bit more subtle and essential in the understanding of the absorption of energy quanta of the field by electrons. Therefore this is explicitly shown.

[TABLE]

When this is substituted in the kinetic equation and the operator $\nabla_{R}$ is replaced by a derivative with respect to the one-dimensional space coordinate $x$ we arrive at a form from which we can go to a diffusion equation.

[TABLE]

The distribution function can be divided in an odd and an even part with respect to $p$.

[TABLE]

Because the field is considered to be uniaxially symmetric the even part of the distribution function only depends on the absolute value of $p$, so that the even part of the distribution function is the distribution function as function of energy only: $f_{e}(p,x,E,T)=f(x,E,T)$. First the kinetic equation is transformed into two equation for positive and negative momentum.

[TABLE]

[TABLE]

The equations are added and subtracted from each other and divided by 2.

[TABLE]

[TABLE]

The identities of the even and odd part of the distribution function can be implemented and the even part is changed to the distribution function as function of energy only.

[TABLE]

[TABLE]

Now if we only consider inelastic impurity scattering we only have an collision integral acting on the odd part of the distribution function. The impurity scattering can only change the momentum of a charge carrier but can not change the energy. When for this collision integral the relaxation time approximation $I[f_{o}(p,x,E,T)]=-f_{o}(p,x,E,T)/\tau_{im}$ is used and the impurity time is considered to be small the time derivative of equation 4.59 can be neglected. Then $f_{o}(p,x,E,T)$ is just a function of the momentum part times f(x,E,T) times $-\tau_{im}$. When this is substituted in equation 4.58 we arrive at the final form of the quantum diffusion equation, where we take $D=v^{2}\tau_{im}$ the diffusion constant.

[TABLE]

A.V. Shytov also studied the energy distribution of electrons in a diffusive, coherent wire. The equation he used to calculate the distribution function is equivalent to that derived above.

4.4 Limit situations for the simple quantum diffusion equation

The quantum diffusion equation 4.32 can be solved analytically for certain limit situations [32]. Therefore it is convenient to express the equation in dimensionless parameters.

[TABLE]

When we introduce the diffusion time for an electron in the wire $\tau_{D}=L^{2}/D$ equation 4.32 becomes

[TABLE]

This differential equation has for the initial and boundary conditions a Fermi distribution

[TABLE]

The limit situations are defined by the ratio of the field frequency $\omega$ and the diffusion time $\tau_{D}$ and the ratio of the photon energy $\hbar\omega$ and the field energy $eV$.

4.4.1 Slow field limit

For $\omega\tau_{D}<<1$ the field oscillates slowly with respect to the time that the electron travels diffusively through the wire. In equation 4.34 the time derivative can be neglected and the solution is obtained by solving the spatial second order differential equation. The solution becomes

[TABLE]

Following the approach in Ref. [32] we take the Fourier transform in energy domain to find the exponent of the finite difference operator which leads to

[TABLE]

Substituting this in the equation for the distribution equation, restoring dimensions and using time averaging $\overline{J_{2n}(2asin(t))}=J^{2}_{n}(a)$ we arrive at the final general expression for the distribution function.

[TABLE]

So we see close resemblance with Tien-Gordon theory where the probability of absorbing $n$ field quanta is also given by squared Besselfunctions. The resemblance with a dc biased wire is also visable in the pre-factors $1-x/L$ and $x/L$, which gives the number of electrons that enter position $x$ from the right and the left reservoir [33].

When the field energy is much larger than the photon energy, $\hbar\omega<<eV$, the asymptotic form of the Bessel function at $x/\omega\approx n>>1$ may be used, which gives

[TABLE]

where $F_{0}(E,x/L)=\frac{1}{\pi}cos^{-1}(\tilde{E})$ for $|\tilde{E}|<1$ with $\tilde{E}=\frac{LE}{xeV}$. For $\tilde{E}<-1$ the occupation is one and for $\tilde{E}>1$ the occupation is zero. Figure 4.3 shows the distribution in slow field, strong signal limit.

4.4.2 Fast field limit

In the fast field limit the diffusion time is much larger than the reciprocal frequency of the field, $\omega\tau_{D}>>1$. This means that equation 4.34 practically becomes time-independent, since the time-derivative is proportional to $1/\omega\tau_{D}$. Averaging equation 4.34 over the field period, leads to an equation for the time-averaged distribution function.

[TABLE]

In the limit $\hbar\omega<<eV$ the finite difference operator $D_{\omega}$ can be replaced by the partial energy derivative $\partial/\partial\epsilon$. This makes equation 4.42 become a Laplace equation in a two-dimensional strip defined by $0<x<1$ and $-\infty<\epsilon<\infty$. This strip can be conformally mapped by the function $w=exp[\pi i(x+i\sqrt{(}2)\epsilon)]$ onto the half-plane $Imw>0$ [32] [34]. The boundary condition on the line $Imw=0$ at zero temperature is set to be $\bar{f}_{\infty}(w)=0$ for $|Rew|<0$ and $\bar{f}_{\infty}(w)=1$ for $|Rew|>0$. The imaginary part of the analytic function gives the solution of this boundary value problem.

[TABLE]

When the original dimensional units are restored the final expression for the time-averaged distribution function is

[TABLE]

In the fast field the energy distribution does not have to go to zero at high energies. The energy gained from the field is not limited by $eV$. Instead an electron has a finite probability of oscillating several times back and forth with the field in the wire before leaving the wire, thereby gaining multiple energy quanta of the field which sum exceeds $eV$. Figure 4.4 shows the electron energy distribution in the fast field, strong signal limit.

4.5 Incorporating inelastic interactions

So far only coherent transport is considered. If the length of the wire is extended in such a way that the diffusion time becomes of the same order as the phase coherence time and energy relaxation time this simple model breaks down. Therefore this model has to be extended to account for electron-electron and electron-phonon interactions. This is done by evaluating the complete Dyson equation 4.5, where we isolate the collision term for the elastic impurity scattering that is treated with a relaxation time approximation in the same way as before.

[TABLE]

By substituting the equation of motion for the Green function of a free particle we obtain again two conjugate equations.

[TABLE]

[TABLE]

As before we are interested in the distribution function, so we concentrate on the lesser Green function by an analytic continuation of the above functions where we concentrate on the self energy part of the functions. The remaining part of the equations in the derivation is similar to the derivation without inelastic interactions.

[TABLE]

[TABLE]

These two equations are subtracted from each other.

[TABLE]

Now the following identities are introduced to gain insight in the derivation [28].

[TABLE]

The terms are arranged so that everything is expressed in commutators and anti-commutators.

[TABLE]

To simplify the calculations we assume the scattering to be local in space, so that the integral operation over $y$ forces the integration variable towards the central space coordinate. Also we can make the assumption of weak interactions, so that we can apply the quasi-particle approximation. Because we also used a gradient expansion of the potential, the first two commutators of the above relation are second order and can be neglected. Basically this means that the density of states of the quasi-particles in the wire is not affected by the vector potential nor the interactions.

[TABLE]

When we now also assume that the scattering is instantaneous, the integral over $\tau$ forces the integration variable towards the second time variable of the self energy in the product of the self energy and the Green function.

[TABLE]

As we assume a slow variation of the Green function induced by the vector potential and we assume the interactions to be weak, we can state that the effect of the self energy on time $t_{1}$ is the same as that at time $t_{1^{\prime}}$. So the self energies at $t_{1}$ and $t_{1^{\prime}}$ can be replaced by a single self energy $\Sigma(x_{1},t_{1};x_{1^{\prime}},t_{1^{\prime}})$.

[TABLE]

By applying the Wigner transformation to this collision term, the product of the self energies with the Green functions can be interpreted as the imaginary in- and out scattering rates $i\hbar\Gamma_{e,h}$ with the electron and hole distribution [35] [36].

[TABLE]

This is again multiplied by $e^{-i(Et-pr)/\hbar}$ and integrated over $t$ and $r$ leading to the final form of the total collision term due to inelastic scattering where this is multiplied by $i$ from the rest of the equation.

[TABLE]

The two parts of the quantum diffusion equation are again connected.

[TABLE]

Same procedure is followed to come to a diffusion equation as for elastic impurity scattering. The equation is divided in an even and odd part, where the elastic impurity scattering only contributes to the even part and the inelastic interactions contribute to the odd part.

[TABLE]

[TABLE]

Taking the same relaxation time approximation $I[f_{o}(p,x,E,T)]=-f_{o}(p,x,E,T)/\tau_{im}$ for the impurity scattering leads to the desired quantum diffusion equation.

[TABLE]

So we see that the quantum diffusion equation 4.32 is extended with a term that controls the in- and outscattering of quasi-particles at energy $E$ due to inelastic collisions. These inelastic collisions could be due to the interaction between two quasi-particles or due to the interaction between a quasi-particle and a phonon. In the next section we will derive expressions for these interactions.

4.6 Inelastic scatterering

The main energy relaxation mechanisms are electron-electron 111the electron is in fact a quasi-particle and electron-phonon scattering and the sum of these contributions give the total interaction term.

[TABLE]

Both collision terms have an inscattering and outscattering term as seen in equation 4.56. A quasi-particle with energy $E$ has an collision term

[TABLE]

The collision terms due to electron-electron scattering and electron-phonon scattering can be calculated independently of each other. First we will tread the interaction between electrons and phonon. Subsequently we look at the interactions between electrons.

4.6.1 Electron-phonon interaction

Let’s first focus on the electron-phonon interactions. To begin some assumptions have to be made. When we only want to consider acoustic phonons with a dispersion relation between energy and wave vector $\epsilon_{k}=\hbar sq$, with $s$ the sound velocity, the phonon temperature $T_{ph}$ has to be small compared to the Debye temperature $T_{D}$. Further the electronic wave functions can be approximated by plane waves, which is justified by the fact that electron-phonon coupling is only relevant for higher energies and from the dispersion relation it is seen that large wave vectors are associated with these energies. Then it is probable that the electronic mean free path is larger than $1/q$. Also the electron-phonon coupling is given by a scalar deformation potential, so only the longitudinal phonons are coupled on the electrons. The matrix element describing the interaction simplifies to $|M(q)|^{2}=|M_{0}|^{2}q/V$, where $|M_{0}|^{2}$ is geometry independent. This only is valid for spherical Fermi surfaces [37].

The transition of an electron to a state with energy $E$ can either be due to the absorption or the emission of a phonon. The same can be said of the transition out of the state with energy $E$. We can define the transition due to absorption by $W^{-}$ and the transition due to emission by $W^{+}$. Further we know that the state from which the particle departes has to be occupied and the state in which the particle arrives has to be unoccupied. The latter is a direct consequence of the fact that we look at fermions and according to the Pauli exclusion principle a state can only be occupied by a single fermion. This leads to the following collision terms [38].

[TABLE]

[TABLE]

Here $n_{ph}$ represents the Bose energy distribution of the phonons, $n_{ph}(E)=(exp(E/kT)-1)^{-1}$. The transition probabilities are given by Fermi’s Golden Rule [38].

[TABLE]

To obtain the collision rate at which an electron with wave vector $k$ emits or absorbs a phonon of energy $E_{|k-k^{\prime}|}$ the equations 4.6.1 and 4.6.1 have to be summed over $k^{\prime}$ with $E(k-k^{\prime})$ fixed. A detailed derivation can be found in Ref. [38].

[TABLE]

[TABLE]

The so called Eliashberg function $\alpha^{2}F(\epsilon)$ is dependent on the coupling between the electrons and phonons. In Ref. [39] this function is determined to be

[TABLE]

Here $|M|^{2}$ is the matrix element depending on the defined deformation potential and $N(0)$ is the electronic density of states at Fermi level. The precise microscopic form of $|M|^{2}$ is dependent on the details of the lattice structure. Therefore in Ref. [39] they present this matrix element in terms of a measurable quantity $\Sigma$ related to the power dissipated to the lattice of volume $V$ by $P=\Sigma VT^{5}$. A detailed form of the electron-phonon interactions and the temperature dependence in disordered conductors can be found in Ref. [40].

4.6.2 Electron-electron interaction

The interaction between quasi-particles is due to the Coulomb potential of the particles. This Coulomb interaction is screened by an effective medium build from all the electrons in the metal. Altshuler et al. showed that multiple scattering events due to disorder in the system reduces the lifetime of the quasi-particle [41]. At zero temperature the lifetime of a particle obeying Fermi statistics in state $|\alpha>$ with energy $\epsilon_{\alpha}$ above Fermi level that interacts with a particle in state $|\gamma>$ with energy $\epsilon_{\gamma}$ directly follows from Fermi’s Golden Rule [42].

[TABLE]

$U$ is the interaction potential from which the states $|\alpha>$ and $|\gamma>$ evolve in the states $|\beta>$ and $|\delta>$. This lifetime has to be averaged over all states having energy $\epsilon$ in order not to single out a give state.

[TABLE]

When the energy of the states $|\gamma>$ are denoted by $\epsilon^{\prime}$ and the energy exchange involved in the scattering is $\omega$, energy conservation leads to energies of the final states $|\beta>$ and $|\delta>$ of $\epsilon-\omega$ and $\epsilon^{\prime}+\omega$. This is depicted in figure 4.5.

Considering all possible initial states $|\gamma>$ leads to integration over $\epsilon^{\prime}$ and $\omega$.

[TABLE]

Now when the requirements of zero temperature and the Fermi statistics are dropped, this approach still holds when we include the occupation numbers of the states in the obtained result 4.72.

[TABLE]

Where

[TABLE]

To complete the collision term for electron-electron interactions it is convenient to let go the notation of Ref. [42] and proceed with the notation used for electron-phonon interactions. We define the kernel $K(\epsilon)$, which follows from $(4\pi)/(\hbar\nu_{0})W^{2}(\omega)$. Further the collision rate can be splitted in the inscattering and outscattering term by multiplying the first part by $f$ and the second part by $1-f$.

[TABLE]

[TABLE]

In Ref. [41] and Ref. [37] the matrix element of the transition in a disorded medium is calculated. Here we will not follow the complete derivation, but directly look at the result for the kernel $K(\epsilon)$.

[TABLE]

The bare Coulomb potential $U_{0}(\textbf{q})$ and the polarizability $\Pi(\textbf{q},\epsilon/\hbar)$ of the electron fluid determines the screened Coulomb potential $U_{\epsilon/\hbar}(\textbf{q})$ effectively experienced by the quasi-particles.

[TABLE]

where

[TABLE]

In a metal the density of states $\nu_{F}$ is so large (order of $10^{47}J^{-1}m^{-3}$) that the polarizability dominates the denominator in the expression of the screened Coulomb potential. Therefore equation 4.78 simplifies to

[TABLE]

and the total kernel becomes

[TABLE]

If we consider a metallic wire with cross-section $S=wt$, where $w$ is the width and $t$ is the thickness of the wire, only the uniform modes in transverse dimensions contribute to $K(\epsilon)$ if the energies $\epsilon$ are smaller than $\hbar D/max(w^{2},t^{2})$. This leads to

[TABLE]

This derivation leads to a difference with the result for the screened Coulomb interactions obtained by Kamanev and Andreev [43]. They found $K(\epsilon)$ to be a factor 2 larger. Experiments showed that the energy dependence of the collision term is accurate, but the intensity is off. A discussion can be found in Ref. [26] and Ref. [44].

4.7 Summary

In this chapter we used the fact that the electrons involved in the ac quantum transport in a diffusive wire can be described as quasi-particles according to the Fermi liquid theory. For coherent transport the energy distribution of the quasi-particles obeys a relative simple quantum diffusion equation. The non-equilibrium in a mesoscopic, diffusive wire induced by a time-dependent field manifests itself in the energy distribution. When the length of the wire is extended, the transport becomes incoherent and the redistribution of energy among the quasi-particles has to be evaluated. For this reason the relative simple quantum diffusion equation is extended with a collision integral accounting for electron-electron and electron-phonon interactions. In the next chapter the model is evaluated using numerical calculation methods.

Chapter 5 Numerical results

5.1 Introduction

The model developed in chapter 4 allows the evaluation of the quasi-particle energy distribution in a mesoscopic wire ac biased with irradiation. For very short wires, where the phase coherence time and energy relaxation time exceed the diffusion time, the transport is fully coherent and the distribution function in the wire is never an equilibrium function. The non-equilibrium description is quite different in the two field limits, $\omega\tau_{D}<<1$ and $\omega\tau_{D}>>1$, as discussed in section 4.4. In the slow field limit ($\omega\tau_{D}<<1$) the quasi-particle energy distribution is varying in time, following the oscillation of the field instantaneously. In the limit $eV>>\hbar\omega$ this shows close resemblance with the dc biased wire and the quasi-particle energy distribution is given by a two step function which varies in time. The fast field limit ($\omega\tau_{D}>>1$) is quite different. In this limit the quasi-particle energy distribution is given by a time-independent multiple step function. For energies $eV>>\hbar\omega$ the steps smooth out and a continuous function is obtained which provides a finite probability of finding a quasi-particle far from the Fermi energy.

To evaluate the slow field regime and fast field regime we can define some ratio $\hbar\omega/eV$ and vary the product $\omega\tau_{D}$. Shytov showed that the crossover from low-frequency behavior to high-frequency behavior occurs at $\omega\tau_{D}\approx 100$ [32]. This is due to the fact that that the quasi-particle energy distribution relaxes at $t\rightarrow\infty$ as $exp(-\mu t)$, where $\mu=\pi^{2}/\tau_{D}$ is the lowest non-zero eigenvalue of the diffusion operator. It is reasonable to assume that the crossover occurs when the relaxation time is of the order of the field period, $2\pi/\omega$. So the crossover is estimated to occur at $\omega\tau_{D}\propto 2\pi^{3}\approx 62$, which is close to 100.

This theoretical research is done in an experimental research group. The strong connection with experimental physics leads to the desire to evaluate the model for realistic situations (THz frequencies and field amplitudes of 1-20 meV), so that when an experimental setup is realized the model can provide the understanding of the experimental results. We apply these conditions in the evaluation of equation 4.60 using numerical calculation methods. The equation is expressed in dimensionless parameters in the same way we did for the discussed limit situations of coherent transport.

[TABLE]

As explained in the previous chapter, the collision term can be neglected for fully coherent transport. The energy of the quasi-particles in the wire is only affected by photon absorption and the diffusive transport itself. This means that the neglect of the collision term is only valid for short wires. To make this somewhat more quantitative, we consider the phase coherence time of a quasi-particle. Two phase breaking mechanisms are distinguished, electron-phonon interaction and electron-electron interaction. The experimental part of the research focuses on aluminum wires with a diffusion coefficient of about 100 cm2s*-1* measured at liquid helium temperatures, so that we first concentrate on this material and temperature. Above temperatures of 1 K the phase breaking mechanism is electron-phonon interaction. The phase coherence time is approximated by [26]:

[TABLE]

Here $E_{F}$ is the Fermi energy, $N(0)$ is the density of states at Fermi energy, $\rho$ is the mass density, $s$ is the speed of sound and $k_{F}$ is the Fermi wave vector. The phase coherence time at 2 K, which can be achieved in a pumped liquid helium cryostat, is approximately 10 ns. This is equivalent to a wire of length $L=\sqrt{D\tau_{\phi}}=10\mu$m. So for wires shorter than this length the transport is coherent. Since this is an approximation we decided to use in our calculations wires of maximum length of 7 $\mu$m, with a diffusion time of 5 ns, to be certain that the transport is coherent. We evaluate the quasi-particle energy distribution for coherent transport at 2 K from the slow field regime to the fast field regime. For the slow field regime $\omega\tau_{D}=1$ we choose a wire of 56 nm and a field frequency of 0.5 THz. In the fast field regime $\omega\tau_{D}=30000$ we take a wire of 7 $\mu$m and a field frequency of 2 THz. A wire of 400 nm and a field frequency of 1 THz makes the evaluation of the intermediate regime $\omega\tau_{D}=100$ possible. We define the ratio $\hbar\omega/eV=0.4$ for all regimes, so the field amplitude varies from 5 meV in the slow field regime to 20 meV in the fast field regime.

For extended wires, the diffusion time can exceed the energy relaxation time, so that the transport is incoherent. In this report we will focus on the fast field regime for incoherent transport. The slow field regime is already quite well understood [4] [6] and thereby the frequency of the field should be extremely low to have a small product $\omega\tau_{D}$, where $\tau_{D}$ should be of the order of $\tau_{E}$ the energy relaxation time.

The effect of electron-phonon interactions is evaluated at a temperature of 2 K. We determined the phase coherence time for electron-phonon interaction to be 10 ns. This will be our reference in defining the ratio between diffusion time and energy relaxation time, since the energy relaxation time is of the same order of magnitude as the coherence time. The intensity of the interaction between quasi-particle and phonon can be calculated from equation 4.68. In Ref. [39] the quantity $\Sigma$ related to the power dissipation is given to be about 1 GWm*-3K-5*. This brings the intensity to about 2 ns*-1meV-3*, which is close to the empirical intensity of 4 ns*-1meV-3* which followed from experiments done by Huard et al. [26].

Below 1 K the situation becomes a bit complicated, since aluminum is no longer a normal metal, but has experienced a phase transition to the superconducting phase. We proceed below 1 K with an undefined material with the same diffusion coefficient of 100 cm2s*-1*, so that we can evaluate the effect of electron-electron interactions on the energy distribution of the quasi-particles. The phase coherence time for electron-electron interactions is approximated with [26]:

[TABLE]

Here $k_{e}e$ is the prefactor in the kernel of equation 4.82 and given by $\left(\sqrt{2D}\pi\hbar^{3/2}N(0)S\right)^{-1}$. A temperature of 500 mK and a cross-section of the wire of 400 nm2 gives an intensity of the interactions of 0.8 ns*-1meV-1/2*. The empirical intensity found by Huard et al. in silver is 0.4 ns*-1meV-1/2*, so this can be used as a realistic value. This leads to a phase coherence time of about 1 ns at 500 mK. The energy relaxation time is of the same order of magnitude.

So the purpose of this chapter is dual. First, we want to investigate the quasi-particle energy distribution for coherent transport and how the slow field regime differs from the fast field regime. Second, we want to investigate the quasi-particle energy distribution for incoherent transport in the fast field regime and how weak interactions are distinguished from strong interactions.

5.2 Calculation method

Numerical calculation principles allow the evaluation of the quantum diffusion equation 5.1 [45]. We use Euler’s method using finite difference approximations for the space, time and energy variables. By iterating the calculation a stable solution for the quasi-particle energy distribution is obtained. This iteration is performed on the time variable, so that every time step $dt$ results in a new function which arises from the old function and the non-time operation part of the equation:

[TABLE]

where

[TABLE]

Here $D_{x,E}$ is the finite difference operator for space and energy, respectively, and $m$ is the number of iteration. For the diffusion equation without inelastic scattering the terms $Iin$ and $Iout$ in equation 5.5 disappear. The finite difference operators are sparse matrices, which means that the percentage of zero elements greatly exceeds the percentage of non-zero elements and their distribution is such that it is advantageous to use this for a more efficient calculations. The MATLAB function sparse provides the possibility to exploit the sparse nature of the operator. What this function does is isolate the non-zero elements, so that only these elements are used in the calculation. The MATLAB code of the simulation program can be found in the appendix.

5.3 Simulation of realistic coherent and incoherent transport situations

5.3.1 Coherent transport in diffusive wires

As explained in the previous section, the crossover from low-frequency behavior to high-frequency behavior occurs at $\omega\tau_{D}=100$. So when we want to evaluate the slow field regime it is sufficient to have a product $\omega\tau_{D}=1$, which is two orders of magnitude below the crossover. As said before this is based on realistic values, but throughout this section we will only work with relative values. The amplitude of the field is such that $\hbar\omega/eV=0.4$. In the slow field regime the energy distribution is highly time-dependent. The time-averaged distribution function is given in figure 5.1 for three different positions in the wire. The full space dependency is shown in the three dimensional figure in appendix C.

The time-dependency in the slow field regime is evaluated by running the simulation during two field periods and plot the normalized occupation at the photon steps. This normalization is performed by taking the value for every iteration on the first and second photon step, determine the maximum value during the iteration process and divide the value determined in every iteration by this maximum value: $|n_{i}|=f_{i}(E+\hbar\omega/2,E+3\hbar\omega/2)/max(f_{i}(E+\hbar\omega/2,E+3\hbar\omega/2))$. The first photon step immediately follows the field, where the second photon step shows a slight delay as shown in the up left picture in figure 5.4.

For the fast field regime we take $\omega\tau_{D}=30000$, two orders of magnitude above the crossover from low-frequency behavior to high-frequency behavior. The amplitude of the field is again defined so that $\hbar\omega/eV=0.4$. The simulation of this situation is time-averaged depicted in figure 5.2 at three different positions in the wire, where we average over a large number of periods to obtain the final time-independent distribution. The full space dependency is shown in the three dimensional figure in appendix C. To evaluate the time-dependence in the fast field regime the normalized value of the occupation in the three photon steps in the distribution is plotted during the evolution of the function. It appears that in the fast field regime the energy distribution indeed becomes time-independent as shown in the up right picture in figure 5.4 and the occupation at the photon energies is maximum when the diffusion time is reached.

The intermediate frequency regime where the crossover occurs from low-frequency behavior to high-frequency behavior is evaluated at $\omega\tau_{D}=100$. The time-averaged distribution function at three different positions is given in figure 5.3. The full space dependency is shown in the three dimensional figure in appendix C. The time-dependency in the distribution function is drastically decreased at the crossover as seen in figure 5.4. When the diffusion time is reached, the occupation at the photon energies is at maximum with a slight oscillatory deviation with the field period.

5.3.2 Incoherent transport in diffusive wires

As seen in the previous section the absorption of field quanta in the short wire, where the transport of quasi-particles is fully coherent, induces a staircase structure in the quasi-particle energy distribution. Now we want to investigate what happens when the wire is extended, so that photon-absorption is no longer the only mechanism that affects the energy distribution, but also interactions between quasi-particles and between quasi-particles and phonons come into play. It appears that these interactions redistribute the quasi-particles with respect to the energy, so that the occupation of the energy levels is changed with respect to the occupation in the coherent situation. The effect of the two phase breaking mechanisms is quite different. We expect that the interactions between quasi-particles cause a smearing in the staircase structure, while the interactions between quasi-particles and phonons cause the annihilation of the photon steps and finally, in the strong interaction limit, leave a Fermi function with the bath temperature. In this section we will limit ourself to the fast field regime which is, as explained in the previous section, the most interesting domain.

Electron-phonon interactions

Since the first experiments are planned to be done at liquid helium temperatures, we first focus on the effect of electron-phonon interactions on the energy distribution of the quasi-particles. We can distinguish different interaction regimes. The weak interaction regime is found for $\tau_{D}\approx\tau_{E}$ and the strong interaction regime is found for $\tau_{D}>>\tau_{E}$.

Let’s first look at the weak interaction regime where the diffusion time is of the order of the energy relaxation time, so we define $\omega\tau_{D}=50000$ and $\tau_{D}\approx\tau_{E}$. What we expect is that the energy gained from the field by a quasi-particle is redistributed, where the photon steps due to the absorption of multiple field quanta are first influenced. The result of this simulation is shown in figure 5.5. The expected disappearance of the photon steps due to the absorption of multiple field quanta is indeed observed and the transition in the first photon step is smoothed. The space dependency is shown in the three dimensional figure in appendix C.

When the length of the wire is increased, the diffusion time becomes much higher than the energy relaxation time. So we enter the strong interaction regime and define $\omega\tau_{D}=10^{7}$ and $\tau_{D}\approx 200\tau_{E}$. In the fast field regime we expect a Fermi function at the bath temperature on every position in the wire. The result, shown in figure 5.6, indeed shows a Fermi function. There is a slight deviation from the Fermi function with the bath temperature.

When we look at the deviation of the calculated energy distributions from the equilibrium function at bath temperature we see what the effect of weak and strong interactions is. For weak interactions the deviation is clearly defined by the photon energy, but the photon step is smoothed. For strong interactions the deviation is no longer defined by the photon energy and the width and height of the peak is small. The height of the peak is in both situations however for energies below Fermi energy somewhat larger. This observed deviation is probably caused by the discretization of the variables and the fact that the interactions are calculated after $n$ iterations, instead for each iteration, to increase calculation speed. It is reasonable to believe that this has no physical meaning, but is just some numerical error which can be solved by solving the equations with a program written in C. This should provide a much higher calculation speed, so that the discretization can be optimized.

Electron-electron interactions

The effect of electron-electron interactions is quite different from the effect of electron-phonon interactions and dominant for lower temperature so we will evaluate this effect at a temperature of 500 mK. Lets first look at the effect of weak interactions when the diffusion time is of the order of the relaxation time, $\tau_{D}\approx\tau_{E}$ and $\omega\tau_{D}=10000$. Figure 5.8 shows that for weak interactions there is some smearing, but the photon steps are still good defined. The full space dependency is depicted in the three dimensional figure in appendix C.

When we increase the length of the wire, the time that the electron spends traveling through the wire increases also and the effect of the interactions becomes more significant. Figure 5.9 shows that the interactions indeed are more relevant for longer wires and the smearing in the photon steps is clearly visible.

By further increasing the length of the wire, the energy relaxation rate becomes dominant with respect to the diffusion time. The energy gained from the electric field is redistributed in a Fermi function with an effective temperature, so that an effective temperature profile arises across the wire analogously to the dc biased macroscopic wire evaluated in the introduction. In figure 5.10 three Fermi functions are given at different positions in the wire for $\omega\tau_{D}=2000000$ and $\tau_{D}\approx 200\tau_{E}$ at 500 mK. Figure 5.11 shows the effective temperature profile across the wire obtained by fitting the energy distribution on every position in the wire to a Fermi function using a least square method. The distribution function on every position in the wire is depicted in the three dimensional figure in appendix C.

Figure 5.12 shows the theoretical prediction for the effective temperature profile given by $T_{e}(x)=\sqrt{T^{2}+\frac{x}{L}\left(1-\frac{x}{L}\right)3/\pi^{2}\left(e/k_{b}\right)^{2}V^{2}}$ when the voltage across the wire is taken to be $V=\hbar\omega/e$. It shows resemblance with the effective temperature profile resulting from the simulation. It seems plausible to say that the observed difference is due to numerical inaccuracy. Another possible reason for this observed deviation could come from the fact that we look at the fast field regime and the position of photon absorption is responsible for the difference. However, a closed statement on this calls for further study.

In figure 5.13 the deviation from the equilibrium function at bath temperature is shown. It appears that the smearing induced for weak interactions, $\tau_{D}\approx\tau_{E}$, causes the smooth deviation of two subtracted Fermi functions at different temperature for strong interactions, $\tau_{D}>>\tau_{E}$.

Chapter 6 Probing the quasi-particles energy distribution

6.1 System description

Obviously it is desirable to be able to obtain experimentally the quasi-particle energy distribution in a diffusive wire subject to high frequency irradiation to test the model. To this end a system is designed which should provide this possibility. An aluminum diffusive wire is connected to large aluminum reservoirs where the quasi-particle can relax to equilibrium.

These equilibrium reservoirs are designed in such a way that they also function as antenna for the high frequency radiation. This means that the thickness of the metal reservoirs should be larger than the penetration depth of the radiation. The penetration depth can be calculated by $\delta=\sqrt{\frac{2}{\omega\mu\sigma}}$ [46]. Here $\omega$ is the radial frequency, $\mu$ is the permeability and $\sigma$ is the conductivity. The penetration depth of radiation with a frequency of 1 THz is in aluminum approximately 82 nm [46]. Based on this number the reservoirs are designed 100 nm thick.

The antenna design is not the subject of this research and since it is a very sophisticated field of science, we did not spend time on calculating field profiles. Instead we used a bow-tie antenna design, so that imperfections in the design matter less due to the broadband character. The bow-tie antenna is designed in such a way that it is self-complementary [47], by designing the triangles of the bow-tie with a 45 degrees angle with respect to the wire [48]. In this way the system is frequency independent and we are not limited by a capacitive effect. The length of the triangular sides of the antenna are chosen in multiples of the wavelength of irradiation. By choosing the frequency of irradiation we take into account that for pronounced photon steps in the energy distribution the photon energy $\hbar\omega$ should exceed the thermal broadening $k_{b}T$.

We use the same probing method as Pothier et al. did for their experiments on the quasi-particle distribution function in a dc biased wire. On top of the aluminum wire, an insulating layer is positioned whereon superconducting probes of niobium are placed perpendicular, so that when polarized light is used the probes are not affected. In this way NIS junctions are formed and from differential conductance measurement the quasi-particle energy distribution can be obtained. The photon energy should not exceed the gap energy of the used superconductor, otherwise cooper pairs will be broken and trouble the probing of the energy distribution. Figure 6.1 shows a schematic overview of the system where two probes are connected on the wire with one connection to set the current through the NIS junction and one connection to measure the voltage drop over the junction. The THz radiation coupled by the antennas causing the ac bias on the aluminum wire is given by $Vcos(\omega t)$.

First measurement were planned at liquid helium temperature (pumped 2 K - unpumped 4.2 K). The phase coherence time can be calculated using equation 5.2 since the dominant phase breaking mechanism is electron-phonon interaction. This approximation is used to define the length of the wires ranging from 500 nm for coherent transport to 100 $\mu$m for fully incoherent transport. When we compare these values with the values stated in the introduction of the previous chapter, we see that we are not measuring in the slow field regime. This has its origin in the fact that we are interested in fields with THz frequencies.

6.2 Fabrication

The fabrication of the samples is not yet proven to be fully successful. The samples to test the NIS junctions however are functioning well enough to conclude that with that recipe the NIS junctions provide the ability to probe the distribution function in a wire driven out of equilibrium. Figure 6.2 shows a scanning electron microscope (SEM) picture of a NIS junction.

The fabrication of the structures described is done in three steps, in which the nanowires, the antenna and the tunneling probes are defined. The different structures are aligned by means of markers. To define the Al nanowires a double resist layer is patternd with electron beam lithography. The Pmma 950k/Pmma 495k resist is developed in MIBK:IPA for 60 seconds and rinsed for 30 seconds in IPA. A 20 nm thick Al film is evaporated at a rate of 1 angstrom/sec at a pressure of 1e-7 mbar, and lift off is done in hot acetone. The procedure for the definition of the antennas is similar to the nanowires, except for a cleaning step prior to the deposition of the antennas, and an Al film which is 90 nm thick instead of 20 nm. To fabricate the tunneling probes the sample is cleaned for six minutes in an Argon plasma, after which it is oxidized for 40 minutes in a pure Oxygen atmosphere of 1 mbar. A 80 nm thick Niobium film is sputtered in situ, at a pressure of 8e-3 mbar and a rate of 1 nm/sec. An etch mask is created using SAL resist and electron beam lithography, using MF-322 for developing. The Niobium is subsequently etched in an SF6/O2 plasma for 5 minutes with end-point detection. The remaining resist is removed in PRS 3000 resist stripper and the sample is cleaned in acetone.

Figure 6.3 shows SEM pictures of two samples with wires of different length. Measurements on most recent samples showed an improvement of performance. The outlook for experimental results becomes more promising.

6.3 Measurements on NIS junctions

6.3.1 Differential conductance of a NIS junction

The energy distribution function of the quasi-particles in the wire can be obtained from conductance measurements of the NIS junctions of the superconducting wires on top of the mesoscopic wire of interest. In appendix D the differential conductance of a NIS junction is derived. The differential conductance is a convolution of the distribution function in the normal wire and the density of states in the superconducting probe. So the unknown distribution function can be obtained by deconvolve the differential conductance with the BCS density of states.

[TABLE]

The differential conductance of the NIS junction is only usefull for probing the energy distribution when the material used for the superconductor is indeed superconducting at liquid helium temperatures. The phase transition of niobium is measured by doing a RT-measurement using a dipstick. The sample is mounted in the vacuum tube of the dipstick with a heating resistance connected. By applying a current to this heating resistance the temperature of the sample is increased from 4.2 K to the desired temperature above the critical temperature. So when simultaneously the resistance of a niobium wire is measured using a four point measurement we can obtain the RT-characteristic. The four point measurement is done by setting a current bias to the niobium wire using a current source and measure the voltage drop across the wire. It appears that the niobium indeed becomes superconducting at the expected critical temperature of 9 K. A result of a RT-measurement is shown in figure 6.4.

To test the NIS junctions, differential conductance measurements are performed for wires (not connected to equilibrium reservoirs) with the superconducting wires on top forming the NIS junction. The measurements are performed at liquid helium temperature (4.2 K) using a dipstick, where again the sample is mounted in a vacuum tube. A current source is used to apply a current bias to the NIS junction. By measuring the voltage drop over the junction in a four points measurement setup, the IV characteristic is obtained. Using a lock-in amplifier the differential conductance of this characteristic is determined. The distribution function can be obtained from such measurements on the NIS junction because of the non-linear IV behavior. Figure 6.5 shows a measurement of the IV-curve and the dI/dV-curve of a NIS junction.

The deconvolution is executed using a steepest descent method. This method is commonly known for the use in minimizing functions [49]. This is used to deduce the distribution function from the differential conductance using equation 6.1. First the effective density of states of the superconductor in the NIS junction is deduced from a $dI/dV$ measurement on a wire in equilibrium, so that fits for the energy gap, tunneling resistance and electron temperature can be implemented. Then an initial distribution function can be chosen in such a way that the calculation converges in a relative short time [50]. The initial distribution function is used to calculate the differential conductance and this is compared to the measured differential conductance for a wire out of equilibrium. If the difference is equal or smaller than the desired precision the deconvolution is completed and the initial distribution is the distribution of the electrons in the wire out of equilibrium. If the difference between the two $dI/dV$’s is larger than the desired precision a new distribution is calculated using the square deviation:

[TABLE]

The occupation probability at each energy in the distribution is incremented by the partial derivative of the square deviation with respect to the occupation factor $f^{k}$ at that energy:

[TABLE]

This calculation is iterated until the desired precision is achieved. The procedure is illustrated in figure 6.6. The deconvolution is executed in MATLAB. The script that is used is given in appendix E.

6.3.2 NIS differential conductance for ac biased wire

The fabrication of samples with wires connected to antennas and superconducting probes on top of the wires is not yet proven to be successful. Therefore the proposed model of chapter 4 cannot yet be verified. We can however calculate what the expected differential conductance measurements will look like when we measure it while driving the wire out of equilibrium with a time-dependent electric field. We do this by taking a calculated distribution function from the model and calculated dI/dV with equation 6.1. This is done for the fast field regime where $\omega\tau_{D}=30000$ and $\hbar\omega/eV=0.4$ at a temperature of 2 K. So we look at a coherent transport situation comparable to that discussed in chapter 5, where we take into account that the photon energy should be below the gap energy of niobium. The distribution function from the model and the calculated differential conductance of the NIS junction on the wire are given in figure 6.7.

The experiments with the described system will proceed and hopefully lead to a satisfying result that can be used to verify the theoretical model.

Chapter 7 Conclusion and discussion

7.1 Conclusion

In this project we studied the ac quantum transport in a quasi-one dimensional, normal metal wire, where the transport is diffusive, connected between equilibrium reservoirs. For coherent transport, where the phase of a charge carrier is preserved, photon absorption and the diffusive character of the transport influence the energy distribution of the quasi-particles inside the wire. When the diffusion time, i.e. the time that a quasi-particle spends in the wire, exceeds the energy relaxation time, the mutual interaction of quasi-particles and the interaction between quasi-particles and phonons causes incoherent transport and influences the energy distribution.

Often scattering theory is used to describe the transport of charge carriers through nano-structures. However, in our situation many processes involving the energy of the charge carriers come into play in the scattering region, i.e. the diffusive wire. Therefore we studied the effect of an ac bias applied to a diffusive wire by looking at the energy distribution of the quasi-particles inside the wire.

Previous work on this subject was independently done by R. Schrijvers [6] and A.V. Shytov [32]. R. Schrijvers approached the situation by applying Tien-Gordon theory to the reservoirs and calculated the energy distribution in the wire with a semi-classical diffusion equation, which is a Boltzmann equation extended for inelastic interactions. He concluded that in this way only the slow field limit without inelastic interactions is adequately described. This is due to the fact that the assumption is made that the path traveled by the charge carrier in the wire is of no influence on the energy distribution. Also Tien-Gordon theory assumes averaging over time, which troubles the correct evaluation of the collision integral.

A.V. Shytov calculated the energy distribution in a diffusive wire where the energy relaxation time exceeds the diffusion time with a quantum diffusion equation. This approach seems valid in all frequency regimes, but only for coherent transport.

Stimulated by these approaches we derived from Green function formalism a quantum diffusion equation equivalent to that of Shytov and we did a second derivation from the full Dyson equation to extend this model to account for inelastic scattering processes. This approach seems successful, as the quasi-particle energy distribution can be calculated in every frequency regime and the limit situations of strong interaction processes provide the correct distribution. This approach is justified by the Landau theory of Fermi liquids. This states that there is an one-to-one correspondence between the states of a non-interacting particle system and the states of an interacting particle system provided that the excitations are near the Fermi level. This means that the excited states of an interacting system are labeled with the same quantum numbers as those of a non-interacting system. The interactions have the effect that the electrons are treated as quasi-particles, particles which are closely related to their environment.

For wires where the energy relaxation time exceeds the diffusion time, the transport is fully elastic and the energy distribution is calculated by a quantum diffusion equation without an inelastic collision term. The distribution function behaves differently in the two field limits, $\omega\tau_{D}<<1$ and $\omega\tau_{D}>>1$. For the slow field, strong signal limit ($\omega\tau_{D}<<1$, $\hbar\omega/eV<<1$) the dc situation is approached and the energy distribution is given by a time-varying two step function. In the fast field, strong signal limit ($\omega\tau_{D}>>1$, $\hbar\omega/eV<<1$) a quasi-particle can oscillate multiple times with the field in the wire before leaving the wire, therefore gaining more energy quanta. This results in a time-independent electron energy distribution which does not go to zero at higher energies.

Numerical simulations for a finite ratio $\hbar\omega/eV$ shows the photon steps in the energy distribution. In the slow field regime the distribution is highly time-dependent. The energy distribution directly follows the field and the photon steps oscillate from zero photon absorption to maximum photon absorption. In the fast field limit the energy distribution becomes completely time-independent. The maximum photon absorption is reached when the diffusion time is exceeded. In the crossover the two effects are both observed. There is a slight oscillation around the maximum photon absorption value of the fast field regime.

When the energy relaxation time becomes comparable to the diffusion time, the transport is no longer coherent and scattering is inelastic. Numerical simulations show how the energy exchange processes of mutual quasi-particles and between quasi-particles and phonons influence the energy distribution in the fast field regime, which is of interest for experimental situations. The interaction between quasi-particles and phonons annihilates first the photon steps in the distribution. In the strong electron-phonon interaction limit the Fermi function at bath temperature is found on every position in the wire.

The interaction between quasi-particles is quite different. It causes a smearing in the photon steps. In the strong electron-electron interaction limit a local equilibrium is reached on every position in the wire. The photon absorption, diffusive transport and the interactions cause an effective temperature profile across the wire. The effective temperature profile is determined by fitting the distribution from the simulation on every position with a Fermi function. The obtained temperature profile deviates from the calculated profile. It is not yet fully understood whether this is caused by a numerical error or that it is caused by some physical effect, like for instance the position in the wire that the photon absorption takes place.

So the complicated interplay between the effect of photon absorption, diffusive transport and inelastic scattering on the quasi-particle energy distribution seems to be accurately described by our model.

7.2 Discussion

The model developed in this work is not yet verified by experiments. The fabrication of the required samples is not yet proven to be completely succesfull due to the failure of fabrication apparatus, however the outlook is promising and the experimental work will be continued by members of this research group.

To calculate the collision term, the interactions are assumed to be instantaneous and local. An experiment can verify this assumption, so it is desirable to proceed with the experimental part of this project. The comparison between experiment and theoretical model has to provide a full insight in the ac quantum transport in diffusive quasi-one dimensional wires and how the non-equilibrium is shown in the quantum statistics of the quasi-particles in the wire.

While our MATLAB code seems to calculate the quasi-particle energy distribution influenced by photon absorption, diffusive transport and inelastic collisions in a correct manner, it also appears that the MATLAB code is not very efficient. A program written in C should in principle work more efficient. This provides the opportunity to optimize the discretization of the variables, so that a more accurate result is obtained.

Acknowledgements

The phrase Tempus fugit is the first thing that comes to mind when I look back at the time spent in this group working on my master thesis project. It was certainly an interesting and challenging time. The project did not go as planned. One of the main goals was to experimentally obtain data on the subject. However, reality seems to play tricks on you when you want to control things. Apparatus failure caused such delay in the fabrication of the samples that I did not get to experience the beauty of experimental success. This pushed me further to theoretical research and I have to say, I didn’t mind. It was really intriguing how the physics on small scale revealed itself to me by doing the math and combining concepts. In addition I learned different experimental techniques in preparation of the big experiment which didn’t come.

I would like to thank prof. Teun Klapwijk for giving me the opportunity of doing this challenging work. It really gave me the chance to evolve in different disciplines. I also would like to express my gratitude towards my daily supervisors, Nathan and Rik. There were times that I was not certain whether everything would work out, but they gave me the confidence that leaded to this satisfying result and educated me in the necessary basics of experimental and theoretical research.

Furthermore, it was just a fun time. I really experienced to be a part of the group. I especially want to thank Nathan, Rik, Eduard and David for the interesting conversations on everything and nothing. Also I want to thank Reinier for his company during the whole year in the huge students office. I want to thank the whole group for this fun time and for their part in my education. I am confident that the gained experience in this group will be of great help in the continuation of my career.

Appendix A Shot noise

Shot noise, the phenomenon that arises because of charge quantization, provides information about the statistics of charge carriers involved in quantum transport. In chapter 2 the mean square fluctuations in the occupation of incident, reflected and transmitted state are derived to be:

[TABLE]

It appears that at zero temperature, when the distribution is given by a step function at chemical potential, the shot noise disappears for full reflectance or full transmittance. At finite temperatures the mean square fluctuations fluctuate like the incident state with occupation $f$.

Now when we proceed with this simplified model of a single incident charge carrier to investigate the fluctuations in the current, we can consider a perfect conductor with three separated channels. One for the incident state, one for the reflected state and one for the transmitted state. The carriers move in one directions with a velocity $v(E)$ dependent on the energy of the charge carrier. For an energy interval $dE$, the incident current is given by $dI_{in}=ev(E)d\rho(E)$. $\rho(E)$ is the energy dependent density of carriers per unit length. It is given by $\rho(E)=n_{in}(E)\nu(E)dE$, where $\nu(E)$ is the density of states per unit length. In a perfect conductor $\nu(E)=(2\pi\hbar v(E))^{-1}$. This leads to $dI_{in}=e(2\pi\hbar)^{-1}n_{in}(E)dE$. For the transmitted and reflected channel the same procedure can be followed leading to $dI_{T}=e(2\pi\hbar)^{-1}n_{T}(E)dE$ and $dI_{R}=e(2\pi\hbar)^{-1}n_{R}(E)dE$. Integrating gives the expressions for the current. When the occupation numbers vary slowly in time, the derivation can be expanded by just taking the occupation not only energy-dependent, but also time-dependent.

[TABLE]

For low frequency fluctuations in the current, these expressions can be Fourier transformed giving the frequency dependent current.

[TABLE]

The fluctuations in current and occupation number are directly related. The current noise power is in the zero frequency limit $S_{II}=e(2\pi\hbar)^{-1}\int S_{nn}(E)dE$. From the current it is seen that the charge carriers arrive at a rate of $dE/(2\pi\hbar)$ in each energy interval. This contributes to the noise with the mean square fluctuations of the relevant state. Therefore $S_{nn}(E)=1/(\pi\hbar)\left\langle\delta n\delta n\right\rangle$. Substitution in the current noise power relations leads to

[TABLE]

So here we see explicitly the earlier found conclusion for the fluctuations that shot noise disappears for full transparency or full reflectance at zero temperature. At finite temperature the distribution function is thermally broadened and the shot noise will not disappear due to thermal fluctuations. When now the system under consideration is extended to a situation where a multi-channel scatterer is placed between two terminals with respectively distribution $f_{L}$ and $f_{R}$ the noise power is given by

[TABLE]

The scale of the energy dependence of the transmission coefficients is usual much bigger than the thermal and bias energy. Therefore these quantities can be taken in equation A.14 at Fermi energy. Then the noise power becomes

[TABLE]

Appendix B MATLAB code of the simulation program

B.1 Script for the simulation of coherent transport

%========================================================================== %CLEAN UP, FUNDAMENTALS %========================================================================== clear all; tic;

%fundamental constants e = 1.602e-19; hbar = 6.63e-34/2/3.141592; %6.6e-16; % kb = 1.38e-23; % 8.6e-5;

%========================================================================== %PARAMETERS %==========================================================================

%general saveall=0; saverepeat=1250;

%time discretization Nt0=7; dt=2pi/Nt0; int_method=’euler’; Nt=5000Nt0;%1000000;

%energy discretization T=2; omega=0.3e1223.141592; V=3e-3; limit_E=2; %limit in multiples of eV hw=hbaromega/e/V; Nw=16; %must be even! dE=hw/Nw; NE=ceil(limit_E/dE); E=-NEdE:dE:(NE-1)*dE+dE/2; NE=length(E);

kT=kb*T/e/V;

%space discretization D_method=’lagrange_2’; Nx=100; x = linspace(0,1,Nx); dx = 1/(Nx-1); %tauD= 1e-9; %z=omega*tauD; z=30000;

fprintf(’homega/eV= %f kT/eV= %f dE=%f Emax= %f omega tau= %f \n’,hw,kT,dE,NE*dE,z);

%========================================================================== %INITIALIZATION %==========================================================================

%INITIAL AND BOUNDARY CONDITIONS Fl=1./(exp((E)/kb/TeV)+1); Fr=1./(exp((E)/kb/TeV)+1); Fold=[ones(Nx-1,1)*Fl; Fr]; Feq=[ones(Nx-1,1)*Fl; Fr]; Fave=[ones(Nx-1,1)*Fl; Fr];

%MATRICES D_x=dx1(Nx,1,dx,D_method); D_xx=dx2(Nx,1,dx,D_method);

D_E=dE1(1,NE,Nw,dE)’; D_EE=dE2(1,NE,Nw,dE)’;

first_step=[ ]; second_step=[]; third_step=[];

%========================================================================== %INTEGRATION %==========================================================================

switch lower(int_method) case ’euler’ %euler for m=1:Nt w=sin(mdt); dF=D_xxFold+2wD_xFoldD_E+w^2FoldD_EE; Fnew=Fold+dFdt/z; Fnew([1 Nx],:)=[Fl;Fr]; Fold=Fnew; Fave=((m-1)Fave+Fold)/m; first_step=[first_step Fnew(round([Nx/2])’,NE/2+Nw/2)]; second_step=[second_step Fnew(round([Nx/2])’,NE/2+3Nw/2)]; third_step=[third_step Fnew(round([Nx/2])’,NE/2+5Nw/2)]; if mod(m,5000)==0 for o=1:Nx for p=1:NE if Fnew(o,p)>1 Fnew(o,p)=1; end if Fnew(o,p)<0 Fnew(o,p)=0; end end end end if mod(m,1000)==0 %saverepeat clc; fprintf(’iteratie %i, time %f’,m,toc);%round(m/Nt0*100) %tic; % hold on; plot(E,Fnew(round([1 Nx/4 Nx/2])’,:)); pause(.2); if saveall save([’Fnew’ num2str(m)],’Fnew’); end end end end

B.2 Functions for the used operators

First spatial derivative

function S=dx1(Nx,NE,h,method)

switch lower(method) case {’lagrange_1’} %lagrange 1st order r=2; w1=[-3 4 -1]; w3=[-1 0 1]; n=1; case ’lagrange_2’ %lagrange r=12; w1=[-25 48 -36 16 -3]; w2=[-3 -10 18 -6 1]; w3=[1 -8 0 8 -1]; n=2; case ’least’ %least squares r=70; w1=[-54 13 40 27 -26]; w2=[-34 3 20 17 -6]; w3=[-2 -1 0 1 2]*7; n=2; end

%the boundaries of the matrix base0=1:NE; v0=ones(1,NE);

%first boundary array1=[]; column1=[]; values1=[]; for k=1:length(w1) array1=[array1 base0]; column1=[column1 base0+(k-1)NE]; values1=[values1 w1(k)v0]; end column1=[column1 NENx-column1+1]; array1=[array1 NENx-array1+1]; values1=[values1 values1];

%second boundary (if needed) array2=[]; column2=[]; values2=[]; if n==2 for k=1:length(w2) array2=[array2 base0+NE]; column2=[column2 base0+(k-1)NE]; values2=[values2 w2(k)v0]; end column2=[column2 NENx-column2+1]; array2=[array2 NENx-array2+1]; values2=[values2 values2]; end

%central part if n==1 base1=1:NE*(Nx-2); v1=ones(1,NE*(Nx-2)); else base1=1:NE*(Nx-4); v1=ones(1,NE*(Nx-4)); end array3=[]; column3=[]; values3=[]; for k=1:length(w3) array3=[array3 base1+n*NE]; column3=[column3 base1+(k-1)*NE]; values3=[values3 w3(k)*v1]; end

%the total matrix S=sparse([array1 array2 array3],[column1 column2 column3],[values1 values2 values3])/r/h;

Second spatial derivative

function S=dx2(Nx,NE,h,method)

switch lower(method) case {’lagrange_1’} %lagrange 1st order r=1; w1=[2 -5 4 -1]; w3=[1 -2 1]; n=1; case ’lagrange_2’ %lagrange r=12; w1=[35 -104 114 -56 11]; w2=[10 -15 -4 14 -6 1]; w3=[-1 16 -30 16 -1]; n=2; case ’least’ %least squares r=14; w1=[9 -15 -2 13 -5]*2; w2=[11 -16 -4 12 -3]; w3=[2 -1 -2 -1 2]*2; n=2; end

%the boundaries of the matrix base0=1:NE; v0=ones(1,NE);

%first boundary array1=[]; column1=[]; values1=[]; for k=1:length(w1) array1=[array1 base0]; column1=[column1 base0+(k-1)NE]; values1=[values1 w1(k)v0]; end column1=[column1 NENx-column1+1]; array1=[array1 NENx-array1+1]; values1=[values1 values1];

%second boundary (if needed) array2=[]; column2=[]; values2=[]; if n==2 for k=1:length(w2) array2=[array2 base0+NE]; column2=[column2 base0+(k-1)NE]; values2=[values2 w2(k)v0]; end column2=[column2 NENx-column2+1]; array2=[array2 NENx-array2+1]; values2=[values2 values2]; end

%central part if n==1 base1=1:NE*(Nx-2); v1=ones(1,NE*(Nx-2)); else base1=1:NE*(Nx-4); v1=ones(1,NE*(Nx-4)); end array3=[]; column3=[]; values3=[]; for k=1:length(w3) array3=[array3 base1+n*NE]; column3=[column3 base1+(k-1)*NE]; values3=[values3 w3(k)*v1]; end

%the total matrix S=sparse([array1 array2 array3],[column1 column2 column3],[values1 values2 values3])/r/h/h;

First energy derivative

function S=dE1(Nx,NE,Nw,h)

%boundaries array_0=[1:Nw NE-[1:Nw]+1]; column_0=[ones(1,Nw) ones(1,Nw)*NE]; values_0=[-ones(1,Nw) ones(1,Nw)];

%central part base=1:NE-Nw; v0=ones(1,NE-Nw);

array_1=[base+Nw base]; column_1=[base base+Nw]; values_1=[-v0 v0];

array=[array_0 array_1]; column=[column_0 column_1]; values=[values_0 values_1];

for k=1:Nx-1 array=[array array_0+NEk array_1+NEk]; column=[column column_0+NEk column_1+NEk]; values=[values values_0 values_1]; end

%the total matrix S=-sparse(array,column,values)/2/h/Nw;

Second energy derivative

function S=dE2(Nx,NE,Nw,h)

%boundaries array_0=[1:2Nw 2:2Nw 1:2Nw]; array_0=[array_0 NE+1-array_0]; column_0=[ones(1,2Nw) 2:2Nw 2Nw+1:4Nw]; column_0=[column_0 NE+1-column_0]; values_0=[-1 ones(1,2Nw-1) -2ones(1,2Nw-1) ones(1,2*Nw)]; values_0=[values_0 values_0];

%central part base=1:NE-2Nw; %4Nw v0=ones(1,NE-2Nw); %4Nw

array_1=[base+1Nw base+1Nw base+1Nw]; %2Nw ; 2Nw ; 2Nw column_1=[base base+1Nw base+2Nw]; %2Nw ; 4Nw values_1=[v0 -2*v0 v0];

array=[array_0 array_1]; column=[column_0 column_1]; values=[values_0 values_1];

for k=1:Nx-1 array=[array array_0+NEk array_1+NEk]; column=[column column_0+NEk column_1+NEk]; values=[values values_0 values_1]; end

%the total matrix S=sparse(array,column,values)/(2Nwh)^2;

B.3 Script for the simulation of incoherent transport

%========================================================================== %CLEAN UP, FUNDAMENTALS %========================================================================== clear all; tic;

%fundamental constants e = 1.602e-19; hbar = 6.63e-34/2/3.141592; %6.6e-16; kb = 1.38e-23; % 8.6e-5;

%========================================================================== %PARAMETERS %==========================================================================

%general saveall=0; saverepeat=1250;

%time discretization Nt0=1.1; dt=2pi/Nt0; Nt=79000Nt0;%20000;

%energy discretization T=0.5; omega=1.6e1223.141592; V=16e-3; limit_E=2; %limit in multiples of eV hw=hbaromega/e/V; Nw=50; %must be even! dE=hw/Nw; NE=ceil(limit_E/dE); E=-NEdE:dE:(NE-1)*dE+dE/2; NE=length(E);

kT=kb*T/e/V;

%interaction parameters int_mech=’eph’; %interaction mechansim, ee for electron-electron, eph for electron phonon

rho=2.7e3; EF=12/V; %J dos=2e47eV; s=6.42e3; kf=1.75e10; D=100e14; %nm^2s^-1 dos=2e47; S=20e-9*20e-9; sigma=1e9;

Ke=(sqrt(2*D)pihbar^(3/2)dosS)^(-1);

NE_ee=NE; E_ee=linspace(10,14,NE_ee)/V; dE_ee=4*e/NE_ee/V;

NE_ph=NE; E_ph=linspace(-2,2,NE_ph); dE_ph=4*e/NE_ph;

kph=sigma/24/zeta(5)/dos/kb^5*(eV)^2dE_ph; %kph=4e12*1e-3/V; %has to be expressed in V

n_ph=(1./(exp(abs(E_ph)/kb/2eV)-1))’;

%space discretization D_method=’lagrange_2’; Nx=20; x = linspace(0,1,Nx); dx = 1/(Nx-1); %tauD= 1e-9; %z=omega*tauD; z=2000000;

fprintf(’homega/eV= %f kT/eV= %f dE=%f Emax= %f omega tau= %f \n’,hw,kT,dE,NE*dE,z);

%========================================================================== %INITIALIZATION %==========================================================================

%INITIAL AND BOUNDARY CONDITIONS Fl=1./(exp((E)/kb/TeV)+1); Fr=1./(exp((E)/kb/TeV)+1); Fold=[ones(Nx-1,1)*Fl; Fr]; Feq=[ones(Nx-1,1)*Fl; Fr]; Fave=[ones(Nx-1,1)*Fl; Fr];

%MATRICES D_x=dx1(Nx,1,dx,D_method); D_xx=dx2(Nx,1,dx,D_method);

D_E=dE1(1,NE,Nw,dE)’; D_EE=dE2(1,NE,Nw,dE)’;

Iin=0; Iout=0; h=400;

first_step=[ ]; second_step=[]; third_step=[];

%========================================================================== %INTEGRATION %==========================================================================

switch lower(int_mech) case ’ee’ %euler for m=1:Nt w=sin(mdt); dF=D_xxFold+2wD_xFoldD_E+w^2FoldD_EE+Iin*(z/omega)-Iout*(z/omega); Fnew=Fold+dFdt/z; Fnew([1 Nx],:)=[Fl;Fr]; Fold=Fnew; Fave=((m-1)Fave+Fold)/m; first_step=[first_step Fnew(round([Nx/2])’,NE/2+Nw/2)]; second_step=[second_step Fnew(round([Nx/2])’,NE/2+3Nw/2)]; third_step=[third_step Fnew(round([Nx/2])’,NE/2+5Nw/2)]; if mod(m,50)==0 for o=1:Nx for p=1:NE if Fnew(o,p)>1 Fnew(o,p)=1; end if Fnew(o,p)<0 Fnew(o,p)=0; end end end end if mod(m,h+1)==0 for o=1:Nx for p=1:NE if Fnew(o,p)>1 Fnew(o,p)=1; end if Fnew(o,p)<0 Fnew(o,p)=0; end end end end if mod(m,h-1)==0 for o=1:Nx for p=1:NE if Fnew(o,p)>1 Fnew(o,p)=1; end if Fnew(o,p)<0 Fnew(o,p)=0; end end end end if mod(m,800)==0 Fnew2=Fnew; Fnew1=Fnew; [q,r]=size(Fnew); g1=diag(ones(r-1,1),1); g2=diag(ones(r-1,1),-1); Fnew1=circshift(Fnew1,[0 1])+[Fnew(:,1), (zeros(NE-1,Nx))’]; b=[]; a=[]; for j=1:Nx; k1=Fnew1(j,:); k2=Fnew2(j,:);

                for i=1:NE/2-1
                    h1=Fnew1(j,:)*g1;
                    h1=h1+[ones(1,1)’, zeros(NE-1,1)’];
                    k1=[k1’ h1’]’;
                    Fnew1(j,:)=h1;
                    h2=Fnew2(j,:)*g2;
                    k2=[h2’ k2’]’;
                    Fnew2(j,:)=h2;
                end
                k=[k2’ k1’]’;
                E2=((E_ee).^(-3/2))’;
                W=(1-fliplr(k’))*Fnew(j,:)’;
                W2=Ke*E2.*W;
                %In=W2’*fliplr(k’)*dE;
                Out=W2’*(1-k)*dE;
                In=fliplr(Out);
                b=[b Out’];
                a=[a In’];
            end
           Iout=(Fnew.*b’);
           Iin=((1-Fnew).*a’);
        end
        if mod(m,1000)==0 %saverepeat
            clc;
            fprintf(’iteratie %i, time %f’,m,toc);%round(m/Nt0*100)
            %tic;
            plot(E,Fnew(round([1 Nx/4 Nx/2])’,:));
            pause(.2);
            if saveall
            save([’Fnew’ num2str(m)],’Fnew’);
            end
        end
    end
case ’eph’
        %euler
    for m=1:Nt
        w=sin(m*dt);
        dF=D_xx*Fold+2*w*D_x*Fold*D_E+w^2*Fold*D_EE+Iin*(z/omega)-Iout*(z/omega);
        Fnew=Fold+dF*dt/z;
        Fnew([1 Nx],:)=[Fl;Fr];
        Fold=Fnew;
        first_step=[first_step Fnew(round([Nx/2])’,NE/2+Nw/2)];
        second_step=[second_step Fnew(round([Nx/2])’,NE/2+3*Nw/2)];
        third_step=[third_step Fnew(round([Nx/2])’,NE/2+5*Nw/2)];
        if mod(m,50)==0
            for o=1:Nx
                for p=1:NE
                    if Fnew(o,p)>1
                        Fnew(o,p)=1;
                    end
                    if Fnew(o,p)<0
                        Fnew(o,p)=0;
                    end
                end
            end
        end
        if mod(m,h+1)==0
            for o=1:Nx
                for p=1:NE
                    if Fnew(o,p)>1
                        Fnew(o,p)=1;
                    end
                    if Fnew(o,p)<0
                        Fnew(o,p)=0;
                    end
                end
            end
        end
        if mod(m,h-1)==0
            for o=1:Nx
                for p=1:NE
                    if Fnew(o,p)>1
                        Fnew(o,p)=1;
                    end
                    if Fnew(o,p)<0
                        Fnew(o,p)=0;
                    end
                end
            end
        end
        Fave=((m-1)*Fave+Fold)/m;
        if mod(m,400)==0
            h=m;
            Fnew2=Fnew;
            Fnew1=Fnew;
            [q,r]=size(Fnew);
            g1=diag(ones(r-1,1),1);
            g2=diag(ones(r-1,1),-1);
            Fnew1=circshift(Fnew1,[0 1])+[Fnew(:,1), (zeros(NE-1,Nx))’];
            b=[];
            a=[];
            for j=1:Nx;
                k1=Fnew1(j,:);
                k2=Fnew2(j,:);
                for i=1:NE/2-1
                    h1=Fnew1(j,:)*g1;
                    h1=h1+[ones(1,1)’, zeros(NE-1,1)’];
                    k1=[k1’ h1’]’;
                    Fnew1(j,:)=h1;
                    h2=Fnew2(j,:)*g2;
                    k2=[h2’ k2’]’;
                    Fnew2(j,:)=h2;
                end
                k=[k2’ k1’]’;
                E2=(E.^2)’;
                W=E2.*(n_ph+heaviside(E_ph)’);
                Out=kph*(1-(k’))*W*dE;

% In=kph*fliplr(k’)WdE; In=fliplr(Out’)’; b=[b Out]; a=[a In]; end Iout=Fnew.*b’; Iin=(1-Fnew).a’; end if mod(m,1000)==0 %saverepeat clc; fprintf(’iteratie %i, time %f’,m,toc);%round(m/Nt0100) %tic; plot(E,Fnew(round([1 Nx/4 Nx/2])’,:)); pause(.2); if saveall save([’Fnew’ num2str(m)],’Fnew’); end end end end

Appendix C Space dependency in the distribution function

Appendix D Differential conductance of a NIS junction

To obtain the electron energy distribution function on a certain position in the wire, superconducting probes are used that are positioned orthogonally on top of the wire with an insulating layer in between. Now by applying a current to this NIS junction we can measure the reciprocal value of the differential conductance. This differential conductance is a convolution of the density of states of the used superconductor and the probed energy distribution function. This can be calculated from Fermi’s golden rule when the tunnel matrix elements are considered nearly constant for all energy states consider in the measurement. Because the energy distribution is probed with a superconductor the tunneling current is elastic, as there is no energy dissipation in a superconductor. So from the golden rule the tunneling rate from the wire to the superconductor and from the superconductor to the wire become respectively

[TABLE]

[TABLE]

Here $e$ is the elementary charge, $R_{t}$ is the tunnel resistance, $n_{x}$ and $n_{p}$ are respectively the density of states in the wire and in the superconductor and $f_{x}$ and $f_{p}$ are the distribution functions in respectively the wire and the superconductor. The current across the junction is calculated from these tunnel rates.

[TABLE]

When the tunnel rates are implemented in the formula above and the density of states of the superconductor is taken to be the BCS density of states $n_{BCS}(E)=\Re(|E|/\sqrt{E^{2}-\Delta^{2}})$ and as the wire is a metal the density of states in the wire is taken to be flat, we arrive at the following expression.

[TABLE]

A variable change of $E\rightarrow E-eV$ and taking the derivative with respect to $V$ leads to an expression for the differential conductance of the NIS junction.

[TABLE]

The distribution function shows up explicitly by integration by part and using the fact that the BCS density of states is even.

[TABLE]

In practice the singular behaviour of the density of states of the superconductor is less sharp than the BCS theory predicts. To avoid problems with this aspect, the effective density of states of the superconductor in the NIS junction can be probed first for an unbiased wire. Since in this situation the energy distribution of the electrons is just a quasi-equilibrium Fermi function, the effective density of states is obtained by a deconvolution of the differential conductance of the NIS junction and the distribution function. This is basically finding a fit for the gap energy $\Delta$, the tunneling resistance $R_{t}$ and the electron temperature $T$.

Appendix E MATLAB script for deconvolution

clear all; dVdI = importdata(’testdatadIdV.txt’); dIdV1 = (1./dVdI); m=length(dIdV1); dIdV = [dIdV1’ dIdV1(m)]’; kb = 8.6e-6;%1.3810^-23; T= 4.2; e = 1.602e-19; NE = length(dIdV); NV = length(dIdV); d = 1.310^(-3);%e; E = linspace(-1010^(-3),10*10^(-3),NE);%e; dE = 20e-3/NE;%e; dV = 20e-3/NV; V = linspace(-1010^(-3),1010^(-3),NV);%e; Rt = 1;%;1.210^4; V2=1e-3;

for i=1:length(V) for j = 1:length(E) f(i,j) = 1/Rt*(exp((E(j)-V(i))/kb/T)./((exp((E(j)-V(i))/kb/T)+1).^2kbT)); end end

Ns=f\dIdV/dE; Ns(1:5)=Ns(6); Ns(246:m+1)=Ns(245);

for i = 1:length(V) dIdV99(i) = 1/Rt.(exp((E-V(i))/kb/T)./((exp((E-V(i))/kb/T)+1).^2kb*T)*Ns)*dE; end

dIdV99(1:2)=dIdV99(3); dIdV99(250:m+1)=dIdV99(249);

Diff_n=diag(ones(NE-1,1),1)-diag(ones(NE-1,1),-1); Diff_n(1:2,1)=[1 ;1]; Diff_n(NE-1:NE,NE)=[1 ;1];

DnBCS=Diff_n*Ns/dE; DnBCS(1:2)=DnBCS(3); DnBCS(251)=DnBCS(250);

E=E’;

nBCS=real(abs(E)./(E.^2-d^2).^(1/2));

fold=1./(exp(E/kb/T)+1); Diff_f=-diag(ones(NE-1,1),1)+diag(ones(NE-1,1),-1); Diff_f(1:2,1)=[1 ;1]; Diff_f(NE-1:NE,NE)=[1 ;-1]; dfold=Diff_f*fold/dV;

[m,n]=size(dIdV);

epsilon=3.9e-14*ones(m,1);

chi2=((dIdV99-dIdV’).^2)’;

while sum(chi2 > epsilon) > 0 fprintf(’q %i \n’,sum(chi2 > epsilon)); fnew=fold+5e0*(DnBCS.chi2.^(1/2)); dfnew=Diff_ffnew/dV; g1=diag(ones(m-1,1),1); g2=diag(ones(m-1,1),-1); dfnew1=dfnew; dfnew2=dfnew; dfnew1=circshift(dfnew1,[0 1])+[dfnew(:,1), (zeros(n-1,m))’]; k1=dfnew1; k2=dfnew2; for i=1:NE/2-1 h1=g1dfnew1; h1=h1+[ones(1,1)’, zeros(NE-1,1)’]’; k1=[k1 h1]; dfnew1=h1; h2=g2dfnew2; k2=[h2 k2]; dfnew2=h2; end dfnew_m=[k2 k1]; q=1/Rt*Ns’dfnew_m; g=dIdV(1)/q(1); q=1/RtNs’dfnew_mg; chi2=((q-dIdV’).^2)’; fold=fnew; end

plot(E,fnew)