{'fraction_non_alphanumeric': 0.04056473829201102, 'fraction_numerical': 0.012493112947658402, 'mean_word_length': 4.0498017667107185, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Green's functions for solving differential equations, in non-boundary value problems in near-field optics and in quantum transport through point contacts

This introduction to Green's functions is based on their role as kernels of differential equations. The procedures to construct solutions to a differential equation with an external source or with an inhomogeneity term are put together to derive the Dyson equation for the Green's function of the inhomogeneous system. Very different areas of physics such as, for example, electrodynamics and quantum transport, can profit from this Green's function formalism. The fundamental homogeneous-medium Green's tensor of electrodynamics is deduced from the field of a dipole. Based upon that a numerical procedure is presented to solve the wave-equation for the near-field in a scattering setup for arbitrary material distributions. The full inhomogeneous system's Green's function is not explicitly needed to get the fields, although it can be obtained by a very similar calculation and in optics can be interpreted as a density of states. It is demonstrated how the transport problem for two open free-electron gas reservoirs with arbitrary coupling can be solved by finding the system's Green's function. In this sense the article is an introduction on Green's functions for treating interaction. A very detailed discussion of the current formula is given on an elementary basis.


I.1. Introduction
Green's function formalisms do not present a cure-all for solving differential equations, because essentially the problem of finding the solution of the differential equation is shifted to that of finding the corresponding Green's function. This can, however, be a simplification and even give access to the solution of a more general class of problems. Starting from ordinary differential equations, this short review is meant to introduce how the corresponding Green's functions are defined and how they are involved in constructing the solutions for different types of differential equations. Although bearing similarities, the Green's function formalism can go beyond perturbation theory. Special emphasize is laid here upon the parallels between homogeneous and inhomogeneous ordinary differential and Green's functions equations.

I.2. Homogeneous equation
The starting point is a homogeneous differential equation which we suppose is exactly solvable, although φ 0 (ξ) will not explicitly be needed. D ξ is some differential operator which may include multiplication with a constant or even another function of ξ. ξ is either a space or time variable. What we need to know is the solution g of the corresponding Green's functions equation g is a tensor-like function of two arguments, D ξ only acting on the first of them. Like g replaces φ on the left, the zero on the right side of (1) is replaced by a δ-distribution in (2).

I.4. Inhomogeneity
Instead of a source term, the differential equation can contain a potential term which we shall call an inhomogeneity.
that is a ξ-dependent function V multiplied with φ. V (ξ) is not included in D ξ , because we assume that it so much complicates the equation that a standard solution is no longer known. The minus sign is a useful convention. We put the inhomogeneity term on the right to make (5) formally look like (4). With V φ playing the role of Q a formal solution is constructed analogously to (3): φ ih (ξ) = φ 0 (ξ) + dξ g(ξ, ξ ) V (ξ ) φ ih (ξ ) (6) This presents an implicit equation for φ ih , the so-called Lippmann-Schwinger equation. We have replaced the differential equation by an integral equation. Inserting the solution of the homogeneous equation φ 0 also for φ ih in the integral on the right side of (6) would give the Born approximation φ ih (ξ) ≈ φ 0 (ξ) + dξ g(ξ, ξ ) V (ξ ) φ 0 (ξ ), which is appropriate if V φ is a small perturbation compared to D ξ φ. However, it is the virtue of the Green's functions method that in contrast to perturbation theory the inhomogeneity need not be a small deviation.
Although not yet providing an explicit solution for φ ih , one can make use of (6) in numerical calculations (see section II and [8]). As (2) is the corresponding Green's function equation with a point source to the homogeneous equation (1), the equation defining the Green's function G for the inhomogeneous case derived from (5) reads: However, while supposing that we know or can easily guess g(ξ, ξ ), there is no hint yet how to calculate G(ξ, ξ ).

I.5. General case
Finally we have to treat the most general case with source and inhomogeneity term: In analogy to (3) and (4) a yet formal solution can be written down using (7): Here we have the freedom to add any solution of (6) to the integral.
But let us start again solving the problem directly from (8) which we rewrite as We compare this with (5a) and construct a solution of the same form as (6). φ q takes the role of φ 0 , because it would satisfy (8a) if the right side were zero. φ ih has to be replaced by φ ihq . Using (3) for φ q we obtain φ ihq (ξ) = φ 0 (ξ) + dξ g(ξ, ξ ) Q(ξ ) (10) could also have been established by a different line of thought. In (8) putting everything that differs from the homogeneous equation (1) on the right side gives If it were not for the inhomegeneity term V φ, the solution would be (3). And if it were not for the source term Q, we could use (6). Although this is, of course, not a correct way to solve non-homogeneous differential equations, we can understand (10) as an ansatz adding these two contributions. The φ 0 -part does not have to be written twice. And φ ihq has to appear instead of φ ih in the integral taken from (6). This still unknown φ ihq leaves the neccessary freedom to somehow counterbalance the Q-term not present in (5) and (6), which justifies (10) as an ansatz. Whichever way it was obtained, (10) is an implicit integral equation for φ ihq as (6) is for φ ih .

I.6. Dyson equation
Now we shall profit from the fact that with (9) we have a second representation of φ ihq .
Insert the expression from (9) for φ ihq both on the left and on the right side in (10): φ ih on the left cancels with φ 0 and one integral term on the right according to (6) and we are left with dξ G(ξ, ξ ) Q(ξ ) = dξ g(ξ, ξ ) Q(ξ ) + dξ dξ g(ξ, ξ ) V (ξ ) G(ξ , ξ ) Q(ξ ) This is an implicit equation for the unknown Green's function G. But as G by definition (7) does not depend on any source term, (11) must be valid for arbitrary Q. The special choice which after carrying out the ξ -integrals and then renaming ξ 0 to ξ becomes G(ξ, ξ ) = g(ξ, ξ ) + dξ g(ξ, ξ ) V (ξ ) G(ξ , ξ ) (11) (11) is the general implicit integral equation for the inhomogeneous system's Green's function G and called Dyson equation.
We could have derived (11) by doing all calculations from (8) on with only a point source δ(ξ − ξ 0 ) instead of Q(ξ). Putting in the solutions (3) and (9) for this case would just have looked a bit awkward. (11) is not an explicit solution for G, but still an implicit equation, and even though it is often written as G = g + gV G one must not forget that there is a convolution-like integration over the inner ξ-index behind the sequence of factors gV G. The integration can either be transformed into a discrete finite sum, and thus (11) into a linear system of equations solvable by a matrix inversion (see section II and [9]), or the convolution can be replaced by a multiplication going to frequency space by a Fourier transformation (see section III and [7]).

I.7. Expansions
There is some subtlety about the φ 0 contribution in (6). (6*) has the trivial solution φ * ih (ξ) = 0 and not necessarily another one. Surely, φ ih (ξ) = 0 satisfies (5), but it is not what we want. φ 0 is the background excitation replacing boundary conditions in this kind of problem, which shall become clear when discussing electrodynamics (see following article). φ 0 acts as a source term which -instead of being put in Q -can be more appropriately and conveniently set as a fixed part of the φ we are looking for. Inserting (6) into itself ever and ever again, φ ih is developed into a series in powers of gV : φ ih = φ 0 + gV φ ih = φ 0 + gV (φ 0 + gV (φ 0 + gV (. . . . . . . . .))) Just for shorthand notation we dropped the arguments and integration variables and in the third line also the integral signs. The resulting series is recognized from perturbation theory summing interactions to zeroth, first, second, etc. order. If there were no contribution φ 0 in φ ih , only the contribution of gV to power infinity would exist with nothing to multiply to at the end. In other words we would have no basis on which to develop φ ih . The Green's functions formalism as we use it aims at solving (6), (10) or (11) in a closed form, not cutting the series in (12) at some finite order. The development (12) is only shown here precisely to demonstrate that interactions are included to all orders as well as to explain the importance of the homogeneous background contribution. The Dyson equation (11) can be expanded in an analogous manner to the Lippmann-Schwinger equation, by the way prooving the equivalence to its complementary form G = g + GV g: (Integration over inner arguments is understood in all contributions to the sums.) Either formally or as a matrix calculation in finite discrete ξ-space, (11) is often solved as Developing the (. . .) −1 -factor in (14) into a geometric series just results in the infinite sums written in (13).
Putting together (12) and (13) we can get an alternative representation to (6) for the solution of the inhomogeneous differential equation: Therefore in the case that should it be easier to get the Green's function G than to solve the implicit equation (6) for φ ih , we see that G can be useful also for treating the equation without an external source Q. Nevertheless, (15) again illustrates that V φ 0 plays the role of the source.

I.8. Conclusions and Outlook for section I
Using Green's functions it has been shown how differential equations can be treated that differ from easily solvable ones by an additional potential term or an arbitrary source term. The formalism as presented here is for open-boundary in contrast to boundary-value problems [3,10]. From constructing the solutions of the differential equations, we also obtained the constituting relation for the inhomogeneous system's Green's function, which characterizes the response to a point source and includes all-order interactions. Examples of applications will be given in the two following sections.

II. NON-BOUNDARY VALUE PROBLEMS IN NEAR-FIELD OPTICS
The fundamental homogeneous-medium Green's tensor of electrodynamics is deduced from the field of a dipole. Based upon that a numerical procedure is presented to solve the wave-equation for the near-field in a scattering setup for arbitrary material distributions.
The full inhomogeneous system's Green's function is not explicitly needed to get the fields, although it can be obtained by a very similar calculation and in optics can be interpreted as a density of states.

II.1. Introduction
The typical problem in nano-optics [11,12] is the situation that some tiny structures are illuminated by an extended source, a plane wave for example, and then one is interested in the field distribution that arises from multiple scattering [13], especially to identify places where the field intensity gets considerably enhanced [14,15]. A theory can be based on Green's functions, however, their implication differs slightly from the standardly taught cases of fixed boundary field values [16] or located sources. Modern optical scanning microscopes make it possible to probe and map directly even different quantities of the near-field [17], such as the electric and magnetic field intensities [18] or the density of states [19]. Applications of tayloring nano-structures with respect to optical properties include resonant particles [20,21] and cavities [22,23], squeezed fields [24], wave guides and their adressing [25,26,27] as well as transmission apertures [23,28] and lithography masks [26]. We here present the Green's functions formalism that forms the bases of a finite-element quite effective numerical algorithm used in current research [29,30,31]. This treatise is also given as an application example of the general methods to solve differential equations with certain perturbations presented in the preceeding paper.

II.2. Problem
The discussion can be reduced to monochromatic light, that is a single frequency ω and thus time dependence e −iωt for the fields. With non-magnetic materials the problem is to find the solution E( r) of the wave equation (1) is a differential equation of the type (5) from section I. There is no source term on the right hand side of (1). If the source were, for example, a dipole located at some point, a source term with its oscillation strength and direction would have to be put on the right as Qδ( r − r 0 ). However, we shall see that a plane-wave source can be and is better included in (1) as it is. For simplicity we shall assume that the background medium, in which objects with different permittivities ε are located ( Fig.1), is vacuum with permittivity ε 0 .
For another embedding medium, its dielectric constant ε b would take the role of ε 0 . In (1) ε( r) of the material distribution designates the dimensionless relative permittivity with respect to vacuum or the background medium. To separate (1) into a part representing a homogeneous differential equation with known solution and an inhomogeneity write its as The correspondances to the quantities of the general formalism given in section I are and φ 0 will become E b . We can write down the solution following section I after having prepared the background Green's function in the next section.

II.3. Background Green's function
In any case we need the Green's function g( r, r ) of the homogeneous problem satisfying with k = ω/c. g will be a 3 ⊗ 3 tensor or matrix here. More commonly [32] the small letter g is used for the scalar function with k = ω/c. The Green's function g from (2) is then named G h with index h for homogeneous. There are several ways to obtain g. One is based on the knowledge that if we have a scalar function Ψ( r) solving Ψ + k 2 Ψ = 0, then F 1 = ∇ × ( aΨ) and F 2 = ∇ × ∇ × ( aΨ) with a constant but arbitrary pivot vector a will both solve the vectorial equation − ∇ × ∇ × F + k 2 F = 0. The tensor g looked for in (2) can be constructed out . We shall not enter into the details of this mathematically slightly precarious approach [33]. A second recipe just mentioned here for completeness is given by the following statement [32]: If g scal ( r, r ) satisfies (3), then is the tensor defined by (2). Of course, because of being in homogeneous space g scal and g effectively are functions of R = r − r alone. 1 is the unit matrix in 3 by 3 cartesian coordinate space and ∇ ⊗ ∇ means building a matrix out of derivatives We shall deduce g from a physical reasoning. From standard electrodynamics [34] one has the electric field of an oscillating dipole It is important to take the exact formula here including retardation in contrast to common near-or far-field approximations. (5) gives the space part, the time dependence is just e −iωt everywhere. To get the Green's tensor, we have to evaluate from (5) what field components in x-, y-and z-direction a dipole p at r oriented along x would produce at r, what components a dipole oriented along y would produce and what components a dipole along z would produce and assemble all these in a matrix. A point dipole is the elementary excitation corresponding to the δ on the right side of (2). The physical meaning of g is to tell us what field any such dipole would have. That is shown formally in the first matrix in (5a). Decomposing any p into its cartesian components, (5) can be rewritten as from which we easily see that the matrix to be multiplied with p to produce E( R) is R without arrow means the absolut value | R| = √ X 2 + Y 2 + Z 2 and X, Y , Z stand for x−x , y − y and z − z , respectively. R ⊗ R is the matrix written with X, Y and Z from (5a).
The terms from (5) vectorially oriented along p cause the diagonal matrix contribution, those stemming from terms with R( R p) the full matrix in (5a). Compared to the above expression g( R) from (4) has a minus sign and misses a factor k 2 /ε 0 . As will be discussed later, the source to put into equation (1) corresponding to an oscillating dipole is not the dipole moment p itself, but −µ 0 ω 2 times p. And because − k 2 The formula (6) fails for r = r or R = 0. g( R) including the case R = 0 can be represented using the principal volume method [32]. In practice, working with finite elements, the value to put for g of its two arguments the same place can be derived from the polarization of a dielectric body. The discussion of g( r, r) is postponed to the next section.

II.4. Solution for the field
The starting point to find a solution for the electric field with the objects present is eq.
(6) of section I, which rewritten in the variables of our problem here reads where E b ( r) is a solution of − ∇ × ∇ × E b ( r) + ω 2 c 2 E b ( r) = 0 or already assumed to be the space part of a linearly polarized plane wave E b ( r) = E 0 e i k r with a fixed amplitude vector E 0 . The time dependence e −iωt can be omitted in E b ( r) as well as in E( r). From (7) a numerical procedure can be deduced if the objects with ε( r) = ε 0 only occupy fractions of space rather small on the scale of the wavelength, not at all principally necessarily much smaller than λ, though. Then we divide them up into finite elements ( Fig.1) of volume ∆v, which we enumerate i = 1, . . . , N and associate permittivities ε( r i ) or perturbances V ( r i ) = − ω 2 c 2 (ε( r i ) − 1) and local fields E( r i ) uniform over ∆v. The linear sizes of the object cells should not exceed about λ/10. They need not at all be placed on a regular grid. And the only reason for demanding small enough objects is not to get too many elements N . The integrand in (7) only exists at places where V does not vanish, and to evaluate the field E( r) also at such places r ∈ { r i }, no values outside the objects appear in the equation. Changing to finite elements we thus get a linear system of equations for the fields in the object cells which can be solved by a matrix inversion. To evaluate the resulting field at any other place, that is outside the objects, the E( r i ) just have to be inserted into the finite-element version of (7): Of course, we already needed g( r i , r j ) with r i = r j to set up the system (8). Let us suppose that there is only a single cell with an ε differing from the background ε 0 . This is placed into a homogeneous field E b . If the cell has the shape of a sphere, the local field throughout its inside is aligned in the direction of E b and its value is E loc = 3 2+ε E b [35,36]. No retardation effects have to be considered here as the size of the cell can in principle be made arbitrarily small. In (8) only keeping the term of the sum with i = j, setting from which follows that The factor 1/3 is also valid for cubic elementary cells, however, other shapes require different depolarization factors [32,35].
The background and the resulting field at each point are already 3-vectors. Nevertheless, in order to solve (8), imagine the Es for the object cells assembled into long or "double" vectors of N times 3 components With further the big 3N × 3N matrix M consisting of 3x3 blocks M can be inverted using the procedure described in the appendix, however, as M is a full matrix and the complete invers M −1 is needed, that is of no advantage and a standard inversion algorithm will do as well.
For the problems of a few small scatterers here by introducing finite elements the implicit integral equation (7) for the electric field has been turned into a linear system of equations that is easily solvable. Even shying this effort, the coarsest, so-called Born approximation consists in replacing E in the integral in (7) or in the sums in (8) or (9) by E b . Keeping the summation over finite elements to estimate the integral we directly get No distinction between places r in or outside the object cells is necessary in (11). The Born approximation only takes into account first-order scattering off every object and thus can only be good for weak scatterers with distances between them rather large on the scale of the wavelength. Producing a clearly different field pattern from the exact solution including all scattering orders, Fig.2 demonstrates that the Born approximation is likely to be insufficient to model near-field optics setups.

II.5. System Green's function and density of states
Adding an arbitrary source term to our original wave equation (1) changes it into If now we know a tensor function G( r, r ) satisfying then obviously would give a special solution of (12). Any solution of (1) could be added. As generally deduced in section I the implicit relation to get G from is the Dyson equation G( r, r ) = g( r, r ) + d 3 r g( r, r ) V ( r ) G( r , r ) G( r, r ) for both arguments r and r covering all space is too much information to display at once and usually much more than what one is interested in. The imaginary part of G( r, r) is proportional to the density of states ρ [31,33]. The deduction of this statement found in quantum mechanics book [37], however, rather argues with a system of energy eigenstates and the variation of the Green's function as well as the density of states with energy. No real r-space is explicitly mentioned. Our interest lies in the spatial dependence of the density of states at fixed light frequency ω.
Even if described in terms of fields, concepts like reactance and work known from electrical circuits may be applied [38,39]. The time average of the work done by the fields is given by With no other imposed fields, charges or currents than an oscillating point dipole, the latter will present the only external current J, which will thus be located as δ( r − r).
If the dipole moment oscillates as p( r , t) = p 0 e −iωt δ( r − r), the corresponding current is . Deducing the wave equation for time harmonic fields (in vacuum for simplicity here) from Maxwells equations with current term from which we see that the source term for the dipole has to be set as Q( r ) = −iωµ 0 J( r ) = −ω 2 µ 0 p 0 e −iωt δ( r − r). The integral (16) reduces to the value of J * E at r. The electric field we get from (14): Inserting J * and E into (16) the time factors cancel as expected for a time average and but for a factor ω 3 µ 0 /2 we get Choosing unit vectors along the coordinate axis for the probe dipole p 0 , (20) will filter out the trace elements of the matrix G( r, r). We associate ρ A motivation for taking the negative imaginary part of G( r, r) as a measure for the presence of modes can also be obtained by comparison to the energy resonance of a forced oscillator [40]. For optimal excitation from the energy point of view -in contrast to amplitude resonance -the force has to be π/2 ahead of the elongation or in phase with the velocity of the oscillator. G( r, r) describes the field caused by backaction of the system at the place of the probe dipole moment (taken as reference phase zero), and therefore −ImG( r, r) is the part that can in a resonant manner further enhance the dipole oscillation. (In reality radiation out of the system will provide strong damping.) We now intend to evaluate a map of G( r, r) on, for example, a horizontal plane. The plane may lie above or below object cells or even cut some. Like the objects, the, of course, finite area of interest on the plane is divided into cells. Just depending on the desired resolution of the map the unit cell length of this mesh may well differ from the cell size chosen to discretize the objects (Fig.1). The list of object cell midpoints { r i } from the last section, which shall be called region A, is extended by all cell midpoints from the map in the plane, which shall be called region B and is now understood to to be included in counting i from 1 to a new N . Analogously to (8) the integral in (15) is replaced by a sum: It does not matter that the map B has a different mesh from ∆v as for the objects A, as  (8), the only difference being that j now also runs over the plane cells in addition to the object cells.
G and g are themselves matrices on contrast to vectors E and E b . G consists of 3x3-blocks G( r j , r k ), g is made of blocks g( r i , r k ) and the 3x3-block at position (i, j) in M given by as it is not M that is singular. Quadrant BA will be needed for multiplication with g in (21a). However there is an even more efficient algorithm to get G that already includes the multiplication by g. It directly calculates G = (1 − gV ) −1 g, which is the compact way to write (21a) as the solution of (15), also denoted G = g + gV G for short. The technical details can be found in appendix B.
Trace components of G( r, r) meaning densities of states for the three polarization directions ( Fig.3) above an optical coral [31] in analogy to a quantum coral [41] have been measured [19] in a so called forbidden-light near-field optical microscope [42]. The sample consists of a stadium arrangement of gold particles on a glass surface. The forbidden-light setup prevents detecting light emitted from the fiber tip that has not passed through surface modes that make up the density of states for this system. Like for antinodal and nodal points in a resonator, more energy can go into the system when the excitation is placed at a point of high density of states than when coupling is bad where the density is low.

II.6. Remarks on the source terms and alternative solutions
In section 4 we saw that it is convenient to start from a solution for the field in the form (7) if the excitation comes, for example, from a background field belonging to a plane wave.
Though the matrix to invert bore a certain similarity to the evaluation of the Green's tensor in section 5, with (8) and (9) we directly calculated the field. In contrast, more adapted to localized sources, there is (14) as a solution of (12). If there is no additional background field to cause any excitation, no solution of the equation (1) with zero right side is to be added as further contribution and (14) is the field distribution to be observed. (14) has to be rewritten in terms of finite elements in order to be used in a numerical calculation. In the same way as the objects the source Q has to be devided into discrete cells or elementary dipoles. To distinguish their L locations from those of the objects we shall enumerate them as ρ α , α = 1, 2, . . . , L. The place r to evaluate the field E may be anywhere outside or inside the objects as well as beside or even at a source location. For the following development the Dyson equation for the Green's tensor is needed in a discretized form for both its variants Having in mind a region where and a resolution with which E( r) is to be evaluated like the discretized plane B from the last section, it would be possible to supply G for all needed combinations of arguments ( r, ρ α ) and calculate E( r) as the single sum from (22a). To weave in the influence of the objects, G would have to be set up as a big matrix like in the last section over all combinations of three regions A, S and B here, the objects, the source and the map. Having calculated G in the A-B-scheme from the last section, one could evaluate E( r) with G( r, r j ) as written in (22b). However, the most efficient way is given in (22c).
G( r i , r j ) is merely needed for r i and r j from the set of object cells, keeping a matrix to be inverted as small as possible, namely of (AA)-type. Choosing r i and r k in (21) in the object set A, instead of in A or B as the equation was originally set up for, we see that (21) presents a closed system of equations for all such G( r i , r k ). Then for (22c) more summations over products with g-functions, which are analytically known for any pair of arguments, can be considered less demanding in computing time than the inversion of large matrices.
In the transformation from (22c) to (22d) after swapping index names i and j in the last sum, G = g + GV g and (14) have been exploited. Using (22d) for E( r) renders an implicit equation for the field in the form (10) or (10*) from section I. E( r) is the equivalent of φ ihq and as stated earlier, we assume that physically there is no background field φ 0 that could initiate an additional field distribution φ ih . Although not very convenient, a plane wave as exciting field could be understood as stemming from a sufficiently long and dense array of Huygens elementary dipole sources Q reasonably far away from the objects. The other way round, for a single dipole source or a number of dipole sources distributed in space the field they would produce at any location r in homogeneous space is the superposition of their individual fields, and putting E b ( r) = L α=1 g( r, ρ α ) Q( ρ α ) the ansatz (7) can be used also for this case.
Like the field E anywhere was obtained as a straight-forward summation once having its values at the places of the object cells, finally an alternative way to the procedure from the last section to get the Green's tensor G shall be given, also requiring only the inversion of a matrix with size the number of object cells. Series expansion is used to rewrite the solution of the Dyson equation: Designating regions the spatial arguments belong to on (23a) we get for G BB : As V does not vanish only in region A, the first index of (1 − gV ) −1 obviously must be A.
(1 − gV ) −1 being ∞ n=0 (gV ) n , for power zero the second region index automatically is the same as the first and all other powers ending with V imply second index A. It is sufficient to set up the matrix 1 − gV as an (AA)-block and invert that. Should one prefer to evaluate a complete G, which differs from the above inverted matrix by a factor g, line (23b) like (22c) shows that it is in principle only necessary to get G from some self-consistent implicit equation in the object region A. G AA can be constructed applying the procedure described in appendix B to a matrix set up as (AA)-block only. The inversion has to be completed in this case, though. Going through the diagonal elements, all lines and columns have to be updated in each step, including the ones above and to the left of as well as the ones the respective diagonal element is in. Writing (23b) as summation over discrete elements ready for use in a calculation then reads: Of course, summations run over all object cells here. No numerical advantage can be drawn out of r = r in (23b"). With the same effort of making G AA it can be used to evaluate maps of G( r, r) as well as plots of G( r, r ) with r fixed or even some function of r.

II.7. Conclusions and Outlook for section II
We have presented a method to solve the problem of scattering of electromagnetic waves off an arbitrary distribution of dielectric objects, that is the exact evaluation of the field, especially in the near zone where higher-order multiple reflections can become important.
Besides the field distribution we have obtained the Green's tensor characterizing the system independently from the form of the excitation. It represents the response function and also the density of states for supported electric fields.
For a methodical introduction we have restricted our considerations to dielectric materials and the electric field. Without magnetic susceptibilities the magnetic field distribution can be calculated once having the electric field inside the objets by a formula like (9) with the magnetic background field and replacing g by a tensor including the conversion from the electric to the magnetic field by taking the rotation [43]. It is further possible to treat non-uniform magnetic permeabilities and even mixed systems with dielectric and magnetic objects [30]. The electric Green's tensor presented above is then paired by a magnetic counterpart and genuine mixed response functions also exist. Whereas the calculation of the field distributions even for mixed systems is quite straight forward, the construction of the Green's tensor is more involved. It lives of the idea of handling one kind of objects first and then considering this setup as the background to include the other kind. There is no approximation or ranking in importance in this procedure.
The discussion here has only considered finite objects in a homogeneous background as well as cartesian coordinates where vector and tensor components have been written out. Cylindrical and spherical coordinates are also commonly used [32,44] and the Green's functins formalism has been developped for layered media [26,32,45]. Besides waveguide applications the use for modelling typical near-field optics experiments, where the microstructures to investigate are prepared on a substrate surface, lies in putting the influence of this surface into a background Green's tensor [33,43] which is then implied the way we used g here.
Details of applications of the Green's functions technique in electrodynamics to more complicated situations as well as beautiful results of corresponding experiments can be found in the given references. This text focussed on calculation techniques and further intended to give an overview of slightly different formal ways to calculate Green's tensors and fields of which either may be optimal for a specific problem.
The inverted matrix (1 − M ) −1 can be represented as a geometric series: Truncating and using the sum from the right side is only possible if the series converges whereas the closed form on the left is valid in any case. In contrast to the infinite sum on the right side of (25) we see that the first step in (24) adds to each matrix element m ab the sum of all products m a1 m 11 m 11 . . . m 11 m 1b . In these at least one pair of indices 1 is squeezed between a and b as in m a1 m 1b , the contribution m ab was already there. In the second step all products In the application from the main text i = 1, . . . , N enumerates the object cells. At position (i, j) in (1 − gV ) −1 when expanded into a series having every possible sequence g( r n 1 , r n 2 )V ( r n 2 )g( r n 2 , r n 3 )V ( r n 3 )g . . . g( r n j−1 , r n j )V ( r n j ) shows that the resulting field at any place (inside or outside the objects) is the interference of the background field and the fields reradiated by all the object dipoles having undergone every possible scattering path between the objects (Fig.1). A complication in the electrodynamics application at this stage is the fact that each matrix element m actually in itself is a 3x3 matrix indicating the effect of three field components at one place onto three field components at a another place. g on the discrete space of the objects can be written blockwise with a scheme of quadruple indices V ( r k ) in the product g( r i , r k )V ( r k ) just multiplies the respective column.
One could use the inversion procedure working off the diagonal elements marked by ovals, requiring 3N steps then. However, the process is equally applicable to 3x3 blocks as marked In the following instructions are given how to calculate the Green's tensor being efficient in the way that finally only G( r, r) of equal arguments on the mesh points of map B need to have the correct values [30,46]. For a start nevertheless consider that G( r i , r k ) with arguments r i and r k from the big set of all object and map cell midpoints will have to be the sum of all products g( r i , r n 1 )V ( r n 1 )g( r n 1 , r n 2 )V ( r n 2 )g( r n 2 , r n 3 ) . . . V ( r n j−1 , r n j )g( r n j , r k ).
r n 1 , r n 2 , . . ., r n j can only be object cells and any sequence of them has to be created, including the empty one with no V giving the term g( r i , r k ). Initiate a matrix (1 − G( r n , r n )V ( r n )) −1 and to all entries from line n + 1 and column n + 1 on add the product given below.
In the following step use the updated entries in the above recipe. In the nth step you need not update entries in line n or above or in column n or to the left of it, because these will not be needed as multiplication factors G( r i , r m ) and G( r m , r k ) for m > n any more. The parts of quadrants AB and BA remaining after these restrictions have to be changed by (26) for every n. Merely diagonal 3x3 blocks have to be done in the BB-quadrant. The first step adds all products consisting of any non-zero number of factors V ( r 1 ) between gs, the second step all products eventually containing V ( r 1 ) and one up to any number of V ( r 2 ), and so on. The procedure is finished after r n has run down the diagonal of the AA-quadrant. All sequences of multiple scattering from the objects are then included.
We could have filled the whole BB-quadrant with initial g BB 3x3-blocks and updated the entire BB-quadrant in each step. Like in the matrix inversion procedure from appendix A we should further have worked down the complete diagonal and done (26) for every r n from region B as well. There will, however, be no additions as V ( r n ) = 0 for r n in B (even if a map cell accidently coincides with an object cell). This argument also reveals why non-diagonal elements in BB do not have to be evaluated. They can never appear as multiplication factors G( r i , r n ) or G( r n , r k ). The BB diagonal is the resulting G( r, r)-map we wanted.
Initiating the matrix by g and introducing the V on treating the respective diagonal element makes the outcome of the procedure directly (1 − gV ) −1 g compared to initiating the matrix with gV and as an intermediate step obtaining (1 − gV ) −1 .

III. ALL-ORDER QUANTUM TRANSPORT
It is demonstrated how the transport problem for two open free-electron gas reservoirs with arbitrary coupling can be solved by finding the system's Green's function. In this sense the article is an introduction on Green's functions for treating interaction. A very detailed discussion of the current formula is given on an elementary basis. Despite formal resemblances the stationary transport situation, however, differs in its nature from introducing coupling between energy levels in a closed system where then the interest lies in modified eigenvalues and eigenstates.

III.1. Introduction
By modern lithography techniques so-called point contacts [47] can be arranged between conductors. These have sufficiently small dimensions such that electronic modes get quantized. However, coupling across such constrictions need not be so weak as to be described by a small tunnel probability, but can be influenced by coherent interference of multiple reflections. Point contacts can be obtained by indenting STM-tips into some material [48], electromigration [49] or the break-junction technique [50]. Whereas for constrictions imposed by gate electrodes to a two-dimensional electron gas in semiconductors [51] one observes quantized conductance values in the sense of fully transmitting or totally switched-off modes, the application in mind behind this work is the type of connection like the single-atom contact, characterized by an ensemble of channels [52], which can also have intermediate transmission amplitudes between zero and one [53]. Even with some fully transmitting modes the contact bears a resistance in the order of the quantum resistances R k = h/e 2 = 26kΩ [54], such that viewing the system as a left and a right side with some interaction is an appropriate picture. Furthermore contributions from several channels just add in the current. Although a transmitting channel in a point contact is the application in mind, this article shows how to set up a general procedure to solve the problem of transport for two open reservoirs with more or less strong coupling between them, and thus how the Green's functions formalism from section I is implied in an area of current research [55,56].

III.2. Green's functions formalism
As a preparation, consider systems like, for example, the bulk material on the left or the right side (Fig.1a) without coupling. For these we suppose that we know the Hamiltonian H 0 and the wave function Ψ 0 at any given energy ω satisfying the Schrödinger equation.
Denoting the Hamiltonians H LL = H RR = H 0 (doubling the index makes sense later), the wave functions Ψ 0 L = Ψ 0 R = Ψ 0 and the time derivative as ∂ τ , the Schrödinger equations for both sides -for the moment just formally put into matrix form -are Corresponding to this differential equation we have the Green's functions equation If we knew g = g LL = g RR , we could immediately also give a solution to (1) with a source However, we are more interested in the solution when an interaction between left and right is present, expressed through coupling Hamiltonians H LR and H RL , and we shall refer to this case by Ψ without upper index. Introducing the coupling, of course, was the motivation for writing the Schrödinger equation in matrix form over the site space consisting of L(left) and R(right). Putting the coupling terms on the right side of the equation, they mimic a source and following (4) the solution can formally be written as obtaining the implicit Lippmann-Schwinger equation for the wave functions. Although for the coupled system we do not expect the eigenvectors of a discretized form to obtain field distributions, for our purposes here (7) will merely serve as a formal step in the derivation of the Green's function. What we are precisely interested in from a physical point of view is the current that will flow between the left and the right side, and not necessarily explicitly evaluating the wave functions.
The Green's function G of the coupled system is inferred from an equation analogous to G like the Hamiltonian is a full matrix in site space. With G, a solution of where the index cq stands for coupling and source. As explained in section I, we could also have set up a Lippmann-Schwinger equation for Ψ cq as Now from inserting (10) into (11) and for Q choosing some δ(τ − τ 0 ) either in the L-or in the R-component the Dyson equation for G is obtained:

III.3. Explicit Green's functions in time and frequency domain
The reservoirs on the left and right being bulk metal, the differential operator of the homogeneous differential equation (1) is given by the free-partical Hamiltonian and the corresponding wave functions are Ψ 0 = e ±i k r−iωτ according to the dispersion relation Looking for a Green's function at given frequency ω, in (2) we can replace H 0 by ω and thus have to solve We need an expression, the derivative of which produces δ(τ − τ ). One easily verifies that (14) is satisfied by either The retarded function g r only exists for τ ≥ τ and the advanced function g a for τ ≤ τ . In (14) for simplicity we restricted ourselves to a single energy and therefore (15) Taking out the factor D/h, we get the dimensionless functions g r (ω) = −i and g a (ω) = i in frequency space. Remark that their deduction here did not consist in calculating a Fourier transformation. The phase factor e −iω(τ −τ ) was already there in (15). The ω-integrals from (16) will be implemented as a useful representation of g r/a still in the time domain. The θ-functions have been deliberately skipped after the =-signs. We shall however see that formally even as double Fourier transform (however with different signs in the exponentials with τ and τ ) h D g r/a (τ, τ ) = dω 1 dω 2 g r/a (ω 1 , ω 2 ) e −iω 1 τ e iω 2 τ with g r/a (ω 1 , ω 2 ) = ∓iδ(ω 1 − ω 2 ), however, because different energies stay independent.
(More or less guessing the Green's functions g was easy in our normal-conducting simplest model here. Generally, it has to be found as the solution to an equation like (14) with the differential operator from the uncoupled system's Schrödinger equation and an elementary perturbation δ. In the superconducting state, for example, the two-particle interaction in the Hamiltonian forces the Green's function to include Andreev reflection [57], and the result is such that g(ω) does not just vanish in the gap of the quasi-particle density of states.) Before we can solve (12) for the coupled system's Green's function G, we have to specify the coupling parts of the Hamiltonian H LR and H RL . Like with the density of states the simplest model will assume that the coupling is energy-independent and described by a constant (real) interaction energy W . The reservoirs are unaltered by transport between them. Due to the applied voltage (Fig.1b) incoming charge carriers (from the left on the right) are led away and outgoing ones (on the left to the right) get replaced. In the contact we do not allow relaxation or other energy-changing processes. We choose the repective Fermi levels as zeros of energy on either side (Fig.1c). An electron going from right to left has to strip off its phase e −iω R τ and aquire e −iω L τ to fit in on the left. An analogous argument holds for transitions from left to right, and therefore the coupling terms are The phase factor gives a time-dependence to H LR and H RL , but ω R − ω L = eV / is indeed independent of energy. For G we make an ansatz like (17) as a two-fold Fourier From (12), which is valid for either advanced or retarded functions, as an example, we pick the upper left component of the 2x2 matrix in LR-space and insert (17), (18) and (19): (Multiple integral signs are skipped from now on.) The integral over τ produces 2πδ(ω L3 + eV − ω R ) and in g( anyway, such that the last term of (20) becomes with t = W D/ . Strictly speaking, if (20) were for the retarded function, in the last term there would be θ(τ − τ ) from g r LL and θ(τ − τ ) from G r RL , and if it were for the advanced function, θ(τ − τ ) and θ(τ − τ ), such that the integral over τ only exists between τ and τ instead of having minus and plus infinity as limits. However, with ω L3 and ω R running over any value, one can argue that even with a finite τ -integral the only remaining contribution stems from ω L3 + eV − ω R = 0. A discussion of time ordering will again appear in section 5.
D is the density of states per frequency interval, dividing by makes it the number of states per energy interval. W is an energy. t can be understood as a dimensionless transmission amplitude. Now we set up the convention that all frequency arguments of g and G are written with respect to the left zero level and for the case that they correspond to the right, that is an R-index, it is understood that eV / is added. (20) has to hold for any τ and τ and from comparing Fourier coefficients we get We could have inserted G RL = g RR H RL G LL and again replaced G LL = g LL + g LL H LR G RL and so on. Instead of an implicit equation for G this would have led to an infinite series (see section I): from σ passes the frequency argument from the g in front to the g behind. And as g can only have two identical frequency arguments, no ω different from the first can ever appear, such that ω L1 = ω L2 = ω and G also effectively is a function of only one frequency argument: . With Green's functions of a single frequency argument (20a) and its analogues for the other three components in LR-space become a simple algebraic equation: (The fact that even the Green's function of the coupled system turns out to be a function of a single frequency argument is a special feature of our simple model for the normal conducting case. In the extension of this model to superconducting reservoirs and transmission processes including Andreev reflection, G becomes a function of two frequency arguments, the second, however, restricted to values differing from the first by an integer multiple of eV / , such that effectively there is one continuous and one discrete frequency parameter [58].) Here, with g r/a (ω) = ∓i, (22) is easily solved and G(ω) comes out independent of frequency, too: The same result could have been obtained from (21) by writing out a few more of the matrix multiplications and using the formula for the geometric series in each element. We shall need two further types of Green's functions. T will be introduced in the next section and g +− and G +− when calculating the current.

III.4. Transfer Green's functions
In (21) there was a sum of products of arbitrary many factors g and σ with outer factors g.
Whereas all contributions to G in (21) began and ended with staying some time in a reservoir, described by g -at L from τ to τ for a start in the last term of (20), for example -each term of T in (24) begins and ends with a transition σ and further contains at least one such hopping across the junction (which g does not). Therefore we call T the transfer Green's function. The same as for G, the relation between T as a function of times and as a function of frequency is given by where ω J = ω if J=L and ω J = ω + eV / if J=R and the same for ω K . (Taking out 2π of D in (16) was a convention. The prefactor of the ω-integral for T follows from consistency requirements. (19) for G was in complete analogy to (16) for g. T with two time arguments, however, has a little different character from σ with just one time parameter.) Alternatively to (24) T could be defined through its link to G σG = T g or Gσ = gT (26) Be careful that replacing one by the other can introduce another internal time integration as, for example, σ(τ ) G(τ, τ ) = dτ T (τ, τ ) g(τ , τ ). (24) and (26) hold for retarded and advanced functions. From (24) it is immediately seen that T like G satisfies a Dyson equation G = g + gσG and T = σ + σgT (27) and even the complementary forms G = g + Gσg and T = σ + T gσ (28) are analogues. In Fourier space, like (22) the T -equation (27) is an algebraic equation and the solution like Inserting g and σ explicitly for our model we get Whereas t is the single hopping amplitude, T LR/RL is a renormalized transfer amplitude.
One may wonder why a model for transport could not have been set up adding amplitudes for transfer processes of all orders, as the interaction (18) seems to be introduced the way it is just in order to result in powers of t. However, g r/a (ω) = ∓i deduced from the Schrödinger equation is decisive for the signs in (30) and (31). One may wonder that multiple reflections are not added as t + t 3 + t 5 + . . . = t 1−t 2 . Processes of different order (Fig.2) are not independent, but interfere. T RL is the transfer amplitude per single electron supplied on the left by the voltage source. But for every electron that goes over to the right in an nth order process (n odd) with weight t 2 there is one that has hopped once more to the right and back (n + 2 order process) and thus is not to be newly supplied, but to be again sent through the junction. The amplitude t is renormalized by 1 + t 2 as denominator. However, such interpretations of quantum mechanical amplitudes are precarious, and the full conversion of t to a transmission probability will be established later.
It is quite instructive to solve (27) in a slightly different way than done in (30). Firstly, for the four components in LR-space we have Inserting these into each other, for example, an equation for T LR alone is obtained: This implicit equation is the basis for calulating the transfer Green's function in more complicated cases than discussed here [58,59], like for example the superconducting junction.
In our model, inserting σ LR = σ RL = t and g More easily than the Dyson equation for the ordinary Green's function G, the one for the transfer Green's function T is illustrated as is done for the LR-component in Fig.3.
(Normally indices are read from right to left such that T LR is considered a transition from right to left, but it does not really matter whether they are interpreted the other way round as in Fig.3. The actual sequence of what is earlier or later in time will be discussed when calculating the current in the next section.) Fig.3 demonstrates the implicitness of the Dyson equation: Any transition from left to right is either a single transfer or an electron hopping to the right and back followed by any process beginning on the left and ending on the right, no matter what happens in between. This last part by definition is the sane as the other side of the equation, namely T LR .

III.5. Calculating the current
From the Heisenberg picture of quantum mechanics we know that the time derivative of a not explicitly time-dependent operator A is given by the commutator with the Hamiltonian [60]: The operator of interest here is the projector on either side of the junction As explained earlier, with the junction coupling left and right together, the solution is not limited to one side, however, the projectors take out the respective part: ρ L and ρ R are proportional to the amount of charge on the left and on the right side. Their time derivatives have equal absolute values, but opposite sign and represent the current.
Putting together (33) and (35) H LR (τ ) and H RL (τ ) denote the time occurence of τ from (18) which has to be considered as still belonging to the Schrödinger picture in our case here. T means time ordering [61] and T anti-time ordering. Ψ 0 in the uncoupled system, of course, also stands for a two-vector with left and right component. Replacing (37) or σ(τ ) = |Ψ 0 (τ ) > W < Ψ 0 (τ )| for short, where bra, ket and W still mean the respective matrices. However, this decomposition might rather be confusing and will not be used, anyway. The translation to the interaction picture is the following: or σ H (τ ) = |Ψ 0 (τ 0 ) >< Ψ 0 (τ )| σ(τ ) |Ψ 0 (τ ) >< Ψ(τ 0 )|. The Heisenberg picture always refers the operator back to the undeveloped state at τ 0 .
Analogously σ c will have to be extended by the time-dependence of the uncoupled states.
The meaning of the time-integral over σ H as an exponential is best explained by writing explicitly: where all arguments τ 1 , . . . , τ j , . . . , τ n have to lie between τ 0 and τ and all products from none to arbitrarily many factors σ H have to be added. With time ordering, furthermore Without T the scheme for the exponential series . . ≤ τ 1r ≤ τ and integration over all time arguments except τ and τ 0 is understood. Like in (39) we mean the sum of all products with arbitrary many inner development factors like Ψ 0 of τ k and Ψ 0 of τ j . Spaces and different font types are just used in (37a) to recognize sequential σ H -parts like (38).
Inner development factors like Ψ 0 of τ k and τ j as well as Ψ 0 (τ ) and the inner Ψ 0 (τ 0 ) have to sum over space or the basis of states and therefore should be thought of as matrices as in (38). We have not yet decided which components outer bra and ket with Ψ 0 (τ 0 ) will project out. All < Ψ 0 (τ 0 )|Ψ 0 (τ 0 ) > are the unit matrix and drop out (Ψ 0 being given by a non-normalizable e-function does not pose a problem here).
We regard the coupled system as having developed out of the uncoupled system, but we are looking for a stationary state. To achieve this, the coupling has to have been turned on infinitely long ago, thus we let τ 0 → −∞. The sequence of time arguments is usually represented on the so-called Keldysh contor (Fig.4). We still have to take the expectation value of (37) or (37a). This means summing over the basis of states for the outer Ψ 0 (τ 0 ). Only being interested in the trace as a result, in a matrix multiplication the order of factor matrices can be changed cyclically. As in (37b) every bra and every ket as well as each operator part between the two | is a matrix we can rotate factors to obtain

The uncoupled basis consists of
however, with no restriction as to which argument τ n or τ m is earlier or later in time. We define this new type of Green's function as (An eventually ill-defined single point τ = τ is irrelevant for later integrations.) As the cases τ > τ and τ < τ are mutually exclusive the function can also be given by g +− (τ, τ ) = g a (τ, τ ) − g r (τ, τ ) (conclusions on g +− (ω) from that are risky, to my opinion, though). We have thus deduced the current formula (like discussed in (20) it is rather irrelevant whether the upper integration limits are set as τ or ∞) or or in short notation. The current only needs G +− of two identical time arguments. We shall need the Fourier representation of g +− (ω), but none such for G +− . Although in our simple model for the normal-conducting contact the current will come out the same for any τ , I(τ ) generally does depend on time (our expectation value is an ensemble, not a time average). For the superconducting case the ac parts [58] like the Josephson current are included in the expression (41), although one is mostly interested in the contribution in I that is independent of τ and gives the dc part.
Using g +− (τ n , τ m ) without a third argument ω in (41) we understood integration over frequency like in (16). But this point requires more care as it has not been taken into account so far whether states are occupied or empty at τ 0 → −∞. Writing out the trace from (41a) in LR-components reveals four contributions It is clear that in contrast to the other appearing σ, in σ c from (36) there is an additional relativ minus sign between the LR-and the RL-component. Obviously the net current is the difference between the current from left to right and the current from right to left. In the first -as well as the second -term in (43) a transition from R to L at τ , the time argument of σ c , is picked out to be counted for the current (see also Fig.4). With g +− at R, the charge carrier is supplied from a state originally located on the right. It is an electron if the energy lies below the right-side Fermi level. Until τ the state has evolved to again be on the right. The originally left behind empty state at the same energy below the Fermi level must have evolved to be on the left at τ such that by the R→L transfer the electron can go into it (Fig.5a). Or you might say that the plus and the minus branch of the Keldysh contor represent two possible parts in the evolution of an original wave function Ψ 0 R , between which there is a non-vanishing matrix element of the operator σ c = [H, ρ L ]. However, still regarding the first line of (43), there is the further possibility that an empty state from above the Fermi level on the right evolves to be at R at τ again, but its left-behind complement (a negative charge) has evolved to be at L. σ LR (τ ) in this case means the transition of an empty state or positively charged particle from right to left (Fig.5b). Although one does not usually introduce the concept of holes with transport in normal conducting metals, it makes sense here to call unoccupied states simply "holes". In this way the model already includes the dual nature of charge carriers needed for the superconducting case. In the normal conducting case states do not change in nature (electron or hole) or energy during their evolution (little zigzags are drawn in Fig.5 only to make the multiple hoppings visible). At the superconducting junction, Andreev reflection can be interpreted as changing an electron into a hole or vice versa and mirroring its energy at the Fermi level [59] (see Fig.9 in section 6). To have a charge carrier at a certain energy level at time τ to make a certain transition, it is important that there was one at the corresponding energy in the original uncoupled system. The coupling may have changed the distribution compared to the occupation in uncoupled bulk reservoirs. And the applied voltage imposes a non-equilibrium situation, anyway. For the evolution of a state from the right as shown in Fig.5 it does not matter whether the state at the respective energy on the left is occupied or empty. If, for example the regarded energy lies below the Fermi levels both left and right, there will be two states evolving as an electron on the plus and a hole on the minus branch of the Keldysh contor, one having originated at R and the other at L. These original, uncoupled and independent states are our basis, especially for calculating an expectation value as trace. They do not interfere. Terms with g +− RR and g +− LL are simply added in (43). Schemes analogous to Fig.5 could be drawn for the terms from the last three lines of (43) as well. The conclusion of the whole argumentation of how to let the original Fermi occupation function for the reservoirs left and right enter the current calculation is that g +− LL has to change sign at the left Fermi energy and g +− RR at the right Fermi energy. Let us note g +− like g r and g a for any bulk reservoir with Fermi level at ω = 0. We shall keep g +− corresponding to g r and g a as in (40) for occupied electron states below the Fermi level and change the sign for empty states above it.
Although it might be practical to use (19), (23), (44), (18) and (36) in (41) to quite directly produce an expression that calculates the current finally as an integral over frequency and in our simple model can even be analytically evaluated, in parallel to [58] we shall use the transfer Green's functions here. The current is best translated into an expression of the transfer functions from the form already resolved into LR-components (43). Furthermore eliminate σ c through σ c,LR = −σ LR and σ c,RL = σ RL . L-and R-indices follow logically from (26), for example G LL σ LR = g LL T LR .
As quantities here are no longer matrices, but simply functions of time or frequency, the leading σ LR and σ RL in the second and third line of (43a) can be moved to the ends of the products. From (41) we remember that their argument is τ , the same as the second argument of the last factor to the right. Then using relations from (27a) and complementary forms the current formula simplifies to The terms with the factors 1 from the brackets have elegantly been made to vanish. Now we use the Fourier representations for all functions in (43b). As all terms follow the same scheme, the second is treated in an exemplary way here (all integrals run from −∞ to +∞): Doing the τ 1 ,τ 2 -and τ 3 -integrals gives (2π) 3 δ(ω 1 − ω 2 − eV / )δ(ω 2 − ω)δ(ω − ω 3 ). Then even the exponentials with τ cancel and the term simplifies to In our case (31) tells us that T r/a LR/RL (ω) are all identical and real, g r LL/RR (ω) = −i is the complex conjugate of g a LL/RR (ω) = i and −g +− LL/RR is the complex conjugate of g +− LL/RR of the same ω as g +− (ω) is purely imaginary, too. Thus the last two terms in (43b) are the complex conjugates of the first two and thus twice the real part of these first two can be taken for I/e. In the superconducting version of the model, where T and g actually are ω-dependent, complex conjugate relations [59] also exist between T s as well as gs, and the current formula can be reduced in the same way. Just to note the quite general formula [58] in short form: For the integrand from the second term in our model we get where the signs refer to ω greater or less than zero and come out reversed for the first term with g +− RR , because there is g a instead of g r . Care has to be taken with the reference point for ω in both terms. This may easily be overlooked in the normal conducting case here in contrast to the superconducting case where T indeed is ω-dependent and like G as a function of only one argument always referred to the same Fermi level (the left, for example).
If we call the argument of g +− LL from the second term ω, the one for g +− RR in the first term is ω + eV / . On an ω (L) -axis, the second term changes sign at zero, however, the first jumps at −eV / (Fig.6). A shift of the integration parameter cannot be made independently for both terms. Thus, for the normal conducting model here In principle the convention is needed, that the T -argument always refers to the left side, but for g general formula like (44) for a single bulk with Fermi level at zero frequency are applied. A more involved situation where choosing integration intervals consistently for all contributing current terms is crucial can be found in [59].
From Fig.6 it is easily seen that that the integral is twice the constant t 2 (1+t 2 ) 2 integrated over an interval of length eV / , and outside that interval contributions cancel. With the other factor 2 from twice the real part and the prefactor 1/2π the result for the current finally is The factor θ = 4t 2 /(1 + t 2 ) 2 is the conductance in units of the quantum conductance or its inverse the resistance in units of the quantum resistance [58]. (The conductance of a channel It may seem a contradiction on the one hand calculating the current from a changing amount of charge on one side and on the other hand saying that missing charges are replaced and superfluous ones led away by the voltage source. A slightly different viewpoint may help to get convinced that the calculated quantity is indeed the current in the stationary, but non-equilibrium system. The crucial point was putting on to evaluate dρ L /dt in state Ψ(t).
Without further ado we could not tell whether this state of the coupled system was occupied or not. The Schrödinger equation (5) set up the left and right material properties as well as the coupling across the junction, however, did not take any account of the effect of the voltage source. Without need to specify real locations for the division, just principally view our structure as consisting of a junction region and leads. The Ψ(t) defined through (5) describes states in the junction region. But think of them as offered by the system and   is identical to the ansatz made by putting together (33) and (34) in (37).

III.6. Comment on the two-level system
At first glance the problem posed by the Schrödinger equations without and with coupling, (1) and (5), especially if we regard a single energy level on each side with corresponding ω L and ω R as in Fig.1c, looks like the two-level system known from quantum mechanics textbooks [62]. By the coupling the two energy levels are shifted. The new eigenstates lie further apart from one another. If the system is initiated in one of the uncoupled states, it will oscillate harmonically between the two levels, and both both the period of the oscillation as well as the maximum transition probability to the other state depend on the energy for example (Fig.8), Ψ 0 L (and Ψ 0 R ) cannot be taken as eigenstates any more, because from i ∂ τ e −iω L τ = ω L e −iω L τ with ω L = ω L + eV 2 it does not follow that i ∂ τ e −i(ω L +eV /2)τ equals ω L e −i(ω L +eV /2)τ . The standard treatment of the two-level system is not applicable to the non-equilibrium situation.

III.7. Superconducting junction
Despite having been the simplest example to introduce the Green's functions scheme, applying the formalism to the normal conducting junction to get out that the current is proportional to the applied voltage was breaking a butterfly on the wheel. Although repeatedly mentioned, fully developing the extension to the superconducting case is beyond the scope of this presentation. But the resulting current-voltage characteristics shall be shown as a plea for the usefulness of the method. A different approach based on matching wave functions [63] leads to identical results, though. Some formula shall be listed, because they are not necessarily written out in complete analogy to the presentation here in [58] and other literature. Working in the quasiparticle picture, in the superconducting case each entry in LR-space of a Green's or transfer function expands into another 2x2 matrix in Nambu space over electrons and holes. There is [64]: g is the same for LL and RR. Different signs refer to the retarded and advanced function.
∆ is half the gap of the superconductor. The small imaginary part η is added to get the correct root besides slightly smoothening singularities. g(τ, τ ) = D/h dω g(ω) e −iωτ e iωτ still holds. Here signs refer to the LR-and RL-direction, respectively. Remark the reversed signs for holes in the exponential compared to electrons. There is no electron-hole conversion during a single hopping. The eh-and he-elements of g introduce so-called Andreev reflection.
(Cooper-pair tunneling [65] is not included.) T is a full matrix in eh-space and has to be set up as D T (τ, τ ) = 1 2π n dω T n (ω) e −iωτ e ineV τ / e iωτ We skipped the distinction of reference ω J/K as in (25), because in practice only one type, e.g. T LR , has to be calculated, as then T RL follows from complex conjugation. n runs over all integers. A Dyson equation like (32) leads to a recursion which connects T n to T n−2 and T n+2 . It can be solved by reasonably truncating the n-range. Then, just also summing over n, the established current formula (45) can be used.
Expressions like (32) have to sum over all possible combinations of e-and h-indices. For (normalized to θ) [reproduced after [58]]. The normal conducting IV is added for comparison (keep ∆ just for units).

example:
T ee LR = σ ee LR + σ ee LR g ee RR σ ee RL g ee LL T ee LR + σ ee LR g eh RR σ hh RL g he LL T ee LR + σ ee LR g ee RR σ ee RL g eh LL T he LR + σ ee LR g eh RR σ hh RL g hh LL T he LR (50) To point out the Andreev reflection in the model, we look at a term gσgσg with alternating e-and h-indices from the development (21). Such a sequence will also be contained as parts in higher-order terms from (21) or (24). Contracting inner time arguments results in dτ 1 dτ 2 g (h)e RR (τ, τ 1 ) σ ee RL (τ 1 ) g eh LL (τ 1 , τ 2 ) σ hh LR (τ 2 ) g h(e) RR (τ 2 , τ ) because the τ 1 and τ 2 -integration give 4π 2 δ(ω 2 − eV / − ω 1 )δ(ω 1 − eV / − ω). In Fig.9 for holes the energy axis has to be reversed as compared to electrons. Therefore the initial ω is negative. For the argument here, it does not matter if the beginning factor in (51) is g he RR or g ee RR . If it were h, then the Andreev reflection marked by the dashed arrow would follow. The last index of the last g RR could also be h instead of e. Whereas in the normal conducting case in multiple reflections the charge carrying particle always came back to the same energy with regard to the Fermi level when coming back to the same side of the junction (Fig.1), in (51) the frequency argument of the first g RR is shifted by twice the applied voltage equivalent with respect to the second g RR (Fig.9). Such shifts enable transport even when the gaps of the left and right side are still overlapping, that is for voltages below 2∆/e (Fig.10). In the IV -curves steps towards lower voltages at fractions of 2∆ are associated to ever higher-order Andreev reflections. Of course, these current contributions become the more prominent the greater the channel transmission. Besides the step heights and positions the curvature on each step is an important signature. The model agrees excellently with experimental results [53]. The interference of transport processes of different order is correctly taken into account.
Simpler models [66] that add up tranfer probabilities proportional to θ order cannot reproduce these characteristics and are only valid for low transmissions θ << 1.

IV. CONCLUSIONS AND OUTLOOK FOR SECTION III
It has been explained on a quite basic level how quantum transport between two reservoirs in a stationary non-equilibrium state can be modelled. The purpose was to present a Green's functions technique for handling coupling in the context of a field of current research interest, namely transport through point contacts. Although besides the general formalism (also see section I) requiring the development of quite some more subject specific mathematical framework such as the transfer Green's functions, the transport through a contact with arbitrary transmission is a suitable example to illustrate the inclusion of interaction up to all orders in the implicit Dyson equations. Calculations have been carried through in every detail for the normal-conducting single junction. Basic formula and results have been given for the superconducting junction. The decomposition into transport channels (eigenmodes) of a point contact can be inferred from its superconducting transport characteristics which can be taken like a PIN-code [53]. The presented Green's functions formalism has a great potential for extension. Systems of two more or less coherently linked junctions can be modelled [59]. Models for transport through molecules [67,68,69] or atomic chains [70] so