text:
# Generalized boundary conditions for the circuit theory of spin transport

## I INTRODUCTION

Spin-transport in hybrid ferromagnet (F)-normal metal (N) systems has been object of extensive investigation since the discovery of the giant magnetoresistance (GMR) effect. Magneto-electronic multiterminal systems lead to novel applications, e. g. magneto recording heads or magnetic-field based sensor devices. Future prospectives of non-volatile electronics motivate also many fundamental and applied research. Transistor-like effects were found in a ferromagnet-normal metal three-terminal device, that depend on the relative orientation of the magnetization of the ferromagnets. The dependence of transport properties with the relative angle between the magnetization directions of the ferromagnets, was addressed experimentally  as well as theoretically in the past. These days non-collinear spin transport attracts an increased attention due to the recent interest on the spin-current induced magnetization torques.  A ferromagnetic single-electron transistor in a three terminal configuration has been realized and studied theoretically.

If a normal metal (N) is attached to a superconductor (S), quasiparticles of different spins are coupled via Andreev reflection on the normal side of the NS interface. A strongly modified density of states (DOS) in the normal metal caused by induced superconducting correlations was found.  This is so-called superconducting proximity effect. In contrast with normal metals, in a ferromagnet (F) the presence of a strong exchange field leads to big differences between the two spin bands. The question of the coexistence of ferromagnetism and superconductivity has been extensively investigated in the last decades. The effect of spin splitting on the superconducting proximity effect was already investigated long time ago. New experimental developments on proximity effect in ballistic ferromagnetic layers have been recently reported  and theoretically confirmed. Results obtained in ferromagnet-superconductor nanocontacts have been explained in terms of bias dependent transparency of the FS interface. Diffusive heterostructures of ferromagnets and superconductors showed unexpected large values of the conductance of the ferromagnetic part. Strong mutual influence between the superconductor and the ferromagnet and long range proximity effects have been proposed to explain these results. Effects related to the interplay between spin accumulation and Andreev reflection has been also investigated. On the other hand, the prediction of an exotic superconducting state formed in a superconductor by the presence of an exchange field $`h`$, was independently reported by Larkin-Ovchinnikov and Fulde-Ferrel, some years ago. Experimental confirmation of these effect in bulk superconductors is still needed. Many investigations have focused on the study of thermodynamic properties of FS multilayers.

Non-equilibrium Keldysh Green’s functions in the quasiclassical approximation, have been extensively used in the past to describe non-equilibrium superconductivity. This quasiclassical theory of superconductivity is based on semiclassical transport equations for quasiparticles. Generally to solve these transport equations is technically difficult, however the final results are relative simple and clear. The so-called “circuit theory of Andreev reflection”, was conceived as a generalization of the Kirchhoff’s theory of electronic circuits, to simplify these equations into a handful of accessible rules. The circuit theory of mesoscopic transport was recently extended to describe transport in non-collinear magnetic structures . Interesting phenomena like spin precession effects on the induced spin accumulation have been found . An important concept in the theory of non-collinear spin transport is the so-called mixing conductance $`G^{}`$  which is closely related to the spin-current induced magnetization torques  and the Gilbert damping in thin ferromagnetic films. Moreover a three terminal spin-transistor device has been proposed in the framework of the circuit theory.  The main advantage of the circuit theory description is that it provides a simple approach based on spin-charge current conservation to calculate the transport properties of (multi-terminal) mesoscopic hybrid systems. This is achieved by mapping the concrete geometry onto a topologically equivalent circuit, represented by finite elements. We note that this very same description allows to calculate a variety of transport properties, such as current statistics, weak localization corrections or transmission eigenvalues. Consequently, it is worthwhile to extend this general framework to the all possible combinations of heterostructures.

So far, the developed circuit theory is suitable to describe electron and spin transport in multiterminal hybrid ferromagnet-normal metal FN or superconductor-normal metal SN systems. Obviously, an extension of the circuit theory to hybrid systems combining ferromagnets, normal metals and superconductors is necessary. This is done in the present paper. We derive a general expression for an arbitrary contact, which, however, turns out to be a little bit unhandy. To obtain manageable expressions we present also approximate results for two special cases. Below we list all our results:

* The general matrix current Eq.(23) through an arbitrary spin-dependent connector. It requires the knowledge of the full scattering matrix (or transmission matrix) as consequence of the non-separable Spin- and Nambu structures. This inhibits the transformation into the normal eigenmodes of the scattering problem. Eq. (23) is therefore mostly of numerical interest.
* The matrix current Eq.(30) for a weakly spin-dependent contact. Here an expansion in terms of normal eigenmodes is possible. Eq. (30) is a spin-dependent correction to the matrix current from Ref. ().
* If the spin-dependent contact is a tunnel barrier, the tunneling matrix current Eq. (34) takes a particular simple and transparent expression. The properties of the contact can be expressed in terms of the spin conductances $`G_{()}`$ and the (complex) mixing conductance $`G^{}`$.

The paper is organized as follows. In section II we illustrate some basic aspects of the circuit theory. In Section III we present the microscopic description of a contact region or connector in the circuit theory. In section IV we calculate the matrix current through a general connector between two metallic regions (nodes) of the circuit. The nodes can be of N-, F- or S-type and the contact is assumed to have a arbitrary magnetic structure (e. g. magnetic tunnels junctions, magnetic interfaces). This expression constitutes the generalized matrix current for the circuit theory. In section V we make use of a perturbation expansion in the spin structure to obtain simplified expressions for two cases. Details of the pertubation expansion are presented in the Appendix. In Section VI we present our conclusions.

## II CIRCUIT THEORY

In the circuit theory the system is split up into reservoirs (voltage sources), connectors (contacts,interfaces) and nodes (low resistance islands/wires) in analogy to classical electric circuits. Both reservoirs and nodes are characterized by $`16N_{ch}\times 16N_{ch}`$ Green’s functions $`\stackrel{ˇ}{G}`$, which are matrices in Keldysh$``$Nambu$``$Spin$``$Channels space, where $`N_{ch}`$ are the number of propagating modes. These Green’s functions play the role of generalized potentials of the circuit theory. The Green’s functions in the reservoirs and nodes are assumed to be isotropic in momentum space. This requires that sufficient elastic scattering is present both in reservoirs and nodes, due to the presence of random scatterers and irregularities in the shape. This justifies the use of the diffusion approximation to describe transport. This assumption is reasonable since hybrid (multi-terminal) devices are quite dirty systems. In general for the stationary case, $`\stackrel{ˇ}{G}`$ depends on space coordinates and energy

$$\stackrel{ˇ}{G}(\stackrel{}{r},\stackrel{}{r}^{},\epsilon )=𝑑t\text{ }\stackrel{ˇ}{G}(\stackrel{}{r},\stackrel{}{r}^{},tt^{})\text{ }\mathrm{exp}\{i\frac{\epsilon }{\mathrm{}}(tt^{})\}$$
(2)
being

$$\stackrel{ˇ}{G}(\stackrel{}{r},\stackrel{}{r}^{},tt^{})=\left[\begin{array}{cc}\widehat{G}^R(\stackrel{}{r},\stackrel{}{r}^{},tt^{})& \widehat{G}^K(\stackrel{}{r},\stackrel{}{r}^{},tt^{})\\ 0& \widehat{G}^A(\stackrel{}{r},\stackrel{}{r}^{},tt^{})\end{array}\right].$$
(3)
$$\widehat{G}^R(\stackrel{}{r},\stackrel{}{r}^{},tt^{})=i\theta (tt^{})\{\widehat{\mathrm{\Psi }}(\stackrel{}{r}),\widehat{\mathrm{\Psi }}^{}(\stackrel{}{r}^{})\}$$
(4)
$$\widehat{G}^A(\stackrel{}{r},\stackrel{}{r}^{},tt^{})=i\theta (t^{}t)\{\widehat{\mathrm{\Psi }}(\stackrel{}{r}),\widehat{\mathrm{\Psi }}^{}(\stackrel{}{r}^{})\}$$
(5)
$$\widehat{G}^K(\stackrel{}{r},\stackrel{}{r}^{},tt^{})=i[\widehat{\mathrm{\Psi }}(\stackrel{}{r}),\widehat{\mathrm{\Psi }}^{}(\stackrel{}{r}^{})].$$
(6)
where $`[..,..]`$ and $`\{..,..\}`$ denotes commutator and anticommutator respectively. Semiclassical and diffusive approximations allow to obtain Green’s function which depends only in one spatial coordinate $`\stackrel{ˇ}{G}(\stackrel{}{r},\epsilon )`$. On the nodes the Green’s functions are also assumed to be spatially homogeneous, depending only on energy $`\stackrel{ˇ}{G}(\epsilon )`$. This requires that the resistance of the node is much smaller than the contacts resistances connecting different nodes, which implies also that the current through the device is controlled by the contacts resistances. Regions with spatially dependent Green’s functions (e.g. diffusive wires) are modelled by an appropriate discretization. For example, a quasi one-dimensional diffusive wire can be represented by a series of tunnel junctions and nodes. Internal dynamics along the wire (finite energy transport, spin-flip, etc.) is included in this description as a leakage current from the node.

In analogy to Ohm’s law in classical electric circuits ($`I=V/R`$), in the present circuit theory it is essential to obtain the “spin-charge” matrix current that flows between two nodes (or between a reservoir and a node) through a contact/connector. This matrix current depends on the properties of the connector (analog to the resistance $`R`$) and on the Green’s functions at both side of the connector (analog to the voltage drop through the connector $`V`$ ). In general the matrix current $`\stackrel{ˇ}{I}(z,\epsilon )`$ is defined in Keldysh$``$Nambu$``$Spin$``$Channels space as

$$\stackrel{ˇ}{I}(z,\epsilon )=\frac{e^2\mathrm{}}{m}𝑑\rho \left(\frac{}{z}\frac{}{z^{}}\right)\stackrel{ˇ}{G}(\stackrel{}{r},\stackrel{}{r}^{},\epsilon )|_{\stackrel{}{r}\text{ }=\text{ }\stackrel{}{r}^{}}$$
(7)
where $`\stackrel{}{r}(\stackrel{}{\rho },z)`$, $`z`$ being along the transport direction and $`\stackrel{}{\rho }`$ being perpendicular to the transport direction. The electric current $`I_e`$ is defined in terms of the matrix current $`\stackrel{ˇ}{I}(z,\epsilon )`$ as

$`I_e(z)={\displaystyle \frac{1}{\text{ }4\text{ }e}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\epsilon }{2\pi \mathrm{}}}\text{Tr}\left\{\widehat{\sigma }_z\widehat{\tau }_3\text{ }\stackrel{ˇ}{I}^\text{K}(z,\epsilon )\right\},`$

where $`\widehat{\stackrel{}{\sigma }}`$, $`\widehat{\tau }`$ are Pauli matrices in spin and Nambu space. $`\stackrel{ˇ}{G}(\stackrel{}{r},\stackrel{}{r}^{};\epsilon )`$ can be expanded in transverse modes in the following way

$`\stackrel{ˇ}{G}(\stackrel{}{r},\stackrel{}{r}^{};\epsilon )`$ $`=`$ $`{\displaystyle \underset{nm,\sigma ,\sigma ^{}}{}}\stackrel{ˇ}{G}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})\times `$ (9)
$`\mathrm{exp}(i\sigma k_n^szi\sigma ^{}k_m^s^{}z^{})\chi _n^s(\rho )\chi _m^s^{}(\rho ^{}).`$
Here $`n`$ ($`m`$) label the propagating modes, $`\sigma ,`$ $`\sigma ^{}=\pm 1`$ are the direction of propagation, $`s,s`$ $`\{,\}`$ are spin indices, $`k_n^s`$ is the longitudinal momentum and $`\chi _n^s(\rho )`$ is the transverse wave function. By using this representation, the current can be written like

$`\stackrel{ˇ}{I}^{s,s^{}}(z,\epsilon )`$ $`=`$ $`i\text{ }e^2{\displaystyle \underset{n\text{ };\text{ }\sigma ,\sigma ^{}}{}}(\sigma \upsilon _n^s+\sigma ^{}\upsilon _m^s^{})\stackrel{ˇ}{G}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,\epsilon )`$ (11)
$`\times {\displaystyle }d\rho \text{ }\chi _n^s(\rho ,z)\chi _m^s^{}(\rho ,z).`$
Here $`\upsilon _{n(m)}^{s(s^{})}=\mathrm{}`$ $`k_{n(m)}^{s(s^{})}/m`$ is the velocity in mode $`n(m)`$ with spin $`s(s^{})`$. The transverse wave functions $`\chi _n^s(\rho ,z)`$ are eigenfunctions of the Spin operator $`S=\mathrm{}`$ $`\widehat{\sigma }_z/2`$ and are normalized in a way that the carry unity flux current. In this case the current reduces to

$$\stackrel{ˇ}{I}^{s,s}(z,\epsilon )=2\text{ }i\text{ }e^2\underset{n\text{ };\text{ }\sigma }{}\sigma \text{ }\upsilon _n^s\text{ }\stackrel{ˇ}{G}_{n\sigma ,n\sigma }^{s,s}(z,\epsilon )$$
(12)
Note that the current $`\stackrel{ˇ}{I}^{s,s}(z,\epsilon )`$ in Eq.(12) is a matrix diagonal in spin space. In the case of a ferromagnetic material Eq.(12) is valid if the magnetization of the ferromagnetic material is parallel to the spin quantization axis. In the case of non-collinear transport where there are two different magnetizations defined in the system $`\stackrel{}{M}_1`$ and $`\stackrel{}{M}_2`$, Eq. (12) has to be properly transformed into the basis of the spin quantization axis

$`\stackrel{ˇ}{I}^{s,s^{}}=U\text{ }\stackrel{ˇ}{I}^{s,s}U^1`$

$`U`$ being the spin rotation matrix that transforms $`\stackrel{}{M}_{1(2)}`$ into the spin quantization axis. Note that in this case the transformed matrix $`\stackrel{ˇ}{I}^{s,s^{}}`$can have general spin structure.

## III ARBITRARY CONNECTOR

The contacts or connectors  between the reservoirs and nodes are described in terms of the scattering formalism by a transfer matrix $`\overline{M}`$. The goal is to express the matrix current given by Eq.(11) in terms of isotropic quasiclassical Green’s functions $`\stackrel{ˇ}{G}_{1(2)}`$ at both sides of the contact and transmission/reflection coefficients that characterize the scattering in the contact region. $`\overline{M}`$ is in general a $`16N_{ch}\times 16N_{ch}`$ matrix in Keldysh$``$Nambu$``$Spin$``$Channels space. Particularly $`\overline{M}`$ is proportional to the unit matrix in Keldysh space and diagonal in Nambu space. This is denoted by $`\left(\text{ }\overline{}\text{ }\right)`$.

The Green’s functions $`\stackrel{ˇ}{G}\stackrel{ˇ}{G}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})`$ is not a continuous function at $`z=z`$. The values of $`\stackrel{ˇ}{G}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})`$ for $`z>z^{}`$ and $`z<z`$ are matched at $`z=z`$ by using the following representation (for details see Ref.)

$`2i\stackrel{ˇ}{G}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})`$ $`=`$ $`{\displaystyle \frac{\stackrel{ˇ}{g}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})}{\sqrt{v_nv_m}}}`$ (14)
$`+{\displaystyle \frac{\sigma \delta _{\sigma \sigma ^{}}\text{ }\delta _{nm}\text{ sign}(zz^{})}{v_n}}\stackrel{ˇ}{1}.`$
where the matrix $`\stackrel{~}{g}(z,z^{})\stackrel{ˇ}{g}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})`$ is continuous at $`z=z`$. Note that the matrix $`\stackrel{ˇ}{g}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}`$ has non-trivial structure in channels space. We call this function ballistic Green’s function. At the contact region $`(z=0)`$ the transfer matrix $`\overline{M}\overline{M}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}`$ connects $`\stackrel{~}{g}(z,z^{})`$ at both sides (left $`z,z=0^{}`$ and right $`z,z=0^+`$)

$$\stackrel{~}{g}_{(2)}=\overline{M}\text{ }\stackrel{~}{g}_{(1)\text{ }}\overline{M}^\text{ }$$
(15)
where $`\stackrel{~}{g}_{(2)}\stackrel{~}{g}(z=0^+,z^{}=0^+)`$ and $`\stackrel{~}{g}_{(1)}\stackrel{~}{g}(z=0^{},z^{}=0^{})`$. Note that in Eq.(15) there is implicit a sumation over Channel and Spin indices.

An important assumption of this theory is the isotropization assumption. The ballistic Green’s function $`\stackrel{~}{g}`$ defined at both sides of the contact region becomes isotropic in momentum space, when departing from the contact region. Such isotropization happens in a region of the order of the elastic mean free path $`l_{\text{mfp}}`$ but much smaller than the spin diffusion length $`l_{\text{sf}}=\sqrt{D\text{ }\tau _{\text{sf}}}`$ and the coherence length of a superconductor $`\xi _\mathrm{\Delta }=\sqrt{D/\mathrm{\Delta }}`$. In the isotropization zone the dominant contribution to the self-energy is the elastic scattering, which implies sufficient disorder (or chaotic scattering) in this region. The self-energy is then $`\stackrel{ˇ}{\mathrm{\Sigma }}_{imp}^{1(2)}=i`$ $`\stackrel{ˇ}{G}_{1(2)}/2`$ $`\tau _{imp}`$, $`\stackrel{ˇ}{G}_{1(2)}`$ being the isotropic quasiclassical Green’s function for the left $`(1)`$ and for the right $`(2)`$ side of the contact region, which is proportional to the unit matrix in Channels space and $`\tau _{imp}`$ is the impurity scattering time. In order to assure that the $`\stackrel{ˇ}{G}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z^{})`$ does not diverge at $`z(z^{})0`$, we have to impose the following conditions

$$\left(\overline{\mathrm{\Sigma }}^z+\stackrel{ˇ}{G}_1\right)\left(\overline{\mathrm{\Sigma }}^z\stackrel{~}{g}_1\right)=0$$
(16)
$$\left(\overline{\mathrm{\Sigma }}^z+\stackrel{~}{g}_1\right)\left(\overline{\mathrm{\Sigma }}^z\stackrel{ˇ}{G}_1\right)=0$$
(17)
$$\left(\overline{\mathrm{\Sigma }}^z\stackrel{ˇ}{G}_2\right)\left(\overline{\mathrm{\Sigma }}^z+\stackrel{~}{g}_2\right)=0$$
(18)
$$\left(\overline{\mathrm{\Sigma }}^z\stackrel{~}{g}_2\right)\left(\overline{\mathrm{\Sigma }}^z+\stackrel{ˇ}{G}_2\right)=0,$$
(19)
$`\overline{\mathrm{\Sigma }}^z=`$$`\sigma `$ $`\delta _{\sigma \sigma ^{}}`$ being the z-Pauli matrix in the direction of propagation “sub-space” $`\left\{\sigma \sigma ^{}\right\}`$ of the Channels space. An important consequence of this approach is that after the isotropization zone, the ballistic Green’s function $`\stackrel{ˇ}{g}_{n\sigma ,m\sigma ^{}}^{s,\text{ }s^{}}(z,z)`$ equals the isotropic quasiclassical Green’s function $`\stackrel{ˇ}{G}_{1(2)}`$.

## IV GENERALIZED MATRIX CURRENT

In previous work, $`\overline{M}`$ and $`\stackrel{ˇ}{G}_{1(2)}`$ commute because no spin structure was considered. Now due to their general spin structure, $`\overline{M}`$ and $`\stackrel{ˇ}{G}`$ do not commute. We multiply Eq. (16) and Eq. (18) by $`\overline{M}`$ from the left and by $`\overline{M}^{}`$ from the right. Adding both resulting equations, using Eq.(15) and using the normalization condition $`\stackrel{ˇ}{G}^2=\stackrel{ˇ}{1}`$,  we get the following expression for $`\stackrel{~}{g}_1`$

$$\stackrel{~}{g}_1=\left(\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M}\right)^1\left[2\stackrel{ˇ}{G}_1+\left(\stackrel{ˇ}{1}\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M}\right)\overline{\mathrm{\Sigma }}^z\right]\text{.}$$
(20)
Note that in Eq.(20) all the structure in channels is contained in $`\overline{M}`$. From Eq.(12), the matrix current $`\stackrel{ˇ}{I}`$ is be expressed in terms of $`\stackrel{~}{g}_{1(2)}`$ as

$$\stackrel{ˇ}{I}(\epsilon )=2\text{ }i\text{ }e^2\underset{n\text{ };\text{ }\sigma }{}\sigma \text{ }\upsilon _n\text{ }\stackrel{ˇ}{G}_{n\sigma ,n\sigma ^{}}^{s,\text{ }s^{}}(\epsilon )=e^2\text{Tr}_{n,\sigma }\left[\overline{\mathrm{\Sigma }}^z\stackrel{~}{g}_1\right].$$
(21)
Note that the Green’s function $`\stackrel{ˇ}{G}_{n\sigma ,n\sigma ^{}}^{s,\text{ }s^{}}`$ is assumed to be spatially homogeneous, depending only on energy $`\stackrel{ˇ}{G}(\epsilon )`$. In this case that current does also not depend on position. By using the cyclic property of the trace, we find finally for the matrix current $`\stackrel{ˇ}{I}(\epsilon )`$

$`\stackrel{ˇ}{I}(\epsilon )`$ $`=`$ $`e^2\text{Tr}_{n,\sigma }\{(\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M})^1(\stackrel{ˇ}{1}\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M})`$ (23)
$`+2(\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M})^1\text{ }\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1\}.`$
Eq.(23) is the most general expression for the current $`\stackrel{ˇ}{I}`$ in terms of isotropic quasiclassical Green’s functions $`\stackrel{ˇ}{G}_{1(2)}`$ at both sides of the contact and transmission/reflection coefficients of the contact region. Once such expression of $`\stackrel{ˇ}{I}(\epsilon )`$ is obtained, we can apply the generalized Kirchhoff’s rules, imposing that the sum of all matrix currents into a node is zero. This completely determines all properties of our circuit. Eq.(23) therefore completes our task to find the generalized boundary condition for the circuit theory. However, in this form a concrete implementation requires the knowledge of the full transfer matrix (or equivalently of the scattering matrix). Usually this information is not available for realistic interfaces and one tries to reduce Eq. (23) to simple expressions by some reasonable assumptions. In the next section we will do this for two special cases. Note, that for a spin-independent interface, Eq. (23) can be expressed by the transmission eigenvalues only, which is a formidable simplification. This transformation is not possible anymore for the spin-dependent contact.

## V TWO SPECIAL CASES

We want to obtain more transparent and clear analytical expressions of Eq.(23). The price to pay for that is loss of generality, since we have to make certain assumptions about the contact. The main difficulty lies in the inversion of the matrix $`\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\stackrel{ˇ}{G}_2\overline{M}`$ in Channel-Spin space. Let us assume that the transfer matrix $`\overline{M}`$ can be split in the following form

$$\overline{M}=\overline{M}_0+\delta \overline{M}=\overline{M}_0(\stackrel{ˇ}{1}+\delta \overline{X}),$$
(24)
where $`\overline{M}_0`$ is a transfer matrix with structure in channel space only and proportional to the unit matrix in spin space, while $`\delta \overline{M}\overline{M}_0\delta \overline{X}`$ includes non-trivial structure in spin space. Note that the matrices $`\overline{M}_0`$, $`\stackrel{ˇ}{G}_{1(2)}`$ commute with each other, whereas in general $`\delta \overline{X}`$ does not commute with $`\overline{M}_0`$ and $`\stackrel{ˇ}{G}_{1(2)}`$. Assuming that $`\delta \overline{X}1`$, we can perform a perturbation expansion of $`\left(\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M}\right)^1`$ in the parameter $`\delta \overline{X}`$ (see Appendix). The matrix $`\overline{M}\left(\overline{M}\right)_{n\sigma ,m\sigma ^{}}^{s;s^{}\text{ }}`$ can now be diagonalized in the basis of eigenmodes $`N(N):`$ $`\left(\overline{M}\right)_{n\sigma ,m\sigma ^{}}^{s;s^{}\text{ }}\left(\overline{M}\right)_{N\sigma ,N\sigma ^{}}^{s;s^{}\text{ }}`$.

To 0th order $`\overline{M}\overline{M}_0`$ commutes with $`\stackrel{ˇ}{G}_{1(2)}`$. In this case Eq. (23) reduces to

$$\stackrel{ˇ}{I}^{(0)}(\epsilon )=e^2\text{Tr}_{N,\sigma }\left[\frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2}\left(2\text{ }\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1+\stackrel{ˇ}{1}\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\right)\right].$$
(25)
The hermitian conjugate matrix $`\overline{Q}_0`$ ($`\overline{Q}_0^{}\overline{Q}_0`$), reads in Channel space

$$\overline{Q}_0=\left[\begin{array}{cc}A& B\\ B^{}& A\end{array}\right]$$
(26)
being $`A`$ real and $`B`$ complex $`N_{ch}\times N_{ch}`$ matrices.

The eigenvalues of $`\overline{Q}_0`$, appear in inverse pairs $`(q_N,q_N^1)`$, and are related with the transmission coefficients $`T_N`$ in the following way

$`A={\displaystyle \frac{q_N+q_N^1}{2}}={\displaystyle \frac{2T_N}{T_N}}`$

being

$`\left|B\right|^2=A^21.`$

By performing the trace over directions of the mode indices $`\sigma (\sigma ^{})`$, the expression for the 0th order current reduces to

$$\stackrel{ˇ}{I}^{(0)}(\epsilon )e^2\underset{N}{}2T_N\frac{[\stackrel{ˇ}{G}_2,\stackrel{ˇ}{G}_1]}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}.$$
(27)
Eq.(27) is the expression for the matrix current obtained in Ref.() for a spin-independent contact.

### A CASE I: WEAK FERROMAGNETIC CONTACT.

Now we concentrate in the first order term $`\stackrel{ˇ}{I}^{(1)}(\epsilon )`$, given by Eq.(A16) of the Appendix, $`\delta \overline{X}\left(\delta \overline{X}\right)_A`$ being anti-symmetric with respect to time-reversal transformation. In particular $`\left(\delta \overline{X}\right)_A`$ is anti-symmetric with respect to transformation over spin space and symmetric with respect to transformation over directional space. In this case $`\left(\delta \overline{X}\right)_A`$ describes the case of a weak ferromagnetic contact $`\left(\delta \overline{X}\right)_A\left(\stackrel{}{M}\text{ }\widehat{\stackrel{}{\sigma }}\text{ }\widehat{\tau }_3\right)`$. The unity vector $`\stackrel{}{M}`$ is in the direction of the magnetization, and $`\widehat{\stackrel{}{\sigma }}`$, $`\widehat{\tau }`$ are Pauli matrices in spin and Nambu space, respectively. On the other hand, the $`\left(\delta \overline{X}\right)_S`$ may describe spin-flip processes due to spin-orbit interaction at the contact. At first order in $`\delta \overline{X}`$ the contribution given by $`\left(\delta \overline{X}\right)_S`$ vanishes. In general, to threat spin-flip at the contact, we need to go to higher orders in $`\delta X`$ $`(\delta X^2,\delta X^3,..)`$. From Eq.(A88), Eq.(A89) and Eq.(A90) (see Appendix), $`\stackrel{ˇ}{I}_1`$ can be written as

$`\stackrel{ˇ}{I}^{(1)}(\epsilon )`$ $`=`$ $`e^2{\displaystyle \underset{N}{}}{\displaystyle \frac{2}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}\times `$ (30)
$`[\{t_2^{}\delta t_2+\delta t_2^{}t_2,\stackrel{ˇ}{G}_2\}+\stackrel{ˇ}{C},\stackrel{ˇ}{G}_1]`$
$`\times {\displaystyle \frac{2}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}`$
where

$$\stackrel{ˇ}{C}=(t_2^{}\delta t_2\delta t_2^{}t_2)(\frac{42T_N}{T_N})4\left(\text{ }t_2^{}\delta t_2+r_1^\text{ }\delta r_1\right).$$
(31)
Eq.(30) is the spin-dependent correction to the matrix current given in Eq.(27). This result constitutes an important step in the applicability of the circuit theory for F$`|`$N$`|`$S systems, because it provides of a simple prescription to describe magnetically active interfaces/contacts.

### B CASE II: TUNNELING BARRIERS

For the case when the contact are tunneling barriers $`(T_N1),`$ is possible to neglect the term $`T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)`$ in the denominators in Eq.(27) and Eq.(30). Keeping terms of order $`T_N`$, $`t_2^{}\delta t_2`$/$`T_N`$, $`t_2^{}\delta t_2`$ and $`r_1^\text{ }\delta r_1`$, we can write the total matrix current $`\stackrel{ˇ}{I}(\epsilon )`$ in a very transparent way like

$`{\displaystyle \frac{\stackrel{ˇ}{I}(\epsilon )}{2\pi \mathrm{}}}`$ $`=`$ $`{\displaystyle \frac{G_{\text{}}}{2}}[\stackrel{ˇ}{G}_2,\stackrel{ˇ}{G}_1]+{\displaystyle \frac{G_{\text{MR }}}{4}}[\{\stackrel{}{M}\text{ }\widehat{\stackrel{}{\sigma }}\text{ }\widehat{\tau }_3,\stackrel{ˇ}{G}_2\},\stackrel{ˇ}{G}_1]`$ (34)
$`+i{\displaystyle \frac{G_\varphi }{2}}[\stackrel{}{M}\text{ }\widehat{\stackrel{}{\sigma }}\text{ }\widehat{\tau }_3,\stackrel{ˇ}{G}_1].`$
where

$`G_\text{T}`$ $`=`$ $`G_Q{\displaystyle \underset{N}{}}T_N`$ (35)
$`G_{\text{MR}}`$ $`=`$ $`G_Q{\displaystyle \underset{N}{}}\delta T_N`$ (36)
$`iG_\varphi /2`$ $`=`$ $`G_Q{\displaystyle \underset{N}{}}(t_2^{}\delta t_2+r_1^{}\delta r_1+`$ (38)
$`\left({\displaystyle \frac{T_N2}{2T_N}}\right)(t_2^{}\delta t_2\delta t_2^\text{ }t_2))`$
$`G_Qe^2/2\pi \mathrm{}`$ being the conductance quantum and where the explicit form of $`\delta \overline{X}=\delta X`$ $`\stackrel{}{M}`$ $`\widehat{\stackrel{}{\sigma }}`$ $`\widehat{\tau }_3`$ is used. For a weak ferromagnetic contact, the elements $`\delta t`$ $`(\delta r)`$ introduced in Eqs. (A54)-(A57) in the Appendix can be expressed in terms of spin dependent amplitudes like

$$\widehat{t}=\left[\begin{array}{cc}t^{}& 0\\ 0& t^{}\end{array}\right]=\left[\begin{array}{cc}t+\delta t& 0\\ 0& t\delta t\end{array}\right]$$
(39)
$$\widehat{r}=\left[\begin{array}{cc}r^{}& 0\\ 0& r^{}\end{array}\right]=\left[\begin{array}{cc}r+\delta r& 0\\ 0& r\delta r\end{array}\right]$$
(40)
From these it easy to see that

$$G_{\text{}}=\underset{N}{}G_{Q\text{ }}T_N=\underset{N}{}G_{Q\text{ }}\frac{T_N^{}+T_N^{}}{2}=\frac{G_\text{T}^{}+G_\text{T}^{}}{2}$$
(41)
$$G_{\text{MR }}=\underset{N}{}G_{Q\text{ }}\delta T_N=\underset{N}{}G_{Q\text{ }}\frac{T_N^{}T_N^{}}{2}=\frac{G_\text{T}^{}G_\text{T}^{}}{2}.$$
(42)
Now we see that $`G_{\text{}}`$is the usual conductance of the contact $`G_{\text{}}=\left(G_\text{T}^{}+G_\text{T}^{}\right)/2`$, and $`G_{\text{MR}}=\left(G_\text{T}^{}G_\text{T}^{}\right)/2`$ accounts for the different conductances for different spin directions, leading to a spin polarized current through the junction.

On the other hand, for the case $`T_N=0`$ (the case of a ferromagnetic insulator contact), Eq.(34) is not zero and reduces to

$$\stackrel{ˇ}{I}=i\frac{G_\varphi }{2}[\stackrel{}{M}\text{ }\widehat{\stackrel{}{\sigma }}\text{ }\widehat{\tau }_3,\stackrel{ˇ}{G}_1].$$
(43)
where $`G_\varphi `$ depends only on reflection amplitudes at the normal metal side

$$iG_\varphi /2=G_{Q\text{ }}\underset{N}{}\frac{\left(r_1^{}\left(r_1^{}\right)^{}r_1^{}\left(r_1^{}\right)^{}\right)_N}{4}.$$
(44)
The physical meaning of this term in this particular case, can be now understood as follows: electrons with different spin directions pick up different phases when reflecting at the ferromagnetic insulator. During this process, ferromagnetic correlations are induced in the normal metal node. In this particular case, the coefficient $`G_\varphi `$ is related to the mixing conductance introduced in  via $`G_\varphi =\text{Im}G^{}`$. Eq.(34) has been recently used as boundary condition current for the circuit theory applied in a MI$`|`$N$`|`$S structure, being MI a magnetic insulator and under conditions of superconducting proximity effect. In Ref.() we shown that in this case, the effect of the conductance $`G_\varphi `$ can be seen as an induced magnetic field which give rise to spin-splitting of the induced “BCS like” density of states $`\stackrel{~}{h}G_\varphi \delta /2G_Q`$, being $`\delta `$ the average level spacing in the normal node. In particular,we also show that for a system composed by two coupled MI$`|`$N$`|`$F trilayer structures the absolute spin-valve effect can be achieved for a finite range of voltages.

## VI CONCLUSIONS

The circuit theory of mesoscopic transport is a systematic way to describe transport in multiterminal hybrid structures in which general rules analog to the Kirchhoff’s rules of classical circuits are used to solve the circuit and compute the transport properties of the system under study. One of these rules imposes that the sum of all “matrix current” into a node must be zero. That is why, the expression of the “spin-charge” matrix current that flows between two nodes (or between a reservoir and a node) through a contact/connector is as essential for the the circuit theory, as the Ohm’s law is for classical electric circuits. In this paper we have generalized the expression of such matrix currents for the case of multiterminal systems that included ferromagnetic and superconducting reservoirs connected through magnetically active contacts to one or several normal nodes.

We have derived the most general expression for this matrix current Eq.(23) in terms of of isotropic quasiclassical Green’s functions $`\stackrel{ˇ}{G}_{1(2)}`$ at both sides of the contact which describe the adjacent reservoirs/nodes and transmission/reflection coefficients that characterize the scattering in the contact region. This expressions should be numerically implemented in order to solve any general arrangement of reservoirs, contact and nodes.

Moreover, we have perform a perturbation expansion in the spin asymmetry of the transfer matrix $`\delta \overline{X}`$ associated to the contact region in order to gain more knowledge and obtain more transparent expressions. We found the expression for the matrix current that describes a weak ferromagnetic contact to first order in $`\delta \overline{X}`$ (Eq.(30)). In order to describe spin-flip processes at the contacts we need to go to higher orders in the asymmetry $`\delta \overline{X}`$. For the case of a tunnel barrier, the tunneling current takes a very simple and clear form (Eq.(34)). This expression is characterized by three conductance parameters $`G_{\text{}}`$,$`G_{\text{MR }}`$and $`G_\varphi `$, which can be expressed in terms of the spin conductances $`G_{()}`$ and the (complex) mixing conductance $`G^{}`$.

We thank discussions with Ya. M. Blanter, A. Brataas, and Gerrit E. W. Bauer. This work was financially supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM).

## A PERTURBATION EXPANSION

We perform a perturbation expansion of $`\left(\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M}\right)^1`$ in terms of the parameter $`\delta \overline{X}`$ as

$$\left(\stackrel{ˇ}{1}+\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M}\right)^1=\frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\stackrel{ˇ}{A}}\frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\stackrel{ˇ}{A}}\delta \overline{M}^{(1)}\frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\stackrel{ˇ}{A}}+..$$
(A2)
where

$$\overline{Q}_0=\overline{M}_0^{}\text{ }\overline{M}_0$$
(A3)
$$\stackrel{ˇ}{A}=\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2$$
(A4)
$$\delta \overline{M}^{(1)}=\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\overline{Q}_0\delta \overline{X}+\stackrel{ˇ}{G}_1\delta \overline{X}^{}\overline{Q}_0\stackrel{ˇ}{G}_2$$
(A5)
To first order in $`\delta \overline{X}`$ an expansion in terms of normal eigenmodes is possible. The matrix $`\overline{M}\left(\overline{M}\right)_{n\sigma ,m\sigma ^{}}^{s;s^{}\text{ }}`$ can be diagonalized in the basis of eigenmodes $`N(N)`$. For each eigenmode $`N`$, the elements of the transfer matrix $`\left(\overline{M}\right)_{\sigma ,\sigma ^{}}^{\text{ }s;s^{}}`$ can be expressed in terms of spin-dependent transmission and reflection amplitudes at each side of the junction as

$$m_{\text{ }+,+}^{s;s^{}}=t_2^{s,s^{}}\underset{\alpha \beta }{}r_1^{s,\alpha }\left(t_1^{\alpha ,\beta }\right)^1r_2^{\beta ,s^{}}$$
(A7)
$$m_{+,}^{s;s^{}\text{ }}=\underset{\alpha }{}r_1^{s,\alpha }\left(t_1^{\alpha ,s^{}}\right)^1$$
(A8)
$$m_{,+}^{\text{ }s;s^{}}=\underset{\alpha }{}\left(t_1^{s,\alpha }\right)^1r_2^{\alpha ,s^{}}$$
(A9)
$$m_,^{s;s^{}\text{ }}=\left(t_1^{s,s^{}}\right)^1$$
(A10)
,where $`s,s(\alpha ,\beta )`$ $`\{,\}`$ are spin labels and $`1(2)`$ denotes left(right) side of the contact region. Usually the transmission and reflection amplitudes are introduced in terms of the scattering matrix $`S`$

$$\overline{S}=\left[\begin{array}{cc}r_2^{s\text{ }s^{^{}}}& t_1^{s\text{ }s^{^{}}}\\ t_2^{s\text{ }s^{^{}}}& r_1^{s\text{ }s^{^{}}}\end{array}\right].$$
(A11)
where in this case $`t^{s\text{ }s^{^{}}}`$ and $`r^{s\text{ }s^{^{}}}`$are matrices in Spin space. The scattering matrix $`\overline{S}`$ and the transfer matrix $`\overline{M}`$ are equivalent descriptions of a contact region. Nevertheless the transfer matrix $`\overline{M}`$ obeys a multiplicative composite rule, whereas the scattering matrix obeys a more complicated composition rule..

By writing in analogous way the other term in Eq.(23) $`\left(\stackrel{ˇ}{1}\stackrel{ˇ}{G}_1\overline{M}^{}\text{ }\stackrel{ˇ}{G}_2\text{ }\overline{M}\right)`$ and after some algebra, we obtain the following expression for the matrix current $`\stackrel{ˇ}{I}`$ to first order in $`\delta \overline{X}`$

$$\stackrel{ˇ}{I}(\epsilon )=\stackrel{ˇ}{I}^{(0)}(\epsilon )+\stackrel{ˇ}{I}^{(1)}(\epsilon )$$
(A13)
where

$$\stackrel{ˇ}{I}^{(0)}(\epsilon )=e^2\text{Tr}_{N,\sigma }\left\{\frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2}\left(\stackrel{ˇ}{1}+2\text{ }\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\right)\right\}$$
(A14)
$`\stackrel{ˇ}{I}^{(1)}(\epsilon )`$ $`=`$ $`2\text{ }e^2\text{Tr}_{N,\sigma }\{{\displaystyle \frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2}}\stackrel{ˇ}{0}{\displaystyle \frac{\stackrel{ˇ}{1}}{\stackrel{ˇ}{1}+\overline{Q}_0\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2}}`$ (A16)
$`\times (\stackrel{ˇ}{1}+\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1)\}`$
with

$$\stackrel{ˇ}{0}=\left(\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\overline{Q}_0\delta \overline{X}+\stackrel{ˇ}{G}_1\delta \overline{X}^{}\overline{Q}_0\stackrel{ˇ}{G}_2\right).$$
(A17)
A general property of the transfer matrix is the flux conservation which can be expressed as

$$\overline{M}^{}\text{ }\overline{\mathrm{\Sigma }}^z\text{ }\overline{M}=\overline{\mathrm{\Sigma }}^z$$
(A18)
Substituting Eq. (24) in Eq.(A18) we get

$$(\stackrel{ˇ}{1}+\delta \overline{X}^{})\text{ }\overline{M}_0^{}\text{ }\overline{\mathrm{\Sigma }}^z\text{ }\overline{M}_0\text{ }(\stackrel{ˇ}{1}+\delta \overline{X})=\overline{\mathrm{\Sigma }}^z$$
(A19)
The matrix $`\overline{M}_{0\text{ }}`$ should also obey flux conservation

$$\overline{M}_0^{}\text{ }\overline{\mathrm{\Sigma }}^z\text{ }\overline{M}_0=\overline{\mathrm{\Sigma }}^z\text{,}$$
(A20)
so we found the following equation for $`\delta \overline{X}`$

$$\delta \overline{X}^{}\overline{\mathrm{\Sigma }}^z+\overline{\mathrm{\Sigma }}^z\text{ }\delta \overline{X}+\delta \overline{X}^{}\overline{\mathrm{\Sigma }}^z\text{ }\delta \overline{X}=0.$$
(A21)
To first order in $`\delta \overline{X}`$ we have

$$\delta \overline{X}^{}=\overline{\mathrm{\Sigma }}^z\text{ }\delta \overline{X}\text{ }\overline{\mathrm{\Sigma }}^z.$$
(A22)
From Eq.(A22) we obtain for the elements of $`\delta \overline{X}`$

$$\delta \overline{X}_{+\text{ }+}^{}=\delta \overline{X}_{+\text{ }+}$$
(A24)
$$\delta \overline{X}_{+\text{ }}^{}=\delta \overline{X}_{\text{ }+}$$
(A25)
$$\delta \overline{X}_\text{ }^{}=\delta \overline{X}_\text{ }.$$
(A26)
We see that $`\delta \overline{X}_{+\text{ }+}`$ and $`\delta \overline{X}_\text{ }`$ are pure complex elements and $`\delta \overline{X}_{+\text{ }}`$ and $`\delta \overline{X}_{\text{ }+}`$ are complex conjugate.

The matrix $`\overline{M}_0`$

$$\overline{M}_0=\left[\begin{array}{cc}m_{+\text{ }+}^0& m_{+\text{ }}^0\\ m_{\text{ }+}^0& m_\text{ }^0\end{array}\right]$$
(A27)
is symmetric under time-reversal transformation, which involves both transformation in Spin and Channels space, under the following prescription

$$\stackrel{~}{\overline{M}_0}=\stackrel{}{\sigma }^y\text{ }\overline{\mathrm{\Sigma }}^x\text{ }\left(\overline{M}_0\right)^{}\text{ }\overline{\mathrm{\Sigma }}^x\text{ }\stackrel{}{\sigma }^y=\overline{M}_0$$
(A28)
being $`\overline{\mathrm{\Sigma }}^x`$ the x-Pauli matrix in directional space and $`\stackrel{}{\sigma }^y`$ the y-Pauli matrix in Spin space. As a result, the components of $`\overline{M}_0`$ fulfill the following relations

$$m_{\text{ }+\text{ }+}^{0\text{ }}=m_{\text{ }\text{ }}^0$$
(A30)
$$m_{\text{ }+\text{ }}^{0\text{ }}=m_{\text{ }\text{ }+}^0$$
(A31)
In general $`\delta \overline{X}`$ can be split in symmetric and anti-symmetric parts with respect to time-reversal

$$\delta \overline{X}=\left(\delta \overline{X}\right)_A+\left(\delta \overline{X}\right)_S$$
(A33)
$$\left(\stackrel{~}{\delta \overline{X}}\right)_A=\stackrel{}{\sigma }^y\overline{\mathrm{\Sigma }}^x\left(\delta \overline{X}\right)^{}\overline{\mathrm{\Sigma }}^x\stackrel{}{\sigma }^y=\text{ }\left(\delta \overline{X}\right)_A$$
(A34)
$$\left(\stackrel{~}{\delta \overline{X}}\right)_S=\stackrel{}{\sigma }^y\overline{\mathrm{\Sigma }}^x\left(\delta \overline{X}\right)^{}\overline{\mathrm{\Sigma }}^x\stackrel{}{\sigma }^y=\left(\delta \overline{X}\right)_S$$
(A35)
Time-reversal transformation over Spin gives always an anti-symmetric contribution $`\stackrel{~}{\delta \overline{X}}=\stackrel{}{\sigma }^y\left(\delta \overline{X}\right)^{}\stackrel{}{\sigma }^y=`$ $`\delta \overline{X}`$, so $`\left(\stackrel{~}{\delta \overline{X}}\right)_{A(S)}`$ will be symmetric (anti-symmetric) with respect to transformation in directional space $`\left(\stackrel{~}{\delta \overline{X}}\right)_{A(S)}=\overline{\mathrm{\Sigma }}^x\left(\delta \overline{X}\right)^{}\overline{\mathrm{\Sigma }}^x=\pm \left(\delta \overline{X}\right)_{A(S)}`$. Taking into account the structure of $`\delta \overline{X}`$ given by the condition $`\delta \overline{X}^{}=\overline{\mathrm{\Sigma }}^z\delta \overline{X}`$ $`\overline{\mathrm{\Sigma }}^z,`$ $`\stackrel{~}{\delta \overline{X}}`$

$$\stackrel{~}{\delta \overline{X}}=\left[\begin{array}{cc}\delta X_\text{ }& \delta X_{+\text{ }}\\ \delta X_{+\text{ }}^{}& \delta X_{+\text{ }+}\end{array}\right].$$
(A37)
we find for $`\left(\stackrel{~}{\delta \overline{X}}\right)_A=\overline{\mathrm{\Sigma }}^x\left(\delta \overline{X}\right)^{}\overline{\mathrm{\Sigma }}^x=+\left(\delta \overline{X}\right)_A`$,

$$\stackrel{~}{\delta \overline{X}}=\left[\begin{array}{cc}\delta X_\text{ }& \delta X_{+\text{ }}\\ \delta X_{+\text{ }}^{}& \delta X_{+\text{ }+}\end{array}\right]=\left[\begin{array}{cc}\delta X_{+\text{ }+}& \delta X_{+\text{ }}\\ \delta X_{+\text{ }}^{}& \delta X_\text{ }\end{array}\right].$$
(A38)
In this case $`\delta X_{+\text{ }+}=\delta X_\text{ }`$ and $`\delta \overline{X}_A`$ is in directional Space

$$\delta \overline{X}_A=\left[\begin{array}{cc}\delta X_{+\text{ }+}& \delta X_{+\text{ }}\\ \delta X_{+\text{ }}^{}& \delta X_{+\text{ }+}\end{array}\right]\overline{\mathrm{\Sigma }}^z,\overline{\mathrm{\Sigma }}^x,\text{ }\overline{\mathrm{\Sigma }}^y.$$
(A39)
On the other hand for $`\delta \overline{X}_S`$, $`\left(\stackrel{~}{\delta \overline{X}}\right)_S=\overline{\mathrm{\Sigma }}^x\left(\delta \overline{X}\right)^{}\overline{\mathrm{\Sigma }}^x=\left(\delta \overline{X}\right)_S`$

$$\stackrel{~}{\delta \overline{X}}=\left[\begin{array}{cc}\delta X_\text{ }& \delta X_{+\text{ }}\\ \delta X_{+\text{ }}^{}& \delta X_{+\text{ }+}\end{array}\right]=\left[\begin{array}{cc}\delta X_{+\text{ }+}& \delta X_{+\text{ }}\\ \delta X_{+\text{ }}^{}& \delta X_\text{ }\end{array}\right].$$
(A40)
Now $`\delta X_{+\text{ }+}=\delta X_\text{ }`$, and $`\delta X_{+\text{ }}`$ must to be zero. As a result $`\delta \overline{X}_S`$ is in directional Space

$$\delta \overline{X}_S=\left[\begin{array}{cc}\delta X_{+\text{ }+}& 0\\ 0& \delta X_{+\text{ }+}\end{array}\right]\overline{1}.$$
(A41)
For the matrix $`\delta \overline{M}=\overline{M}_0\delta \overline{X}`$ we have

$`\delta M`$ $`=`$ $`\left[\begin{array}{cc}\delta m_{+\text{ }+}& \delta m_{+\text{ }}\\ \delta m_{\text{ }+}& \delta m_\text{ }\end{array}\right]`$ (A44)
$`=`$ $`\left[\begin{array}{cc}_\sigma m_{+\text{ }\sigma }\delta \overline{X}_{\sigma \text{ }+}& _\sigma m_{+\text{ }\sigma }\delta \overline{X}_{\sigma \text{ }}\\ \begin{array}{c}\\ _\sigma m_{\text{ }\sigma }\delta \overline{X}_{\sigma \text{ }+}\end{array}& \begin{array}{c}\\ _\sigma m_{\text{ }\sigma }\delta \overline{X}_{\sigma \text{ }}\end{array}\end{array}\right]`$ (A51)
where $`\sigma \{+,\}`$. Analogy to $`\overline{M}`$, the matrix $`\overline{S}`$ can be written

$$\overline{S}=\overline{S}_0+\delta \overline{S}$$
(A52)
where $`\overline{S}_0`$ corresponds to the spin-independent part and $`\delta \overline{S}`$ accounts for the spin structure of $`\overline{S}`$. Now we assume that $`\delta \overline{S}`$ is a small deviation with respect to $`\overline{S}_0`$ $`(\delta \overline{S}\overline{S}_0).`$ In this case the elements of $`\overline{S}`$ are

$$r_1^{s\text{ }s^{^{}}}=r_1+\delta r_1$$
(A54)
$$t_1^{s\text{ }s^{^{}}}=t_1+\delta t_1$$
(A55)
$$r_2^{s\text{ }s^{^{}}}=r_2+\delta r_2$$
(A56)
$$t_2^{s\text{ }s^{^{}}}=t_2+\delta t_2.$$
(A57)
The quantities $`t_{1(2)}\left(r_{1(2)}\right)`$ are spin-independent transmission(reflection) amplitudes, whereas $`\delta t_{1(2)}^{s\text{ }s^{^{}}}\left(\delta r_{1(2)}^{s\text{ }s^{^{}}}\right)`$ are the spin-dependent correction to the total transmission(reflection) $`t_1^{s\text{ }s^{^{}}}\left(r_2^{s\text{ }s^{^{}}}\right)`$ amplitudes. The elements of $`\overline{S}`$ can be expressed in terms of the elements of $`\overline{M}`$ like

$$r_2^{s\text{ }s^{^{}}}=(m_{\text{ }\text{ }}^{s;s^{}})^1\text{ }m_{\text{ }\text{ }+}^{s;s^{}}$$
(A59)
$$t_2^{s\text{ }s^{^{}}}=m_{+\text{ }+}^{\text{ }s;s^{}}m_{+\text{ }}^{\text{ }s;s^{}}(m_\text{ }^{\text{ }s;s^{}})^1\text{ }m_{\text{ }+}^{\text{ }s;s^{}}$$
(A60)
$$r_1^{s\text{ }s^{^{}}}=m_{+\text{ }}^{\text{ }s;s^{}}(m_{\text{ }\text{ }^{}}^{s;s})^1$$
(A61)
$$t_1^{s\text{ }s^{^{}}}=(m_\text{ }^{\text{ }s;s^{}})^1.\text{ }$$
(A62)
By substituting Eq.(A44) and Eqs.(A54-A57), using analog relations between the elements of $`\overline{S}_0`$ and $`\overline{M}_0`$ and by expanding $`(m_{s;s^{}}^{\text{ }\text{ }})^1`$ up to order $`\delta m_\text{ }`$, we obtain

$$\delta r_2^{s\text{ }s^{^{}}}=r_2\text{ }\delta X_{+\text{ }+}^{\text{ }s\text{ }s^{^{}}}\delta X_{\text{ }+}^{\text{ }s\text{ }s^{}}+(r_2\text{ }\delta X_{+\text{ }}^{\text{ }s\text{ }s^{^{}}}\delta X_\text{ }^{\text{ }s\text{ }s^{^{}}})\text{ }r_2$$
(A64)
$$\delta t_2^{s\text{ }s^{^{}}}=t_2(\text{ }\delta X_{+\text{ }+}^{\text{ }s\text{ }s^{^{}}}+r_2\text{ }\delta X_{+\text{ }}^{\text{ }s\text{ }s^{^{}}})$$
(A65)
$$\delta r_1^{s\text{ }s^{^{}}}=t_2\text{ }t_1\text{ }\delta X_{+\text{ }}^{\text{ }s\text{ }s^{^{}}}$$
(A66)
$$\delta t_1^{s\text{ }s^{^{}}}=t_1(r_2\text{ }\delta X_{+\text{ }}^{\text{ }s\text{ }s^{^{}}}\delta X_\text{ }^{\text{ }s\text{ }s^{^{}}}).$$
(A67)
So we see that the elements of $`\delta \overline{X}`$ can be expressed in terms of the transmission and reflection amplitudes $`\delta t`$ $`(\delta r)`$. Note that in Eq.(A64)-Eq.(A67) we explicitly include the spin indices $`s,s`$ of the elements $`\delta \overline{X}`$ and $`\delta t`$ $`(\delta r)`$ to emphasize that both $`\delta \overline{X}`$ and $`\delta t`$ $`(\delta r)`$, contain information about the spin structure of the contact region.

Now we consider the case of $`\delta \overline{X}_A`$ which corresponds to a weak ferromagnetic contact $`\left(\left(\delta \overline{X}\right)_A\left(\stackrel{}{M}\text{ }\widehat{\stackrel{}{\sigma }}\text{ }\widehat{\tau }_3\right)\right)`$. To evaluate the trace of Eq. (A16) we can re-write the expression in the following way

$`\stackrel{ˇ}{I}^{(1)}(\epsilon )`$ $`=`$ $`e^2\text{Tr}_{N,\sigma }\{2{\displaystyle \frac{T_N}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}\times `$ (A71)
$`\left(\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1+\overline{Q}_0^1\right)\left(\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\overline{Q}_0\delta \overline{X}+\stackrel{ˇ}{G}_1\delta \overline{X}^{}\overline{Q}_0\stackrel{ˇ}{G}_2\right)`$
$`\left(\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1+\overline{Q}_0^1\right)\left(\stackrel{ˇ}{1}+\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1\right)`$
$`\times {\displaystyle \frac{T_N}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}\}`$
The central term can be separate in four different terms as

$$\left(\begin{array}{c}\overline{Q}_0\delta \overline{X}\text{ }\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1+\stackrel{ˇ}{G}_2\delta \overline{X}^{}\overline{Q}_0\stackrel{ˇ}{G}_1+\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\text{ }\delta \overline{X}\text{ }\overline{Q}_0^1\\ +\stackrel{ˇ}{G}_1\text{ }\overline{Q}_0^1\delta \overline{X}^{}\text{ }\stackrel{ˇ}{G}_2\end{array}\right)+$$
(A73)
$$\left(\begin{array}{c}\overline{Q}_0\delta \overline{X}\overline{Q}_0^1+\stackrel{ˇ}{G}_2\delta \overline{X}^{}\stackrel{ˇ}{G}_2+\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\text{ }\delta \overline{X}\text{ }\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1\\ +\stackrel{ˇ}{G}_1\text{ }\overline{Q}_0^1\delta \overline{X}^{}\text{ }\overline{Q}_0\text{ }\stackrel{ˇ}{G}_1\end{array}\right)+$$
(A74)
$$\left(\begin{array}{c}\overline{Q}_0\delta \overline{X}\text{ }\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1+\stackrel{ˇ}{G}_2\delta \overline{X}^{}\overline{Q}_0\stackrel{ˇ}{G}_1+\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\text{ }\delta \overline{X}\text{ }\overline{Q}_0^1\\ +\stackrel{ˇ}{G}_1\text{ }\overline{Q}_0^1\delta \overline{X}^{}\text{ }\stackrel{ˇ}{G}_2\end{array}\right)\left(\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1\right)+$$
(A75)
$$\left(\begin{array}{c}\overline{Q}_0\delta \overline{X}\overline{Q}_0^1+\stackrel{ˇ}{G}_2\delta \overline{X}^{}\stackrel{ˇ}{G}_2+\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\text{ }\delta \overline{X}\text{ }\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1\\ +\stackrel{ˇ}{G}_1\text{ }\overline{Q}_0^1\delta \overline{X}^{}\text{ }\overline{Q}_0\text{ }\stackrel{ˇ}{G}_1\end{array}\right)\left(\overline{\mathrm{\Sigma }}^z\text{ }\stackrel{ˇ}{G}_1\right).$$
(A76)
The trace over “direction” indices $`\sigma (\sigma ^{})`$ gives for the first term Eq.(A73)

$`2Re[B^\text{ }\delta X_{12}]\text{ }\stackrel{ˇ}{G}_2\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_2\text{ }(2Re[B^{}\delta X_{12}])`$ (A77)
$`+\stackrel{ˇ}{G}_2\text{ }2Re[B^\text{ }\delta X_{12}]\text{ }\stackrel{ˇ}{G}_1\stackrel{ˇ}{G}_1\text{ }2Re[B^{}\delta X_{12}]\text{ }\stackrel{ˇ}{G}_2.`$ (A78)
Note that $`\stackrel{ˇ}{G}_{1(2)}`$ and $`\delta X_{12}`$ are still matrices in “Keldysh-Nambu-spin” space and do not commute with each other. $`B^{}`$ can be written in terms of reflection and transmission amplitudes like

$$B^{}=2m_{+\text{ }+}m_{+\text{ }}=\frac{2}{T_N}r_2.$$
(A79)
The spin-dependent corrections to the transmission and reflection probabilities are defined like

$$\delta T_{1(2)}=t_{1(2)\text{ }}^\text{ }\delta t_{1(2)}+\delta t_{1(2)}^\text{ }t_{1(2)}$$
(A80)
$$\delta R_{1(2)}=r_{1(2)\text{ }}^{}\delta r_{1(2)}+\delta r_{1(2)}^\text{ }r_{1(2)}.$$
(A81)
For this quantities we find

$$\frac{\delta T_1}{T_1}\frac{\delta T_2}{T_2}\frac{\delta R_1}{T_1}\frac{\delta R_2}{T_2}=2Re[r_2\text{ }\delta X_{+\text{ }}^{\text{ }s\text{ }s^{^{}}}].$$
(A82)
so

$$2Re[B^{}\delta X_{+\text{ }}^{\text{ }s\text{ }s^{}}]=\frac{2}{T_N}(2Re[r_2\delta X_{+\text{ }}^{\text{ }s\text{ }s^{}}])=\frac{2}{T_N}\frac{\delta T_N}{T_N}.$$
(A83)
At first order in $`\delta \overline{X}`$, the contribution given by Eq.(A83) is zero for in the case of $`\left(\delta \overline{X}\right)_S`$. Finally The contribution to this first term for the $`\stackrel{ˇ}{I}_1`$ current in terms of $`T_N,`$ $`\delta T_N`$ and $`\stackrel{ˇ}{G}_{1(2)}`$ is

$`\stackrel{ˇ}{I}^{(1)}(\epsilon )`$ $`=`$ $`e^2{\displaystyle \underset{N}{}}{\displaystyle \frac{2}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}[\{\delta \stackrel{ˇ}{T}_N,\stackrel{ˇ}{G}_2\},\stackrel{ˇ}{G}_1]`$ (A85)
$`\times {\displaystyle \frac{2}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}+..`$
For the second term (Eq.A74), using the property that the trace is cyclic, it is easy to see that Tr$`{}_{\sigma }{}^{}\left\{\delta \overline{X}\right\}`$Tr$`{}_{\sigma }{}^{}\left\{\delta \overline{X}^{}\right\}`$Tr$`{}_{\sigma }{}^{}\left\{\overline{Q}_0\delta \overline{X}\overline{Q}_0^1\right\}`$

$``$Tr$`{}_{\sigma }{}^{}\left\{\overline{Q}_0^1\delta \overline{X}^{}\overline{Q}_0\right\}=0`$, so the second term gives no contribution to $`\stackrel{ˇ}{I}_1`$.

The third and fourth terms include the factor $`\overline{\mathrm{\Sigma }}^z`$ $`\stackrel{ˇ}{G}_1`$ in the right hand side. The presence of the matrix $`\overline{\mathrm{\Sigma }}^z`$ will change the structure of this term in Channels space with respect to Eqs.A73 and A74. Using that $`[\delta \overline{X},\stackrel{ˇ}{G}_2]=0`$, $`\stackrel{ˇ}{G}_2`$ being the Green’s function of a (normal/ferromagnetic) reservoir, finally $`\stackrel{ˇ}{I}_1`$ gives

$`\stackrel{ˇ}{I}^{(1)}(\epsilon )`$ $`=`$ $`e^2\text{ }{\displaystyle \underset{N}{}}{\displaystyle \frac{2}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}\times `$ (A88)
$`[\{\delta \stackrel{ˇ}{T}_N,\stackrel{ˇ}{G}_2\}+\left(42T_N\right)\text{ }\delta \stackrel{ˇ}{\mathrm{\Theta }}_N4\text{ }\delta \stackrel{ˇ}{\mathrm{\Xi }}_N\text{ },\stackrel{ˇ}{G}_1]`$
$`\times {\displaystyle \frac{2}{4+T_N(\{\stackrel{ˇ}{G}_1,\stackrel{ˇ}{G}_2\}2)}}`$
being $`\delta \stackrel{ˇ}{\mathrm{\Xi }}_N=T_N`$ $`\delta X_{+\text{ }+,N}`$ and $`\delta \stackrel{ˇ}{\mathrm{\Theta }}_N=2\left[\delta X_{+\text{ }+}+iIm\left(r_2\text{ }\delta X_{+\text{ }}\right)\right]_N`$ both pure imaginary quantities. Finally inverting Eq.(A64)-Eq.(A67) we can express also the quantities $`\delta \stackrel{ˇ}{\mathrm{\Xi }}_N`$ and $`\delta \stackrel{ˇ}{\mathrm{\Theta }}_N`$ in terms of spin-independent transmission(reflection) amplitudes $`t_{1(2)}\left(r_{1(2)}\right)`$,and its the spin-dependent corrections $`\delta t_{1(2)}\left(\delta r_{1(2)}\right)`$

$$\delta \stackrel{ˇ}{\mathrm{\Xi }}_N=t_2^{}\delta t_2+r_1^{}\delta r_1$$
(A89)
$$\delta \stackrel{ˇ}{\mathrm{\Theta }}_N=\frac{t_2^{}\delta t_2\delta t_2^\text{ }t_2}{T_N}$$
(A90)


format: text

text:
# ENERGY AND MOMENTUM OF A STATIONARY BEAM OF LIGHT

## I Introduction: <br>

Since Einstein proposed the General Theory of Relativity, relativists have not been able to agree upon a definition of the energy-momentum distribution associated with the gravitational field (see Ref. debate  and references therein). In 1990, H. BondiBondi  argued that General Relativity does not permit a non-localizable form of energy, so, in principle, we should expect to be able to find an acceptable definition.

Energy and momentum density are usually defined by a second rank tensor $`T_i^k`$. The conservation of energy and momentum are described by the requirement that the tensor’s divergence is zero. However, in General Relativity, the partial derivative in the usual conservation equation $`T_{i,k}^k=0`$ is replaced by a covariant derivative. $`T_i^k`$ then represents the energy and momentum of matter and all non-gravitational fields and no longer satisfies $`T_{i,k}^k=0`$. A contribution from the gravitational field must be added to obtain an energy-momentum expression with zero divergence. Einstein first obtained such an expression and many others such as Landau and Lifshitz, Papapetrou, and Weinberg gave similar prescriptions (see Ref. definitions ). The expressions they gave are called energy-momentum complexes because they can be expressed as a combination of $`T_i^k`$ and a pseudotensor, which is interpreted to represent the energy and momentum of the gravitational field. These complexes have been heavily criticized because they are non-tensorial, i.e. they are coordinate dependent. For the Einstein, Landau-Lifshitz, Papapetrou, and Weinberg (later ELLPW) energy-momentum complexes, one gets physically meaningful results only in ”Cartesian coordinates” (see Ref. Cartesian ; Moller ). Because of this drawback, many others, including Møller Moller , KomarKomar , and PenrosePenrose , have proposed coordinate independent definitions. Each of these, however, has its own drawbacks (for further discussion see Ref. problems ).

Recently, the energy momentum complexes of ELLPW have shown promise and may prove to be more useful than the others. Much work has been produced in the past decade showing the these complexes give meaningful results for many well known metrics ( Ref. results ; VirKN ; VirFlat ). Also, Aguirregabiria et al.Aguirr  showed that they coincide for all Kerr-Schild class metrics which include the Schwarzschild, Reissner-Nordström, Kerr-Newman, Vaidya, Dybney et al., Kinnersley, Bonnor-Vaidya, and Vaidya-Patel spacetimes (for references see in spacetimes ).

In this Letter, we will calculate, using the energy momentum complexes of ELLPW, the energy and momentum densities of the metric given by BonnorBonnor  describing a stationary beam of light. First we will show that the metric is of Kerr-Schild class, and then we will use the results given by Aguirregabiria et. al. to compute its energy and momentum densities. Finally, we will discuss the physical implications of our calculations.

We use the convention that Latin indices take values from $`0`$ to $`3`$ and Greek indices take values from $`1`$ to $`3`$, and take units where $`G=1`$ and $`c=1`$.

## II Kerr-Schild Class Metrics

A metric is said to be of Kerr-Schild class if it can be written in the form

$$g_{ab}=\eta _{ab}2Vl_al_b$$
(1)
where $`\eta `$ is the Minkowski metric, $`V`$ is a scalar field, and $`l`$ is a null, geodesic, and shear free vector field in the Minkowski spacetime. These three properties of $`l`$ are respectively expressed as

$`\eta ^{ab}l_al_b`$ $`=`$ $`0`$ (2)
$`\eta ^{ab}l_{i,a}l_b`$ $`=`$ $`0`$ (3)
$`(l_{a,b}+l_{b,a})l_{,c}^a\eta ^{bc}(l_{,a}^a)^2`$ $`=`$ $`0`$ (4)
Aguirregabiria et. al.Aguirr  showed that for general Kerr-Schild class metrics the ELLPW energy-momentum complexes coincide up to raising and lowering of indices by the Minkowski metric. Later, VirbhadraKSV  clarified that the vector field $`l`$ need only be null and geodesic for the ELLPW complexes to coincide. These energy momentum complexes for any Kerr-Schild class metric are given by

$`\mathrm{\Theta }_i^k`$ $`=`$ $`\eta _{ij}L^{jk}`$ (5)
$`L^{ik}=\mathrm{\Sigma }^{ik}`$ $`=`$ $`W^{ik}={\displaystyle \frac{1}{16\pi }}\mathrm{\Lambda }_{,lm}^{iklm}`$ (6)
where $`\mathrm{\Theta }`$, $`L`$, $`\mathrm{\Sigma }`$, and $`W`$ are the Einstein, Landau-Lifshitz, Papepetrou, and Weinberg energy momentum complexes respectively and

$$\mathrm{\Lambda }^{ikpq}=2V(\eta ^{ik}l^pl^q+\eta ^{pq}l^il^k\eta ^{ip}l^kl^q\eta ^{kq}l^il^p).$$
(7)
$`\mathrm{\Theta }_0^0`$, $`L^{00}`$, $`\mathrm{\Sigma }^{00}`$, and $`W^{00}`$ represent total energy density and $`\mathrm{\Theta }_\alpha ^0`$, $`L^{0\alpha }`$, $`\mathrm{\Sigma }^{0\alpha }`$, and $`W^{0\alpha }`$ represent momentum density in the $`x^\alpha `$ direction.

## III The Bonnor Metric

The Bonnor metric describing a stationary beam of light flowing in the $`Z`$ direction is given by the line elementBonnor

$$ds^2=dx^2dy^2(1m)dz^22mdzdt+(1+m)dt^2$$
(8)
where $`m`$ is a function of $`x`$ and $`y`$,

$`^2m`$ $`=`$ $`16\pi \rho `$ (9)
$`\rho =T_3^3`$ $`=`$ $`T_3^0=T_0^3=T_0^0`$ (10)
and $`T_a^b`$ is the energy momentum tensor.

The Bonnor metric can be rewritten as

$$ds^2=dt^2dx^2dy^2dz^2m(dtdz)^2.$$
(11)
This is in form required of a Kerr-Schild class metric with

$`V`$ $`=`$ $`m/2`$ (12)
$`l_t`$ $`=`$ $`1`$ (13)
$`l_z`$ $`=`$ $`1`$ (14)
Both components of $`l`$ are constant so $`l`$ is trivially geodesic and shear free. It can be easily shown to be null which proves that the Bonnor metric is of Kerr-Schild class.

We can then easily compute the energy momentum complexes of ELLPW for the metric. We note that, in general $`\mathrm{\Lambda }`$ exhibits the following symmetries

$$\mathrm{\Lambda }^{ikpq}=\mathrm{\Lambda }^{kiqp}=\mathrm{\Lambda }^{pkiq}$$
(15)
Then, we calculate the required, non-vanishing components of $`\mathrm{\Lambda }`$ to be

$`\mathrm{\Lambda }^{0011}`$ $`=`$ $`\mathrm{\Lambda }^{0022}=m`$ (16)
$`\mathrm{\Lambda }^{0033}`$ $`=`$ $`2m`$ (17)
$`\mathrm{\Lambda }^{0311}`$ $`=`$ $`\mathrm{\Lambda }^{0322}=m`$ (18)
$`\mathrm{\Lambda }^{0303}`$ $`=`$ $`2m`$ (19)
From here we calculate the energy-momentum complexes. The non-vanishing components are, as expected, total energy density and momentum density in the $`z`$ direction. These both are given by

$$L^{00}=L^{03}=\frac{1}{16\pi }(m_{,xx}+m_{,yy})$$
(20)
Then, recalling Equations (9) and (10)

$$L^{00}=L^{03}=\rho =T_0^0=T_0^3$$
(21)
## IV Discussion

This result is what would be expected from purely physical arguments. The energy and momentum densities calculated from the complexes of ELLPW coincide and are equal to the energy and momentum density components of $`T_i^k`$.

This simple result supports the Cooperstock hypothesisCooperstock  which states that energy is localized to the region where the energy-momentum tensor is non-vanishing. This hypothesis would imply that there is no energy-momentum contribution from the ”vacuum” regions of spacetime. If true, the hypothesis would have broad implications. For instance, a quantum theory of gravity requires that the gravitons which make up the field carry energy. Also, the hypothesis suggests that gravitational waves have no energy and that current attempts to detect these waves by detecting their energy will be fruitless.

In summary, we have calculated the energy and momentum for a stationary beam of light and found that the results are physically meaningful and support the Cooperstock hypothesis. Certainly the questions of energy and momentum localization in General Relativity as well as the Cooperstock hypothesis are far from resolved, but the Bonnor metric provides one more example where the energy momentum complexes of ELLPW provide, for a physical metric, meaningful results that support the hypothesis.

###### Acknowledgements.

I am grateful to K.S. Virbhadra for his guidance, to Duke University for its support, and also to the referee for helpful suggestions.


text:
# Phase transition in a non-conserving driven diffusive system

## Abstract

An asymmetric exclusion process comprising positive particles, negative particles and vacancies is introduced. The model is defined on a ring and the dynamics does not conserve the number of particles. We solve the steady state exactly and show that it can exhibit a continuous phase transition in which the density of vacancies decreases to zero. The model has no absorbing state and furnishes an example of a one-dimensional phase transition in a homogeneous non-conserving system which does not belong to the absorbing state universality classes.

One dimensional driven diffusive systems have been a subject of extensive study in recent years . A prototypical model for studying these systems is the asymmetric exclusion process in which particles move stochastically with a preferred direction and hard-core exclusion . In this model the local dynamics conserves the particles. Asymmetric exclusion processes with a single species of particles do not exhibit phase transitions on a ring geometry. On the other hand open systems in which the particle number is not conserved at the boundaries can exhibit a variety of phase transitions . Generalisations of these models to more than one species of particles have shown that phase transitions and long-range order do exist in these systems even on a ring geometry .

Generally speaking one-dimensional models that exhibit phase transitions are either (a) of the asymmetric-exclusion-process type with a drive and conserving bulk dynamics or (b) have non-conserving bulk dynamics with one or more absorbing states. The latter case is related to the directed percolation or other absorbing state universality classes . It is therefore of interest to find a one-dimensional model where the bulk dynamics is not conserving that exhibits a phase transition unrelated to the absorbing state universality classes.

In this letter we introduce such a model and solve for its steady state exactly. The model is defined on a ring of $`L`$ sites where each site can be occupied by either a positive ($`+`$) particle, a negative ($``$) particle or a vacancy ($`0`$). The model evolves through the following conserving rates:

$$+\mathrm{\hspace{0.33em}0}\stackrel{1}{}\mathrm{\hspace{0.17em}0}+0\stackrel{1}{}\mathrm{\hspace{0.17em}0}+\underset{q}{\overset{1}{}}+$$
(1)
augmented by the following non-conserving rates

$$+\mathrm{\hspace{0.17em}0}\underset{1}{\overset{w}{}}\mathrm{\hspace{0.17em}0\hspace{0.17em}0}0\underset{1}{\overset{w}{}}\mathrm{\hspace{0.17em}0\hspace{0.17em}0}.$$
(2)
Thus, the model generalizes the model of Derrida et al  and Arndt et al  through the introduction of the process of creation and annihilation of particles. Note that any configuration can be reached from any other except states with no vacancies which are not dynamically accessible.

The rate $`w`$ controls the density of vacancies, denoted by $`\theta `$, in the system. For large $`w`$ one expects the density of vacancies to be high and for small $`w`$ it is expected to be low. We will show that for $`q<1`$ the dependence of $`\theta `$ on $`w`$ is not smooth and that for small enough $`w`$, $`\theta `$ is zero in the thermodynamic limit. Thus the model exhibits a phase transition from a “fluid” phase with a finite $`\theta `$ to a “maximal current” phase where $`\theta `$ is zero (the nomenclature will be explained below). The transition occurs at a critical value $`w_c>0`$. The phase transition is found to be continuous with

$$\theta |ww_c|^\beta ,$$
(3)
where $`\beta =1`$. For $`q>1`$, on the other hand, the system is always strongly phase separated  with a single zero and two extensive pure domains of positive and negative particles. Here there is no phase transition as $`w`$ is varied.

A powerful technique for solving the steady-states of asymmetric exclusion processes is the matrix product ansatz . This involves representing the steady-state weights as the trace of a product of matrices which depends on the microscopic configuration. The matrices then obey certain algebraic rules which are derived from the dynamics of the model. This technique has yielded the exact behaviour at the various phase transitions in many asymmetric exclusion models. The transitions solved within the matrix product have been found to be robust for a large class of systems . However, in most models solved so far using the matrix product, particle numbers are conserved in the bulk (two exceptions are  and ). In this work we employ the matrix product technique to obtain an exact solution for the steady state and phase transitions of the model defined by (1,2). It turns out that despite the non-conserving dynamics of the present model we can use the same matrices that have been previously used to solve models with conserving dynamics.

We now proceed to outline the matrix product solution for the steady state. The steady state weight of a configuration $`𝒞`$ is represented by the trace of a product of matrices:

$$W_L(𝒞)=\mathrm{Tr}\underset{i=1}{\overset{L}{}}\left[\delta _{\tau _i,+}D+\delta _{\tau _i,}E+\delta _{\tau _i,\mathrm{\hspace{0.17em}0}}A\right],$$
(4)
where $`L`$ is the system size and $`\tau _j=+,,0`$ if site $`j`$ is occupied by a $`+,`$ or $`0`$ respectively. That is, a matrix $`D`$ ($`E`$) represents a positive (negative) particle and $`A`$ represents a vacancy. It is easy to show using the technique of  that $`D,E`$ and $`A`$ should satisfy

$$DEqED=D+E;DA=AE=A;AA=wA$$
(5)
to give the correct steady-state weights. These equations can be satisfied by choosing

$$A=w|VV|,$$
(6)
where $`V|V=1`$, $`D|V=|V`$ and $`V|E=V|`$. Then $`D`$, $`E`$ and $`V|`$ are identical to the matrices and vectors studied in  where they are used to solve the single-species partially asymmetric exclusion process with open boundaries.

We now wish to calculate the normalisation, i.e. the partition function, which is given by the sum of the weights (4) of all accessible configurations:

$$𝒵_L=\mathrm{Tr}\left[(A+D+E)^L(D+E)^L\right]$$
(7)
Using the form (6) of $`A`$ we first write the sum of the weights of configurations with exactly $`M`$ vacancies on a lattice of size $`L`$:

$$𝒵_{L,M}=w^M\underset{\mu =1}{\overset{M}{}}\underset{n_\mu =0}{\overset{\mathrm{}}{}}V|C^{n_\mu }|V\delta _{_\mu n_\mu ,LM}$$
(8)
where $`C=D+E`$. The sum over each $`n_\mu `$ corresponds to the possible number of particles between two consecutive zeros. The delta function enforces the constraint that the total number of particles equals $`LM`$ and the factor $`w^M`$ arises from the $`M`$ zeros. This form neglects the degeneracy in placing a given configuration $`\{n_\mu \}`$ on a ring geometry, which is bounded from above by $`L`$. It is straightforward to check using bounds on the true partition function that this does not affect any of the results present here.

In order to calculate the partition function $`𝒵_L=_{M=1}^L𝒵_{L,M}`$ it is convenient to replace the delta function by a contour integral. Taking the limit of the sum over $`M`$ to infinity yields

$`𝒵_L`$ $`=`$ $`{\displaystyle \underset{M=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{dz}{2\pi i}\frac{w^M}{z^{LM+1}}\underset{\mu =1}{\overset{M}{}}\underset{n_\mu =0}{\overset{\mathrm{}}{}}z^{n_\mu }V|C^{n_\mu }|V}`$ (9)
$`=`$ $`{\displaystyle \frac{dz}{2\pi iz^{L+1}}\left[\frac{zwU(z)}{1zwU(z)}\right]},`$ (10)
where $`U(z)`$ is defined as

$$U(z)=\underset{n=0}{\overset{\mathrm{}}{}}z^nG_n\text{with}G_nV|C^n|V.$$
(11)
The weight $`G_n`$ has been studied before. It is the normalisation sum of the single species partially asymmetric exclusion model on an open lattice of size $`n`$ and with particle injection rate $`1`$ at the left boundary and removal rate $`1`$ at the right boundary .

Below, we shall consider the density of vacancies given by

$$\theta =\underset{L\mathrm{}}{lim}\frac{\overline{M}}{L}=\underset{L\mathrm{}}{lim}\frac{w}{L}\frac{\mathrm{ln}𝒵_L}{w},$$
(12)
where $`\overline{M}`$ is the average number of vacancies in the system. We shall also consider the current $`J_L`$ of positive particles (which is equal to that of the negative particles). We define the current in the same way as for the conserving system although here the density is not conserved. Thus taking into account the inaccessibility of configurations with no vacancies we find

$`J_L`$ $`=`$ $`{\displaystyle \frac{1}{𝒵_L}}\mathrm{Tr}\left[\left(DA+DEqED\right)\left(C+A\right)^{L2}\left(DEqED\right)C^{L2}\right]`$ (13)
$`=`$ $`{\displaystyle \frac{𝒵_{L1}}{𝒵_L}},`$
where we have used the algebraic rules (5) and the form (7) of $`𝒵_L`$.

We now discuss the two distinct cases $`q<1`$ and $`q>1`$ separately:

(i) The case of $`q<1`$

Here the normalisation $`𝒵_L`$ can be evaluated for large $`L`$ from the integral (10) by the saddle point method. This method amounts to working in the grand canonical ensemble. The term in the square brackets of (10) is then just the grand canonical partition function and $`𝒵_L(z^{})^L`$ where $`z^{}`$ is the saddle point value of $`z`$, i.e. the fugacity. Thus, in the thermodynamic limit, we see from (13) that $`J_Lz^{}J`$ which identifies the particle current in the system as the fugacity. Also, we see from (12) that the density of vacancies is given by

$$\theta =w\frac{\mathrm{ln}z^{}}{w}.$$
(14)
Thus one may identify $`\theta `$ as the order parameter and $`\mathrm{ln}z^{}=\mathrm{ln}J`$ as the analogue of the free energy density.

The saddle point equation for the integral (10) may be written as

$$L=\frac{1}{1wzU(z)}\left[\frac{zU^{}(z)}{U(z)}+wzU(z)\right].$$
(15)
For $`L\mathrm{}`$ this equation is satisfied either by $`1wzU(z)𝒪(1/L)`$ or by the term inside the brackets being of order $`L`$. To analyse which of the two scenarios pertains one needs to consider the properties of the increasing function $`U(z)`$. Note from (11) that the convergence of $`U(z)`$ is determined by the large $`n`$ form of $`G_n`$. For $`q<1`$ this quantity is known  to behave for large $`n`$ as

$$G_na\frac{K^n}{n^{3/2}}$$
(16)
where

$$K=\frac{1q}{4}\text{and}a=\frac{4}{\sqrt{\pi }}\left[\underset{i=1}{\overset{\mathrm{}}{}}\frac{(1q^i)^3}{(1+q^i)^4}\right].$$
(17)
From (16) one deduces that $`U(z)`$ converges for $`zK`$ and $`U^{}(z)`$ diverges as $`zK`$. Thus, when $`1wzU(z)=0`$ for some $`zK`$ we solve the saddle-point equation (15) by choosing $`z^{}`$ so that $`1wz^{}U(z^{})𝒪(1/L)`$. On the other hand when $`1wzU(z)>0`$ for all $`zK`$, the divergence must come from the term in the square brackets of equation (15) and we need $`z^{}=K(1𝒪(1/L^2))`$. This can be deduced by noting that $`U^{}(z)|Kz|^{1/2}`$ for $`zK`$. Thus, in the thermodynamic limit $`L\mathrm{}`$, $`z^{}`$ increases to $`K`$ as $`w`$ is decreased to $`w_c`$, given by

$$w_c=\frac{1}{KU(K)}.$$
(18)
Any further decrease in $`w`$ leaves the value of $`z^{}`$ unchanged. Using equation (14) one can see that for $`w>w_c`$, $`\theta `$ decreases as $`w`$ decreases while for $`w<w_c`$, $`\theta =0`$. Using $`U(z)U(K)|Kz|^{1/2}`$ it is easy to show by expanding (15) that $`\theta |ww_c|`$ as $`ww_c`$, recovering (3).

Since $`z^{}`$ is the current of particles and it saturates at $`K=(1q)/4`$ for $`w<w_c`$ we refer to this phase as the maximal current phase. We refer to the phase $`w>w_c`$ as the fluid since the typical configuration is a disordered arrangement of $`+,`$ and $`0`$s.

(ii) The case of $`q>1`$

For $`q>1`$ it is known  that the $`n`$ dependence of $`G_n`$ is given by

$$G_n\left(\frac{q}{q1}\right)^nq^{n^2/4}.$$
(19)
In this case $`U(z)`$ diverges for all $`z>0`$ and we have to impose cutoffs on the sums in equation (9). That is, we write

$$𝒵_L=\underset{M=1}{\overset{L}{}}\frac{dz}{2\pi i}\frac{w^M}{z^{LM+1}}\left(\underset{n=0}{\overset{LM}{}}z^nG_n\right)^M.$$
(20)
From (19) it is clear that the dominant contribution to the sum is when $`M=1`$ and $`n=L1`$. Then $`J_L=𝒵_{L1}/𝒵_L(q1)q^{L/21}`$. Thus, the current of particles is exponentially small in the system size and the density of vacancies is zero. This corresponds to a strongly phase separated state with a single vacancy followed by an extensive block of positive particles (to the right of the vacancy) followed by an extensive block of negative particles (to the left of the vacancy).

Finally, we discuss the transition between $`q<1`$ and $`q>1`$. Using equations (16) and (17) it is easy to show that as $`q1`$ from below, $`KU(K)`$ tends to zero. Thus using equation (18) we see that $`w_c`$ diverges as $`q1`$. Therefore at the $`q=1`$ transition the system changes from a maximal current phase $`J=(1q)/4`$ for $`q<1`$ to a strongly phase separated state with $`J=𝒪(q^{L/2})`$ for $`q>1`$. In both of these phases the density of vacancies is zero. Therefore the transition has no relation to the non-conserving dynamics and is instead related to the reversal of the bias .

It is interesting to make a comparison between the transition from the fluid to the maximal current phase and the denaturation transition in DNA where the two strands of the molecule unbind at a certain temperature. In models of this transition  one assigns a Boltzmann weight $`w`$ to bound base pairs and a weight $`K^n/n^c`$ to unbound segments of length $`n`$. Here $`K`$ is a non-universal constant whereas $`c`$ is a universal exponent depending on the dimension and self-avoidance properties of the unbound DNA. For example, using a random walk model yields $`c=3/2`$ in three dimensions . As temperature is raised and $`w`$ decreases there is an unbinding transition where the fraction of bound base pairs, $`\theta `$, vanishes. The DNA models and the present model can be related by identifying vacancies with bound pairs and blocks of particles with unbound DNA loops. Then $`G_n`$ corresponds to the weight of DNA loop of length $`n`$ and the grand canonical partition functions of the two systems are the same .

The difference between the systems is that in the DNA the exponent $`c=3/2`$ is a result of the fact that the one-dimensional molecule is embedded in three dimensions. However, for the model considered here the exponent $`3/2`$ arises from truly one-dimensional phenomena related to the currents flowing in the system. This can be seen from the relation between $`G_n`$ and the current of the single-species partially asymmetric exclusion process. More explicitly, for such a process the current in a system of size $`n`$ is given by $`G_{n1}/G_n`$. Thus, the $`n`$ dependence of $`G_n`$ is related to the $`n`$ dependence of the current in a single species system of size $`n`$. In the context of the two species model we have studied here such a system corresponds to an uninterrupted block of particles, of length $`n`$, bounded between two vacancies. This picture has recently been used to study the conditions for phase separation in conserving models .

The model can easily be generalised to contain two additional parameters $`\alpha `$ and $`\beta `$ by modifying the rates (1) to read

$$+\mathrm{\hspace{0.33em}0}\stackrel{\beta }{}\mathrm{\hspace{0.17em}0}+0\stackrel{\alpha }{}\mathrm{\hspace{0.17em}0}+\underset{q}{\overset{1}{}}+$$
(21)
and the rates (2) to

$$+\mathrm{\hspace{0.17em}0}\underset{1}{\overset{\beta w}{}}\mathrm{\hspace{0.17em}0\hspace{0.17em}0}0\underset{1}{\overset{\alpha w}{}}\mathrm{\hspace{0.17em}0\hspace{0.17em}0}.$$
(22)
Then one can take $`A=w|VW|`$ where $`\beta D|V=|V`$ and $`\alpha W|E=W|`$. This generalisation should allow for a richer phase diagram than that presented here.

Finally, we point out that the transition from the fluid to the maximal current phase is lost if particles are not conserved inside a block consisting only of particles. For example, a process where $`+0`$ regardless of its neighbours destroys the maximal current phase.

The support of the Israeli Science Foundation is gratefully acknowledged. MRE thanks the Einstein Center for support during a visit to the Weizmann Institute.