format: text

text:
# On the perturbative expansion of the magnetization in the out-of-equilibrium Kondo model.

## I Introduction

In recent years, much experimental and theoretical work has been devoted to exploring the properties of so-called ‘quantum dots’ Aleiner et al. (2002). These are mesoscopic devices in which a ‘dot’ containing a small number of electrons is isolated from two macroscopic leads (denoted ‘left’ and ‘right’) by potential barriers, through which tunneling processes take place. Experimentally, these systems are small devices fabricated using a two-dimensional electron gas Goldhaber-Gordon et al. (1998); Cronenwett et al. (1998); Schmid et al. (1998); van der Wiel et al. (2000) or carbon nanotubes Nygard et al. (2000). The number of electrons on the dot, $`N`$, is controlled by a gate voltage $`V_\mathrm{g}`$. From the experimental point of view, one is primarily interested in the current flowing through the dot as a function of $`V_\mathrm{g}`$ and of the potential difference $`V`$ between the two leads. When tunneling is weak enough, the Coulomb blockade phenomenon appears Aleiner et al. (2002): the conductance through the dot is essentially zero except in the vicinity of certain special values of $`V_\mathrm{g}`$, where the energy difference between the ground states of the dot with $`N`$ and $`N+1`$ electrons vanishes. At these points, conductance peaks are observed.

A simple model for such a system is the Anderson model, where the localized level represents the dot and the hopping term describes its hybridization with the leads. In this paper, we shall restrict ourselves to the regime where the occupation of the dot is not fluctuating, and where naïve application of Coulomb blockade ideas would predict a strongly suppressed conductance. However, when $`N`$ is odd there is exactly one unpaired spin, which is coupled to the leads via a Kondo interaction. In this case, in the linear response regime, spin physics opens up a new transmission channel via the Kondo effect. Although the occupancy of the dot remains fixed, spin-flip interactions permit the formation of strong dot-lead hybridization for temperatures $`T<T_K`$, where $`T_K`$ is the Kondo temperature. As the temperature approaches zero, this leads to unitary limit conductance ($`G=G_02e^2/h`$) via a sharp resonance at the Fermi surface—the Abrikosov-Suhl resonance. This effect was predicted in the context of quantum dots fourteen years ago Glazman and Raikh (1988); Ng and Lee (1988), and was recently observed in a series of experiments Goldhaber-Gordon et al. (1998); Cronenwett et al. (1998); Schmid et al. (1998); van der Wiel et al. (2000); Nygard et al. (2000).

By contrast, in the large-voltage regime, full non-equilibrium calculations are required, and much less is known than close to equilibrium. The out-of-equilibrium Anderson and Kondo models have been studied by various methods in the last decade. Much of this work Meir et al. (1993); Wingreen and Meir (1994); Rosch et al. (2001) has concentrated on the non-crossing approximation (NCA) approach, which adopts a slave-boson description of the problem, and then renders it tractable by neglecting certain vertex corrections. This yields a picture in which the Kondo resonance in the density of states is both split and broadened as the voltage $`V`$ between the leads is increased. Recent work Konik et al.  has explored a new approach, where one attempts to use the Bethe Ansatz results for the Anderson model to construct a Landauer-type picture of transport through the dot. This approach also involves approximation, when one comes to construct the ‘in’ and ‘out’ scattering states from the dressed excitations of the model. Another thread Hershfield et al. (1992); Oguri (2001) has involved studying the Anderson model via perturbation theory in the on-site Coulomb repulsion $`U`$. While these works provide approximate information on the behavior of the current-voltage characteristic, they shed little light on the nature of the many-body state of the system when $`V>T_K`$.

In particular, a basic question recently debated Coleman et al. (2001); Rosch et al. (2001) is whether the Kondo problem has a strong coupling regime at low temperature and high voltage. In discussing this point, a previous paper Coleman et al. (2001) used a second order perturbative expression for the magnetic susceptibility in the out-of-equilibrium steady state induced by $`V`$. This putative result was however incorrect, even at order 0 in $`J`$, the strength of the Kondo couplings; in this note, we correct this result at order 0, discuss the related physics and present a systematic method for calculating higher-order corrections in $`J`$. Our main conclusion is that, even at order 0 in $`J`$, the Keldysh function of the spin deviates from its equilibrium value and therefore the steady state magnetization of the dot is not given by $`M_{\mathrm{eq}}\frac{1}{2}\mathrm{tanh}(B/2T)`$. Rather, it must be computed by solving a transport equation. We emphasize that this issue is not directly related to the so called “decoherence time” Rosch et al. (2001) but is a basic point about perturbation theory to be addressed before discussing out-of-equilibrium RG equations and the existence of a strong coupling regime at large voltage. In particular, it has ramifications for other physical quantities, e.g. the current and the current-current correlation function in a magnetic field (see section III.3.2).

The paper is organized as follows. In Section II, we present the model, our main result (the magnetization at order 0 given by Eq. 4), and the associated physical discussion. In Section III, we present a detailed pedagogical derivation of (4) using the Keldysh method. Finally, in Appendix D, we present more details about the computation presented in ref. Coleman et al. (2001), and explain why it was incorrect.

## II Results and discussion

Our starting point is the Hamiltonian of the two-lead Kondo model. For a discussion of the modeling of the quantum dot, and for a derivation of this model from the Anderson model via an out-of-equilibrium Schrieffer-Wolff transformation, we refer the reader to the literature, in particular Kaminski et al. (2000) and references therein. The Hamiltonian is:

$`H=`$ $`{\displaystyle \underset{\alpha 𝐤\sigma }{}}\left(\epsilon _𝐤\mu _\alpha \right)c_{\alpha 𝐤\sigma }^{}c_{\alpha 𝐤\sigma }+H_{\mathrm{refl}}+H_{\mathrm{trans}}B_sS_zB_c{\displaystyle \underset{\alpha 𝐤\sigma }{}}\sigma c_{\alpha 𝐤\sigma }^{}c_{\alpha 𝐤\sigma },`$ (1)
$`H_{\mathrm{refl}}=`$ $`J_R{\displaystyle \underset{𝐤,𝐤^{},\sigma ,\sigma ^{}}{}}\left(c_{R𝐤\sigma }^{}\stackrel{}{\sigma }_{\sigma \sigma ^{}}c_{R𝐤^{}\sigma ^{}}\right)\stackrel{}{S}+(RL),`$
$`H_{\mathrm{trans}}=`$ $`J_{RL}{\displaystyle \underset{𝐤,𝐤^{},\sigma ,\sigma ^{}}{}}\left(c_{R𝐤\sigma }^{}\stackrel{}{\sigma }_{\sigma \sigma ^{}}c_{L𝐤^{}\sigma ^{}}\right)\stackrel{}{S}+(RL),`$
where $`c_{\alpha 𝐤\sigma }^{}`$ creates an electron in lead $`\alpha \{L,R\}`$ with momentum $`𝐤`$ and spin $`\sigma `$, and $`J_L`$, $`J_R`$ and $`J_{LR}=(J_{RL})^{}`$ are Kondo coupling constants between the electrons and the spin of the dot $`\stackrel{}{S}`$. The first term in $`H`$ describes the electrons in the leads, with $`\epsilon _𝐤`$ being the bare energy of an electron of momentum $`𝐤`$ at zero voltage (the same for each lead) and $`\mu _\alpha `$ is the potential in the lead $`\alpha `$. Each lead consists of a free electron gas with a density of states $`\rho (ϵ)`$ of bandwidth $`D`$: ultimately, we will be interested in the result in the large bandwidth ($`D\mathrm{}`$) limit, but the computations are first performed for finite $`D`$. We will make the physical assumption that the leads are in thermal equilibrium at a temperature $`T`$. The voltage is applied by taking the chemical potentials of the two leads to be different, $`\mu _L\mu _R=V`$. $`H_{\mathrm{refl}}`$ describes regular Kondo processes, where an electron from a given lead is spin-flip scattered back into the same lead; $`H_{\mathrm{trans}}`$ describes ‘spin-flip cotunneling’, where an electron from one lead is spin-flip scattered into the other lead. If the model (1) is derived from the Anderson model, one finds that $`J_RJ_L=|J_{RL}|^2`$. In this paper, however, we relax this relation between the coupling constants and treat them as independent parameters. Finally, the last terms represent the coupling to the magnetic field. We allow two different magnetic fields, $`B_s`$ for the spin and $`B_c`$ for the lead electrons; this permits the calculation of the total and local spin susceptibilities within the same computation.

In this paper, we shall be interested in the values taken by the following quantities in the non-equilibrium steady state induced by the voltage $`V`$: the dot magnetization $`M_{\mathrm{dot}}`$, the magnetization of the leads $`M_{\mathrm{leads}}`$, the total magnetization $`M_{\mathrm{tot}}`$, and the total and local susceptibilities. These are given by:

$`M_{\mathrm{dot}}(B_s,B_c)`$ $`=S_z,`$ (2a)
$`M_{\mathrm{leads}}(B_s,B_c)`$ $`={\displaystyle \underset{\alpha 𝐤\sigma }{}}\sigma c_{\alpha 𝐤\sigma }^{}c_{\alpha 𝐤\sigma },`$ (2b)
$`M_{\mathrm{tot}}(B)`$ $`=M_{\mathrm{dot}}(B,B)+M_{\mathrm{leads}}(B,B)M_{\mathrm{Pauli}},`$ (2c)
$`\chi _{\mathrm{tot}}`$ $`={\displaystyle \frac{M_{\mathrm{tot}}(B)}{B}}|_{B=0},`$ (2d)
$`\chi _{\mathrm{loc}}`$ $`={\displaystyle \frac{M_{\mathrm{dot}}(B_s,B_c=0)}{B_s}}|_{B_s=0},`$ (2e)
where angle-brackets $`\mathrm{}`$ denote an expectation value taken in the steady (i.e. long-time) state of the system. $`M_{\mathrm{Pauli}}`$ is simply the Pauli paramagnetic contribution from the lead electrons which would be present even in the absence of the impurity, and which we therefore exclude from $`M_{\mathrm{tot}}`$. We consider the perturbative expansions of these steady state quantities; more precisely, we define:

$$J_R\theta _RJJ_L\theta _LJJ_{RL}\theta _{RL}J$$
(3)
and we let $`J`$ go to zero while keeping the coefficients $`\theta _R`$, $`\theta _L`$ and $`\theta _{RL}`$ fixed. In the following, the expression “order $`n`$” refers to the order $`n`$ of this expansion in $`J`$.

Our main result is the order $`0`$ term of the perturbative expansion of the magnetization:

$$M_{\mathrm{tot}}(B)=M_{\mathrm{dot}}(B,B_c)+O(J)=\frac{1}{2}\mathrm{tanh}\left(\frac{B}{2T}\right)\frac{\phi \left(\frac{B}{T}\right)\left(1+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}\right)}{\frac{1}{2}\left[\phi \left(\frac{B+V}{T}\right)+\phi \left(\frac{BV}{T}\right)\right]+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}\phi \left(\frac{B}{T}\right)}+O(J),$$
(4)
where $`\phi `$ is defined by:

$$\phi (x)\frac{x}{\mathrm{tanh}\left(\frac{x}{2}\right)}.$$
(5)
(At this order, $`M_{\mathrm{dot}}(B,B_c)`$ does not depend on $`B_c`$.) As a result, the magnetic susceptibility at order 0 is:

$$\chi (T,V)=\frac{1}{4T}\frac{1+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}}{\frac{1}{2}\phi \left(\frac{V}{T}\right)+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}}+O(J),$$
(6)
and in particular for $`|V|\mathrm{}`$

$$\chi (T,V)\frac{1+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}}{2|V|}.$$
(7)
A striking feature of (4) is that the magnetization at order 0 in $`J`$ is not given by the equilibrium expression $`M_{\mathrm{eq}}=\frac{1}{2}\mathrm{tanh}(B/2T)`$. This may seem surprising: since the couplings to the leads are relaxed to zero, why should we not find $`M_{\mathrm{eq}}`$, the magnetization of a free spin? Physically, the answer is that $`M_{\mathrm{eq}}`$ is not the “magnetization of a free spin”, but rather the magnetization of a spin weakly coupled to an equilibrium bath. On the other hand, Eq. (4) gives the magnetization of a spin weakly coupled, via the Hamiltonian (1), to two leads with different chemical potentials. It is simply an out-of-equilibrium extension of the Curie law. In particular, we emphasize that the finite susceptibility at $`T=0`$ cannot be interpreted as a renormalization (or a screening) of the spin due to the voltage.

At long times, the state of the spin is completely determined by the properties of these leads, and hence so is its distribution function, which describes the population of its two states as a function of temperature, voltage, and magnetic field. Since the whole system is not in equilibrium, this steady state is not described by the Gibbs distribution; in particular, the fluctuation-dissipation theorem (FDT) need not hold. Hence the magnetization need not be (and is not) $`M_{\mathrm{eq}}`$ at order 0 in $`J`$. Rather it should be computed by solving a transport equation in the steady state regime, i.e. a quantum Boltzmann equation. (At dominant order this procedure is equivalent to using a semi-classical master equation; see Appendix A.)

The crucial point is that $`J_R,J_L,J_{RL}`$ are relaxed to zero, assuming that they are still bigger than the coupling of the spin to any other thermal bath. If we were to take into account such a coupling (denoted by $`J_0`$) then the result would cross over to $`M_{\mathrm{eq}}`$ when $`J_R,J_L,J_{RL}J_0`$ (all the couplings going to zero while maintaining fixed ratios). In fact, (4) implies that the equilibrium value is only an upper bound: $`0M_{\mathrm{tot}}(J=0)/M_{\mathrm{eq}}1`$, which follows from the convexity of $`\phi `$. This bound is saturated only in equilibrium, i.e. for $`\theta _{RL}0`$ or $`V0`$.

Moreover, as expected on physical grounds, the non-equilibrium result is much less universal than the equilibrium one. In equilibrium, the magnetization of a spin weakly coupled to a bath depends neither on the properties of the bath except the temperature $`T`$ nor on the form of the couplings between the spin and the bath. Neither of these statements holds true for the non-equilibrium Kondo model: $`M_{\text{dot}}`$ (at order 0) depends not only on $`T`$ but also on the voltage $`V`$ and on a ratio involving the parameters $`\theta _R`$, $`\theta _L`$, and $`\theta _{RL}`$. Thus $`M_{\text{dot}}`$ is perturbative in $`J`$, but not in the three variables $`J_R,J_L,J_{RL}`$, and in particular it is not analytic in $`J_R,J_L,J_{RL}`$ around $`(0,0,0)`$. One can find such a dependence on the ratio of couplings even in a simple free model (a single level coupled to two leads), as illustrated in Appendix B. Furthermore, it should be remembered that we computed $`M_{\text{dot}}`$ with free leads: introducing interactions in the leads would change the function $`\phi `$ (the important quantity being the electron-electron bubble in the leads). For example, we expect a different result to hold for a spin coupled to Luttinger liquids, even at order 0.

It should be noted that the result (4) gives a non-trivial expression for the magnetization even at zero temperature. Since $`\varphi (x)|x|`$ as $`x\mathrm{}`$, we find that

$$M_{\mathrm{tot}}(B,V)|_{T=0}=\frac{B}{2}\left(\frac{1+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}}{\frac{1}{2}\left(|B+V|+|BV|\right)+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}|B|}\right)+O(J).$$
(8)
We then have two cases (we can take $`B,V>0`$): for $`B>V>0`$, we obtain the equilibrium result $`M_{\mathrm{tot}}(B,V)=\text{sgn}(B)/2+O(J)`$, but for $`0<B<V`$, the magnetization is a still a function of $`B/V`$:

$$M_{\mathrm{tot}}|_{T=0}\frac{B}{2}\left(\frac{1+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}}{|V|+\frac{\theta _R^2+\theta _L^2}{2\theta _{RL}^2}|B|}\right),$$
(9)
in agreement with (7) in the limit $`|V||B|`$. Physically, there are two sources of energy available to flip the spin: the thermal fluctuations of both baths (represented by $`T`$) and the fact that an electron can go from L to R and give to the spin an energy of order $`V`$. If we decrease $`B`$ from high values ($`BT,V`$), the spin is locked until $`B`$ reaches the largest of these energy scales. Thus the magnetization at zero temperature is expected to saturate only for $`B>V`$. Similarly, the susceptibility is in general expected to behave as $`\chi 1/E`$, where $`E`$ is the largest energy available to flip the spin. The fact that our result is still non-trivial at zero temperature implies that it could be seen in numerical computations, such as the DMRG approach of Cazalilla and Marston Cazalilla and Marston .

## III Derivation of Eq. (4)

The purpose of this section is to present the derivation of Eq. (4) and more generally the procedure for obtaining the perturbative expansion of physical quantities in the steady state. It is organized as follows: in paragraph III.1, we first give an overview of the derivation; full details are given in the following sections (III.2 and III.3) and in the appendices, including a presentation of the Keldysh method.

### III.1 Overview

When doing perturbation theory in the steady state, there are two important small couplings: $`J`$, the strength of the Kondo couplings, and $`s`$, the small regulator that appears in Green’s functions such as the “spin” retarded function (defined in section III.2.3):

$$R_B^0(\omega )=\frac{1}{\omega B+is}.$$
(10)
The scale $`s`$ should be thought of as being due to the coupling to an auxiliary thermal bath, which in the physical system would be the substrate. If the impurity is coupled to such a bath whose density of states is $`\rho _2`$ by a coupling $`g`$, one finds $`sg^2\rho _2`$ (in the large bandwidth limit). In the physical quantum dot system, $`Js`$, meaning that the correct order of limits to take is $`s0`$ followed by $`J0`$, as pointed out in section II.

However, ‘straightforward’ perturbation theory in $`J`$ takes the limit in the opposite order: one first expands in $`J`$ while keeping $`s`$ finite, and only then takes $`s0`$ term by term in the perturbation series. In equilibrium, these two limits commute, but out of equilibrium they do not. This is explicitly shown in section III.3.1; the out-of-equilibrium Keldysh Green’s function is not analytic around $`(s,J)=(0,0)`$. This non-analytic behaviour is quite generic in situations where an impurity is coupled to several leads, and may be seen even in a simple free model (see Appendix B). The signature of that non-commutativity is that ‘straightforward’ perturbation theory fails: its terms (starting at $`O(J^2)`$) exhibit divergences of a $`1/s`$ form. These are similar to the infrared divergences in equilibrium perturbation theory that signal an incorrect choice of reference state. As shown below, $`1/s`$ divergences at order $`J^2`$ signal an incorrect choice for the Keldysh Green’s function at order $`0`$, or equivalently of the distribution function which describes the non-thermal population of the two levels of the spin. Consequently, contrary to the claims of Coleman et al. (2001), there is no possibility of regulating these $`1/s`$ divergences order by order in $`J`$. We discuss the putative regulation procedure of Coleman et al. (2001) in Appendix D, and explain why it is incorrect.

The solution is to begin with the Dyson equation in the steady state, considered as a functional equation for the full Green’s function $`G`$, using the skeleton self-energy diagrams (see Appendix E for further details). After taking the $`s0`$ limit in this equation, a perturbative expansion for $`G`$ may be inserted into it, and a solution obtained order by order in $`J`$. Solving the Keldysh component of this Dyson equation is equivalent to solving the quantum Boltzmann equation; as shown below, the solution is non-thermal even at order $`0`$ in $`J`$. Another procedure for solving the problem of $`1/s`$ divergences was proposed in Wingreen and Sivan (1996) in the context of the $`U=\mathrm{}`$ out-of-equilibrium Anderson model: its method was to choose the zeroth-order spin Keldysh function (see section III.2.3) in order precisely to cancel the divergences at order $`J^2`$. This is completely equivalent to our approach, as explained in section III.3.2.

Finally, we note that one may give a simple semiclassical derivation of the out-of-equilibrium result (4) based on a master equation. This derivation is due to L. Glazman and A. Kaminski Glazman and Kaminski , and is presented in Appendix A. The success of such a semiclassical approach (at this lowest order) is related to the fact that one can compute the equilibrium magnetization of a free quantum spin using a classical Ising model. In fact, at this order, the semi-classical master equation is strictly equivalent to the Keldysh component of the steady-state Dyson equation, so the two apparently disparate derivations yield the same result.

### III.2 Technical preliminaries: the Keldysh method

Let us now turn to the technical details. In the following, for simplicity, we will write some equations for a generic fermionic field $`\psi `$, which will be specialized afterwards to the fields representing the electrons and the spin.

#### III.2.1 Generalities

The basic idea of the non-equilibrium Keldysh method Schwinger (1961); Keldysh (1965); Rammer and Smith (1986) consists in taking the system at an initial time $`t=0`$ in an initial state described by a density matrix $`\rho _0`$ and letting the system relax, using the Hamiltonian evolution given by $`H`$, to a long-time regime. In order to ensure that the system relaxes, it may be necessary to add some additional coupling terms to the Hamiltonian, in particular to break conservation laws: see the discussion in Appendix E. Depending on the system, the long-time regime can be an equilibrium state, a non-equilibrium steady state, a non-time-translation-invariant steady state, or even an aging regime (in glassy systems). In the quantum dot problem, we assume that the system reaches at finite voltage $`V`$ a non-equilibrium steady state, in which we want to compute physical quantities.

A “Keldysh” average of any quantity $`A`$ is defined by

$$A(t)_𝒦e^{iHt}Ae^{iHt}_0\text{Tr}\left(\rho _0e^{iHt}Ae^{iHt}\right),$$
(11)
where $`\mathrm{}_0`$ is the average taken using the initial density matrix of the system. The steady-state average is given by:

$$A=\underset{t\mathrm{}}{lim}A(t)_𝒦.$$

Using the usual representation of the evolution operator $`e^{iHt}`$ as a $`T`$-ordered exponential in the interaction picture (and the anti-$`T`$-ordered one for $`e^{iHt}`$), one can obtain an expansion in the coupling constant $`J`$. Following the usual conventions, it is convenient to keep track of the two exponentials using a closed time contour, running from 0 to $`+\mathrm{}`$ and back to 0 Keldysh (1965): we denote by $`+`$ the upper contour (from 0 to $`\mathrm{}`$), which arises from expanding $`e^{iHt}`$, and by $``$ the lower contour. In accordance with this notation, we define the four Green’s functions:

$`G_\psi ^{++}(t,t^{})`$ $`iT\psi (t)\psi ^{}(t^{})_𝒦,`$ $`G_\psi ^+(t,t^{})`$ $`i\psi (t)\psi ^{}(t^{})_𝒦,`$
$`G_\psi ^{}(t,t^{})`$ $`i\stackrel{~}{T}\psi (t)\psi ^{}(t^{})_𝒦,`$ $`G_\psi ^+(t,t^{})`$ $`i\psi ^{}(t^{})\psi (t)_𝒦.`$ (12)
Here $`T`$ is the time-ordering operator, $`\stackrel{~}{T}`$ the anti-time-ordering operator, and $`\psi `$ is any fermionic field. The two indices of the matrix $`G`$ will be called “indices in Keldysh space”. In the following, the equations (III.2.1) will be summarized with the notation:

$$G_\psi (t,t^{})=i\psi (t)\psi ^{}(t^{})_𝒦$$
(13)
and $`G`$ will always denote a $`2\times 2`$ Keldysh matrix.

The Green’s functions (III.2.1) are not independent but can be expressed as functions of the retarded, advanced and Keldysh Green’s functions defined respectively by <sup>1</sup><sup>1</sup>1Throughout this paper, $`A`$ is the advanced function and should not be confused with a spectral function.:

$`R_\psi (t,t^{})`$ $`i\theta (tt^{})\{\psi (t),\psi ^{}(t^{})\}_𝒦,`$
$`A_\psi (t,t^{})`$ $`i\theta (t^{}t)\{\psi (t),\psi ^{}(t^{})\}_𝒦,`$
$`K_\psi (t,t^{})`$ $`i[\psi (t),\psi ^{}(t^{})]_𝒦,`$ (14)
as is shown from the transformation Rammer and Smith (1986):

$$G_\psi \left(\begin{array}{cc}G_\psi ^{++}& G_\psi ^+\\ G_\psi ^+& G_\psi ^{}\end{array}\right)\stackrel{~}{G}_\psi \left(\begin{array}{cc}R_\psi & K_\psi \\ 0& A_\psi \end{array}\right)\stackrel{~}{G}_\psi =\frac{1}{2}\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right)G_\psi \left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).$$
(15)
We will denote the first set of Green’s functions (III.2.1) “the $`\pm `$ basis”, and the second (III.2.1) the “Larkin-Ovchinnikov (LO) basis”.

There are further relations between these Green’s functions. In general, $`R_\psi (t,t^{})=\left(A_\psi (t^{},t)\right)^{}`$. Moreover, in equilibrium FDT reads:

$$K_\psi (\omega )=h_{\mathrm{eq}}(\omega )\left(A_\psi (\omega )R_\psi (\omega )\right)h_{\mathrm{eq}}(\omega )\mathrm{tanh}\left(\frac{\omega }{2T}\right)$$
(16)
thus the retarded Green’s function is the only remaining independent Green’s function, and so contains all the information about the state of the system. Out of equilibrium, however, $`K_\psi (\omega )`$ and $`R_\psi (\omega )`$ should be independently determined.

In order to write down the diagrammatic expansion, we also define the “bare” counterparts of the Green’s functions in (III.2.1, III.2.1) by:

$$G_{\psi 0}(t,t^{})i\psi _{\mathrm{int}}(t)\psi _{\mathrm{int}}^{}(t^{})_0,$$

where the field $`\psi _{\mathrm{int}}`$ is in the interaction picture. From the expansions of the evolution operators, we obtain the diagrammatic expansion, provided that $`\mathrm{}_0`$ satisfies Wick’s theorem. For the problem at hand, we take as the initial condition the density matrix of the model with $`J_R=J_L=J_{RL}=0`$ and for the spin a finite $`s`$ (arising from the coupling of the spin to a thermal bath; see section III.1 above). Going to the long time limit, every function becomes a function of the difference of the times, and transforming to Fourier space we obtain standard Feynman rules.

#### III.2.2 Lead electrons

Let us begin with the lead electrons. We care here only about the local Green’s functions, so we will drop the spatial indices. We denote by $`(𝐆_𝐜)_{\alpha \sigma ,\beta \sigma ^{}}(t,t^{})`$ the Green’s function describing the creation in lead $`\beta `$ of an electron with spin $`\sigma ^{}`$ at time $`t^{}`$, and a corresponding annihilation in lead $`\alpha `$ of an electron with spin $`\sigma `$ at time $`t`$:

$$(𝐆_𝐜)_{\alpha \sigma ,\beta \sigma ^{}}(t,t^{})=ic_{\alpha \sigma }(t)c_{\beta \sigma ^{}}^{}(t^{})_𝒦.$$
(17)
(Here, the indices $`\alpha ,\beta \{L,R\}`$, while $`\sigma ,\sigma ^{}=+,`$, and the local electron operators are defined by $`c_{\alpha \sigma }_𝐤c_{\alpha 𝐤\sigma }`$.) The Green’s function $`𝐆_𝐜`$ is thus a priori a $`4\times 4`$ matrix (in lead-spin space), whose entries are themselves $`2\times 2`$ Keldysh matrices. In the following, a bold notation with subscript $`c`$ always designates such a $`4\times 4`$ matrix. The bare Green’s function $`𝐆_{\mathrm{𝐜𝟎}}`$ is however diagonal. The bare (diagonal) density of states is :

$$(𝝆_{𝒄\mathrm{𝟎}})_{\alpha \sigma ,\beta \sigma ^{}}=\rho (\omega +\sigma B_c\mu _\alpha )\delta _{\alpha \beta }\delta _{\sigma \sigma ^{}}$$

(the energy levels are the same is both leads, but shifted by the Zeeman energy and the voltage). We assume that the baths, being much bigger than the impurity, are permanently in thermal equilibrium so that the bare electrons’ functions read:

$`𝐑_{\mathrm{𝐜𝟎}}(\omega )`$ $`={\displaystyle 𝑑ϵ\frac{𝝆_{𝒄\mathrm{𝟎}}(ϵ)}{\omega ϵ+i0^+}},`$ (18a)
$`𝐊_{\mathrm{𝐜𝟎}}(\omega )`$ $`=2i\pi 𝐡_{\mathrm{𝐜𝟎}}(\omega )𝝆_{𝒄\mathrm{𝟎}}(\omega ),`$ (18b)
$`(𝐡_{\mathrm{𝐜𝟎}})_{\alpha \sigma ,\beta \sigma ^{}}`$ $`h_{\mathrm{eq}}\left(\omega \mu _\alpha \right)\delta _{\alpha \beta }\delta _{\sigma \sigma ^{}},`$ (18c)
where the first two equations are matricial, and $`\mu _\alpha `$ is the potential of the lead $`\alpha =R,L`$. The voltage difference between the leads is given by $`V=\mu _L\mu _R`$.

#### III.2.3 Spin

Since the spin operator is not appropriate for diagrammatic computations (it does not satisfy Wick’s theorem), we represent the spin 1/2 by three Majorana fermions $`\eta ^a,a\{x,y,z\}`$ which satisfy Wick’s theorem and the relations

$`S^a`$ $`={\displaystyle \frac{i}{2}}ϵ^{abc}\eta ^b\eta ^c,`$ (19a)
$`\left(\eta ^a\right)^{}`$ $`=\eta ^a,`$ (19b)
$`\{\eta ^a,\eta ^b\}`$ $`=\delta ^{ab}.`$ (19c)
Using (19), one can easily show that $`\stackrel{}{S}`$ satisfies the correct commutation relations and that $`\stackrel{}{S}^2=3/4`$. Note that this last constraint is automatically satisfied, unlike in the case of a Dirac fermion representation, where a Lagrange multiplier would have been required to fix the magnitude of the spin: the Majorana representation therefore makes the computation simpler. In this paper, we consider only spin 1/2, but our computations could be extended to higher spin, provided that one used another representation for $`\stackrel{}{S}`$.

Let us now discuss the propagators $`G_{ab}`$ ($`a,b\{x,y,z\}`$) of the $`\eta `$, in the presence of a magnetic field $`B_s`$ along the $`z`$ direction. The general form is:

$$𝐆=\left(\begin{array}{ccc}G_{xx}& G_{xy}& 0\\ G_{yx}& G_{yy}& 0\\ 0& 0& G_{zz}\end{array}\right),G_{ab}(t,t^{})i\eta ^a(t)\eta ^b(t^{})_𝒦,$$
(20)
where the elements are $`2\times 2`$ Keldysh matrices. Indeed $`G_{xz}=G_{yz}=0`$ by symmetry. To prove this, note that the Hamiltonian is invariant under a $`\pi `$ rotation around the $`z`$ axis, which is implemented by $`(\eta ^x,\eta ^y,\eta ^z)(\eta ^x,\eta ^y,\eta ^z)`$ (and the corresponding rotation for the $`c`$ electrons). Furthermore, making a $`\pi /2`$ rotation around the $`z`$ axis, implemented by $`(\eta ^x,\eta ^y,\eta ^z)(\eta ^y,\eta ^x,\eta ^z)`$, we obtain the relations

$$G_{xy}=G_{yx},G_{xx}=G_{yy}.$$
(21)
The Hamiltonian is also invariant under a $`\pi `$ rotation around the $`x`$ axis ($`(\eta ^x,\eta ^y,\eta ^z)(\eta ^x,\eta ^y,\eta ^z)`$) together with a change of sign of the magnetic fields $`B_s`$ and $`B_c`$, and hence $`G_{xy}`$ and $`G_{xx}`$ are respectively odd and even in the magnetic field. In particular, for $`B_s=B_c=0`$, $`G_{xy}=0`$. The foregoing arguments apply to the full propagator and to the free propagator $`G^0`$ (computed with only the magnetic field $`B_s`$).

A different basis is also useful: defining the Dirac fermion $`f`$ and its Green’s function $`G_B`$ by

$`f`$ $`{\displaystyle \frac{\eta ^xi\eta ^y}{\sqrt{2}}},`$ (22)
$`G_B(t)`$ $`=if(t)f^{}(0)_𝒦,`$ (23)
we have the relations:

$`G_{xx}(\omega )`$ $`=G_{yy}(\omega )={\displaystyle \frac{1}{2}}\left(G_B(\omega )+G_B(\omega )\right),`$ (24a)
$`G_{xy}(\omega )`$ $`=G_{yx}(\omega )={\displaystyle \frac{i}{2}}\left(G_B(\omega )G_B(\omega )\right).`$ (24b)
In the $`(f,\eta ^z)`$ basis, the propagator is diagonal, so it is more convenient for example to write the Dyson equation, whereas the original basis $`(\eta ^x,\eta ^y,\eta ^z)`$ is more convenient for the diagrammatics. The bare propagators in the $`(f,\eta ^z)`$ basis are given by:

$`R_z^0(\omega )`$ $`={\displaystyle \frac{1}{\omega +is}},`$ (25a)
$`K_z^0(\omega )`$ $`=2ih_{\mathrm{eq}}(\omega ){\displaystyle \frac{s}{\omega ^2+s^2}},`$ (25b)
$`R_B^0(\omega )`$ $`={\displaystyle \frac{1}{\omega B+is}},`$ (25c)
$`K_B^0(\omega )`$ $`=2ih_{\mathrm{eq}}(\omega ){\displaystyle \frac{s}{(\omega B)^2+s^2}},`$ (25d)
where $`s`$ is a small regulator which, as discussed in section III.1 above, should be thought of as the width due to coupling to an auxiliary thermal bath.

#### III.2.4 Vertex factors and Dyson equations

The vertex factors can be extracted simply from the Hamiltonian, and the Feynman rules are summarized in Fig. 1 (in the $`\pm `$ basis). Note that we have oriented the Majorana fermion lines, despite the operator property that $`(\eta ^{})=\eta `$; a proof that the lines are orientable is given in Appendix C.

We now derive the expression for the Dyson equation describing the $`\eta `$ fields. In general a Dyson equation reads

$$𝐆^1(\omega )=𝐆_{\mathrm{𝟎}}^{}{}_{}{}^{1}(\omega )𝚺(\omega ),$$
(26)
where the inversion has to be taken in the tensor product of the $`x,y,z`$ space and the Keldysh space (thus with $`6\times 6`$ matrices). The free propagator is given by (25). In the $`x,y,z`$ basis, the self energy can be written as:

$$𝚺(\omega )=\left(\begin{array}{ccc}\mathrm{\Sigma }_d(\omega )& i\mathrm{\Sigma }_a(\omega )& 0\\ i\mathrm{\Sigma }_a(\omega )& \mathrm{\Sigma }_d(\omega )& 0\\ 0& 0& \mathrm{\Sigma }_z(\omega )\end{array}\right)\mathrm{\Sigma }_\mu (\omega )=\left(\begin{array}{cc}\mathrm{\Sigma }_\mu ^R(\omega )& \mathrm{\Sigma }_\mu ^K(\omega )\\ 0& \mathrm{\Sigma }_\mu ^A(\omega )\end{array}\right)\mu \{a,d,z\}.$$
(27)
Using the definitions (24), we find:

$$G_{\pm B}=\left(\begin{array}{cc}\frac{1}{\omega B\mathrm{\Sigma }_{\pm B}^R(\omega )+is}& \frac{\mathrm{\Sigma }_{\pm B}^K(\omega )+2ish_{\mathrm{eq}}(\omega )}{\left|\omega B\mathrm{\Sigma }_{\pm B}^R(\omega )+is\right|^2}\\ \\ 0& \frac{1}{\omega B\mathrm{\Sigma }_{\pm B}^A(\omega )is}\end{array}\right)\mathrm{\Sigma }_{\pm B}(\omega )\mathrm{\Sigma }_d(\omega )\pm \mathrm{\Sigma }_a(\omega ).$$
(28)
### III.3 Perturbation expansion in the steady state

After these preliminaries, let us now come back to the perturbative expansion in the steady state and develop on the points introduced in III.1. In paragraph III.3.1 we show that the limits $`J0`$ and $`s0`$ do not commute; in paragraph III.3.2 we solve that problem and present the derivation of Eq. (4) together with an algorithm for computing higher order terms in the perturbative expansion.

#### III.3.1 Non-commutativity of limits

We are interested in the steady-state values of the magnetizations (2a) and (2b), which can be expressed in the following way:

$`M_{\text{dot}}`$ $`={\displaystyle \frac{i}{4}}{\displaystyle \frac{d\omega }{2\pi }\left(K_B(\omega )K_B(\omega )\right)},`$ (29a)
$`M_{\text{leads}}`$ $`={\displaystyle \frac{i}{2}}{\displaystyle \frac{d\omega }{2\pi }\underset{\genfrac{}{}{0.0pt}{}{\alpha =R,L}{u=+,}}{}(1)^u(𝐊_𝐜)_{\alpha u,\alpha u}},`$ (29b)
Indeed, the dot magnetization is given by:

$$M_{\text{dot}}=S^z=\frac{i}{2}[\eta ^x,\eta ^y]=\frac{1}{2}K_{xy}(t=0).$$

$`M_{\text{leads}}`$ is derived analogously.

We begin with the expression for the dot magnetization derived from (29a) using the Dyson equation (28):

$$M_{\text{dot}}=\frac{1}{4\pi }𝑑\omega \text{Im}\left(\frac{h_{\mathrm{eq}}(\omega )}{\omega B\mathrm{\Sigma }_B^R(\omega )+is}\right)\frac{d\omega }{8i\pi }\frac{\mathrm{\Sigma }_B^K(\omega )h_{\mathrm{eq}}(\omega )[\mathrm{\Sigma }_B^A(\omega )\mathrm{\Sigma }_B^R(\omega )]}{\left|\omega B\mathrm{\Sigma }_B^R(\omega )+is\right|^2}(BB),$$
(30)
where the last term denotes an antisymmetrization in $`B`$. The first term in (30) can be transformed into a sum over Matsubara frequencies and thus we see that the limits $`J0`$ and $`s0`$ commute in that term. However, the second term is more interesting. Firstly, it vanishes in equilibrium since the numerator cancels, as required by the FDT. Moreover, using (for $`ϵ0`$)

$$\frac{1}{\left|\omega B+iϵ\right|^2}\frac{1}{ϵ}\delta \left(\omega B\right),$$
(31)
we see that in that second term these limits do not commute:

* $`J0`$ and then $`s0`$. We first take $`\mathrm{\Sigma }_B^R=0`$ in the denominator. We find that this second term in (30) gives a diverging term at second order in $`J`$ proportional to $`J_{RL}^2/s`$. More precisely, using the expression for the self-energy at second order derived in Appendix F, and (31) with $`ϵ=s`$, we find

$$M_{\text{dot}}\frac{J_{RL}^2}{s}\mathrm{tanh}\left(\frac{B}{2T}\right)[\phi \left(\frac{B}{T}\right)\frac{1}{2}(\phi \left(\frac{B+V}{T}\right)+\phi \left(\frac{BV}{T}\right)\left)\right]$$

up to a finite (i.e. not diverging as $`s0`$) term of order $`J^2`$ and to $`O(J^3)`$ terms ($`\phi `$ is defined in (5)).
* $`s0`$ and then $`J0`$. Using (31) with $`ϵ=\text{Im}\mathrm{\Sigma }_B^R(B)`$ (we work at dominant order in $`J`$), we find:

$$M_{\text{dot}}=\frac{1}{4}\left(\frac{\mathrm{\Sigma }_B^K(B)}{\mathrm{\Sigma }_B^A(B)\mathrm{\Sigma }_B^R(B)}\right)(BB)+O(J)$$
(32)
where the self-energy has to be expanded at order 2 in $`J`$ (the first term in (30) cancels a part of the second term).

The second limit is the physical one (as explained in section II) and it gives a formula (32) for the magnetization at order 0. However, to make use of this formula one needs to know the Keldysh ($`\mathrm{\Sigma }_B^K(B)`$) and spectral ($`\mathrm{\Sigma }_B^A(B)\mathrm{\Sigma }_B^R(B)`$) parts of the self-energy. The leading terms of these self-energies are of order $`2`$ in $`J`$, but because the ratio of them is taken, they determine the magnetization at order $`0`$. The crucial point is that (32) is in fact an implicit equation for $`h_B(B)`$, which appears on the left-hand side because $`M_{\mathrm{dot}}`$ is defined in terms of it by (29), and on the right-hand side because the self-energies at $`O(J^2)`$ depend on the Keldysh Green’s function at $`O(1)`$, and hence on $`h_B`$. The most elegant way to capture this ‘feedback’ effect is to reformulate the problem in terms of the perturbative expansion of the steady-state Green’s function in powers of $`J`$, starting from the Dyson equation in which the $`s0`$ limit has already been taken. We now describe this method in more detail, and specify the procedure for calculating the steady-state Green’s functions perturbatively in $`J`$ to arbitrary order.

#### III.3.2 Perturbative expansion of steady-state quantities

Let us now give a general method for computing the perturbative expansion of the Green’s functions, and use it to derive Eq. (4). First, we reformulate slightly the diagrammatic expansion in terms of the full Green’s function $`G_\psi `$ and of skeleton diagrams, and we explicitly derive the steady-state Dyson equation as a functional equation for $`G_\psi `$ (again, for simplicity, we write some equations for a generic field $`\psi `$). The “generic” Dyson equation reads:

$$G_\psi (J,s,\omega )=\left(G_{\psi 0}^1(\omega )\mathrm{\Sigma }_{\text{skel}}[G_\psi ](J,s,\omega )\right)^1,$$
(33)
where $`\mathrm{\Sigma }_{\text{skel}}`$ is a functional of $`G_\psi `$ defined by the skeleton expansion of the self-energy (see Appendix E for a derivation). We first take the $`s0`$ limit in (33) and then solve it order by order in $`J`$. It is in principle sufficient to solve (33) for $`G_\psi `$ since the magnetization can be directly extracted from a Green’s function, and other physical quantities (e.g. currents) are given by their skeleton expansions.

We first derive the explicit form of the Dyson equation for all fields. Defining

$$\sigma _\psi ^K\frac{\mathrm{\Sigma }_\psi ^K}{2i\pi },\sigma _\psi \frac{\mathrm{\Sigma }_\psi ^A\mathrm{\Sigma }_\psi ^R}{2i\pi },$$
(34)
the full set of Dyson equations for the electrons and the Majorana fields can be rewritten explicitly as:

$`𝐑_𝐜(\omega )`$ $`={\displaystyle \frac{1}{𝐑_{\mathrm{𝐜𝟎}}^1𝚺_𝑹(\omega )}},`$ (35)
$`𝐊_𝐜(\omega )`$ $`=\left(1𝐑_{\mathrm{𝐜𝟎}}𝚺_𝑹\right)^1𝐊_{\mathrm{𝐜𝟎}}\left(1𝚺_𝑨𝐀_{\mathrm{𝐜𝟎}}\right)^1+\left(𝐑_{\mathrm{𝐜𝟎}}^1𝚺_𝑹\right)^1𝚺_𝑲\left(𝐀_{\mathrm{𝐜𝟎}}^1𝚺_𝑨\right)^1`$ (36)
$`R_B(\omega )`$ $`={\displaystyle \frac{1}{\omega B\mathrm{\Sigma }_B^R(\omega )+i0^+}},`$ (37)
$`K_B(\omega )`$ $`=h_B(\omega )\left(A_B(\omega )R_B(\omega )\right),`$ (38)
$`h_B(\omega )`$ $`{\displaystyle \frac{\sigma _B^K(\omega )}{\sigma _B(\omega )}},`$ (39)
$`R_z(\omega )`$ $`={\displaystyle \frac{1}{\omega \mathrm{\Sigma }_z^R(\omega )+i0^+}},`$ (40)
$`K_z(\omega )`$ $`=h_z(\omega )\left(A_z(\omega )R_z(\omega )\right),`$ (41)
$`h_z(\omega )`$ $`{\displaystyle \frac{\sigma _z^K(\omega )}{\sigma _z(\omega )}},`$ (42)
where the bold symbols are $`4\times 4`$ matrices (in lead-spin space). For completeness, we have also written the definitions of the $`h`$ functions. We have three blocks of equations, for the electrons, the $`f`$ field and the $`\eta ^z`$ field respectively. Within these blocks, we have an equation for the retarded function (35,37,40), an equation for the Keldysh function (36,38,41), and for the Majorana fermions the definition of the $`h`$ function (39,42). We define $`h_B`$ with (39) rather than with (38) since the spectral density is a delta peak at order 0, whereas the self-energy is a smooth function.

The spin and the lead electrons appear on a different footing: the order 0 part of the electronic Keldysh function is given by $`𝐊_{\mathrm{𝐜𝟎}}`$ whereas the order 0 parts of $`h_B`$ and $`h_z`$ must be computed using equations (39,42). Note that to obtain the order $`n`$ part of $`h_B`$, one needs to compute the self-energies at order $`n+2`$. The method of obtaining the perturbative expansion in $`J`$ order by order is as follows (denoting by $`f^{(n)}`$ the order $`n`$ part of any function $`f`$):

1. Assume that we have the expansion of all functions to order $`n1`$.
2. Since $`\mathrm{\Sigma }_{\text{skel}}`$ is at least of order 1, using (35,37,40), we compute $`R_B^{(n)}`$, $`R_z^{(n)}`$, $`𝐑_{𝐜}^{}{}_{}{}^{(n)}`$ and $`𝐊_{𝐜}^{}{}_{}{}^{(n)}`$.
3. Since $`\mathrm{\Sigma }_{\text{skel}}^K`$ and $`\mathrm{\Sigma }_{\text{skel}}^{\prime \prime }`$ begin at order 2, we compute the $`\sigma ^K`$ and $`\sigma `$ to order $`n+2`$, as functions of the unknowns $`h_B^{(n)}`$ and $`h_z^{(n)}`$.
4. We then obtain closed equations for $`h_B^{(n)}`$ and $`h_z^{(n)}`$ from (39,42) (expanded to order $`n`$).

Thus the order 0 part of the impurity magnetization is given by (from (29a)):

$$M_{\mathrm{dot}}=\frac{1}{2}h_B^{(0)}(B).$$
(43)
The order 0 parts of $`R_B`$ and $`K_B`$ are

$$R_B^{(0)}=\frac{1}{\omega B+i0^+},K_B^{(0)}=2i\pi h_B^{(0)}(B)\delta (\omega B),$$
(44)
and the bare Green’s functions of the electrons are given by (18). We compute $`h_B(\omega =B)`$ at order 0 by expanding the self-energies at second order (given by the diagram of Fig. 2), and then solving (39) for $`h_B(B)`$. Finally we find (4) in the large bandwidth limit ($`D\mathrm{}`$). The computation is presented in detail in Appendix F.

This completes the computation of the Green’s function to order 0.

Note that $`𝐊_𝐜`$ is thermal at order 0, which expresses the fact that the leads are in thermal equilibrium. Therefore the leads’ contribution to the total magnetization at order 0 is given by the Pauli term, which was explicitly excluded from $`M_{\text{tot}}`$; hence $`M_{\text{tot}}=M_{\text{dot}}`$ at this order, as claimed in (4). Moreover, we have not written explicitly the full forms of the functions $`h_B(\omega )`$ and $`h_z(\omega )`$ to this order, since they are not required in the calculation of the zeroth-order Keldysh functions: $`h_B(B)`$ and $`h_z(0)=0`$ are sufficient, since the spectral density (at this order) is a delta peak. The full functions would, however, be needed to compute at second order; the function $`h_B(\omega )`$ can easily be extracted from Appendix F, and the calculation of $`h_z(\omega )`$ proceeds along similar lines.

Our earlier interpretation of the $`1/s`$ divergences in ‘straightforward’ perturbation theory is borne out by this result. As stated above, the $`1/s`$ divergences result from an incorrect choice of zeroth-order distribution function $`h`$. If we insert the $`O(1)`$ part of (39) into (30) in place of $`h_{\mathrm{eq}}`$, we see that the divergences are cancelled at order $`J^2`$, since to this order we are now using the correct long-time distribution function. This shows that our method and that of Wingreen and Sivan (1996) are equivalent.

It is important to note that these corrections to the zeroth-order terms in perturbation series are in no way restricted to quantities such as the magnetization. On the contrary, since what we have really calculated is the correction to the zeroth-order Keldysh Green’s function, they manifest themselves in many quantities. As an example, we may consider the current-current correlator: the leading terms in this quantity are of order $`J^2`$, and are calculated by inserting the zeroth-order Green’s functions into the skeleton diagrams shown schematically in Fig. 3.

The result is that

$`I(\omega )I(\omega )`$ $`=`$ $`{\displaystyle \frac{\pi \left(\rho _0J_{LR}\right)^2}{2\mathrm{}^2}}\{\left[h_B^{(0)}(B)\right]^2p(\omega ,V)h_B^{(0)}(B)[p(\omega +B,V)p(\omega B,V)]`$ (45)
$`+[p(\omega ,V)+p(\omega B,V)+p(\omega +B,V)]\},`$
where the function $`p(x,y)`$ is defined by

$$p(x,y)=\frac{xy}{e^{\beta (xy)}1}+\frac{x+y}{e^{\beta (x+y)}1}.$$
(46)
The result (45) is a function of $`h_B^{(0)}(B)`$, and is therefore clearly sensitive to the corrections made to the zeroth-order Green’s function, as expected on the basis of the discussion above.

## IV Conclusion

The purpose of this note has been to present the expression for the order 0 magnetization in the Kondo model out of equilibrium (Eq. 4), and a systematic procedure for obtaining higher-order corrections to this result. The result we obtain may seem surprising, in the sense that even at order 0 in $`J`$ it does not coincide with the equilibrium expression $`\frac{1}{2}\mathrm{tanh}(B/2T)`$. Indeed the out-of-equilibrium distribution function which describes the population of the two levels of a weakly coupled spin is in general not thermal but must be computed by solving a transport equation: it is determined by the steady state into which the voltage difference forces the system. Moreover, this distribution function also enters the computation of other physical quantities (e.g. currents) and their perturbative expansions therefore exhibit similar phenomena.

Finally, we stress that the issue discussed in this note is not directly related to the so-called “decoherence time” issue. Answering the question of strong coupling at finite voltage requires computations at higher orders in $`J`$. It is clear from the above, however, that the behavior of the out-of-equilibrium perturbative expansion will be markedly different from that of its equilibrium counterpart. For example, the $`O(J)`$ contribution to the Keldysh Green’s function of the Majorana fermions depends on the $`O(J^3)`$ contribution to the self-energies. These, however, contain terms that diverge like $`\mathrm{ln}D`$ in the $`D\mathrm{}`$ limit, and so such logarithmic divergences may be expected to appear at $`O(J)`$ in some of the Green’s functions, and therefore in physical properties such as the magnetization. Indeed, this phenomenon has been reported recently Rosch et al. ; the full interpretation of this striking departure from equilibrium behavior merits further work.

###### Acknowledgements.

This work has been supported by the Center for Materials Theory at Rutgers University, by the EPSRC (UK) under grant number GR/M70476, and by the U.S. Department of Energy under grant number DOE grant DE- FG02 - 00ER45790. We thank L. I. Glazman and A. Kaminski for useful discussions and for the master equation derivation presented in Appendix A. We also thank E. Abrahams, N. Andrei, G. Biroli, R. Chitra, G. Kotliar, D. Langreth, Y. Meir, B. Mouzykantskii, A. Ruckenstein and A. Schofield for useful discussions. Our thanks are also due to the first author of ref. Coleman et al. (2001), P. Coleman, for numerous discussions, though we note that he is in disagreement with the corrected analysis presented above Coleman and Mao , and continues to support the perturbative result presented in Eq. (2) of ref. Coleman et al. (2001). During the completion of this paper, we learned that the same conclusion (Eq. 4) has been independently obtained by A. Rosch, P. Wölfle and J. Kroha. Part of this work was completed during our stay at the Aspen Center for Physics in July 2001.

## Appendix A A semi-classical derivation of (4)

The magnetization at order 0 in $`J`$ (4) can also be derived using a semiclassical master equation approach Glazman and Kaminski . By symmetry, the reduced density matrix of the spin is diagonal, so we consider the spin as classical and characterize its state by the probability of its being up, $`P_{}`$, or down, $`P_{}`$. The spin dynamics is governed by a master equation:

$`{\displaystyle \frac{dP_{}(t)}{dt}}`$ $`=\mathrm{\Gamma }_{}P_{}(t)\mathrm{\Gamma }_{}P_{}(t),`$ (47a)
$`{\displaystyle \frac{dP_{}(t)}{dt}}`$ $`=\mathrm{\Gamma }_{}P_{}(t)\mathrm{\Gamma }_{}P_{}(t),`$ (47b)
where $`\mathrm{\Gamma }_{AB}`$ is the rate of the spin-flip process $`AB`$ induced by the Kondo terms. At second order in perturbation theory, these rates are given by:

$`\mathrm{\Gamma }_{}`$ $`{\displaystyle \underset{\alpha ,\beta \{R,L\}}{}}J_{\alpha \beta }^2{\displaystyle 𝑑ϵ_1𝑑ϵ_2\rho (ϵ_1B_c\mu _\alpha )\rho (ϵ_2+B_c\mu _\beta )n_F(ϵ_1\mu _\alpha )\left(1n_F(ϵ_2\mu _\beta )\right)\delta \left(ϵ_1ϵ_2+B_s\right)}`$ (48a)
$`\mathrm{\Gamma }_{}`$ $`{\displaystyle \underset{\alpha ,\beta \{R,L\}}{}}J_{\alpha \beta }^2{\displaystyle 𝑑ϵ_1𝑑ϵ_2\rho (ϵ_1+B_c\mu _\alpha )\rho (ϵ_2B_c\mu _\beta )n_F(ϵ_1\mu _\alpha )\left(1n_F(ϵ_2\mu _\beta )\right)\delta \left(ϵ_1ϵ_2B_s\right)}`$ (48b)
Taking the large bandwidth ($`D\mathrm{}`$) limit and integrating, the rates become:

$`\mathrm{\Gamma }_{}`$ $`=\left(J_R^2+J_L^2\right){\displaystyle \frac{B_s}{e^{B_s/T}1}}+J_{RL}^2\left({\displaystyle \frac{B_sV}{e^{(B_sV)/T}1}}+{\displaystyle \frac{B_s+V}{e^{(B_s+V)/T}1}}\right),`$ (49a)
$`\mathrm{\Gamma }_{}`$ $`=\left(J_R^2+J_L^2\right){\displaystyle \frac{B_s}{e^{B_s/T}1}}+J_{RL}^2\left({\displaystyle \frac{B_sV}{e^{(B_s+V)/T}1}}+{\displaystyle \frac{B_s+V}{e^{(B_sV)/T}1}}\right).`$ (49b)
The steady state value of the probability is given by equating the left hand sides of (47) to 0, and the magnetization is given by $`M_{\text{dot}}=(P_{}P_{})/2`$, thus

$$M_{\text{dot}}=\frac{1}{2}\left(\frac{\mathrm{\Gamma }_{}\mathrm{\Gamma }_{}}{\mathrm{\Gamma }_{}+\mathrm{\Gamma }_{}}\right),$$
(50)
which leads to Eq. (4).

One can see that the Dyson equation in steady state, at second order in $`J`$, maps exactly to the master equation, although they appear to have different transients. So at order 0, this computation is just a reformulation of the one presented above.

## Appendix B A free model

In this appendix, we recall Caroli et al. (1971); Oguri (2001) the solution of a simple free model in an out-of-equilibrium steady state regime, which displays a result very similar to (4). We consider a free level coupled to some reservoirs of free electrons via hopping terms. The Hamiltonian is given by:

$$H=\underset{\alpha 𝐤}{}\epsilon _𝐤c_{\alpha 𝐤}^{}c_{\alpha 𝐤}+\underset{\alpha }{}g_\alpha \left(c_\alpha ^{}d+d^{}c_\alpha \right)+ϵ_dd^{}d,$$
(51)
where the local reservoir-electron operator $`c_\alpha `$ is defined by $`c_\alpha =c_\alpha (𝐱=\mathrm{𝟎})\mathrm{\Omega }^{1/2}_𝐤c_{\alpha 𝐤}`$. ($`\mathrm{\Omega }`$ is the usual normalization factor related to the volume of the reservoirs.) As before, we assume that the reservoirs are thermalized with the same temperature as each other, but with different chemical potentials $`\mu _\alpha `$. In particular, we are interested in a model with three reservoirs: 1 and 2 are the leads (at different potentials) and 3 an additional thermal bath (with $`\mu _3=0`$) to which the level is coupled. Since it is a Gaussian model, we can simply solve the Dyson equation:

$$G^1=G_0^1\mathrm{\Sigma }.$$
(52)
We use the $`(d,c_1,c_2,c_3)`$ basis. The inverse bare Green’s function is given by:

$$G_0^1(\omega )=\left(\begin{array}{cccc}\omega ϵ_d& 0& 0& 0\\ 0& G_1^1& 0& 0\\ 0& 0& G_2^1& 0\\ 0& 0& 0& G_3^1\end{array}\right),$$
(53)
where $`G_\alpha `$ is the bare Green’s function of the reservoir $`\alpha `$, given in the large bandwidth limit by (in the LO basis):

$$G_\alpha ^1(\omega )=\frac{1}{\pi \rho _0}\left(\begin{array}{cc}i& 2ih_{\mathrm{eq}}(\omega \mu _\alpha )\\ 0& i\end{array}\right),$$
(54)
where $`\rho _0`$ is the density of states, which we take to be the same for each reservoir. The self-energy is given by:

$$\mathrm{\Sigma }(\omega )=\left(\begin{array}{cccc}0& g_1& g_2& g_3\\ g_1& 0& 0& 0\\ g_2& 0& 0& 0\\ g_3& 0& 0& 0\end{array}\right).$$
(55)
Solving (52), we find the occupation of the dot in the steady state:

$$n(ϵ_d)d^{}d.$$
(56)
In the limit where $`g_1,g_2,g_3`$ go to zero in fixed ratios, we obtain:

$$n(ϵ_d)=\frac{g_1^2n_F(ϵ_d\mu _1)+g_2^2n_F(ϵ_d\mu _2)+g_3^2n_F(ϵ_d)}{g_1^2+g_2^2+g_3^2}$$
(57)
where $`n_F`$ is the Fermi function. The properties of this result are similar to those of (4): the occupation of the dot, in the limit of zero couplings, depends on the ratios of these couplings and is not in general given by the Fermi function. Moreover, if $`g_3g_1,g_2`$, we recover the equilibrium result since the physics is dominated by the thermal bath 3. If, on the other hand, we take $`g_3g_1,g_2`$, we find a non-equilibrium result since the occupation at order 0 in perturbation theory in the $`g`$’s is determined by the leads 1 and 2.

## Appendix C Orientability of Majorana fermion lines

In this appendix, we demonstrate that we can treat the Majorana lines as oriented in the diagrams. To show this, it is simplest to take the lines to represent the bare zero-field Majorana functions, and to treat both the magnetic field and the Kondo interaction as vertices. We have three species of Majorana fermion, $`\{\eta ^x,\eta ^y,\eta ^z\}`$; formally, we may represent each of these as the sum of the creation and annihilation operators of a Dirac fermion:

$$\eta ^a=\frac{1}{\sqrt{2}}\left(f_a+f_a^{}\right);$$
(58)
the dual operator $`f_af_a^{}`$ decouples from the problem and may be ignored. In the absence of a magnetic field and interactions, all Majorana Green’s functions are diagonal, so the only Green’s functions that occur are:

$`G_a`$ $`=`$ $`i\eta ^a(t)\eta ^a(0)_𝒦`$ (59)
$`=`$ $`{\displaystyle \frac{i}{2}}f_a(t)f_a^{}(0)+f_a^{}(t)f_a(0)_𝒦`$
$`=`$ $`\left({\displaystyle \frac{i}{2}}f_a(t)f_a^{}(0)_𝒦\right)+\text{p.h.},`$
where the notation “p.h.” stands for “particle-hole”, i.e. $`ff^{}`$. We thus see that the bare Majorana Green’s function may be written simply as the sum of two bare Dirac fermion Green’s functions of opposite orientations.

It is easy to see that the vertex factors at the magnetic field and Kondo vertices do not depend on whether the Dirac fermion lines are incoming or outgoing. In the case of the magnetic field, for example, the interaction term is

$`H_{\mathrm{mag}}`$ $`=`$ $`BS^z`$ (60)
$`=`$ $`iB\eta ^x\eta ^y`$
$`=`$ $`{\displaystyle \frac{iB}{2}}\left(f_x+f_x^{}\right)\left(f_y+f_y^{}\right),`$
so we see that the vertex factor is the same irrespective of the orientation of the two $`f`$-lines; a parallel argument may be given for the Kondo interaction vertex. Hence, each diagram consists of a sum of $`2^N`$ diagrams that differ only in the orientations of their $`N`$ Dirac fermion lines.

But these orientational differences do not alter the value of the diagram, since the bare Green’s function of the $`f`$-fermion is particle-hole symmetric, and hence (59) corresponds to the Green’s function for a single orientation of the Dirac fermion line, with the prefactor $`1/2`$ removed. Hence we may represent the Majorana Green’s functions in all diagrams using oriented lines.

## Appendix D Erratum to Coleman et al. (2001)

In a previous paper Coleman et al. (2001), a formula (Eq. 2) was proposed for the (second order) perturbative expansion of the magnetic susceptibility in the out-of-equilibrium steady state which was of the form:

$$\chi (T,V)=\frac{1}{4T}+O(J).$$
(61)
In this appendix, we briefly rediscuss its derivation and explain why it is incorrect. The method used was straightforward Keldysh perturbation theory to second order with $`s`$ finite in which the bare Keldysh function of the Majorana spin was taken to be thermal. The $`1/s`$ divergences were regulated using a “point splitting” procedure: since they occur due to the coincidence of two poles in the integrals, one splits these poles on the real axis to a distance $`\delta `$ and sends $`\delta `$ to 0 at the end of the computation, term by term. The justification given in footnote 20 of Coleman et al. (2001) is however incorrect: endowing the Majorana fermions on the dot with a fictitious dispersion does not lead to this prescription, since divergences reappear when the bandwidth is sent to zero. The “point splitting” regularization prescription can however be described physically as follows: let the applied magnetic fields $`B_s`$ and $`B_c`$ oscillate slowly at a frequency $`\omega _0`$; compute the total magnetization at frequency $`\omega _0`$ up to second order in $`J`$; then take the limit $`\omega _00`$ in the coefficient of each power of $`J`$.

However, this regulation procedure is based on an interchange of the order of the limits $`\omega _00`$ and $`J0`$. We want to calculate the static magnetic response, and consequently wish to take $`\omega _00`$ before $`J0`$; but in fact the technique used in Coleman et al. (2001) does the opposite: it expands in $`J`$ ($`J0`$) before taking the $`\omega _00`$ limit. It is simple to show that these limits do not commute; this can be seen explicitly from the form of the second term of (30):

$$\text{Re}𝑑\omega \frac{g(\omega )}{\left(\omega +\omega _0+i\gamma \right)\left(\omega \omega _0i\gamma \right)},$$
(62)
where $`g`$ is some function and $`\gamma \mathrm{\Sigma }^{\prime \prime }(B)`$ (the limit $`s0`$ having already been taken). If we take the limits $`J0`$ (i.e. $`\gamma 0`$) and then $`\omega _00`$ we obtain the result of the point splitting prescription of Coleman et al. (2001):

$$𝑑\omega \left(\frac{g^{}(\omega )}{\omega }\right),$$
(63)
whereas if we take them in the opposite order, $`\omega _00`$ and then $`J0`$, we recover (4).

## Appendix E The steady state Dyson equation

In this appendix, we present a formal derivation of the steady state Dyson equation (33). First of all, let us emphasise that our calculation is carried out in the steady state regime, that is in the long-time limit after the switching on of the interactions. We assume that, in this limit, the system evolves to a time-independent steady state under the time evolution described by its Hamiltonian (as required by the Keldysh method: see section III.2). Strictly speaking, this is not possible, since the Hamiltonian (1) conserves the total magnetization of the system and that conservation law prevents the magnetization of the system from relaxing. However, this conservation law is not physical (we have omitted, for example, spin-orbit terms in the leads); therefore, to allow the system to relax to its steady state, we introduce a coupling $`g`$ that breaks the conservation laws. As a specific example, one could consider an anisotropic ($`J_xJ_yJ_z`$) Kondo model. As this extra coupling is relaxed to zero, the transient time taken to reach the steady state diverges but we make the assumption that the values of physical quantities in the steady state are smooth functions of $`g`$. Therefore, once we have taken the $`t\mathrm{}`$ limit, we can set $`g=0`$ in the equation which determines the steady state.

The details of the derivation are as follows:

1. Using the perturbative expansion, we establish the Dyson equation at finite times:

$$𝑑u\left(G_{\psi 0}^1(t,u)\mathrm{\Sigma }_{\text{skel}}[G_\psi ](t,u)\right)G_\psi (u,t^{})=\delta (tt^{})1,$$
(64)
where $`\mathrm{\Sigma }_{\text{skel}}`$ is a functional of $`G_\psi `$ defined by the skeleton expansion of the self-energy. The product should be understood as a matrix product in the LO basis. $`G_\psi `$ is a function of two times and of $`J,g`$. To obtain this equation, we write the Dyson equation in the finite time diagrammatic expansion outlined in section III.2.4 Rammer and Smith (1986), and use the definition of the skeleton diagrams Nozières (1964); Abrikosov et al. (1963). In (64), the times $`t,t^{}`$ and $`u`$ run from $`\mathrm{}`$ to $`\mathrm{}`$, and the couplings $`J`$ are time dependent: $`J_\alpha (t)=J_\alpha \theta (t)`$, i.e. we switch on the interaction suddenly at $`t=0`$.
2. The assumption that the system relaxes to a non-equilibrium time translation invariant steady state (the coupling $`g`$ to the relaxation bath is finite) is transcribed mathematically as the existence of the limit

$$\underset{\genfrac{}{}{0.0pt}{}{t,t^{}\mathrm{}}{tt^{}=\tau }}{lim}G(t,t^{})=G(\tau )$$
(65)
(we shall denote it with the same function name). This assumption is not trivial. In particular, in Coleman et al. , the existence of steady oscillating states has been suggested using a large-$`N`$ slave boson treatment. We exclude them here on physical grounds, since we do not expect such states to appear in the regime where the perturbation theory is applicable anyway (i.e. at high temperatures or at high magnetic fields in the Kondo problem). We can thus take the long time limit, and Fourier transform the Green’s functions and self-energy in (64) to obtain:

$$G_\psi (J,s,g,\omega )=\left(G_{\psi 0}^1(\omega )\mathrm{\Sigma }_{\text{skel}}[G_\psi ](J,s,g,\omega )\right)^1.$$
(66)
3. Taking the $`g0`$ limit as discussed above, we finally obtain the Dyson equation (33):

$$G_\psi (J,s,\omega )=\left(G_{\psi 0}^1(\omega )\mathrm{\Sigma }_{\text{skel}}[G_\psi ](J,s,\omega )\right)^1.$$
(67)

## Appendix F Computation of the self-energy diagram

In this appendix, we present the computation of $`h_B(B)`$ at order 0. The computation is in three steps. First, we compute the self-energy diagrams at second order as a function of a “generic” diagram. (This simplifies the problem by separating the Majorana and spin indices from the Keldysh structure.) Second, we compute this generic diagram. Finally, we solve the resulting implicit equation for $`h_B(B)`$. Throughout this appendix, we streamline our notation by omitting the temperature $`T`$; in the final formulas, therefore, $`B`$ and $`V`$ should be replaced by $`B/T`$ and $`V/T`$ respectively.

The first part of the computation reduces the spin and lead indices, and thus expresses the self energy diagrams as functions of the “generic” diagram presented in Fig. 4, where we allow any potentials $`V_1`$ and $`V_2`$ for the electrons and any field $`B`$ for the internal Majorana line.

The main formula is:

$$\mathrm{\Sigma }_B(\omega )=\underset{\genfrac{}{}{0.0pt}{}{\alpha ,\beta \{R,L\}}{\sigma =+,}}{}\frac{|J_{\alpha \beta }|^2}{4}\left(f(\omega ,B_s,V_\alpha ,\sigma ,V_\beta ,\sigma )+f(\omega ,0,V_\alpha ,1,V_\beta ,1)\right).$$
(68)
In this expression, the Keldysh structure is implicit and by convention, the Majorana line is an $`\eta ^z`$ line when $`B_s=0`$ and an $`f`$ line otherwise. To establish (68), we compute the spin and lead indices of $`\mathrm{\Sigma }_{xx}`$ and $`\mathrm{\Sigma }_{xy}`$ which are given by the diagrams of Fig. 5 (with the Feynman rules given in section III.2), and we use (28).

We now compute the Keldysh structure of the generic diagram $`f`$.

We have:

$`f_+(\omega )`$ $`=i\pi {\displaystyle 𝑑x𝑑y\rho (x+\sigma _1B_cV_1)\rho (y+\sigma _2B_cV_2)n_F(xV_1)\stackrel{ˇ}{n}_F(yV_2)(1+\stackrel{~}{h}(B))\delta (\omega +yxB)},`$
$`f_+(\omega )`$ $`=i\pi {\displaystyle 𝑑x𝑑y\rho (x+\sigma _1B_cV_1)\rho (y+\sigma _2B_cV_2)\stackrel{ˇ}{n}_F(xV_1)n_F(yV_2)(1\stackrel{~}{h}(B))\delta (\omega +yxB)},`$ (69)
where $`\stackrel{ˇ}{n}_F(x)n_F(x)`$, and the function $`\stackrel{~}{h}(\omega )`$ is defined by

$$\stackrel{~}{h}(\omega )=\{\begin{array}{cc}h_B(\omega )\hfill & \omega 0,\hfill \\ h_z(0)\hfill & \omega =0.\hfill \end{array}$$
(70)
We then use the relations Rammer and Smith (1986)

$`f^K`$ $`=\left(f_++f_+\right),`$ (71)
$`f^Af^R`$ $`=\left(f_+f_+\right).`$ (72)
We can perform the integrals in the large bandwidth limit $`D\mathrm{}`$ using

$$𝑑xn_F(x+A)\stackrel{ˇ}{n}_F(x)=\frac{A}{e^A1},$$

and replacing the densities $`\rho `$ by their finite value $`\rho _0`$ (we can take the limit $`D\mathrm{}`$ under the integral). Using the definition $`V=V_2V_1`$, we find:

$`{\displaystyle \frac{f^K}{2i\pi }}`$ $`={\displaystyle \frac{\rho _0^2}{2}}(\omega B+V)\left(1{\displaystyle \frac{\stackrel{~}{h}(B)}{\mathrm{tanh}\frac{\omega B+V}{2}}}\right),`$ (73)
$`{\displaystyle \frac{\left(f^Af^R\right)}{2i\pi }}`$ $`={\displaystyle \frac{\rho _0^2}{2}}(\omega B+V)\left({\displaystyle \frac{1}{\mathrm{tanh}\frac{\omega B+V}{2}}}\stackrel{~}{h}(B)\right).`$ (74)
Using (68), and introducing

$$\phi (x)\frac{x}{\mathrm{tanh}\frac{x}{2}},$$
(75)
we now have:

$`\sigma _B(\omega )=`$ $`{\displaystyle \frac{|\rho _0J_{RL}|^2}{4}}\left[h_B(B)(\omega B+V)h_z(0)(\omega +V)+\phi (\omega B+V)+\phi (\omega +V)+(VV)\right]`$
$`+{\displaystyle \underset{\alpha =R,L}{}}{\displaystyle \frac{|\rho _0J_\alpha |^2}{4}}\left[h_B(B)(\omega B)h_z(0)(\omega )+\phi (\omega B)+\phi (\omega )\right],`$ (76)
$`\sigma _B^K(\omega )=`$ $`{\displaystyle \frac{|\rho _0J_{RL}|^2}{4}}\left[(\omega B+V)\left(1{\displaystyle \frac{h_B(B)}{\mathrm{tanh}\frac{\omega B+V}{2}}}\right)+(\omega +V)\left(1{\displaystyle \frac{h_z(0)}{\mathrm{tanh}\frac{\omega +V}{2}}}\right)+(VV)\right]`$
$`{\displaystyle \underset{\alpha =R,L}{}}{\displaystyle \frac{|\rho _0J_\alpha |^2}{4}}\left[(\omega B)\left(1{\displaystyle \frac{h_B(B)}{\mathrm{tanh}\frac{\omega B}{2}}}\right)+\omega \left(1{\displaystyle \frac{h_z(0)}{\mathrm{tanh}\frac{\omega }{2}}}\right)\right].`$ (77)
Using the definition of the Majorana Green’s function $`G_{ab}^K(t)`$,

$$G_{ab}^K(t)=i[\eta _a(t),\eta _b(0)],$$
(78)
we have $`G_{ab}^K(t)=G_{ba}^K(t)`$. In the frequency representation, this reads $`G_{ab}^K(\omega )=G_{ba}^K(\omega )`$ and hence we infer that $`G_{zz}^K(\omega )`$ is odd in frequency. Consequently, given the form $`G_{zz}^K=h_z(\omega )\delta (\omega )`$ it is clear that

$$h_z(0)=0.$$

Using this result together with (39) and denoting $`x=h_B(B)`$, we have:

$$x=\frac{|J_{RL}|^2\left(2x\phi (V)+2B\right)+(J_R^2+J_L^2)\left(2x+B\right)}{|J_{RL}|^2\left(2\phi (V)+\phi (B+V)+\phi (BV)\right)+(J_R^2+J_L^2)\left(2+\phi (B)\right)}.$$

Solving for $`x`$ and substituting into (43) gives Eq. (4) of the text.


text:
# A limit on the detectability of the energy scale of inflation

## Abstract

We show that the polarization of the cosmic microwave background can be used to detect gravity waves from inflation if the energy scale of inflation is above $`3.2\times 10^{15}`$ GeV. These gravity waves generate polarization patterns with a curl, whereas (to first order in perturbation theory) density perturbations do not. The limiting “noise” arises from the second–order generation of curl from density perturbations, or rather residuals from its subtraction. We calculate optimal sky coverage and detectability limits as a function of detector sensitivity and observing time.

Few ideas have had greater impact in cosmology than that of inflationguth81 ; albrecht82 ; linde82 . Inflation makes four predictions, three of which provide very good descriptions of data: the mean curvature of space is vanishingly close to zero, the power spectrum of initial density perturbations is nearly scale–invariant, and the perturbations follow a Gaussian distribution. As the data have improved substantially (e.g., halverson01 ; netterfield01 ; lee01 ) they have agreed well with inflation, whereas all competing models for explaining the large–scale structure in the Universe have been ruled out (e.g., pen97 ; allen97 ; albrecht97 ).

We must note though that these three predictions are all fairly generic111For example, a much more impressive success would be prediction of features in the power spectrum followed by their detection.. Further, although existing models for the formation of structure have been ruled out, there is no proof of inflation’s unique ability to lead to our Universe. Indeed, alternatives are being invented khoury01 .

The fourth (and yet untested) prediction may therefore play a crucial role in distinguishing inflation from other possible early Universe scenarios. Inflation inevitably leads to a nearly scale–invariant spectrum of gravitational waves, which are tensor perturbations to the spatial metric. Detection of this gravitational–wave background might allow discrimination between competing scenarios (e.g. khoury01 ), and different inflationary models (e.g. dodelson97 ).

The amplitude of the power spectrum of tensor perturbations to the metric is directly proportional to the energy scale of inflation. One can use a determination of the tensor contribution to CMB temperature anisotropy, here parameterized by the quadruple variance, to determine this energy scale turner95 :

$$V_{}^{1/4}/m_{Pl}=1.2Q_T^2^{1/4}=3.0\times 10^3r^{1/4}$$
(1)
where $`rQ_T^2/Q_S^2`$, $`S`$ stands for scalar (density) perturbation and $`Q_S^24\times 10^{11}`$ from observations222The * subscript on $`V`$ means that it is evaluated when the scale factor is about $`e^{50}`$ times smaller than it is at the end of inflation, when the relevant perturbations for large–scale CMB fluctuations are exiting the Horizon. There is a slight dependence on $`\mathrm{\Omega }_\mathrm{\Lambda }`$ knox95 , here taken to be $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$. Other cosmological parameters assumed throughout are $`h=0.65`$, $`\mathrm{\Omega }_b=0.05`$, $`\mathrm{\Omega }_m=0.35`$ and $`\delta _H=4.2\times 10^5`$.. Without detection of gravitational waves, the energy scale of inflation remains uncertain by at least 12 orders of magnitude. Pinning down this energy scale could be crucial to understanding how inflation arises in a fundamental theory of physics.

In Fig. 1 we show the angular power spectrum of CMB temperature perturbations contributed by scalar perturbations and by tensor perturbations with $`r=10^3`$. By determining the total CMB temperature power spectrum we can determine or limit the energy scale of inflation, based on the presence or absence of extra power at low $`l`$. The scalar temperature perturbations inevitably limit our ability to detect the tensor temperature perturbations to those cases with $`r>r_{\mathrm{lim}}=0.13`$. knox94 333Here and throughout, the detection limit is set at 3.3 times the standard deviation. For a normal distribution this means that one can make a detection at 95% confidence for 95% of possible realizations.

In kamionkowski97 ; seljak97  it was pointed out that tensor perturbations result in CMB polarization patterns with a curl, whereas scalar perturbations do not. By analogy with electromagnetism, these modes are called “B-modes”, and the curl–free modes are called “E-modes”. This was an exciting development, because the new signature of tensor perturbations could not be confused with scalar perturbations. Figure 1 also shows the power spectrum of this B-mode with the amplitude it would have if $`V_{}^{1/4}=6.4\times 10^{15}`$ GeV, corresponding to $`r=10^3`$. It is a very weak signal, even at its peak $`10^7`$ times less than the power spectrum of the temperature anisotropy. Ignoring any contaminating source of curl–mode polarization the detectability limit is solely a matter of sufficient sensitivity to the CMB polarization. For a detector of sensitivity $`s`$ uniformly observing the entire sky for time $`t`$:

$$r_{\mathrm{lim}}=10^2\left(\frac{s}{\mu K\sqrt{\mathrm{sec}}}\right)^2\left(\frac{t}{1\mathrm{year}}\right)^1.$$
(2)
Equation 2 kamionkowski98  does not take into account a contaminating source of B mode that arises from the lensing of the E mode by density perturbations along lines–of–sight between the observer and the last–scattering surface zaldarriaga98 . This scalar contribution to the B-mode power spectrum is shown in Fig. 1. As we will see below this contamination sets the detectability limit at $`r_{\mathrm{lim}}=6\times 10^4`$, similar to what was found in lewis01 .

This lensing contaminant can be cleaned from the maps as discussed in hu01 ; hu02 . Here we show that there are limits to how well these maps can be cleaned. Our main result is that in the no–noise limit the tensor perturbations can be detected only if $`V_{}^{1/4}>3.2\times 10^{15}\mathrm{GeV}`$. Interestingly, this is below the mass scale for the gauge coupling–constant convergence at $`10^{16}\mathrm{GeV}`$ that occurs in minimal extensions to the standard model of particle physics (e.g. langacker91 ).

The Stokes parameters, $`I`$, $`Q`$ and $`U`$ are related to the unlensed Stokes parameters (denoted with a tilde) by

$`I(\stackrel{}{\theta })=\stackrel{~}{I}(\stackrel{}{\theta }+\delta \stackrel{}{\theta })`$ $`Q(\stackrel{}{\theta })=\stackrel{~}{Q}(\stackrel{}{\theta }+\delta \stackrel{}{\theta })`$ $`U(\stackrel{}{\theta })=\stackrel{~}{U}(\stackrel{}{\theta }+\delta \stackrel{}{\theta }).`$ (3)
The deflection angle, $`\delta \stackrel{}{\theta }`$, is the tangential gradient of the projected gravitational potential,

$$\varphi (\stackrel{}{\theta })=2𝑑r\frac{(rr_s)}{rr_s}\mathrm{\Psi }(r\widehat{n},r)$$
(4)
where $`r`$ is the coordinate distance along our past light cone, $`s`$ denotes the CMB last–scattering surface and $`\widehat{n}`$ is the unit vector in the $`\stackrel{}{\theta }`$ direction.

The effect on the B–mode power spectrum is zaldarriaga98

$$C_{Bl}=C_{\stackrel{~}{B}l}+\underset{l^{}}{}W_l^l^{}C_{\stackrel{~}{E}l^{}}$$
(5)
where

$`W_l^l^{}={\displaystyle \frac{l^3}{4}}{\displaystyle _0^\pi }\theta d\theta \times \{`$
$`\sigma _2^2(\theta )\left[J_0(l\theta )J_2(l^{}\theta ){\displaystyle \frac{1}{2}}J_4(l\theta )\left[J_2(l^{}\theta )+J_6(l^{}\theta )\right]\right]+`$
$`\sigma _0^2(\theta )[J_4(l\theta )J_4(l^{}\theta )+J_0(l^{}\theta )J_0(l\theta )]\}.`$ (6)
The functions $`\sigma _2^2(\theta )`$ and $`\sigma _0^2(\theta )`$ depend on the statistical properties of the displacement potential $`\varphi `$ and are given by

$`\sigma _0^2(\theta )`$ $`=`$ $`{\displaystyle \frac{ldl}{2\pi }l^2C_l^\varphi (1J_0(l\theta ))}`$
$`\sigma _2^2(\theta )`$ $`=`$ $`{\displaystyle \frac{ldl}{2\pi }l^2C_l^\varphi J_2(l\theta )}`$ (7)
where

$$C_l^\varphi \varphi _{lm}\varphi _{lm}$$
(8)
and $`\varphi _{lm}`$ is the spherical–harmonic transform of $`\varphi (\stackrel{}{\theta })`$.

If we have a means of determining $`\varphi `$, and therefore $`\delta \stackrel{}{\theta }`$, we can reconstruct the unlensed maps from the lensed maps by use of Eq. 3. This procedure cleans out the lensing–induced B–mode. The $`C_l^\varphi `$ above can either be interpreted as the power spectrum of the lensing potential (in the case of uncleaned maps in which case we will call it $`C_l^{\varphi (S)}`$) or the power spectrum of the lensing potential residuals (in the case of lensing–cleaned maps). In the uncleaned case hu00  (in the Limber approximation valid at small scales):

$$C_l^{\varphi (S)}=\frac{8\pi ^2}{l^3}H_0^2_0^{r_s}𝑑rr\left[\frac{(rr_s)}{rr_s}\right]^2\mathrm{\Delta }_\mathrm{\Phi }^2(k,r)|_{k=l/r}$$
(9)
where $`r_s`$ is the comoving distance to the last-scattering surface and $`\mathrm{\Delta }_\mathrm{\Phi }^2(k,r)k^3/(2\pi )^2P_\mathrm{\Phi }(k,r)`$ where $`P_\mathrm{\Phi }(k,r)`$ is the power spectrum of the gravitational potential at the time corresponding to coordinate distance $`r`$ on our past light cone. We plot $`C_l^{\varphi (S)}`$ in Fig. 2

We consider the $`\varphi `$ reconstruction procedure given in hu02  which exploits the fact that lensing leads to a mode–mode coupling with expectation value proportional to $`\varphi `$. Their estimator is a minimum–variance (MV) average over pairs of map modes with $`ll^{}`$. Since we know the statistics of the signal, $`C_l^\varphi `$, we can Wiener filter (WF) the MV estimate and further reduce the errors in the reconstruction. The error in the WF estimate of $`\varphi _{lm}`$ has variance:

$$C_l^{\varphi (WF)}=\frac{C_l^{\varphi (S)}C_l^{\varphi (MV)}}{C_l^{\varphi (S)}+C_l^{\varphi (MV)}}.$$
(10)
The Wiener–filtering is important in the low signal–to–noise regime $`(C_l^{\varphi (S)}/C_l^{\varphi (MV)}1)`$ where the variance of the reconstructed $`\varphi `$ is $`C_l^{\varphi (S)}`$ instead of $`C_l^{\varphi (MV)}`$.

Cleaning can greatly reduce the amplitude of $`C_l^\varphi `$ as is shown in Fig. 2 where the angular power spectrum of the residual $`\varphi `$ is shown for two different experiments. The first is the ‘reference’ experiment of hu02  which makes temperature maps with weight-per-solid angle of $`w=(1\mu \mathrm{K}\mathrm{arcmin})^2`$ and $`Q`$ and $`U`$ maps each with half this weight, all with $`7^{}`$ (full-width-half-maximum) resolution. The second is the no–noise, perfect angular resolution case. One can see that high signal–to–noise reconstruction of $`\varphi `$ is possible.

Because the cleaning is not perfect, there is still lensing-induced scalar B-mode in the cleaned maps. In the no-noise limit the residual scalar B-mode power spectrum, $`C_{Bl}^{WF}`$, has been reduced by a factor of 10 from the uncleaned amplitude, as shown in Fig. 1.

Writing the tensor–induced B-mode power spectrum as $`C_{Bl}^T=rC_{Bl}^{T,1}`$ the resulting error in $`r`$, $`\sigma _r`$, is given by

$$\frac{1}{\sigma _r^2}=f_{\mathrm{sky}}\underset{l}{}\frac{2l+1}{2}\left(C_{Bl}^{T,1}\right)^2\left(C_{Bl}^{WF}+w^1e^{l^2\sigma _b^2}\right)^2$$
(11)
where the $`Q`$ and $`U`$ maps cover $`f_{\mathrm{sky}}`$ of the sky with uniform weight per solid angle $`w`$ and have been smoothed by a circular Gaussian beam profile with full-width-half-maixmum = $`\sqrt{8\mathrm{ln}2}\sigma _b`$. For $`r_{\mathrm{lim}}=3.3\sigma _r`$ we find for noise–free maps, but with no cleaning done, $`r_{\mathrm{lim}}=6\times 10^4`$. For the ‘reference’ experiment of hu01  we find $`r_{\mathrm{lim}}=4\times 10^4`$ after cleaning with the MV estimator and $`1.8\times 10^4`$ with the WF estimator. Finally, for noise–free maps with WF cleaning $`r_{\mathrm{lim}}=6\times 10^5`$.

All of the above results are for $`f_{\mathrm{sky}}=1`$. Given a detector (or array of detectors) with total sensitivity $`s`$ and an amount of observing time $`t`$ there is an optimal amount of sky to cover. Too much sky coverage means the maps are noisy and the lensing contaminant can not be removed accurately enough, too little sky hampers the statistical subtraction of the residual scalar B-mode power spectrum. The optimal sky coverage is given by $`f_{\mathrm{sky}}^{\mathrm{opt}}=10^8(\mu \mathrm{K})^2t/s^2`$. In Fig. 3 we plot $`r_{\mathrm{lim}}`$ as a function of $`s^2/t`$ assuming optimal sky coverage (constrained to $`f_{\mathrm{sky}}1`$). The detectability limit is not strongly sensitive to deviations from optimal sky coverage, which is fortunate since other considerations can influence sky coverage choice as well. Note that $`s^2/t=10^8(\mu \mathrm{K})^2`$ could be achieved by observing for a year with an array of 30,000 detectors each with sensitivity $`100\mu \mathrm{K}\sqrt{\mathrm{sec}}`$. Ten times more detectors (or integration time) are required to reach $`s^2/t=10^9(\mu \mathrm{K})^2`$ and approach the smallest possible $`r_{\mathrm{lim}}`$. We have ignored aliasing of the scalar E–mode wich can be important for $`f_{\mathrm{sky}}1`$, as studied in lewis01 . Aliasing considerations increase $`f_{\mathrm{sky}}^{\mathrm{opt}}`$ and $`r_{\mathrm{lim}}`$. For $`s^2/t=10^7(\mu \mathrm{K})^2`$ and $`f_{\mathrm{sky}}=1`$, $`r_{\mathrm{lim}}=7\times 10^4`$. Therefore, including effects of aliasing would not increase $`r_{\mathrm{lim}}`$ by more than 40% over our results as shown in Fig. 3.

Since we are interested in determining a fundamental limit to the detectability of tensor perturbations, undeterred by the virtual cost of imaginary experiments, we explore one further method of cleaning out the lensing contaminant. Imagine we reconstruct the projected gravitational potential or the density field on our past light cone out to some limiting redshift, $`z_{\mathrm{lim}}`$444Perhaps this is done with a POTENT–likebertschinger89  use of peculiar velocity measurements, or analysis of the statistical distribution of apparent galaxy shapes.. Then we could directly calculate the lensing generated by the matter at $`z<z_{\mathrm{lim}}`$. Subtracting this lensing contaminant from CMB maps leaves only the contribution from $`z>z_{\mathrm{lim}}`$. The effect of this remaining lensing can be calculated by altering the lower limit of integration in Eq. 9 from 0 to $`r(z_{\mathrm{lim}})`$. To improve upon the noise–free self–cleaned maps result of $`r_{\mathrm{lim}}=6\times 10^5`$ one needs to take $`z_{\mathrm{lim}}\begin{array}{c}>\hfill \\ \hfill \end{array}7`$. We know of no way of obtaining this information at such high $`z`$.

Our analysis ignores polarized emission from galactic and extragalactic sources. Multi–frequency observations can be used to clean out these signals based on their distinct spectral shapes. However, even in the no–noise limit these cannot be cleaned perfectly since their spectral shapes are not perfectly well known and vary spatially (e.g. tegmark00b ). Our estimate of the detectability limit should be viewed as a lower limit.

How likely is it that the energy scale of inflation is high enough to be determined? It is conceivable that $`V_{}^{1/4}`$ is as low as $`10^3\mathrm{GeV}`$knox93 . Given even odds between $`10^3`$ and $`10^{19}`$ GeV the chances are small. On the other hand gauge coupling unification gives us a hint that something interesting may be occurring at $`10^{16}`$ GeV. If inflation has anything to do with grand unification, or physics at higher energy scales, then the chances are good.

The likelihood of detection may also be addressed by better determination of the scalar spectrum. If we assume particular functional forms for $`V(\varphi )`$ and a scalar perturbation spectrum with power–spectral index, $`n`$, near the scale–invariant value of unity, then we can give constraints on the range of possible values of $`r`$. Here we follow the nomenclature and calculations of dodelson97 . Exponential potentials have $`r=5(1n)`$ and “small–field polynomial” potentials ($`V(\varphi )=\mathrm{\Lambda }^4\left[1\left(\varphi /\mu \right)^p\right]`$) with $`p=2`$ have $`r=10\mathrm{exp}\left(50\left(1n\right)\right)`$. Still other classes of models (hybrid inflation and the small–field polynomial potentials with $`p>2`$) leave little or no relationship between $`r`$ and $`n`$; $`r`$ can take on detectable values or vanishingly small ones. Perhaps the happiest case is the simplest one: polynomial potentials with $`V(\varphi )=\mathrm{\Lambda }^4(\varphi /\mu )^p`$. With $`p2`$ all have $`r>0.1`$, well above our detectability limit.

In sum, we have calculated how large the amplitude of tensor perturbations must be in order for their effect on the CMB to be distinguishable from those of lensing. By reconstructing the lensing potential, using the mode–mode coupling induced by lensing, one can reduce the lensing contamination but not eliminate it. The residual lensing signal prevents detection of tensor perturbations if $`V_{}<(3.2\times 10^{15}\mathrm{GeV})^4`$. If $`V_{}`$ is slightly larger then an ambitious observational program may succeed in determining the energy scale of inflation—possibly a key step towards finding inflation a comfortable home in a fundamental theory of physics.

We thank A. Albrecht, D. Chung, N. Dalal and J. Kiskis for useful conversations. We used CMBfast seljak96 .


text:
# 1 Introduction

## 1 Introduction

Based on max-plus algebra, a class of discrete event processes can be described by linear recursive equations. Such a system exhibits asymptotic periodic behavior, and its performance can be evaluated by calculating the eigenvalue of the system matrix in max-plus algebra. On the other hand, control techniques based on switching among different controllers have been explored extensively in recent years, where they have been shown to achieve better dynamic performance.

This paper proposes a new model for a class of switched discrete event systems. Such a model consists of a finite set of discrete event subsystems, and a switching law that orchestrates the switching among them. We show that, the switched system can be transformed into a non-switched system, and under certain conditions, the switched system exhibits asymptotic periodic behavior, and its performance can be evaluated by calculating the eigenvalue of certain matrix in max-plus algebra.

## 2  Preliminaries

Denote

$$R_e=R\{\mathrm{}\}$$

$$ϵ=\mathrm{}$$

and for any $`x,yR_e`$, define

$$xy=\mathrm{max}\{x,y\}$$

$$xy=x+y$$

A matrix $`A=(a_{ij})R_e^{n\times n}`$ is said to be irreducible if $`i,j,(i_1=i,i_2,\mathrm{}\mathrm{},i_{k1},i_k=j)`$, s.t. $`a_{ii_2}+a_{i_2i_3}+\mathrm{}\mathrm{}+a_{i_{k1}j}>\mathrm{}`$.

For any matrices $`A,BR_e^{n\times n}`$, define

$$AB=(a_{ij}b_{ij})$$

$$AB=(\underset{k=1}{\overset{n}{}}(a_{ik}b_{kj}))=:AB$$

Given any matrix $`AR_e^{n\times n}`$, the corresponding directed graph (digraph) is a graph with $`n`$ nodes, and there is a directed arc from node $`j`$ to node $`i`$ with weight $`a_{ij}`$ if and only if $`a_{ij}\mathrm{}`$.

In a digraph, a circuit is a directed path that starts and ends at the same node. In a circuit, the sum of the weights of all its arcs divided by the number of arcs is called the mean weight. The circuit with the maximal mean weight in a digraph is called the critical circuit.

A zero vector is a vector with all its entries equal to $`\mathrm{}`$.

For an irreducible matrix $`AR_e^{n\times n}`$, if there exist a real number $`\lambda `$ and a nonzero vector $`hR_e^{n\times 1}`$ such that $`Ah=\lambda h`$. Then $`\lambda `$ and $`h`$ are called the eigenvalue, eigenvector of $`A`$, respectively.

Lemma 1

For an irreducible matrix $`AR_e^{n\times n}`$, there is a unique eigenvalue $`\lambda `$, and it equals the mean weight of the critical circuit of its corresponding digraph.

Lemma 2

For an irreducible matrix $`AR_e^{n\times n}`$, there exist positive integers $`k_0`$ and $`d`$ such that

$$A^{k+d}=\lambda ^dA^k,kk_0$$

$`d`$ is called the period order of $`A`$.

Lemma 3

For an irreducible matrix $`AR_e^{n\times n}`$, suppose its eigenvalue is $`\lambda `$, and its period order is $`d`$. Then there exists a positive integer $`k_0`$ such that the solution of

$$X(k+1)=AX(k)$$

satisfies

$$X(k+d)=\lambda ^dX(k),kk_0$$

This shows that the system will exhibit periodic behavior asymptotically. The mean period is exactly equal to the eigenvalue of $`A`$. Hence, the eigenvalue of $`A`$ is an important performance index of the system.

## 3  Switched Systems

For notational simplicity, we first discuss switching between two subsystems. That is, the switched system is governed by

$$X(k+1)=A_iX(k)$$
(1)
where $`A_iR_e^{n\times n}`$, and the switching law is

$$i=\{\begin{array}{cc}1\hfill & k\text{ even}\hfill \\ 2\hfill & k\text{ odd}\hfill \end{array}$$

Namely

$$X(1)=A_1X(0)$$

$$X(2)=A_2X(1)$$

$$X(3)=A_1X(2)$$

$$X(4)=A_2X(3)$$

$$\mathrm{}\mathrm{}$$

That is

$$X(2)=A_2X(1)=A_2A_1X(0)$$

$$X(4)=A_2X(3)=A_2A_1X(2)$$

$$\mathrm{}\mathrm{}$$

Let

$$Y(k)=X(2k)$$
(2)
Then

$$Y(1)=A_2A_1Y(0)$$

$$Y(2)=A_2A_1Y(1)$$

$$\mathrm{}\mathrm{}$$

$$Y(k+1)=A_2A_1Y(k)$$
(3)
In this way, we tramsform a switched system into a non-switched system. Thus, the following problem naturally arises: Suppose $`A_1`$ and $`A_2`$ are irreducible matrices, is their product $`A_2A_1`$ still irreducible?

The answer is NO in general case. Consider the two irreducible matrices

$$A_1=\left[\begin{array}{ccc}ϵ& 1& ϵ\\ ϵ& ϵ& 1\\ 1& ϵ& ϵ\end{array}\right]A_2=\left[\begin{array}{ccc}ϵ& ϵ& 1\\ 1& ϵ& ϵ\\ ϵ& 1& ϵ\end{array}\right]$$

Then, their product is

$$A_2A_1=\left[\begin{array}{ccc}2& ϵ& ϵ\\ ϵ& 2& ϵ\\ ϵ& ϵ& 2\end{array}\right]$$

Clearly, $`A_2A_1`$ is reducible. However, if every main diagonal entry of $`A_1`$ (or $`A_2`$) is not the null element $`ϵ`$, then the answer to the question above is YES.

Theorem 1

Suppose $`A,BR_e^{n\times n}`$ are irreducible matrices, with all the main diagonal entries of $`A`$ (or $`B`$) not equal to $`ϵ`$. Then, $`AB`$ is irreducible, too.

Proof: Without loss of generality, suppose all the main diagonal entries of $`A`$ are not equal to $`ϵ`$. Then, for any $`1s,tn`$, $`(AB)_{st}ϵ`$ whenever $`(B)_{st}ϵ`$. Moreover, since $`B`$ is irreducible, by definition, $`AB`$ is irreducible, too.

Theorem 2

Suppose $`A_1,A_2R_e^{n\times n}`$ are irreducible matrices, with all the main diagonal entries of $`A_1`$ (or $`A_2`$) not equal to $`ϵ`$. Then, there exist positive number $`\lambda `$, positive integers $`d`$ and $`k_0`$, such that the switched system (1) satisfies

$$X(k+d)=\lambda ^dX(k),kk_0$$

Proof: By the transformation (2), the switched system (1) can be transformed into a non-switched system (3). That is

$$Y(k+1)=A_2A_1Y(k)$$

By Theorem 1, $`A_2A_1`$ is irreducible. Hence, by Lemma 3 and by the transformation (2), we get the result.

Example 1

Consider the two irreducible matrices

$$A=\left[\begin{array}{ccc}2& ϵ& 3\\ 6& 2& ϵ\\ ϵ& 4& 3\end{array}\right]B=\left[\begin{array}{ccc}ϵ& 3& ϵ\\ ϵ& ϵ& 2\\ 4& ϵ& ϵ\end{array}\right]$$

It is easy to see that

$$\lambda (A)=\frac{13}{3},\lambda (B)=3$$

Moreover

$$AB=\left[\begin{array}{ccc}7& 5& ϵ\\ ϵ& 9& 4\\ 7& ϵ& 6\end{array}\right]$$

is also irreducible, and

$$\lambda (AB)=9$$

Note that

$$\lambda (AB)>\lambda (A)+\lambda (B)$$

But this inequality is not always true in general case.

Example 2

Consider the two irreducible matrices

$$A=\left[\begin{array}{ccc}10& 1& ϵ\\ ϵ& 1& 1\\ 1& ϵ& 1\end{array}\right]B=\left[\begin{array}{ccc}1& 1& ϵ\\ ϵ& 1& 1\\ 1& ϵ& 10\end{array}\right]$$

It is easy to see that

$$\lambda (A)=10,\lambda (B)=10$$

Moreover

$$AB=\left[\begin{array}{ccc}11& 11& 2\\ 2& 2& 11\\ 2& 2& 11\end{array}\right]$$

is also irreducible, and

$$\lambda (AB)=11$$

Hence

$$\lambda (AB)<\lambda (A)+\lambda (B)$$

## 4  Some Extensions

More complicated switching laws can be accommodated for performance evaluation. Suppose $`A_iR_e^{n\times n},i=1,2,\mathrm{}\mathrm{},m`$ are irreducible matrices, with all their main diagonal entries not equal to $`ϵ`$. This switched system is governed by

$$X(k+1)=A_iX(k)$$
(4)
with switching law

$$i=\{\begin{array}{cc}1\hfill & k=0,1,2,\mathrm{}\mathrm{},k_1mod(K)\hfill \\ 2\hfill & k=k_1+1,k_1+2,\mathrm{}\mathrm{},k_2mod(K)\hfill \\ 3\hfill & k=k_2+1,k_2+2,\mathrm{}\mathrm{},k_3mod(K)\hfill \\ \mathrm{}\hfill & \mathrm{}\hfill \\ m\hfill & k=k_{m1}+1,k_{m1}+2,\mathrm{}\mathrm{},k_mmod(K)\hfill \end{array}$$

where $`K=k_m+1`$.

In this case, the transformed system is

$$Y(k+1)=A_m^{k_mk_{m1}}\mathrm{}\mathrm{}A_3^{k_3k_2}A_2^{k_2k_1}A_1^{k_1+1}Y(k)$$

and

$$Y(k)=X(Kk)$$

Similar asymptotic periodic properties can be established as follows.

Theorem 3

Suppose $`A_iR_e^{n\times n},i=1,2,\mathrm{}\mathrm{},m`$ are irreducible matrices, with all their main diagonal entries not equal to $`ϵ`$. Then, for any positive integers $`l_1,i=1,2,\mathrm{}\mathrm{},m`$, $`A_m^{l_m}\mathrm{}\mathrm{}A_3^{l_3}A_2^{l_2}A_1^{l_1}`$ is irreducible, too.

Theorem 4

Suppose $`A_iR_e^{n\times n},i=1,2,\mathrm{}\mathrm{},m`$ are irreducible matrices, with all their main diagonal entries not equal to $`ϵ`$. Then, there exist positive number $`\lambda `$, positive integers $`d`$ and $`k_0`$, such that the switched system (4) satisfies

$$X(k+d)=\lambda ^dX(k),kk_0$$

Note that even if the matrix $`A`$ is irreducible, its power $`A^l`$ can be reducible for some integer $`l`$. For example, let

$$A=\left[\begin{array}{ccc}ϵ& 1& ϵ\\ ϵ& ϵ& 1\\ 1& ϵ& ϵ\end{array}\right]$$

Then

$$A^3=\left[\begin{array}{ccc}3& ϵ& ϵ\\ ϵ& 3& ϵ\\ ϵ& ϵ& 3\end{array}\right]$$

which is reducible. This is why we assume that all the main diagonal entries are not equal to $`ϵ`$.

## 5  Future Research

Two issues are under investigation.

1. what is the necessary and sufficient condition for the product of some matrices to be irreducible? In some cases, even if each individual matrix is reducible, their product can still be irreducible. For example

$$A=\left[\begin{array}{cc}ϵ& 1\\ ϵ& 1\end{array}\right]B=\left[\begin{array}{cc}ϵ& ϵ\\ 1& 1\end{array}\right]$$

This issue is important in performance evaluation of switched discrete event systems.

2. how is the eigenvalue of the product of some matrices related to the the eigenvalue of each individual matrix? The eigenvalue of the product of some matrices represents the asymptotic mean period of the switched system, thereby plays an important role in performance evaluation.

Yet another interesting research direction is to study the asymptotic behavior of general 2-D discrete-event systems

$$X(m+1,n+1)=A_1X(m+1,n)A_2X(m,n+1)A_3X(m,n)$$

with the boundary condition

$$X(m,0)=X_{m0},X(0,n)=X_{0n},m,n=0,1,2,\mathrm{}\mathrm{}$$

Under what conditions does the system exhibit periodic behavior (with respect to $`m,n`$) asymptotically? and how to evaluate its asymptotic performance?

A popular model for 2-D systems is the so-called Roesser model

$$\left[\begin{array}{c}X^h(i+1,j)\\ X^v(i,j+1)\end{array}\right]=\left[\begin{array}{cc}A_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right]\left[\begin{array}{c}X^h(i,j)\\ X^v(i,j)\end{array}\right]$$

with the boundary condition

$$X^h(0,j)=X_j^h,X^v(i,0)=X_i^v,i,j=0,1,2,\mathrm{}\mathrm{}$$

How the system (in the max-plus algebra sense) evolves asymptotically, and how to evaluate its asymptotic performance are the subjects of current research.