Are current discontinuities in molecular devices experimentally observable?

An ongoing debate in the first-principles description of conduction in molecular devices concerns the correct definition of current in the presence of non-local potentials. If the physical current density ${\bf j}=(-ie\hbar/2m)(\Psi^* \nabla \Psi- \Psi \nabla \Psi^*)$ is not locally conserved but can be re-adjusted by a non-local term, which current should be regarded as real? We prove that the extended Maxwell equations by Aharonov-Bohm give the e.m.\ field generated by such currents without any ambiguity. For an oscillating dipole we show that the radiated electrical field has a longitudinal component proportional to $ \omega \hat{P}$, where $\hat{P}$ is the anomalous moment $\int \hat{I}(\mathbf{x})\mathbf{x} d^3x$ and $\hat{I}$ is the space-dependent part of the anomaly $I=\partial_t \rho+\nabla \cdot \mathbf{j}$. In the case of a stationary current in a molecular device, a failure of local current conservation causes a"missing field"effect that can be experimentally observable, especially if its entity depends on the total current.


I. INTRODUCTION
There are different opinions in the literature concerning the possible existence and physical interpretation of quantum systems for which the electric charge is not locally conserved, i.e., the continuity condition ∂ t ρ + ∇ • j = 0 is not valid everywhere.
Wang et al. [1,2] point out explicitly that in first-principles calculations of the current in molecular devices based on single-particle non-equilibrium Green's function (NEGF) and density-functional theory (DFT) the usual quantum-mechanical current j = (−ie /2m)(Ψ * ∇Ψ − Ψ∇Ψ * ) does not satisfy the continuity condition and must be complemented by a non-local term.The recourse to the NEGF-DFT formalism is necessary in order to take into account the contribution of the inner atomic shells.In principle a full quantum field theory of the system, if it can be formulated in a standard way including all the internal electrons, would yield a conserved bare current; in practice, such a formulation does not exist for realistic systems [3], and even if it could be achieved, it cannot be excluded that it needs to be renormalized, giving rise to quantum anomalies (effective breaking of symmetries and local conservation properties of the bare theory, like for the ABJ anomaly in condensed matter [4,5]).
An appealing feature of the approach by Wang et al. is that the additional current they propose is exactly the same as that originating from the extended Maxwell equations of Aharonov-Bohm without any consideration of the microscopic/quantum aspects.These equations represent essentially the only possible covariant extension of Maxwell's theory which is compatible with the established standard phenomenology of classical electromagnetism and QED and is also applicable to currents that are not locally conserved [6][7][8][9][10][11][12][13][14][15][16][17].
In a series of papers starting in 2108, Jensen, Garner and collaborators have analysed in depth the concept of density current in quantum transport, applying it to specific molecules and obtaining results compatible with the approach by Wang et al.Ref. [18] sets the general theoretical framework and notes that besides legitimate physical reasons for the non-conserving character of local currents, there are technical problems related to the choice of the basis in the first-principles calculations.The discussion is illustrated by simulating elastic and inelastic local currents in a benzenedithiol junction.Simulations show that the local flux does not necessarily follow molecular bonds, with significant part of the flux going "through space".
In Ref. [19], the current density is investigated in saturated chains of alkanes, silanes and germanes.The current density is defined in this context as where G < ij (E) are the matrix elements of the lesser Green function in a standard nonorthogonal basis {ψ i }.The authors show that an enlargement of the eigenbasis for the ab initio calculations does not substantially improve the conservation of current density in this case (while it does lead to an improvement in other cases, notably for graphene ribbons [20]).They then apply the recipe by Wang et al. computing the secondary currents via a Poisson equation.In Refs.[21] and [22] linear carbon wires are considered.Detailed plots of the integrated local current density are given, as compared to the (constant) total current.There are also wave equations in quantum mechanics which do not derive from a microscopic theory but are proposed as effective models with several important applications, in which the current is not locally conserved [23][24][25][26][27][28][29][30][31].It is important in our opinion to develop a formalism that allows to compute the electromagnetic field generated by local currents also in those cases.
In this work, after recalling recent progress in the theory and numerical solutions of the extended Maxwell-Aharonov-Bohm equations, we derive in Sect.II the corresponding wave equations for the electric and magnetic field.These equations allow to compute E and B directly from the physical sources ρ and j, without any reference to the scalar field S which in the traditional formulation has the role to restore local conservation through an additional (or "secondary") charge density proportional to ∂ t ρ and an additional current density proportional to ∇S.The general solution for localized sources is written in the form of retarded integrals.In Sect.II A we compute for the first time the electric and magnetic dipole radiation in far-field approximation.We point out that in general a longitudinal electric radiation field E L can be present; this is one of the main predictions of the extended theory for the case of an oscillating current that is not locally conserved.On the other hand we note that the magnetic field B is simply proportional to the curl ∇ × j, like in Maxwell's theory.Therefore in the limit of stationary currents, in the anomalous case in which ∇ • j is not zero everywhere (presence of charge sinks/sources), B is insensitive to the secondary current ∇S, but can reveal the discontinuities in the current simply because to the "missing links" of current correspond missing contributions to B in the Biot-Savart formula.Such effects are expected to be small, but detectable with accurate experiments, as briefly discussed in Sect.III.Finally, in Sect.IV we present some numerical solutions for the case of stationary currents, obtained not through the direct wave equations derived in this work, but through the double-retarded integrals written in [14,32].In this way it is possible to display explicitly the contributions of the auxiliary field S, confirming that such contributions cancel out and the only consequence for B is the missing field effect.
Our final message can be summarized as follows: microscopic models for the computation of the current density in molecular devices can and should continue to improve their performance and their precision without worrying about local conservation of j.The extended Maxwell equations allow in any case a reliable and efficient computation of the resulting fields.This investigation may lead to interesting discoveries in those cases where the e.m. fields generated by molecular currents are strong enough to play a significant role.

II. GENERAL EQUATIONS FOR E AND B AND THEIR RADIATIVE SOLUTION
In this section we write the general wave equations for the electric and magnetic field in extended electrodynamics, in a form in which the auxiliary field S is completely eliminated, and we find their radiative solution.
The extended Maxwell equations in the Aharonov-Bohm theory are [13,14] ∇ where CGS units have been employed and S is an auxiliary field whose source is the extra- The field S is clearly zero in the pure Maxwell theory, which requires strict local conservation of the current.
It is possible to interpret the extended equations (2), ( 5), (6) as involving some additional or "secondary" sources, namely a secondary charge density proportional to ∂ t S and a secondary current density proportional to ∇S.Including these additional charge and current densities gives total densities which satisfy the continuity equation.Note that the solution of ( 6) for a localized extra-source yields S as a retarded integral: . Therefore S is not localized in the region where the physical sources are present.For this reason we have called the secondary charge and current in our previous work "cloud charge" and "cloud current" and we have evaluated them in some specific cases (see also the numerical simulations in Sect.IV of this work).
The secondary current, in particular, coincides with the additional current predicted by the Landauer-Büttiker theory for systems with quantum transport in which local conservation of the current fails [1,2].
It is possible, however, to write wave equations for E and B in which the field S and the secondary charge and current are completely absent.This is in some sense reassuring, because it implies that the physical fields only depend on the localized, physical sources, and that there is no reason to regard the secondary currents as real and to care, for instance, about their dissipation properties.(See also an alternative proof of the independence of the electric and magnetic fields on the scalar source in Appendix B.) Although the extended Aharonov-Bohm theory has only a limited gauge invariance, it is possible to define potentials φ and A with the usual relations to E and B, namely and to write their wave equations as where once again the role of the secondary charge and current is evident.The solutions of ( 7) and ( 8) are By taking the gradient of ( 7) and adding the time derivative of (8) one obtains an equation for E without S, while taking the curl of ( 8) one obtains an equation for B without S: These equations have also been derived in [16], starting from the extended field equations (2), ( 5), (6) and using some vector calculus identities.
For localized sources the corresponding solutions are We recall that the current density j can be in general written as the sum of an "irrotational" part, which is the gradient of a scalar field and has zero curl, plus a "solenoidal" part, which is the curl of a vector field and has zero divergence: According to eq. ( 12), the irrotational part has no influence on the magnetic field.It follows in particular that interruptions in the current inevitably cause a "missing field" effect, as discussed in Sects.III, IV. (Because the secondary current ∇S which restores current conservation in the extended equations is purely irrotational.)On the other hand, any change in the solenoidal part (like e.g. in [33]) does not affect current conservation but is reflected in a change in the field.
One can further operate on the expressions (11), (12) considering that ∂j ∂t and that where the surface integral is zero because there is no charge on the external surface of the volume of the source, and the change from ∇ to −∇ in the integral in the second line was done because the function on which it operates is of argument x − x .
We have in a similar manner The expressions of the fields thus reduce to These expressions show that the EM fields are obtained from the potentials evaluated without considering their source terms depending on S.

A. Radiative solution in the dipole approximation
For a pure temporal Fourier mode: j (x, t) = j (x) exp (−iωt), ρ (x, t) = ρ (x) exp (−iωt), the general solution for the corresponding scalar potential, without the source depending on S, is of the form φ (x, t) = φ (x) exp (−iωt), with where k = ω/c.For a distant observation point and wavelength large as compared with the source dimension, the expression is approximated by with r = |x|, and n = x/r.
Even if the charge is not conserved locally, it is conserved globally, so that ρ (x ) d 3 x = 0, and thus where p is the Fourier amplitude of the usual charge dipole.
The solution to the corresponding Fourier amplitude of the vector potential is, with the same approximations, Using the identity where in the second line it was used that there is no current leaving or entering the volume of the source, and also allowing for charge non-conservation, so that we have where the dipolar moment of the Fourier amplitude of the extra source I was defined.
The Fourier amplitude of the radiative electric field is thus Analogously, Transforming back to the time domain the EM fields are given by and The first term of E vanishes when it is multiplied by n and therefore represents the transverse component.The longitudinal component is just E • n, i.e.
To fix the ideas, suppose that the moment P of the extra-current is directed along the z-axis.
Then the component E L at a fixed distance r is seen to be maximum on the z-axis and zero in the x-y plane.The opposite happens with the transverse component.
A simple formal example of oscillating extra-current has been introduced in [34].Consider a point-like charge q which oscillates between the positions (0, 0, −a) and (0, 0, a), without a corresponding current: The More realistically we can suppose that if a local violation of charge conservation occurs, this will only concern a small fraction η of the oscillating charge, while the rest of the charge will have a corresponding current and consequently will not generate any longitudinal radiation field.In that case, the ratio between E L and E T will approximately be equal to η.
In [34] we performed a numerical simulation of a source of this kind with a = 10 −7 cm, ω = 2π • 10 9 Hz, θ = π/4 (θ is the 3D polar coordinate), and we obtained that for a totally anomalous source (i.e.η = 1, the entire charge oscillates without a current) the longitudinal field is much larger that the transverse field that would be generated by a corresponding regular source at the same position.The field was computed at a relatively small distance (r varied between 3λ and 13λ), therefore the result cannot be directly compared with eq. ( 23), which holds for larger distances; still it confirms the presence of a longitudinal field, because even in the near-field range it is impossible to obtain such big longitudinal components if the source satisfies local conservation.
In conclusion, if oscillating currents exist that violate even partially the continuity condition, a sizable longitudinal electric radiation should be generated which is obviously incompatible with the standard Maxwell equations.

III. POSSIBLE EXPERIMENTAL SIGNATURES OF A MISSING B AND A RADIATIVE E L
Let us first consider a possible experimental observation of the missing field effect for stationary currents.As discussed in [32] and confirmed in Sects.II and III of this work, if in a linear conductor there are regions in which a violation of continuity occurs and a fraction η of the current does not flow as "physical current ρv" but as "secondary current ∇S", then the magnetic field generated by the secondary current is zero, as a consequence of eq.(10).
To detect the missing field, we proposed in [32] to measure at a fixed distance r the field of a normal conductor carrying a current i, and then the field of the anomalous conductor carrying the same current.A differential measurement with three wires was devised, in order to minimize errors.The scheme was especially suited for brief current pulses and for the case when the supposed anomalous conductor is a superconductor, which can be driven in/out a normal state by changing its temperature.
Another possible technique, which employs stationary currents and one single wire, is based on the measurement of the ratio B/i at a fixed distance in dependence on i (supposing that a variable external bias allows to change i).This would work if the ratio η depends on i, in which case B/i is expected to change, in open violation of the Maxwell equations.
In fact, the current patterns and local discontinuities observed in simulations of molecular devices depend in general on the total current.
Possible errors could originate from slight changes in the spatial distribution of the current at different i, if the conductor has a radius that cannot be disregarded compared to the distance r.This radius depends on how "elementary" the conductor is (for example, a single carbon nanowire vs. a bundle of nanowires); in turn, that depends on how much current is needed for the measurement, and thus indirectly on r (because of the size and sensitivity of the detector).
For example, suppose that a single carbon nanowire can carry a current of 10 −10 A and the detector can be placed at a distance of 10 −4 m, very large compared to the radius of the wire.The field would then be of the order of 10 −12 T, i.e. accessible to a sensitive SQUID.
If the ratio η is of the order of 10 −2 and the error on the current is negligible, the SQUID is required to detect a variation in the field of 10 −14 T.
Turning now to electric fields, the possible generation and detection of a longitudinal e.m. radiation has rarely been explored in the past decades [35][36][37][38][39][40].A recent preliminary experiment and its relation to the extended electromagnetic theory has been described in [16], including a discussion of error sources.A distinctive feature of longitudinal electric radiation would be its ability to penetrate thin layers of good conductors much easier than a transverse radiation.Even in favourable cases, however, a predominant transverse component would be present, causing interference and noise.If future developments can lead to a clean selective detection of the longitudinal component, the possible technological applications would be manifold.
There is also much to do, of course, concerning the design of efficient antennas.First, one would need to identify materials with local non-conservation that can support a sufficiently large current at high frequency.Then the most appropriate antenna geometry should be studied.For this purpose, the general equations developed in this paper constitute a firm starting point.

IV. NUMERICAL SOLUTIONS OF THE EXTENDED EQUATIONS WITH THE POTENTIALS
Consider again the equations for the potentials, written in the form where I(t, x) is an assigned extra-current.By solving eqs.( 25) and ( 26) we obtain for φ and A expressions which contain "double-retarded integrals" of the extra-current I.In [32] we gave some numerical computations of these integrals for a slowly varying source, as recalled below.In this section we consider the stationary case and the resulting expressions are much simpler.For the numerical integration we employ a 6D Monte Carlo, which can be quite time consuming but is straightforward and does not require any analytical approximations.
The idea is to take, instead of a source with a slow temporal variation, a truly stationary source.Then the retarded integrals become simple space integrals.We recall that the typical time scale of the source in the work [32] was τ 10 −5 s, chosen because (1) it corresponds to the proposed experimental conditions, being the characteristic discharge time of the circuit; (2) it allows to disregard certain phenomena occurring only at high frequency, like the intervention of stray capacitance and temporary charge accumulation.
In [32] we considered a case of local non-conservation involving a point-like sink, where current partially disappears, and a point-like source where the current reappears.More precisely, the sizes of sink and source were given by a regulator ε of the order of 10 −7 cm.The regulator was applied after the first 3D integration (analytical) in d 3 z.Then the second numerical integration was made (in d 3 y, over the extended secondary cloud current decreasing like 1/|y| 3 ), using the command NIntegrate of Mathematica and checking the results by comparison with a Monte Carlo integration in 3D.
The disadvantage of that procedure is the complication at the formal-algebraic level.
The expressions obtained after the first integration are quite bulky.At the same time, the assumption of point-like sources is limiting, because when the failure of local conservation occurs in a quantum wavefunction, the regions where ∂ t ρ + ∇ • j = 0 are not pointlike but extended, with a shape more similar to that of a couple of disks.
Therefore we model here the extra-current source with a double Gaussian (Fig. 1), and the method actually applies to any geometrical shape: The auxiliary anomalous vector potential written as double integral of the extra-current is When we compute the magnetic field at the position x = (r, 0, 0), its only non-zero component is B 2 : The derivative of the first term is The other derivatives are similar, and in total we obtain, also considering that x 3 = 0 in this particular configuration: This expression is simpler and more direct than those obtained in [32] and suitable for a 6D Monte Carlo algorithm which first generates random values of z in a quite narrow range, corresponding to the support of the Gaussian (28), i.e., inside the primary source.The main computational difficulty is (as before, actually) that for y it is necessary to sample a much wider region.Typically for each value of z we need to generate tens of values of y in the sampling.
While we sample the integrand in order to evaluate the total integral, we also make a map of the function resulting from the partial integration in z.For this purpose we divide the 3D cube of integration in y (with variable side 2R y ) into 100 cells along each direction.
The contribution to the magnetic field from each cell can be interpreted as being generated by a secondary current, or "cloud" current, proportional to the gradient of the field S. As seen from Tables I and II and from Figures 2-5, the regions with the largest secondary current density are those between the current source and sink and close to them, but the large cloud that extends around them and whose density decreases slowly with distance also gives a relevant contribution to the total integral; in the end, the latter contribution exactly cancels that of the localized sources.In order to evaluate the far contributions we sample the integrand over cubes of y of increasing size, excluding each time a cubic core equal to the previous cube.This is necessary because the number of sampling points cannot be increased beyond a practical limit of the order or 10 10 , and if we would sample the cores together with the periphery, the cores would produce too much noise.All the physical parameters have been chosen of the same magnitude order as in [32].In summary, we have shown on general grounds that for any charge density ρ and current density j given by a microscopic model, the generated e.m. field is given by eqs.( 9) and (10), no matter if the continuity equation is satisfied or not.
In the case of oscillating charges and currents, the radiated longitudinal electric field is given by eq. ( 17) or (22).For example, if a fraction η of a charge q oscillating along the z-axis range R y (units 10  The field is computed at a distance r = 10 −4 cm from the origin, on the x 1 axis.The shape of the current sink and source is spherical (ε = D = 0.5 • 10 −7 cm).The integration range in z 1 and z 2 is R = 2 • 10 −7 cm, in z 3 is R 3 = 5 • 10 −7 cm.In the y range [800 − 1100] • 10 −7 cm, which contains the point where B s is computed, the integration is further split into two parts in order to reduce the noise due to the factor 1/|x − y|.As seen from the last row of the table, the total anomalous field is zero within errors ("missing field" effect).
range R y (units 10  I, but for a sink and source with the shape of a disk/ellipsoid (ε = 0.5 • 10 −7 cm, D = 2.5 • 10 −7 cm).The integration range in z 1 and z 2 changes accordingly: here R = 10 • 10 −7 cm.The total anomalous field is again zero within errors; the single contributions differ significantly from the case of spherical sink and source only at small R y range.The data for a wider disk (D = 10 • 10 −7 cm, not shown here) reveal a similar behavior: in the first two y ranges we have respectively contributions 0.185 and 0.143, with no substantial differences in the other ranges.In the 6D Monte Carlo integration, the 3D region of the y variable has been subdivided into 100 3 cells; for each cell, the contribution of sampling points falling into it is displayed.The numerical values on the color scale must be normalized to compute the field, and are therefore not meaningful as absolute value, but their sign and relative values are of interest.In this figure we have y 3 = 0, i.e., we are observing the cloud on the plane equidistant from source and sink.All contributions are positive, since we are between the sources (compare Fig. 3).(b) Here we have y 1 = 0, therefore we are observing the cloud in the plane that cuts sink and source.Note that sink and source both give positive contributions on the inner side, and negative contributions on the outer side.From Tab.I we deduce that the total contribution of the region shown in this figure is positive.
The equation for the magnetic field (10) indicates that it is possible in principle to observe experimentally if a fraction η of a stationary current flows in a conductor with microscopic spatial "interruptions" affecting a volume that is a fraction χ of the total volume: in that case, the average Biot-Savart magnetic field generated will be (1 − ηχ)B 0 , where B 0 is the  field generated by the same total current flowing in a conductor where ∇ • j = 0 everywhere.
The microscopic models cited in this work allow in principle to estimate the fractions η and χ.On the experimental side, a possible strategy to prove some anomalies (though not to disprove them) consists of looking for changes in the ratio B/i in dependence on i.We have briefly discussed the possible uncertainties associated with this measurement.
Finally we note that the control of magnetic field patterns in molecular devices may be interesting for applications to high-density memory storage consisting of single molecular magnets (see [44] and refs.).It has been further suggested that NMR-type experiments can be performed in order to observe the spatial fluctuations of magnetic fields generated by DC current flow [45].
Appendix A: Equations in SI units In SI units the equations for the potentials and the auxiliary field S read with ε 0 and µ 0 the permittivity and permeability of free space, respectively.
The EM fields are expressed as For a localized extra-source I, the solution for S is S (x, t) = µ 0 4π The wave equations for the electric and magnetic fields are: For localized sources the corresponding solutions are Appendix B: Alternative proof of the independence of the EM fields on the scalar source.
The contribution to the potentials given by the scalar source is φ (x, t) = − 1 4πc ∂S (x , t ) ∂t where it is assumed that at infinity no scalar field exists (no incoming scalar field from large distances in the remote past, and local additional sources acting during finite time intervals).Since the scalar field is assumed to originate from a localized source that was not acting in the infinite past, the surface integrals at infinity are zero for finite (x, t).This is so because the contributions from |x | → ∞ to the fields at finite values of (x, t), come from the scalar at t → −∞, which is zero by hypothesis.
moment P in this case has the only non-zero component | Pz | = 2ωqa.If we assign the physical parameters q, a and ω it is straightforward to compute the longitudinal far field and we have along the z-axis | ÊL | = 2ω 2 qa/(c 2 r) (in SI units: | ÊL | = µ 0 ω 2 qa/(2πr)).The transverse far field is vanishing.

FIG. 1 :
FIG. 1: Parameters and location of the Gaussian extra-sources employed in the numerical simulations.(Orange: current sink.Blue: current source.)

FIG. 2 :
FIG. 2: (a) Contributions of the "cloud" of secondary current to the anomalous field B s for spherical sink and source (see Tab. I), in the y range [0 − 10] • 10 −7 cm.In the 6D Monte Carlo integration, the 3D region of the y variable has been subdivided into 100 3 cells; for each cell, the contribution of sampling points falling into it is displayed.The numerical values on the color scale must be normalized to compute the field, and are therefore not meaningful as absolute value, but their sign and relative values are of interest.In this figure we have y 3 = 0, i.e., we are observing the cloud on the plane equidistant from source and sink.All contributions are positive, since we are between the sources (compare Fig.3).(b) Here we have y 1 = 0, therefore we are observing the cloud in the plane that cuts sink and source.Note that sink and source both give positive contributions on the inner side, and negative contributions on the outer side.From Tab.I we deduce that the total contribution of the region shown in this figure is positive.

FIG. 3 :
FIG. 3: Same sources and cutting plane as in Fig. 2, (b), but in the y range [10 − 50] • 10 −7 cm.(The central region [0 − 10] is not sampled because it would give a strong noise, compared to the larger region.)Note that the pattern of Fig. 2, (b) is confirmed concerning the inner/outer regions with positive/negative contributions respectively.

FIG. 4 :
FIG. 4: This figure can be compared with Fig. 2, (a), with the only difference that sink and source are here Gaussian disks/ellypsoids of diameter D = 10 • 10 −7 cm instead of D = 0.5 • 10 −7 cm.(D corresponds to the σ of the Gaussian density distribution.)As in Fig. 2, (a), the contributions shown lie on the plane between source and sink, with y 3 = 0.

TABLE I :
Contributions of the secondary current ∇S, in various regions of the y-space, to the ratio B s /B 0 .B s is the anomalous field generated by the secondary current.B 0 is the Biot-Savart field that the same total current would generate if flowing as a primary local current through the gap of length 2a.

TABLE II :
Same as in Tab.