text:
# Spin transport in diffusive superconductors

(Received 28 January 2004)

## Abstract

We employ the Keldysh formalism in the quasiclassical approximation to study transport in a diffusive superconductor. The resulting $`4\times 4`$ transport equations describe the flow of charge and energy as well as the corresponding flow of spin and spin energy. Spin-flip scattering due to magnetic impurities is included. We find that the spin-flip length is renormalized in the superconducting case and propose an experimental system to measure the spin-accumulation in a superconductor.

Manipulation of spin-polarized currents can be used to study fundamental transport processes and might also provide new functionality in electronic devices. In ferromagnets (F), the current is spin-polarized due to the spin-dependent density of states and the spin-dependent scattering potentials. In contrast, in $`s`$-wave superconductors (S), electrons with spin up and spin down and opposite momentum form Cooper pairs with no net spin. Nano-scale superconductors therefore display strikingly different properties when driven out of equilibrium by spin transport than by charge transport.

Most of the recent activities on the transport properties of F/S junctions have studied effects caused by the physical properties on the F side of the junction. The zero spin Cooper pairs prevent spin-polarized electrons to flow into S. Consequently, a spin polarized current from F injected into S can result in nonequilibrium spin accumulation near the F/S interface. The competition between electron-hole correlations and spin accumulation on the F side has recently attracted considerable interestfs . Possible influence of the ferromagnetic order parameter on the superconductor has received less attention. Singlet pairing does not allow a spin accumulation in the superconductor. Consequently, spin accumulation can reduce the superconducting gap and change the transport properties both for transport via quasiparticles and for the supercurrent. Experimentally, spin transport in diffusive S has recently been studied Gu et al. (2002). Here, the reduced quasiparticle penetration due to spin accumulation results in loss of spin memory which can be measured as a decreased magnetoresistance.

Although the theory of nonequilibrium superconductivity is widely used and developed, it has to the best of our knowledge not been completely generalized to study spin transport. In this work we thus use the Keldysh formalism and the quasiclassical approximation Belzig et al. (1999); Kopnin (2001); Rammer and Smith (1986) to rigorously obtain a set of equations describing the transport of charge and energy in a diffusive weak coupling S, as well as the transport of spin. This will describe the penetration of spins into S and the associated suppression of the superconducting order parameter. Our description of the transport properties will be based on a 4$`\times `$4 matrix equation formalism to include spin accumulation as well as electron-hole correlations. Spin-flip scattering from magnetic impurities is included as the dominant spin relaxation process inside the superconductor. We find that the spin-flip length is renormalized in the BCS case, and propose an experimental system to measure the properties resulting from the superconducting correlations. Many, but not all, experimental systems involving spin-transport in superconductors are in the elastic transport regime Tserkovnyak and Brataas (2002), which is considered here. Complementary studies based on the Boltzmann equation for spin-transport by quasiparticles in the inelastic transport regime have recently been published Mae . Note that spin-injection is qualitatively different in these opposite transport regimes due to the strong energy-dependence of quasiparticle flow in superconductors Mae .

Let us now outline the derivation of our main results. We use natural units so that $`\mathrm{}=k_\text{B}=1`$, and the electron charge is $`e=|e|`$. To describe the out-of-equilibrium electron-hole correlations as well as spin accumulation, we define the Keldysh Green’s function as

$$\widehat{G}_{ij}^\text{K}(1,2)=\underset{k}{}(\mathrm{i})\left(\widehat{\rho }_3\right)_{ik}[\left(\psi (1)\right)_k,\left(\psi ^{}(2)\right)_j]_{},$$
(1)
where $`\psi =[\psi _{},\psi _{},\psi _{}^{},\psi _{}^{}]^\text{T}`$ is a four-vector and $`\psi ^{}`$ the corresponding adjoint vector. The matrix $`\widehat{\rho }_3`$ is the third Pauli matrix generalized to 4$`\times `$4 space, $`\widehat{\rho }_3=\mathrm{diag}(1,1,1,1)`$. The coordinates are $`1=(𝐫_1,t_1)`$ and $`2=(𝐫_2,t_2)`$. Similarly, we define $`4\times 4`$ retarded and advanced Green’s functions ($`\widehat{G}^\text{R},\widehat{G}^\text{A}`$) in spin- and particle-hole space. $`4\times 4`$ matrices are denoted by a “hat” superscript. A compact notation can be obtained by construction of an 8$`\times `$8 matrix in the Keldysh space (denoted by a “check” superscript) Rammer and Smith (1986).

The quasiclassical Green’s function is defined by $`\stackrel{ˇ}{g}(𝐑,T,𝐩_\text{F},E)=\frac{\mathrm{i}}{\pi }d\xi _𝐩\stackrel{ˇ}{G}(𝐑,T,𝐩,E)`$. This function is determined by the Eilenberger equation which in the mixed representation for a stationary state can be written

$$[E\widehat{\rho }_3+\mathrm{i}\frac{𝐩}{m}\widehat{\mathbf{}}e\varphi \widehat{1}\widehat{\mathrm{\Delta }}\stackrel{ˇ}{\sigma },\stackrel{ˇ}{g}]_{}=0.$$
(2)
Here $`\widehat{\mathbf{}}=\widehat{1}\mathrm{i}e𝐀\widehat{\rho }_3`$ is the gauge invariant derivative, $`\widehat{1}`$ is the 4$`\times `$4 unit matrix, $`\varphi `$ is the electromagnetic scalar potential, $`\widehat{\mathrm{\Delta }}`$ contains the superconducting gap and $`\stackrel{ˇ}{\sigma }`$ is the self-energy due to elastic impurity scattering and spin-flip scattering by magnetic impurities in quasiclassical approximation. In the case of strong impurity scattering (dirty limit) transport is diffusive. Expansion of the quasiclassical Green’s function in spherical harmonics then gives the Usadel equations. The symmetries and normalization of the Green’s function allows for a parameterization of the quasiclassical, retarded component Belzig et al. (1999)

$$\widehat{g}_s^\text{R}=\left(\begin{array}{cc}\overline{1}\mathrm{cosh}(\theta )& \mathrm{i}\overline{\tau }_2\mathrm{sinh}(\theta )\mathrm{e}^{\mathrm{i}\chi }\\ \mathrm{i}\overline{\tau }_2\mathrm{sinh}(\theta )\mathrm{e}^{\mathrm{i}\chi }& \overline{1}\mathrm{cosh}(\theta )\end{array}\right),$$
(3)
where $`\overline{1}`$ is the 2$`\times `$2 unit matrix, $`\overline{\tau }_2`$ is the second Pauli matrix and $`\theta `$ and $`\chi `$ are position and energy dependent functions. We assume colinear magnetizations along the $`z`$-axis and $`s`$-wave singlet superconducting state. We choose a gauge where the superconducting order parameter $`\mathrm{\Delta }`$ is real and positive, and then the supercurrent is contained in the electromagnetic vector potential $`𝐀`$ and the chemical potential of the Cooper pairs is included in $`\varphi `$. Inspection of the self-consistency relation for $`\mathrm{\Delta }`$ reveals that $`\chi =0,\pi `$ depending on the boundary conditions. This ansatz simplifies the calculations considerably. The advanced Green’s function is related to the retarded through $`\widehat{g}^\text{A}=\left[\widehat{\rho }_3\widehat{g}^\text{R}\widehat{\rho }_3\right]^{}`$. Because of normalization, the Keldysh Green’s function can be expressed as $`\widehat{g}^\text{K}=\widehat{g}^\text{R}\widehat{h}\widehat{h}\widehat{g}^\text{A}`$ where $`\widehat{h}`$ is a diagonal distribution function matrix.

We will now consider a stationary state. A kinetic equation can be derived from the Usadel equations if we include Keldysh components. The important quantities are the physical particle and energy currents (including particles and holes), which we will denote by $`𝐣_\text{T}`$ and $`𝐣_\text{L}`$ respectively, with the corresponding distribution functions carrying the same indices, $`h_\text{T}`$ and $`h_\text{L}`$ Belzig et al. (1999). The physical spin current is denoted $`𝐣_{\text{TS}}`$ and the spin energy current $`𝐣_{\text{LS}}`$, with distribution functions $`h_{\text{TS}}`$ and $`h_{\text{LS}}`$. The spin-resolved distribution functions can be expressed by the particle distribution function as $`h_{\text{TS(LS)}}=(f_{}(E)f_{}(E))/2(+)(f_{}(E)f_{}(E))/2`$. The current components $`𝐣_\text{T}`$ etc. are spectral quantities, and the total charge current is given as an integral $`𝐣_{\text{charge}}(𝐫,t)=|e|N_0_{\mathrm{}}^{\mathrm{}}dE𝐣_\text{T}(𝐫,t,E)`$, and the spin current is obtained by a similar integral of $`𝐣_{\text{TS}}`$. Energy current is given by $`𝐣_{\text{energy}}(𝐫,t)=|e|N_0_{\mathrm{}}^{\mathrm{}}dEE𝐣_\text{L}(𝐫,t,E)`$, and the difference in energy current carried by opposite spins by a similar integral of $`𝐣_{\text{LS}}`$.

The equilibrium solutions for the distribution functions are $`h_{\text{L,0}}=\mathrm{tanh}\left(\beta E/2\right)`$ and $`h_{\text{T,0}}=h_{\text{LS,0}}=h_{\text{TS,0}}=0`$. We derive kinetic equations and find,

$`𝐣_\text{L}`$ $`=0,`$ (4a)
$`𝐣_\text{T}`$ $`=2|\mathrm{\Delta }|\alpha _{\text{TT}}h_\text{T},`$ (4b)
$`𝐣_{\text{LS}}`$ $`=\left(2|\mathrm{\Delta }|\alpha _{\text{TT}}+{\displaystyle \frac{1}{\tau _{\text{sf}}}}\alpha _{\text{LSLS}}\right)h_{\text{LS}},`$ (4c)
$`𝐣_{\text{TS}}`$ $`={\displaystyle \frac{1}{\tau _{\text{sf}}}}\alpha _{\text{TSTS}}h_{\text{TS}}.`$ (4d)
The right-hand side terms represent renormalized scattering because of superconductivity:

$`\alpha _{\text{TT}}=`$ $`\text{Im}\left\{\mathrm{sinh}(\theta )\right\},`$ (5a)
$`\alpha _{\text{LSLS}}=`$ $`\left(\text{Re}\left\{\mathrm{cosh}(\theta )\right\}\right)^2\left(\text{Im}\left\{\mathrm{sinh}(\theta )\right\}\right)^2,`$ (5b)
$`\alpha _{\text{TSTS}}=`$ $`\left(\text{Re}\left\{\mathrm{cosh}(\theta )\right\}\right)^2+\left(\text{Re}\left\{\mathrm{sinh}(\theta )\right\}\right)^2.`$ (5c)
The $`|\mathrm{\Delta }|\alpha _{\text{TT}}`$ terms on the right hand side in (4) are due to conversion of quasiparticle current into supercurrent, and the $`\alpha _{\text{LSLS}}/\tau _{\text{sf}},\alpha _{\text{TSTS}}/\tau _{\text{sf}}`$ terms are due to spin-flips. The spin-flip time in the normal state is $`\tau _{\text{sf}}^1=8\pi n_{\text{sf}}N_0S(S+1)|v_{\text{sf}}|^2/3`$, where $`n_{\text{sf}}`$ is the magnetic impurity density, $`N_0`$ the density of states at the Fermi level, $`S`$ the impurity spin quantum number and $`v_{\text{sf}}`$ is the Fourier transformed spin-flip impurity potential. We assume isotropic scattering. Our definition of $`\tau _{\text{sf}}`$ differs from the usual spin-flip lifetime by a renormalization factor $`4/3`$. This definition reproduces the diffusion equation with a spin-flip length $`l_{\text{sf}}^{\text{(N)}}=\sqrt{D\tau _{\text{sf}}}`$ in the normal state. Thus there is a difference between the spin-flip lifetime measured in e.g. electron spin resonance and spin-flip transport time.

We introduce generalized energy-dependent diffusion coefficients

$`D_\text{L}=`$ $`D\left[\left(\text{Re}\left\{\mathrm{cosh}(\theta )\right\}\right)^2\left(\text{Re}\left\{\mathrm{sinh}(\theta )\right\}\right)^2\right],`$ (6a)
$`D_\text{T}=`$ $`D\left[\left(\text{Re}\left\{\mathrm{cosh}(\theta )\right\}\right)^2+\left(\text{Im}\left\{\mathrm{sinh}(\theta )\right\}\right)^2\right],`$ (6b)
where $`D=\tau v_\text{F}^2/3`$ is the diffusion constant. The currents can then be expressed as

$`𝐣_\text{L}=`$ $`D_\text{L}h_\text{L}+\text{Im}\left\{𝐣_\text{E}\right\}h_\text{T},`$ (7a)
$`𝐣_\text{T}=`$ $`D_\text{T}h_\text{T}+\text{Im}\left\{𝐣_\text{E}\right\}h_\text{L},`$ (7b)
$`𝐣_{\text{LS}}=`$ $`D_\text{T}h_{\text{LS}}+\text{Im}\left\{𝐣_\text{E}\right\}h_{\text{TS}},`$ (7c)
$`𝐣_{\text{TS}}=`$ $`D_\text{L}h_{\text{TS}}+\text{Im}\left\{𝐣_\text{E}\right\}h_{\text{LS}}.`$ (7d)
Here we have defined the spectral supercurrent as $`𝐣_E=D(\chi 2e𝐀)\mathrm{sinh}^2(\theta )`$. The self-consistency relation is

$`\mathrm{\Delta }(𝐫)=`$ $`{\displaystyle \frac{1}{2}}\mathrm{sgn}(\mathrm{\Delta }_0)N_0\lambda {\displaystyle _{\mathrm{}}^{\mathrm{}}}dE\mathrm{sinh}\left(\theta \right)h_\text{L},`$ (8)
where the factor $`\mathrm{sgn}(\mathrm{\Delta }_0)`$ is determined from the boundary condition to give the correct sign and $`\lambda `$ is the interaction parameter. The complex part of this equation is neglected as a consequence of charge conservation Schmid (1981).

The functions $`\theta `$ and $`\chi `$ are determined by the retarded components of the Usadel equation. We obtain

$`𝐣_E=`$ $`0,`$ (9)
$`D\left[^2\theta {\displaystyle \frac{1}{2}}(\chi 2e𝐀)^2\mathrm{sinh}(2\theta )\right]=`$ $`2\mathrm{i}E\mathrm{sinh}(\theta )`$
$`2\mathrm{i}\mathrm{cosh}(\theta )|\mathrm{\Delta }|+{\displaystyle \frac{3}{4}}`$ $`{\displaystyle \frac{1}{\tau _{\text{sf}}}}\mathrm{sinh}(2\theta ),`$ (10)
where Eq. (9) implies that the spectral supercurrent is conserved. In addition we have the following symmetry conditions, $`\theta ^{}(E)=\theta (E),\chi ^{}(E)=\chi (E)`$. Equations (4)-(Spin transport in diffusive superconductors) determine all transport properties of S.

In general, in a hybrid F/S system, the superconductor cannot be described as in terms of BCS-formulas close to the F/S interface due to the proximity effect. Nevertheless, to gain insight into the physics implied by the above-mentioned formulas let us now consider the limit of a homogeneous BCS superconductor, and select $`\chi =0`$. This is relevant for the proposed experiment below. For energies $`|E|<|\mathrm{\Delta }|`$ $`\alpha _{\text{TT}}=\mathrm{\Delta }/\sqrt{\mathrm{\Delta }^2E^2}`$ and the spin-flip renormalization factors are $`\alpha _{\text{TSTS}}=0,\alpha _{\text{LSLS}}=\mathrm{\Delta }^2/(\mathrm{\Delta }^2E^2)`$. The generalized diffusion constant $`D_\text{L}=0`$ while $`D_\text{T}=D\mathrm{\Delta }^2/(\mathrm{\Delta }^2E^2)`$. From Eq. (7a) this means that there is no energy current carried by quasiparticles with energy $`|E|<|\mathrm{\Delta }|`$. Gap scattering for quasiparticle energies below the superconducting gap corresponds to a transformation of the charge current ($`𝐣_\text{T}`$) into supercurrent. Such scattering is not possible for the physical spin current ($`𝐣_{\text{TS}}`$). Consequently, in the absence of spin-flip scattering the quasiparticle spin-current into the superconductor vanishes for $`|E|<|\mathrm{\Delta }|`$ since $`D_\text{L}=0`$ in Eq. (7d) and $`\alpha _{\text{TSTS}}=0`$ in the kinetic equation (4d). Note that this result relies on the fact that there are different effective diffusion coefficients for charge current ($`D_\text{T}`$) and for spin current ($`D_\text{L}`$). We also observe that the term $`\alpha _{\text{LSLS}}`$ is negative below the gap, acting as a source of spin energy.

Above the gap $`(|E|>|\mathrm{\Delta }|)`$ the factor $`\alpha _{\text{TT}}`$ vanishes while $`\alpha _{\text{LSLS}}=E^2/(E^2\mathrm{\Delta }^2),\alpha _{\text{TSTS}}=(E^2+\mathrm{\Delta }^2)/(E^2\mathrm{\Delta }^2)`$. For the generalized diffusion coefficients we find that $`D_\text{L}=D`$ and $`D_\text{T}=DE^2/(E^2\mathrm{\Delta }^2)`$. Now consider the kinetic equations in the BCS case. A charge current carried by quasiparticles with energy $`|E|>|\mathrm{\Delta }|`$ can propagate into S. For quasiparticles at $`|E|>|\mathrm{\Delta }|`$ we see that there is no renormalization for the spin-energy diffusion length in Eq. (4c), whereas the spin diffusion length in Equation (4d) has an energy dependent renormalization factor which diverges for energies $`|E|=|\mathrm{\Delta }|`$ causing massive spin-flip scattering.

We will now apply this formalism to study spin diffusion, and demonstrate the significance of the renormalization of the spin diffusion length. Experimental studies of spin accumulation and spin injection has recently been performed spi  in metallic spin valves. The spin accumulation in the physically different inelastic regime for a superconductor in this experimental system has also been calculated theoretically Mae . We will consider the simplified geometry shown in Figure 1, where there is no charge transport in the superconductor, and calculate the spin accumulation signal in the elastic regime.

The F1/N/F2 systems act as a spin-battery which is connected via a tunnel junction to the superconductor. A voltage bias between F1 and F2 induces a spin accumulation that can flow into S. The superconducting wire is connected to an S reservoir in equilibrium BCS state by a good metallic contact at distance $`L^{(\text{S})}`$ from the N/S interface. On top of the S wire there is a ferromagnet connected by tunnel barrier which upon switching of the magnetization direction acts as a detector for the spin signal. Measurement of the relative voltage of this electrode between parallel and antiparallel (with respect to the top F reservoir) magnetization gives $`\mathrm{\Delta }\mu =\mu ^{(\text{P})}\mu ^{(\text{AP})}`$ which describes the difference between electrochemical potential of spin-up and spin-down quasiparticles located a distance $`L`$ from the N/S interface. This quantity can be calculated $`\mathrm{\Delta }\mu =_{\mathrm{}}^{\mathrm{}}dEP^{(\text{D})}h_{\text{TS}}(L,E)`$, where $`P^{(\text{D})}`$ is the spin polarization of the tunnel barrier between S and the F detector. We assume a homogeneous order parameter and BCS spectral properties in the S wire since there are tunnel barriers between the N,F and S elements and perturbation from current and spin-flip is weak.

We can express the difference between the spin-up and spin-down distribution functions in N close to S as $`\mathrm{\Delta }f^{\text{(N)}}f_{}^{\text{(N)}}f_{}^{\text{(N)}}=P^{(\text{FN})}(f(EeV/2)f(E+eV/2))`$, where $`P^{(\text{FN})}=\left(G_{\text{maj}}G_{\text{min}}\right)/\left(G_{\text{maj}}+G_{\text{min}}\right)`$ is the spin polarization between the F reservoirs and N, $`f(E\pm eV/2)`$ is the Fermi-Dirac distributions in the F reservoirs and $`G_{\text{maj(min)}}`$ is the conductance of majority (minority) spin electrons from ferromagnetic reservoir to the middle of N. There is thus no charge current or supercurrent anywhere in S, however there may be a spin-current. Equation (7d) states that there is no spin-current for energies below the gap, thus for these energies the N/S interface is effectively insulating. Since the S wire is connected to a reservoir in the other end for $`|E|<\mathrm{\Delta }`$ the spin distribution function equals the equilibrium value $`h_{\text{TS}}=0`$. We solve the TS kinetic equation (4d) for energies $`|E|>\mathrm{\Delta }`$. This equation reduces to a diffusion equation with renormalized spin-flip length $`l_{\text{sf}}^{(\text{S})}(E)=l_{\text{sf}}\sqrt{(E^2\mathrm{\Delta }^2)/(E^2+\mathrm{\Delta }^2)}`$, where $`l_{\text{sf}}=\sqrt{D\tau _{\text{sf}}}`$ is the normal state spin-flip length. The boundary condition at the S reservoir is that the distribution function attains the equilibrium value, and at the S/N interface we match at each energy the tunnel spin current to the spin current inside S, $`|e|N_0𝐣_{\text{TS}}`$. We assume that $`L^{(\text{S})}/l_{\text{sf}}^{(\text{S})}1`$ which is a relevant physical situation.

The position and energy dependent solution $`h_{\text{TS}}`$ is substituted into the expression for the measured difference in electrochemical potential for parallel and antiparallel configuration, and we obtain

$`\mathrm{\Delta }\mu =2P^{(\text{D})}{\displaystyle _\mathrm{\Delta }^{\mathrm{}}}dE\mathrm{\Delta }f^{\text{(N)}}\mathrm{e}^{L/l_{\text{sf}}^{(\text{S})}}{\displaystyle \frac{R_{\text{sf}}^{(\text{S})}}{R_{\text{sf}}^{(\text{S})}+R^{(\text{I})}}},`$ (11)
where $`R^{(\text{I})}(E)=1/(|T|^2N_{\text{BCS}}(E)N_0)`$ is the resistance of the N/S tunnel barrier, $`|T|`$ is the tunneling matrix element, $`N_{\text{BCS}}(E)`$ is the BCS density of states, $`R_{\text{sf}}^{(\text{S})}(E)=l_{\text{sf}}^{(\text{S})}(E)\rho /A`$ is the resistance of the S wire within a spin-flip length and $`\rho `$ is the resistivity of the material in S when in the normal state ($`T>T_c`$). This result can be understood as follows. The spin accumulation close to the tunnel interface is exponentially attenuated by spin-flip scattering in S. The spin signal is also decreased by the tunnel resistance, and since spin current is strongly suppressed for energies $`|E|<\mathrm{\Delta }`$ only quasiparticles with energies higher than the gap contribute. The effective total resistance is a series of the tunnel interface resistance with the resistance of S within one spin-flip length.

We will now consider some simplified limits for the quantity $`\mathrm{\Delta }\mu `$ defined above. In the normal state where $`\mathrm{\Delta }0`$ we find that $`\mathrm{\Delta }\mu /eV=2P^{\text{(D)}}P^{\text{(FN)}}R_{\text{sf}}^{\text{(S)}}\mathrm{exp}(L/l_{\text{sf}})/(R_{\text{sf}}^{\text{(S)}}+R^{\text{(I)}})`$ where $`R_{\text{sf}}^{\text{(S)}}`$ and $`R^{\text{(I)}}`$ assume their normal state (energy independent) values. At $`k_\text{B}T\mathrm{\Delta }`$ the signal measured by $`\mathrm{\Delta }\mu `$ vanishes when the bias is lower than the energy gap $`eV<\mathrm{\Delta }`$ since spin current is suppressed for quasiparticles below the gap. For higher bias, $`eV>\mathrm{\Delta }`$, and at zero temperature when the bulk resistance dominates, $`R_{\text{sf}}^{(\text{S})}R^{(\text{I})},`$ an approximate solution is $`\mathrm{\Delta }\mu =2P^{\text{(D)}}P^{\text{(FN)}}\mathrm{\Delta }\mathrm{e}^{L/l_{\text{sf}}}\{\mathrm{e}^{Lr^2/2l_{\text{sf}}}/r\mathrm{e}^{L/2l_{\text{sf}}}+\sqrt{\pi L/2l_{\text{sf}}}(\mathrm{erf}[r\sqrt{L/2l_{\text{sf}}}]\mathrm{erf}[\sqrt{L/2l_{\text{sf}}}])\},`$ where $`r=2\mathrm{\Delta }/eV`$. In this case the relation between the energy gap and the bias determines the magnitude of the spin signal, and the exponential decrease of the signal.

The temperature dependence of $`\mathrm{\Delta }\mu `$ in the general case is given by decrease from a constant value above $`T_c`$ as the temperature approaches zero. An example of this behavior is shown in Figure 2. Here we have used the approximate temperature dependence $`\mathrm{\Delta }=1.76T_c\mathrm{tanh}(1.74\sqrt{T_c/T1})`$.

Our calculations show that the spin signal decreases due to superconducting correlations. For a large energy gap the spin accumulation vanishes completely at low temperatures. These effects can be explained by suppressed subgap spin-current and massive spin-flip at energies close to the gap because of the superconducting correlations.

In conclusion, we have presented a formalism to describe elastic spin transport in superconductors with spin-flip scattering. We find different effective diffusion coefficients for charge- and spin-current. The spin-flip length is renormalized in the superconducting case, and at energies close to the gap there is massive spin-flip. As an illustration we compute the difference in electrochemical potential due to spin accumulation in an experiment sensitive to the renormalization of spin-flip length.

###### Acknowledgements.

This work was supported in part by The Research Council of Norway, NANOMAT Grants No. 158518/431 and 158547/431, RTN Spintronics, the Swiss NSF and the NCCR Nanocience.


format: text

text:
# Front-End electronics and integration of ATLAS pixel modules

## 1 Introduction

The Pixel Detector of the ATLAS experiment at the LHC, which is currently under construction at CERN, Geneva, is the crucial part for good secondary vertex resolution and high b-tagging capability. Therefore high spatial resolution and fast read-out in a multi-hit environment is needed and leads to the use of a hybrid silicon pixel detector. This detector will consist of three barrel layers and three disks in each forward direction to have at least three space points up to $`|\eta |=2.5`$ for particles with a transversal momentum greater than $`500`$ MeV .

The smallest detector unit will be a module made of one silicon pixel sensor and sixteen readout chips connected with high density bump bonding techniques to the sensor. Each module consists of about $`50,000`$ pixel cells with the size of $`50\mu m400\mu m`$ to reach the required spatial resolution of $`12\mu m`$ in $`r\varphi `$-direction. The main requirement for the ATLAS pixel modules is high radiation tolerance of all components in a harsch radiation environment close to the interaction point. The integrated design fluence is $`10^{15}cm^2`$ 1 MeV neutron equivalent coming predominantly from charged hadrons over ten years for the outer barrel layers and disks and five years for the innermost layer. Secondly the total need of about $`1,800`$ modules requires very good testability of all components before assembly and high fault tolerance in long term operation as they are not supposed to be exchanged during the whole foreseen ten years lifetime of the ATLAS experiment. The material budget for the whole pixel detector is very strict to affect the later detector parts as less as possible. This leads to a total amount of material less than $`1.2\%`$ of one radiation length per module not including further services and support structures.

## 2 Module concept

A cross-section of an ATLAS pixel module can be seen in figure 1. The module basically consists of a so called bare module which meant the sixteen readout chips bump bonded to the silicon pixel sensor. The size of the module is roughly given by the size of the sensor $`18.663.0`$mm2. The readout chips with are placed in two rows of eight. The front-end chips are slightly longer than the sensor width to be able to reach the wire-bond pads of the chips. The interconnection techniques using fine pitch bump bonding is either done with Pb/Sn by IZM111Institut für Zuverlässigkeit und Mikrointegration, Berlin, Germany. or with Indium by AMS222Alenia Marconi Systems, Roma, Italy.. The total active area of one module is $`16.460.8`$mm2 with $`46080`$ pixel.

For routing the signals to the MCC  and further to off-module electronics and providing the analog and digital supply powers to the readout chips a two layer high density interconnect (flexible capton cable) is used. Also the flexible hybrid carries the MCC chip which performs control and clock operations. It is glued to backside of the sensor and has roughly the size of the sensor tile. To connect the flex hybrid with the Front-End chips and with the MCC standard wire bond techniques are used; altogether around 350 wire-bonds are needed for a complete module. The module is glued chips down to the mechanical support which also serves as cooling structure. The whole power consumption of about $`4`$W per module has to be dissipated through the readout chips to allow an operation of the module at $`6^{}`$C.

## 3 Front-End electronics

Each FE-chip contains of 2880 pixel cells arranged in 18 columns times 160 rows. Every pixel cell consists a charge sensitive amplifier, a discriminator for zero suppression and a digital readout logic to transport the hit information to the periphery. The signal charge can be determined by measuring the width of the discriminator output signal. Hits are temporarily stored in one of the 64 buffers per column pair located at the chip periphery until a trigger signal selects them for readout. Hits with no corresponding trigger signal are discarded after the programmable trigger latency of up to 256 clock cycles of 25 ns. All hits corresponding to trigger are sent through one serial LVDS link to module controller chip which builds full modules events containing hit information of all 16 FE chips on a module.

Figure 2 shows a block diagram of most important elements in the pixel cell. Charge deposited in the sensor is brought to the pixel through the bump bond pad which connects to the input of a charge sensitive amplifier. An inverting folded cascode amplifier is fed back by a capacitor of $`C_f=6.5`$fF integrated into the lower metal layers of the bump bond pad. The feedback capacitor is discharged by a constant current so that a triangular pulse shape is obtained. As consequence of this particular pulse shape, the width of the discriminator output signal is nearly proportional to the deposited charge. This signal is called ’time over threshold’ or ToT, is measured in units of the bunch crossing clocks (25 ns) and provides an analogue information of evry hit with a resolution of 4-6 bit. The feedback current $`I_f`$ is set globally with an 8 bit on-chip bias DAC so that a compromise between a good ToT resolution and small dead time in the pixel can be found. The constant feedback current generation is part of a circuit which compensates for detector leakage current. This must be sourced by the electronics because the pixel sensor is dc-coupled to readout chip. The compensation circuit can cope with more than 100 nA detector leakage current per pixel.

The threshold generation of each discriminator is done within every single pixel. A small threshold dispersion between the 2880 pixel is crucial in order to achieve low threshold settings without an increased noise hit rate. Therefore a number of 7 individually programmable trim bits in each pixel covers a wide range of thresholds between $`1,000e^{}`$ and $`10,000e^{}`$. These local trim bits affects both side of a differential pair amplifier leading to a linear threshold vs. DAC behaviour in the central part of the covered threshold range. Furthermore a globally set 5 bit DAC, so called GDAC, allows to move the threshold for the whole chip without losing the triming between the individual pixel.

Every pixel contains a calibration circuit used to inject known charges into the input node. This circuit generates a voltage step which is connected to one of two injection capacitors for different charge ranges.

The output signal of the discriminator in every pixel can be disabled with a local mask bit. A fast OR of a programmable subset of pixel is available for testing purposes and for self-triggered operation of the chip by using this OR signal after delaying it as a trigger. The 14 configuration bits used in every pixel are stored in SEU tolerant DICE cells  which can be written and read back by means of a shift register running vertically through the column. The layout of this cell has been significantly changed with respect to the previous chip (FE-I1) in order to improve the SEU tolerance. For the globally set DACs a basic triple redundancy scheme with majority readback is used.

The basic philosophy of the readout is to associate the hit to unique bunch crossing by recording a time stamp when the rising edge of the discriminator occurs. The time of the falling edge is memorized as well in order to calculate the ToT as the difference of the two values in units of the bunch crossing clock. The two 8 bit RAM cells used for this purpose are classical static memory cells.

The readout scheme can be divided into four elementary tasks running on the chip in parallel. First time stamps generated by an 8 bit Gray counter and clocked with bunch crossing clock are stored for the rising and falling edge of the discriminators in the pixel. After receiving the falling edge a hit flag is set and the pixel is ready to be processed. Secondly as soon as hits are flagged to the Column Control Logic the uppermost hit pixel is requested to send its rising and falling edge time stamps together with its ID down the column where is information is stored in a free location of the End-of-Column (EoC) buffer pool. The pixel is then cleared and the scan continues the search for other hit pixel. Thirdly the hit stays in the EoC buffer until the trigger latency has elapsed. The leading edge time stamp of the hits in the buffer are therefore permanently compared to the actual time stamp minus the fix but programmable latency. When the comparison is true and a trigger signal is present the hit is flagged as valid for readout, it is discarded otherwise. Incoming triggers are counted on the chip, the trigger number is stored together with the flag in the EoC buffer so that several trigger do not lead to confusion. A list of pending triggers is kept in a FIFO. Lastly a Readout Controller initiates the serial readout of hit data as soon as pending triggers are present in the Trigger FIFO, the EoC buffers are searched for valid data with correct trigger numbers. The column and row address and the ToT of these hits are serialized and send to the MCC. The Readout Controller adds a start-of-event and an end-of-event word together with error and status bits to the data stream.

The recent chip (FE-I2) has been designed in a quarter micron technology with 6 metal layers using radiation tolerant layout rules . Tests of the first batches received back from the vendor from May 2003 on showed that the design is fully functional except one problem dealing with race condition in the control block leading to the need to operate the chip with a reduced digital supply voltage of 1.6 V instead of the nominal 2.0 V. But this problem has been fixed via a re-distribution of the clock signal by changing 2 metal layers slightly in the back-end processing (FE-I2.1).

## 4 Results of prototype modules

Several ten modules have been built with older chip generation (FE-I1) in order to qualify the bumping process and the whole module production chain. As an example of the overall good quality of these modules one can see in figure 3 the hitmap and spectrum obtained with one module illuminated by an Am241-source from the top using the self-trigger possibility of the chip. Only less than 10 out of $`46,080`$ pixel don’t see source hits showing the excellent bump quality of this module which is true for almost all modules. For every channel the ToT response was calibrated individually using the integrated charge generation circuit. Afterwards this calibration was applied to the source measurement in order to obtain the shown spectrum which is a sum over all pixel. Since no clustering of hit data was applied the spectrum is relatively broad but one can clearly see the main $`60`$keV photo-peak roughly at the expected value of $`16,600e^{}`$. More results obtained with FE-I1 modules can be found in .

Recently first modules with the newer chip generation (FE-I2.1) has been built. Figure 4 and 6 show the threshold distribution and noise measured for a whole module after the trimming of the individual pixel to a threshold of roughly $`3,100e^{}`$. All these measurements were done with internal charge injection circuit. One can see that the reached overall threshold dispersion is only $`40e^{}`$ which has to be compared to value of $`100e^{}`$ reached by the older chip generation. Note further that there is no single pixel with an threshold lower than $`2,900e^{}`$. This shows the high tuning capability of this chip which is important to reach small thresholds on the whole module without any extra noisy pixel especially after irradiation.

The measured mean noise of the module is with $`150e^{}`$ for the standard pixel much better than the specification of $`300e^{}`$, even the longer and so called ganged pixel located in the inter-chip regions which show due to the higher input capacity a higher noise of $`170e^{}`$ and $`270e^{}`$ resp. are well below this specification.

The most demanding requirement of LHC is the time correct hit detection within one bunch crossing cycle (25 ns) together with the strict power budget of less than 4 W for the whole 16 chip module. Therefore not the real threshold is important but the threshold taken only hits into account which arrive within in time slot of 20 ns after the injection. Using a special delay circuit inside the MCC measurements with the needed accuracy are possible and in figure 6 the results are shown for the same module tuned to a threshold of $`3,100e^{}`$. Because the timing is very sensitive to the input capacity of the amplifier the in-time-threshold is slightly different for the different kind of pixel. But overall a in-time-threshold of $`4,200e^{}`$ is reached and also the special pixel don’t show a threshold higher than $`4,600e^{}`$ meaning that only an extra charge of $`1,100e^{}`$ and $`1,500e^{}`$ resp. is needed to meet the timing requirement of ATLAS. This results are very much improved with respect to the previous chip generation where the in-time-threshold was in the order of $`5,000`$-$`6,000e^{}`$ by a tuned threshold of $`3,000e^{}`$. This exceeds specification of $`5,000e^{}`$ given by the fact that after irradiation damage only half of the initial charge of $`20,000e^{}`$ is available.

More extensive studies have been done with these modules, including beam test measurements and measurements during and after proton irradiation up to 100 MRad or $`210^{15}`$neqcm-2 well above the design fluence of 50 MRad and $`110^{15}`$neqcm-2 respectively. All these measurement confirmed the good quality of the chip and the whole module, in particular it was shown that these modules meet all requirements and are able to operate in the environment of ATLAS and LHC.

## 5 Conclusions

The recent ATLAS pixel Frond-End chip generation has been successfully produced, assembled to full size modules and tested. Improvements concerning the timing behaviour and chip threshold tuning capability with respect to the previous chip generation has been obtained. All further results including beam tests and irradiation confirmed the overall good quality of the chip reaching all requirements. By building more than 100 modules with high quality a lot of experience has been gained with the fine pitch bump bonding process for pixel module using two different vendors and technologies showing that the target of building $`2,000`$ modules needed for ATLAS in the next two years is feasible.


text:
# Wetting phenomenon in the liquid-vapor phase coexistence of a partially miscible Lennard-Jones binary mixtures

## I Introduction

Wetting phenomena occurs very often in daily life and has a strong interdisciplinary character. It is of great relevance for fundamental areas of knowledge such as physics koshevnikov97 ; physics , chemistry chemistry , biology biology0 , and several other applied sciences metallurgy  as well as technology technology0 . The occurrence of wetting phenomenon is usually associated with the existence of three phases, at least one of which is liquid and no more than one phase is solid. In most practical situations the solid phase is wet by the liquid phase and the disfavored phase is the vapor. This is expressed in terms of surface free energies as, $`\gamma _{SL}<\gamma _{SV}`$. Although the vapor phase is disfavored, a solid-vapor coexistence can be possible and in such case it is said that the liquid ”partially wets” the solid. However, as the liquid phase is further adsorbed by the solid, it is possible that the contact of the vapor with the solid be excluded by the creation of a liquid layer between them. In such a situation one says that the liquid “completely wets” the solid. The surface free energies involved in the wetting process are related by Young’s rule, $`\gamma _{SV}=\gamma _{SL}+\gamma _{LV}\mathrm{cos}\theta `$, with $`\theta `$ the contact angle. For partial wetting, $`0<\theta <90^0`$, and for complete wetting, $`\theta =0`$. Usually the presence of a wall in most of the theoretical as well as experimental studies complicates a detailed and precise investigation of the interfacial structures. This is due to the fact that the atomic interactions between the wall and the fluid components play an important role in the wetting phenomena. Fortunately, binary liquid mixtures offer a great opportunity to investigate in detail wetting phenomena since they involve only fluid phases in coexistence. An important number of experimental studies of interfacial wetting have studied only a small subset of possible binary mixtures koshevnikov97 ; kwon82 ; houessou85 ; kahlweit88 ; bonn92 ; kellay93 . Wetting and prewetting phenomena in which one of the phases is solid has also been studied by means of Monte Carlo simulations using a 3D lattice-gas model nicolaides89 . A 12-6 LJ potential modeling argon and a 9-3 LJ potential modeling a CO2 covered solid wall have been chosen to study the wetting transitions by means of Monte Carlo simulations finn89 . However, to the best of our knowledge, there are not numerical simulations that explicitly study the interfacial behavior and the associated wetting phenomena in fluid mixtures in absence of a wall. In this paper we apply a well established methodology using MD simulations chapela ; alejandre ; mecke ; nosotros ; holc93  to study interfacial properties and surface phenomena in a partially miscible fluid mixture. We consider a model binary LJ fluid mixture in which the attractive part of the A-B interactions is weaker than the A-A and B-B interactions. By studying the properties of the density profiles of the fluid phases as a function of temperature and miscibility parameter, we show clear evidence that the vapor phase spontaneously excludes the liquid A-liquid B interface at and above the wetting transition temperature, $`T_\text{w}^{}`$. That is, the vapor phase fully wets the liquid A-liquid B phases. This phenomenon occurs when the miscibility parameter $`\alpha `$ is in the range $`0<\alpha 0.5`$ and in the temperature region $`T_W^{}(\alpha )<T^{}<T_{\text{cons}}^{}(\alpha )`$. Here, $`T^{}=(k_BT)/ϵ`$, is the reduced temperature, $`T_\mathrm{w}^{}(\alpha )`$, is the wetting transition temperature (WTT), and $`T_{\text{cons}}^{}(\alpha )`$, is the consolute temperature of the mixture. By analyzing the structural properties of the mixture as a function of $`\alpha `$ we estimate the wetting phase diagram $`T_\mathrm{w}^{}(\alpha )`$. A further quantitative analysis, for $`\alpha =0.25`$, of the surface free energy $`\gamma (T^{})`$ and the adsorption, $`\mathrm{\Gamma }(T^{})`$, at the LL interface indicates that the wetting transition is of first order.

The layout of the remainder of this paper is as follows. In section II we introduce the model of the symmetric binary mixture, in section III we explain the details of the simulations, then in section IV we present and discuss the representative results of the extensive MD simulations. Finally, we end with the conclusions in section V.

## II The model

The model binary mixture studied in this paper consists of fluids A and B made up of spherical molecules of the same size, $`\sigma _{AA}=\sigma _{BB}`$, and at concentrations of 50% each. The interaction between molecules of the same type is through the classical 12-6 LJ potential. However, the interaction between molecules of fluids A and B is given by the LJ potential,

$$u_{ij}(r)=\{\begin{array}{cc}4ϵ_{ij}\left[\left(\frac{\sigma _{ij}}{r_{ij}}\right)^{12}\left(\frac{\sigma _{ij}}{r_{ij}}\right)^6\right],\hfill & \text{if}rR_c=3\sigma _{AA},\hfill \\ 0,\hfill & \text{if}r>R_c=3\sigma _{AA},\hfill \end{array}$$
(1)
where the mixing rule is defined by

$`\sigma _{AB}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\sigma _{AA}+\sigma _{BB}),`$ (2)
$`ϵ_{AB}`$ $`=`$ $`\alpha _{AB}ϵ_{AA},`$ (3)
with $`ϵ_{AA}=ϵ_{BB}`$, and $`\alpha _{AB}`$ the parameter that controls the miscibility of the two fluids. For the sake of simplicity from now on we will use the short notation, $`\alpha _{AB}=\alpha `$. Notice that when, $`\alpha =0`$, we obtain two independent single LJ fluids while in the opposite case, $`\alpha =1`$, the system reduces to a single LJ fluid. By choosing, $`0<\alpha <1`$, the attractive part of the A-B interactions becomes weaker than that of the AA and BB interactions, and then, the liquid phases are immiscible in a wide range of temperatures. Thus, one can obtain the coexistence of three fluid phases: liquid A-liquid B and the vapor.

## III Details of the simulations

We have carried out extensive MD simulations to investigate the structural properties of this model binary mixture as a function of $`\alpha `$, varying this quantity in steps, $`\delta \alpha =0.05`$, in the interval $`0.2\alpha 0.5`$. For $`\alpha =0.25`$ we quantitatively related the interfacial properties with the corresponding phase diagram properties. In all the simulations we applied periodic boundary conditions along the $`x,y`$ and $`z`$ directions. At the lowest temperature, $`T^{}=0.65`$, the simulations were initiated from a configuration where the molecules of type A and B form two contiguous FCC crystals. At higher temperatures we take as the initial configuration the final configuration of the previous temperature. The initial velocities of the molecules were assigned from a Maxwell-Boltzmann distribution. The equations of motion were integrated using a leap frog algorithm with a time step size, $`\delta t^{}=0.005`$. This corresponds to $`1.1\times 10^5`$ nanoseconds in the scale of argon. At each time step iteration we monitor the temperature of the system, by means of the equipartition theorem, and rescale the linear momentum of the molecules to keep the temperature constant. This method of rescaling the linear momenta is known as the isokinetic thermostat. To check that this isokinetic thermostat produces meaningful results we have also carried out some MD simulations using the Nosé-Hoover thermostat. In figure 1 we show the kinetic energy distributions as the dynamics of the system evolves at, $`T^{}=0.90`$, applying both methods. As one would expect both thermostats yield a Gaussian distribution of kinetic energies. The only difference is the width of the distributions. The isokinetic thermostat produced a distribution with a variance that is one order of magnitude smaller than the variance of the distribution obtained with the Nosé-Hoover thermostat. As a further check of the isokinetic thermostat we also calculated and monitored some thermodynamic and surface quantities applying both thermostats. The results of this comparison test showed full consistency. Since MD simulations using the Nosé-Hoover thermostat are computationally more demanding we used the isokinetic thermostat in all the simulations reported in this paper. Furthermore, to check the stability of the interfaces as well as the distribution of the species, we simulated the system for as long as 55 ns in the scale of Argon. Thermodynamic quantities and interfacial properties of interest were measured by averaging over the last million of time-step iterations. To minimize correlations between measurements we calculated thermodynamic, structural and surface quantities every 50 time steps. We also investigated the role of finite size effects for the value of the miscibility parameter, $`\alpha =0.25`$. To this end we carried out MD simulations with three system sizes, N=$`1728,\mathrm{\hspace{0.25em}4096}`$ and 6144 molecules. We found that for all the quantities studied here, simulations with $`N=4096`$ provided reliable results. Therefore, most of the simulations were carried out with $`N=4096`$ molecules. The discussion of these results will be presented, where appropriate, in the next section. On the other hand, interfacial properties are sensitive to the cross section area of the simulational box that is parallel to the interfaces, as discussed in previous MD simulations of the LV interface of a single LJ fluid holc93 ; li95 . These authors arrived to the conclusion that a reliable value of the cross section area of the computational box should be at least $`(8\sigma )^2`$. So, to be on the safe side, in the present simulations we have considered a computational box with a cross section area, $`Lx\times Ly=(9\sigma _{AA})^2`$. The length, $`L_z`$, of the simulational parallelepiped was adjusted such that the average density of the system laid somewhere inside the LV coexistence curve. In this way one readily gets the liquid-vapor phase coexistence. The average densities of the simulated systems were in the range $`0.2\rho ^{}0.4`$, where the reduced density defined as, $`\rho ^{}=\rho \sigma _{\mathrm{AA}}^3`$. In the following section we present, analyze and discuss the results of the thermodynamic and interfacial properties calculated from our MD simulations.

## IV Results and discussion

### IV.1 Structural properties and phase diagrams

We performed extensive MD simulations for mixtures with $`0.2\alpha 0.65`$ changing this parameter in steps of, $`\delta \alpha =0.05`$. We studied the density profiles, $`\rho ^{}(z)`$, of the liquid-vapor coexistence in different temperature regions. From an analysis of $`\rho ^{}(z)`$, as function of $`T^{}`$, and, $`\alpha `$, we estimated the $`T_\mathrm{W}`$ versus $`\alpha `$ phase diagram. In what follows we present some representative results for, $`\rho ^{}(z)`$, when $`\alpha =0.25`$ and 0.30. Then, for the particular value, $`\alpha =0.25`$, we present a quantitative analysis of the $`T^{}`$ versus $`\rho ^{}`$, and $`T^{}`$ versus $`P^{}`$ phase diagrams. We also give a brief description of the procedure we followed to locate the coexistence courve and the –$`\lambda `$ line– mixing-demixing line.

Once the system reached equilibrium we calculated the structural properties of the system from simulations with $`N=4096`$ molecules. In figure 2 we show the density profiles at the relatively low temperatures, $`T^{}=0.65`$ and $`0.75`$, when $`\alpha =0.25`$. In this temperature region the liquid-vapor equilibrium structure of the mixture consists of a liquid A-liquid B interface in coexistence with the vapor phase. As temperature increases and reaches the region $`0.80T^{}1.25`$, however, this fluid phase structure rearranges in such a way that, the vapor phase spontaneously sets in between the Liquid A-Liquid B phases, as shown in figure 3. This structure remains stable during all the time of the simulation, about 55 ns in the scale of argon. This is a clear evidence that the vapor phase wets the liquid phases. In figure 4 we show similar results for $`\alpha =0.3`$. There one sees that at, $`T^{}=0.82`$, the structure of the system is such that there is a liquid A-liquid B interface. Nonetheless, at higher temperatures, for instance, $`T^{}=0.9`$, the vapor phase wets the liquid phases. This behavior of $`\rho ^{}(z)`$ for $`\alpha =0.3`$, suggest that $`T_\mathrm{W}^{}(\alpha =0.25)<T_\mathrm{W}^{}(\alpha =0.3)`$. In fact, comparing the structures plotted in figures 3 and 4, one should note that the mixture with $`\alpha =0.25`$ already wets at, $`T^{}=0.83`$. Following a systematic analysis of the structure of density profiles as function of temperatures for all the values of $`\alpha `$ we estimated the wetting transition temperatures $`T_\mathrm{W}^{}(\alpha )`$. We found that $`T_\mathrm{W}^{}`$ increases monotonically as a function of $`\alpha `$, whenever $`0<\alpha 0.5`$. For higher values, $`0.5<\alpha <1.0`$, this wetting phenomenon does not occur. The results are summarized in the wetting phase diagram in figure 5. We believe that the reason for which the system wets below $`\alpha =0.5`$, and no longer does above this value is due to the equal size of the molecules of type A and B and because $`\epsilon _{AA}=\epsilon _{BB}`$.

Now we will try to relate these density profiles structure and wetting phenomenon with the properties of the phase diagram of the mixture. To this end we will discuss in detail the properties of the corresponding phase diagram for $`\alpha =0.25`$. To begin it is important to remind that the present model binary mixture corresponds to the type III in the classification of Scott and Konynenburg scott70 . The phase diagram properties of this kind of mixtures has been quantitatively studied by Wilding et. al. wilding98 , using a square well potential for the intermolecular interactions. They showed that for a strong immiscible binary mixture a tricritical point exists. This means that the $`\lambda `$ line meets the LV coexistence curve at the critical point. In figure 6 we show the $`T^{}`$ versus $`\rho ^{}`$ phase diagram of the mixture obtained from extensive MD simulations with $`\alpha =0.25`$. Since for $`\alpha =0`$ and 1 our model reduces to single LJ fluids, we have also included in the same figure, for comparison, the phase diagram of a single LJ fluid. Both phase diagrams were calculated simulating a system with $`N=4096`$ molecules and using a shifted intermolecular potential with a cutoff of $`3\sigma `$. Therefore, the LV critical temperature of the single fluid became $`T_C^{LJ}1.2`$. One should note that, the mixture critical density is higher than that of the single LJ fluid, and its critical temperature shifts upwards. This shift in $`\rho _c^{}`$ and $`T_c^{}`$ occurs due to the fact that the less miscible the mixture is, –smaller $`\alpha `$–, the larger the temperature range of immiscibility. In addition, the $`\lambda `$ line appears to touch the LV critical point, and the critical point becomes tricritical. These results agree well with those obtained recently for square well strong immiscible binary mixtures wilding98 .

To check how sensitive are the phase diagram properties of the mixture, for $`\alpha =0.25`$, to the number of molecules in the simulations, we also calculated the phase boundaries simulating a system with $`N=1728`$ molecules. In fig. 7 we show the results of this finite size analysis. For a mixture with $`N=1728`$ the density of the liquid phase decreases while the density of the vapor phase increases. We also observe that at low temperatures the results are system size independent. However, as temperature increases, in particular, close to the critical point, there are differences in the phase boundaries obtained using $`N=1728`$ and 4096 molecules. As one approaches the critical point it turns out more difficult to determine the coexistence densities with a system with $`N=1728`$ molecules. This is so because the difference in the coexistence densities becomes smaller, and the number of particles in the system is not sufficient to give rise to bulk fluid phases. Nevertheless, for a system with $`N=4096`$ this is not the case and we indeed obtained the liquid and vapor fluid phases. On the other hand, it is known that near to the critical point the fluctuations of the density are strong and the vapor and liquid densities are not well defined. This fact complicates the location of the critical point. To circumvent this difficulty we proceed as in reference binder  and calculated the total density distributions in a system with $`N=4096`$, at several temperatures around $`T_c`$. The simulational box was divided in several slabs, parallel to the interface, of width between $`\sigma `$ and $`3\sigma `$ . The density of particles was calculated in each slab every 50 time steps of the MD simulations. A block average histogram of the densities is obtained every $`5\times 10^4`$ time steps. The resulting total density distribution, $`\rho _A+\rho _B`$, was calculated after averaging over 20 blocks. The result is presented in figure 8 . At temperatures slightly below $`T_c`$, we obtained density distributions that show two maxima. The low density maximum corresponds to the vapor phase and the higher maximum corresponds to the liquid density. Nonetheless, at temperatures above $`T_c`$, the density distributions show only one maximum. In Fig 8 we show the total density distributions at $`T^{}=1.1,\mathrm{\hspace{0.25em}1.15}`$ and $`1.2`$. The low density maximum is higher because the volume of the vapor phase is larger than the volume of the liquid phases. Thus, the vapor density appears with a higher frequency in the histograms.

To locate of mixing-demixing transition temperatures, –$`\lambda `$ line– for $`\alpha =0.25`$ and $`N=4096`$, we followed a similar procedure as that described for the location of the LV critical point. However, in the determination of the $`\lambda `$ line, we only consider the density distribution of one of the species, $`\rho _A`$ or $`\rho _B`$. The reason is that when the system is in the demixing region, the density distribution of one of the species shows two maxima. The low density peak corresponds to the poor fluid phase and the high density peak correspond to the rich fluid phase. Nonetheless, when the system is in the mixing region the fluid phases become homogeneous and therefore one should observe only one peak in the density distribution. Therefore, to locate the mixing-demixing points we looked at the transition from the two peak structure density distribution to one peak density distribution. This analysis was done as a function of the total density of the system and at three different temperatures. The results are shown in figure 9.

Another way of locating the $`\lambda `$ line is calculating the pressure versus temperature phase diagram, shown in figure 10, for a system with $`N=4096`$ molecules. The pressure was calculated as the average of the pressure tensor component perpendicular to the interface via the virial formula rowlin-capillary . There we also included for comparison, the results for a single LJ fluid with the same number of molecules. We found that the LV phase boundary is located right at $`T^{}1.25`$. At higher temperatures there are two branches that were obtained by approaching the mixing-demixing boundary from both sides of the $`\lambda `$ line. These branches enclose a region that is narrower in size, as compared to the region obtained in the $`\rho ^{}`$ versus $`T^{}`$ phase diagram. These results suggest that the calculation of the $`P^{}`$ versus $`T^{}`$ phase diagram gives more precise way to locate the mixing-demixing line. Again, we did find evidence that suggest that the $`\lambda `$ line meets the LV coexistence line at the LV critical point. These results give a strong support to the existence of a tricritical point wilding98 .

### IV.2 Interfacial properties

Now we turn to the discussion of the interfacial properties and surface phenomenon of the mixture at the liquid-vapor phase coexistence. We carried out a quantitative analysis of these properties for $`\alpha =0.25`$. To evaluate the wetting transition temperature with precision we calculated the interfacial free energies as a function of temperature. To this purpose we use the well known formula,

$$\gamma =_{\mathrm{bulk}_1}^{\mathrm{bulk}_2}\left(P_n(z)P_t(z)\right)𝑑z,$$
(4)
where the integrations were carried out up to the middle of the corresponding bulk phases. The tangential and normal pressure profiles were calculated using the definition of the Irving-Kirkwood pressure tensor rowlin-capillary . For a planar interface these pressure profiles are given by the formula rowlin-capillary ; rao .

$`P_n(z)`$ $`=`$ $`\rho (z)k_BT`$
$``$ $`{\displaystyle \frac{1}{2A}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{z_{ij}^2u_{ij}^{}(r_{ij})}{r_{ij}|z_{ij}|}}\theta \left({\displaystyle \frac{zz_i}{z_{ij}}}\right)\theta \left({\displaystyle \frac{z_jz}{z_{ij}}}\right),`$
$`P_t(z)`$ $`=`$ $`\rho (z)k_BT`$
$``$ $`{\displaystyle \frac{1}{4A}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{[x_{ij}^2+y_{ij}^2]u_{ij}^{}(r_{ij})}{r_{ij}|z_{ij}|}}\theta \left({\displaystyle \frac{zz_i}{z_{ij}}}\right)\theta \left({\displaystyle \frac{z_jz}{z_{ij}}}\right).`$
According to Young’s rule the difference,

$$\mathrm{\Delta }(T)=2\gamma _{\mathrm{LV}}\gamma _{\mathrm{LL}},$$
(7)
must be zero at the WTT. So, the wetting by the vapor phase occurs when the free energy difference, $`\mathrm{\Delta }`$, becomes negative. This quantity is plotted as a function of reduced temperature in the inset of figure 11 for $`\alpha =0.25`$ and $`N=4096`$. A linear extrapolation of the data indicates that the wetting occurs at about $`T_\mathrm{W}^{}(\alpha =0.25)=0.80`$. Notice that due to the planar geometry of the interfaces there is no contact angle and the three surface tensions are independent and considered separately. To improve the accuracy of $`T_\mathrm{W}^{}(\alpha =0.25)`$, one needs to carry out even more demanding simulations. This is due to the fact that the interfacial tension always shows relatively large fluctuations. The situation complicates even more when the simulations are performed at temperatures very close to the WTT. A second independent way to estimate the WTT and to figure out the nature of the WT, we calculated the adsorption of molecules at the LL interfaces shown in figures 2 . This is done using the formula,

$$\mathrm{\Gamma }=_{\mathrm{bulkA}}^{\mathrm{bulkB}}\left(\rho (z)\rho _{\mathrm{bulk}}\right)𝑑z,$$
(8)
The results of the calculations for $`\alpha =0.25`$ and $`N=4096`$ and $`6148`$ are plotted in figure 12. One sees that $`\mathrm{\Gamma }(T^{},\alpha =0.25)`$ decreases monotonically in the temperature range, $`0.75T^{}0.78`$, and it consistently shows negative values. This is due to the inhomogeneity at the LL interface, since the density there, is much smaller than the density of the liquid bulk phases. Nonetheless, as $`T^{}0.80`$ from below the adsorption jumps from a finite negative value up to minus infinite since the vapor wets the LL interface at $`T^{}(\alpha =0.25)0.80`$. Note that the WTT shifts slightly towards higher temperatures as the number of molecules in the system increases from $`N=4096`$ up to $`N=6148`$. This jump or discontinuity  is a strong indication that the WT is of first order. As expected the closer the temperature approach to $`T_W^{}(\alpha )`$ the stronger the fluctuations in $`\mathrm{\Gamma }(T^{},\alpha )`$. In figure 12 the closest approach to the WTT, in reduced temperature, was $`\frac{\delta T}{T_W}=6.25\times 10^3`$. One would expect that the wetting transition is of first order even for other values of the miscibility parameter, $`0<\alpha 0.5`$. This is so, since the behavior of $`\rho ^{}(z,\alpha )`$ is similar to that of $`\alpha =0.25`$, whose surface properties were studied in detail.

## V Conclusions

We carried out extensive MD simulations to study the LV phase coexistence, the structural properties and interfacial phenomena of a partially miscible symmetrical LJ binary mixture. By analyzing the density profiles as a function of temperature and miscibility parameter we estimated the wetting phase diagram, $`T_\mathrm{W}^{}`$ versus $`\alpha `$. The wetting of the vapor phase happens whenever $`0<\alpha 0.5`$. We also found that $`T_\mathrm{W}^{}(\alpha )`$, monotonically increases up to $`\alpha =0.5`$. For other values of $`\alpha `$, this wetting phenomenon does not occur. In addition, we also studied quantitatively the $`T^{}`$ versus $`\rho ^{}`$ and $`P^{}`$ versus $`T^{}`$ phase diagrams for $`\alpha =0.25`$. The results indicate that the former phase diagram shows a similar topology as that obtained for a square-well potential mixture estimated by means of mean-field theory and Monte Carlo simulations wilding98 . An analysis of the behavior of the adsorption of particles, at the LL interface, as a function of $`T^{}`$, for $`\alpha =0.25`$, led to the conclusion that the WT is of first order. These results should be valid for a family of mixtures of the type III, in the classification of Scott and Konynemburg. The phase diagrams discussed here together with the wetting phenomenon are an explicit quantitative demonstration of the scenario suggested some time ago based on a microscopic expression for the Hamaker constant getta93 . To the best of our knowledge this is the first time that this wetting phenomenon is quantitatively studied by means of MD simulations in binary LJ fluid mixtures in absence of a wall.. We would like to point out that this wetting phenomena agrees well with some experimental studies carried out with fluid binary mixtures. Finally, the results reported in this paper provide with a more complete understanding of the surface phenomena in partially miscible fluid binary mixtures.

## VI Acknowledgments

EDH would like to thank the DAAD-Germany for financial support during the summer of 1999. Partial support from CONACYT contracts L0080-E9607 (EDH) and 25298-E (GRS) as well as from DGAPA-UNAM contracts IN110103, and IX107704, is also acknowledge.


text:
# Extinction due to amorphous carbon grains in red quasars from the SDSS

## 1 Introduction

It is well known that the spectra of Seyfert galaxies – particularly Seyfert 2 galaxies - are considerably modified by circumnuclear extinction. It is also believed that Narrow Line Seyfert galaxies are on average redder than typical Seyfert 1 galaxies because of reddening due to the presence of dust (e.g. Kuraszkiewicz & Wilkes 2000, Constantin & Shields 2003). Evidence is now accumulating that a significant fraction of bright AGN, including high redshift quasars, are affected by dust.

The dust is seen directly in emission (IR continuum of dust and emission lines of accompanying molecular gas; e.g. Priddey et al. 2003, Omont et al. 1996, Kuraszkiewicz et al. 2003, Andreani, Francheschini & Granato 1999). The presence of dust may also be detected indirectly because of its effect on the quasar continuum in the optical/UV band. Recent surveys show the widespread presence of red quasars in quasar samples, and reddening by dust is the most plausible explanation for their spectral shape (e.g. FIRST quasar sample, White et al. 2003; SDSS quasar survey, Richards et al. 2003; 2MASS red quasars, Marble et al. 2003).

Knowledge of dust properties in the close vicinity of an active nucleus is thus important both for determining intrinsic AGN spectra and understanding the conditions in the gas surrounding the nucleus.

It has already been argued that the properties of the dust in such extreme conditions need not be similar to the dust properties in the interstellar medium of our Galaxy. It has been found that the extinction curve of Prévot et al. (1984) determined for the Small Magellanic Cloud (SMC) roughly applies to AGN dust (e.g. Czerny, Loska & Szczerba 1991, Laor & Draine 1993, Pitman et al. 2000, Kuhn et al. 2001, Richards et al. 2003).

In this paper we use the composite spectra of quasars from the Sloan Digital Sky Survey (SDSS) prepared by Richards et al. (2003) to determine an observational quasar extinction curve. We then compare this curve with model curves and argue that the circumnuclear regions of AGN are filled with amorphous carbon grains. We provide a simple empirical formula for this new extinction law.

## 2 Extinction curve from red SDSS quasars

Results from the SDSS were used by Richards et al. (2003) to create 6 composite quasar spectra (see their Fig. 7), selected according to the adopted ranges of the (g-i) color excess. Composite 1, made from 770 objects, contains the bluest quasars in the sample. Their optical/UV slope is well fitted by a power law, $`F_\nu \nu ^\alpha `$, with energy index $`\alpha =0.25`$. Composites 1 to 4 each consist of 770 objects. Composite 5 consists of 211 objects well representing a red quasar population; only objects with large color excesses were included, but the few reddest objects were rejected in order to preserve the completeness of the sample.

We assume that the bluest composite (composite 1) is essentially unaffected by extinction (we will return to this point in the Discussion). We then assume that other composites intrinsically represent the same spectrum as composite 1, except affected by dust.

We can therefore directly construct the observed quasar extinction curve for each of the composites from 2 to 5 vs. composite 1 by using its usual definition,

$$X(\lambda )=\frac{A_{\lambda V}}{E_{BV}}=\frac{A_\lambda A_V}{A_BA_V},$$
(1)
where

$$A_\lambda =2.5\mathrm{log}(F_\lambda ^1/F_\lambda ^i).$$
(2)
Here $`F_\lambda ^1`$ is the composite 1 spectrum and $`F_\lambda ^i`$ is one of the redder composites. The values $`A_V`$ and $`A_B`$ are calculated at the corresponding visual (5500 Å) and blue (4400 Å) wavelengths.

The curves obtained from composites 2 and 2n were very noisy since the differences between the spectra were small, so we neglected them. We next averaged the curves obtained from composites 3, 4 and 5 in order to reduce errors, thus creating a single curve.

The resulting extinction curve in the range covered by the SDSS data is shown in Fig. 1. The curve is monotonic with wavelength, without any trace of the 2200 Å absorption feature. It thus resembles the observationally determined extinction curve for dust in the Small Magellanic Cloud (Prévot et al. 1984) more than the Seaton (1979) extinction curve that characterizes the interstellar dust in our galaxy. Accordingly, if the composite spectra are indeed shaped by circumnuclear extinction, the dust composition close to quasar nuclei must be non-standard.

In order to analyze the quasar extinction curve properties, we compare it with several theoretical extinction curves.

## 3 Theoretical extinction curves

We construct the model extinction curves for the dust surrounding an active nucleus following the approach of Czerny et al. (1995).

We assume that the nucleus is a source of radiation with a spectral shape given by the sum of two components: a power law component extending from IR to X-rays, and an accretion disc component dominating the optical/UV part. The power law component is spherically symmetric while the disc component intensity depends on the cosine of the inclination angle to the disc symmetry axis.

The innermost part of an AGN is surrounded by a spherical dust shell. The inner radius of the dusty shell is determined by the temperature of dust evaporation. The outer radius of the shell is assumed to be at 300 pc, and the adopted density profile is given in the form of a power law, $`n(r)r^{0.75}`$, the favored one by models of the far-IR spectra of AGN (e.g. Loska, Szczerba & Czerny 1993, Granato, Danese & Franceschini 1997). The results are not very sensitive to these parameters so long as the spatial distribution is rather flat.

We determine the theoretical extinction curve by solving for the radiation transfer through the dust surrounding the envelope, determining the dust-reddened radiation spectrum, and comparing this spectrum to the intrinsic spectrum of the source, using Eqs. 1 and 2. The radiative transfer is solved for as described by Loska et al. (1993). The temperature of the grains of a given size and chemical composition is calculated from the energy balance. We assume a maximum carbon grain temperature of 1500 K and maximum silicate dust temperature of 1000 K, and we allow for selective grain evaporation close to the inner edge of the dusty cloud. The radiative transfer is solved for assuming that the dust is optically thin, i.e. we neglect multiple scattering effects.

Optical properties of astronomical silicate and carbon grains are taken from Draine (see Laor & Draine 1993), and the amorphous carbon properties were calculated using Mie theory and the optical constants of Martin (see Rouleau & Martin 1991). Some modifications were introduced for large grains/short wavelengths (see Czerny et al. 1995 for details). The distribution of grain sizes was assumed to be in a power law form, $`n(a)a^p`$, with minimum and maximum grain sizes given by $`a_{min}`$ and $`a_{max}`$. The standard values of these parameters in the ISM are: $`p=3.5`$, $`a_{min}=0.005\mu `$m and $`a_{max}=0.25\mu `$m (Mathis, Rumpl & Nordsieck 1977). We generally tried to minimize departure from these values.

The extinction curves representing various chemical compositions of dust grains but otherwise based on the dust parameters $`p,a_{min}`$, and $`a_{max}`$ of the ISM are shown in Fig 1. The observed extinction curve roughly follows at long wavelengths the extinction curve of dust consisting of AC amorphous carbon grains. The departure is seen only at short wavelengths above $`1500`$ Å.

We can find a better representation of the observed curve assuming AC composition but allowing for a change in the minimum and maximum grain sizes. In Fig. 1 we show such a curve with $`a_{min}=0.016\mu `$m and $`a_{max}=0.12\mu `$m.

We also analyzed the properties of dust of other chemical compositions, but they did not mimic the quasar extinction curve equally well. In particular, no other standard or non-standard extinction curve could reasonably trace the observational extinction curve in the $`.2\mu `$m to $`.5\mu `$m range. Furthermore, we do not expect large amounts of silicates in AGN surroundings; silicate emission features expected in the IR are weak or absent in most AGN (Roche et al. 1991), and efficient formation of CO molecules offers a natural explanation for the depletion of the oxygen necessary for silicate formation (Rawlings & Wiliams 1989). Any significant contribution from graphite also resulted in the presence of a 2200 Å peak and no traces of such a peak are in the data. We thus argue that AC amorphous carbon grains best describe the observational extinction curve based on the Richards et al. (2003) data because there are significant errors when attempting to model the curve based on any other dust grains.

## 4 Discussion

The lack of the 2200 Å absorption feature due to graphite grains characteristic of interstellar extinction was already noted in several papers concerning AGN (e.g. Pitman et al. 2000, Kuhn et al. 2001, Richards et al. 2003) and therefore an SMC extinction curve seems to apply better than the Seaton ISM extinction curve. Our results based on SDSS quasars from Richards et al. (2003) support this view.

Our best model of the observational quasar extinction curve is made from amorphous AC carbon grains. We assumed a standard power law size distribution of grains, $`p=3.5`$, but better agreement with the observational curve required a narrower size range (minimum size $`0.016\mu `$m, maximum size $`0.12\mu `$m) than the canonical one of Mathis et al. (1977), 0.005 and 0.25$`\mu `$m.

This extinction curve is slightly flatter than the SMC extinction curve at short wavelengths. The reddening law corresponding to it can be well approximated by the following empirical formula

$$\frac{A_\lambda }{E_{BV}}=1.36+13\mathrm{log}(1/\lambda )[\mu m],$$
(3)
in the range of $`1/\lambda `$ between 1.5 and 8.5 $`\mu `$m<sup>-1</sup>.

Our analysis was based on an average of several observational curves derived from Richards et al. (2003) quasar data. We derived extinction curves assuming that their composite 1 was the intrinsic spectrum of all the composites. We were then able to construct extinction curves based on composite 1 vs. composites 2 through 5 and composite 1n vs. composites 2n through 5n. We ended up disregarding the composite 1 vs. composite 2 and composite 1n vs. composite 2n extinction curves because they were too noisy. We then averaged the extinction curves based on composite 1 vs. composites 3 through 5. We selected the composite 1 vs. composite 5 extinction curve over the composite 1n vs. composite 5n curve because the former has higher signal-to-noise. The latter was corrected for redshift and magnitude differences, but we found that the two curves were quite similar with the latter having more noise. We did not include the 2n through 5n subgroups in our average because these quasars were already included in the 2 through 5 composites. We found that the composite 1 vs. composites 3 through 5 and composite 1n vs. composites 3n through 5n created a small range of possible extinction curves. We believe that the true observational extinction curve lies somewhere within this range, and so our final observational extinction curve is an average of the curves that are contained in this range. We find that AC amorphous carbon grains are able to model the extinction curves within this range as well as they model the average extinction curve. We also find that even the lower extinction curves in this range are modeled with difficulty when taking into account other types of dust grains.

Our analysis used the composite spectra. Any such spectra, including those of Richards et al. (2003), contain systematical errors which are difficult to assign. The errors come partially from the calibration uncertainty, and partially from the systematic differences in the objects constituting the short wavelength and the long wavelength parts of the final spectrum. Therefore, determination of the errors of the derived extinction curve is rather difficult. Although the conclusion about the dominating role of amorphous carbon seems to be firm, the grain size range needs confirmation.

This dust composition is not surprising taking into account that the conditions around an active nucleus are quite similar to conditions in the envelopes of asymptotic giant branch (AGB) stars. Radiation flux, pressure and densities may be comparable (see e.g. Elvis, Marengo and Karovska 2002), particularly taking into account that the cool gas in AGN is possibly in the form of clouds. This is indicated by the presence of the Broad Line Region and Narrow Line Region implied by AGN spectra. Extinction in carbon-rich AGB stars is also dominated by amorphous carbon, although in some stars with additional silicon carbide (e.g. Piovan, Tantalo & Chiosi 2003).

The 2200 Å feature seen in the interstellar medium is caused by single graphite and/or clusters of graphite grains (see e.g. Draine & Lee 1984, Andersen et al. 2003). Laboratory data indicate that graphite grains form from amorphous carbon grains as a result of exposure to strong UV radiation (Ogmen & Duley 1988; see also Dorschner & Henning 1995). AGN environments provide intensive UV flux, much more than the relatively cool environments of AGB stars (their temperature is of the order of 2000 - 4000 K), yet graphite grains apparently do not form in AGN surroundings. This may support the view of Elvis et al. (2002) that the dust forms in an outflowing wind. The dust forms and is exposed to strong UV radiation only for a time period of the order of years to hundreds of years, most probably not enough time to form the regular crystal structures of graphite. Alternatively, dust may form relatively far from the nucleus. However, time delay measurements between the optical and IR band indicate hot dust distances not larger than $`1`$ pc for bright objects (e.g. Sitko et al. 1993, Ulrich, Maraschi & Urry 1997).

The composite 1 spectrum is most probably only weakly reddened, if at all. We estimate this in the following way. We assume that the amorphous carbon grain extinction curve also well represents the reddening within the composite 1 objects. We allow for an arbitrary amount of reddening and compare the resulting spectra with the standard accretion disc spectra of Shakura & Sunyaev (1973). Such models represent bright and blue quasar spectra quite well (see e.g. Koratkar & Blaes 1999 for the discussion of the quasar composite spectrum of Francis et al. 1991).

We can model the composite 1 spectrum assuming any value of $`E_{BV}`$ between 0 and $`0.05`$ by changing the mass and accretion rate of our model. Reddening larger than $`E_{BV}=0.05`$ is excluded since it leads to an intrinsic spectrum harder than the asymptotic energy index value of 0.33 predicted by the disc model. These two extreme cases of $`E_{BV}=0`$ and 0.05 are shown in Fig. 2. For all intermediate reddening values the model well represents the overall continuum. However, there is evidently some flattening of the observed spectrum towards IR (see Fig. 2). This effect is well known in AGN although its nature is not clear. It might be caused by the contribution of some starlight in the host galaxy (clearly important for less luminous AGN although not necessarily for bright quasars, see e.g. Elvis et al. 1994), the illumination of the outer disc by the radiation generated in its innermost region, or perhaps some contribution of non-thermal emission to the overall spectrum.

The good correspondence between the dereddened spectrum of the bright quasars with a simple accretion disc model indicates that in these sources the production of X-ray emission does not modify an optically thick flow. In contrast to much fainter Seyfert 1 galaxies, X-ray emission in quasars must therefore make up only a small fraction of the bolometric luminosity. Studies of broad band $`\alpha _{ox}`$ index trends with luminosity indeed seem to support such a view (e.g. Bechtold et al. 2003).

## 5 Conclusions

Analysis of SDSS quasar composites shows that quasar spectra are affected by circumnuclear dust composed primarily of AC amorphous carbon grains. The derived observational extinction curve is similar to SMC extinction but flatter at the shortest wavelengths. Eq. 3 provides a simple empirical formula for the derived quasar extinction law.

## Acknowledgements

This work was supported in part by grants 2P03D 003 22 (BCz & ZL) of the Polish State Committee for Scientific Research (KBN). JL acknowledges financial support from the Institute of Physics, Polish Academy of Sciences.

This paper has been processed by the authors using the Blackwell Scientific Publications style file.


text:
# Dynamical interplay between local connectivity and global competition in a networked population

\[

## Abstract

We show, both numerically and analytically, that the consequences of ‘wiring up’ a competitive population depend quite dramatically on the interplay between the local connectivity and the global resources. With modest global resources, adding small amounts of local connectivity yields an increasingly heterogeneous population. With substantial global resources, high-performing yet reasonably homogenous collective states emerge instead.

PACS numbers: 02.50.Le, 87.23.Kg, 89.65.Ef, 05.40.2a

\]

Is ‘getting connected’ a good or bad thing? How does increased access to both global (i.e. public) information and local (i.e. private) information, affect both the population as a whole and its individual members? Thinking in terms of future technologies, what are the possible benefits and dangers of introducing communication channels between collections of intelligent devices, microsensors, semi-autonomous robots, nanocomputers, and even biological micro-organisms such as bacteria ? These questions are relevant to a wide range of computational, technological, biological and socio-economic systems. However such questions are very difficult to answer since these systems are simply ‘too complex’ .

Here we examine these questions using a minimal generic model which mimics a group of selfish, rational individuals or components, referred to simply as ‘agents’. These agents repeatedly compete for a global resource by making inductive decisions about future group behavior. The population and its agents may be biological (e.g. a population of cellular organisms competing for nutrients), computational (e.g. a grid of software modules competing for processing time), mechanical (e.g. a constellation of sensors or devices competing for communications bandwidth or operating power) or social (e.g. a population of companies competing for business). We uncover a rich interplay between the global competition for resources and the local inter-agent connectivity. In a population with modest resources, low levels of inter-connectivity increase the disparity between successful and unsuccessful agents. By contrast in a higher-resource population with low inter-connectivity, high-performing collective states can spontaneously arise in which nearly all agents are reasonably successful. At high levels of inter-connectivity, the overall population becomes fairer (i.e. smaller spread in success-rates) but less efficient (i.e. smaller mean success-rate) irrespective of the global resource level.

Our B-A-R (Binary Agent Resource) model is a minimal generic model of a competitive networked population, which is based on Arthur’s so-called El Farol Problem  in which a population of agents repeatedly compete for a limited resource (i.e. seating in a crowded bar) . Despite its everday human setting, the El Farol Problem embodies key generic characteristics of Complex Systems : feedback and adaptation at the macroscopic and/or microscopic level, many interacting parts, non-stationarity, evolution, coupling with the environment, and observed dynamics which depend upon the particular realization of the system. Although the El Farol Problem is more general than, say, the Minority Game  in that it allows for arbitrary levels of global resource and is non-binary, the agents still only have access to global information. Most biological, informational and socio-economic systems have at least some degree of underlying connectivity among agents  leading to exchange of local information. It is the resulting interplay of network connectivity and competition for global resources that we wish to explore and understand. Therefore in contrast to previous works  our B-A-R model features network connectivity and a general global resource level $`L`$. Here $`L`$, which is not announced to the agents, may represent the optimal load capacity of a regional power grid, public utility, or electronic microcircuit; the available space in a given urban region, public area or popular bar; the data capacity of a given communications link in an informational or biological system; the vehicular or passenger capacity on a given road or transport system.

At each timestep $`t`$, each agent decides whether to access resource $`L`$ (i.e. action $`+1`$ which might correspond to an electronic component or device deciding to draw power, a computer sending a data-packet down a given route, a surfer accessing a particular website, a commuter taking a given road to work) or not to access this resource (i.e. action $`1`$). The two global outcomes at each timestep, ‘resource over-used’ and ‘resource not over-used’, are denoted as ‘0’ and ‘1’. If the number of agents $`n_{+1}[t]`$ choosing action $`+1`$ exceeds $`L`$ (i.e. resource over-used and hence global outcome ‘0’) then the $`n_1[t]`$ abstaining agents win. By contrast if $`n_{+1}[t]L`$ (i.e. resource not over-used and hence global outcome ‘1’) then these $`n_{+1}[t]`$ agents win. Each agent decides its actions in light of (i) global information which takes the form of the history of the $`m`$ most recent global outcomes, and (ii) local information obtained via network connections. Adaptation is introduced by randomly assigning $`S`$ strategies to each agent. Each of the $`2^{2^m}`$ possible strategies is a bitstring of length $`2^m`$ defining an action ($`+1`$ or $`1`$) for each of the $`2^m`$ possible global outcome histories. For example, $`m=2`$ has $`2^2=4`$ possible global outcome histories: $`00`$, $`01`$, $`10`$ and $`11`$. Consequently there are $`2^{2^{m=2}}=16`$ possible strategies. Strategies which predicted the winning (losing) action at a given timestep, are assigned (deducted) one point. At each timestep, each agent compares the score of his own best-scoring strategy (or strategies) with the highest-scoring strategy (or strategies) among the agents to whom he is connected. The agent adopts the action of whichever strategy is highest-scoring overall, using a coin-toss to break any ties. For simplicity we here assume a random network, where the connection between any two agents exists with a probability $`p`$. We emphasize that any of the above model assumptions can easily be generalized: the numerical results are reasonably robust.

Figure 1 shows the mean success-rate per agent, averaged over many numerical realizations of the underlying connections and strategy distributions (left axis: solid curve) and the success-rate distribution within the population for a typical run (right axis: histograms) as a function of the inter-connectivity $`p`$ in a modest-resource population. We have chosen $`L=(N1)/2`$, which is representative of modest resources since there are more losers than winners. We focus on the small $`m`$ or ‘crowded’ regime  since there are relatively few strategies compared to the number of agents and hence many agents will simultaneously be competing to win with the same strategy. The mean success-rate decreases rapidly as $`p`$ increases (left axis: solid curve) before saturating at $`p=p_{sat}0.1`$. As $`m`$ increases, $`p_{sat}`$ increases – full results for the uncrowded, large $`m`$ regime will be discussed elsewhere. The success-rate distribution for a typical run (right axis: histogram) at $`p=0`$ exhibits a definite ‘class structure’ in terms of success. Dramatic changes then arise as $`p`$ increases. The spread in the success-rate – shown explicitly in the inset and by vertical arrows in Fig. 1 – indicates a rapid increase in the population’s heterogeneity in the range $`0p0.02`$. The success-rate distribution becomes almost continuous, washing out the $`p=0`$ class structure. Above $`p0.02`$, there is a rapid drop in the proportion of highly-successful agents. Further increasing $`p`$ leads to a decrease of the spread in success-rates. High levels of inter-connectivity therefore provide increased fairness (i.e. small spread in success-rate) but decreased efficiency (i.e. small mean success-rate).

Figure 2 compares our analytical results using the Crowd-Anticrowd theory  to the numerical results for the mean success-rate per agent. As will be explained in more detail elsewhere, the mean success-rate per agent can be obtained from the standard deviation of the number of agents making a given choice (e.g. $`+1`$), by averaging this quantity over the attractors of the system (i.e. averaging over the Eulerian Trail  or subsets of it). The mean success-rate can hence be expressed  in terms of the mean number or Crowd of agents using the $`K`$’th highest-scoring strategy (i.e. $`n_K^{\mathrm{net}}`$) and the mean number or Anticrowd of agents using strategy $`\overline{K}`$ which is anticorrelated to $`K`$ (i.e. $`n_{\overline{K}}^{\mathrm{net}}`$ where $`\overline{K}=2^{m+1}+1K`$). Explicitly

$$n_K^{\mathrm{net}}=n_K+n_Kn_K$$
(1)
where $`n_K`$ is the number of agents who would have used strategy $`K`$ in the absence of the network because of the initial random strategy allocation:

$$n_K=N\left(\left[1\frac{(K1)}{2^{m+1}}\right]^S\left[1\frac{K}{2^{m+1}}\right]^S\right).$$
(2)
In short, $`n_K`$ is the intrinsic Crowd-size due to crowding in the global strategy space. $`n_K`$ is a sum over all agents who use strategy $`K`$ as a direct result of a network connection:

$$n_K=\left[\underset{J>K}{}n_J\right]\left[(1p)^{_{G<K}n_G}\right]\left[1(1p)^{n_K}\right].$$
(3)
Hence $`n_K`$ represents the increase in Crowd-size due to the local connectivity. By contrast, $`n_K`$ is a sum over all agents who would have used $`K`$ in the absence of a network, but who will now use a better strategy as a result of their network connections:

$$n_K=n_K\left[1(1p)^{_{G<K}n_G}\right].$$
(4)
Hence $`n_K`$ represents the decrease in Crowd-size due to the local connectivity. The resulting analytic expressions for the mean success-rate are in excellent agreement with numerical results (see Fig. 2) implying that the population’s dynamical evolution is indeed governed by crowds resulting from the interplay between local connectivity and global competition.

Figure 3 shows the mean success-rate per agent (thin solid lines) as a function of the inter-connectivity $`p`$ in higher-resource populations, for a range of $`L`$ values and small $`m`$. If the global resource level $`L`$ exceeds a critical amount $`L_{crit}=3N/4`$, the mean success-rate per agent can exhibit a maximum at small but finite inter-connectivity, for the following reason. When $`L>L_{crit}`$ and $`p=0`$, the population inhabits a ‘frozen’ state in the sense that the global outcome is persistently 1, i.e. $`\mathrm{}111111`$. With $`S=2`$ as in Fig. 3, each agent has probability $`1/4`$ of being assigned two strategies which both define action $`1`$ following a string of $`m`$ ‘1’s. Thus approximately $`1/4`$ of the population always lose and $`3/4`$ of the population always win, leading to $`L_{crit}=3N/475`$. The global resource is therefore under-used at $`p=0`$ by $`\mathrm{\Delta }L(L3N/4)`$ at each timestep. Therefore increasing $`p`$ away from $`p=0`$ will benefit some of the less successful agents by connecting them up to successful agents, thereby giving them access to $`\mathrm{\Delta }L`$, until a $`p`$ value is reached where the connectivity is sufficient to break the outcome series of 1’s. These run-averaged results can be better understood by considering run-specific results for $`L=100`$ as an illustrative case (scattered circles). Specific runs appear to aggregate into groups having similar temporal dynamics and success-rate distributions – these groups correspond to different dynamical states of the system. The increase in the mean success-rate at low $`p`$, corresponds to the system following an Efficient State E which has an outcome series of 1’s as discussed above. The success-rate distribution in State E (inset) is characterized by groups of persistent winners and losers. As $`p`$ increases further, the outcome series for high $`L`$ tends to move through a set of Intermediate States which combine a low spread in success-rate with a high mean (e.g. Intermediate States A and B). In Intermediate State A, the least-successful agents (which in State E had zero success) now have a success-rate which exceeds the average success-rate in the high connectivity limit. In the high $`p`$ limit, the outcome series tends toward the period-4 Eulerian Trail  given by $`\mathrm{}00110011\mathrm{}`$. The corresponding number of agents taking action $`+1`$ follows the pattern $`\mathrm{}N,N,N/2,N/2,N,N,N/2,N/2\mathrm{}`$. The resulting Fair State F corresponds to all agents having access to the best performing strategies, either by being assigned that strategy or by being connected to another agent with that strategy. The resulting system is fair but inefficient, having a small spread but also a small mean success-rate $`0.25`$.

Figure 4 shows the mean success-rate per agent (thin solid lines) as a function of the inter-connectivity $`p`$ in a high-resource population ($`L=95`$) for various history bit-string lengths $`m`$. The Intermediate States, characterized by a reasonably high mean success-rate and a reasonably small spread (see Fig. 3 inset), have an increasingly dominant effect on the system’s behavior as $`m`$ increases. For a particular value of $`m`$, increasing $`p`$ moves the system toward cycles of increasing period and hence the fractions of $`1`$’s and $`0`$’s in the output series become more equal. The resulting State F, which represents an increasingly noisy version of the Eulerian Trail as $`m`$ increases, is fair but inefficient (see Fig. 3 inset).

We have derived analytic expressions (thick solid lines in Figs. 3 and 4) which describe well the three dynamical regimes exhibited by the numerical results. The upper analytic branch at low $`p`$, describing the efficient State E, is given by:

$$\frac{3}{4}+\frac{1}{4}\left[1(1p)^{\frac{3N}{4}}\right].$$
(5)
The middle analytic branch at intermediate $`p`$, describing the Intermediate States, is given by:

$$\frac{3}{4}\frac{9}{32}(1p)^{\frac{7N}{16}}\frac{1}{32}(1p)^{\frac{15N}{16}}+\frac{1}{8}(1p)^{\frac{3N}{4}}.$$
(6)
The lower analytic branch at high $`p`$, describing the fair State F, is given by:

$$\frac{1}{4}+\frac{9}{128}(1p)^{\frac{7N}{16}}+\frac{1}{128}(1p)^{\frac{15N}{16}}+\frac{1}{16}(1p)^{\frac{3N}{4}}.$$
(7)
The outcome series of 1’s at low $`p`$ which yield State E, can persist up to $`p_{crit}1(14\mathrm{\Delta }L/N)^{4/3N}`$. For $`L=80,90,95`$ and $`100`$, this yields $`p_{crit}0.002,0.011,0.019`$ and $`0.042`$ respectively, which are also all in excellent agreement with the numerical results. A full analysis of the dynamics of, and switching between, these regimes will be given elsewhere.

In conclusion, we have reported a rich dynamical interplay between local connectivity and global competition in a generic networked population. Apart from the intrinsic interest regarding functionality in complex system networks, our results suggest that the internal connectivity and global resources in such systems can be engineered in order to enhance performance.

PMH acknowledges support from the Research Grants Council of the Hong Kong SAR Government (grant CUHK4241/01P).


text:
# Time Reversal and Exceptional Points

## A The phase $`\xi `$

The product

$$U_0^{}U_1=U^{}(\varphi _0,\tau _0)U(\varphi _1,\tau _1)$$
(A1)
is equal to

$$U_0^{}U_1=\left(\begin{array}{cc}\mathrm{cos}\varphi _0\mathrm{cos}\varphi _1+\mathrm{sin}\varphi _0\mathrm{sin}\varphi _1e^{i(\tau _0\tau _1)}& \mathrm{cos}\varphi _1\mathrm{sin}\varphi _0e^{i\tau _0}\mathrm{cos}\varphi _0\mathrm{sin}\varphi _1e^{i\tau _1}\\ \mathrm{cos}\varphi _1\mathrm{sin}\varphi _0e^{i\tau _0}+\mathrm{cos}\varphi _0\mathrm{sin}\varphi _1e^{i\tau _1}& \mathrm{cos}\varphi _0\mathrm{cos}\varphi _1+\mathrm{sin}\varphi _0\mathrm{sin}\varphi _1e^{i(\tau _0\tau _1)}\end{array}\right).$$
(A2)
We introduce the phase of the (1,1)-element, viz.

$$\gamma =\mathrm{arg}\left(\mathrm{cos}\varphi _0\mathrm{cos}\varphi _1+\mathrm{sin}\varphi _0\mathrm{sin}\varphi _1e^{i(\tau _0\tau _1)}\right),$$
(A3)
and write (A1) in the form

$$U_0^{}U_1=\left(\begin{array}{cc}\mathrm{cos}\beta & \mathrm{sin}\beta e^{i\xi }\\ \mathrm{sin}\beta e^{i\xi }& \mathrm{cos}\beta \end{array}\right)z(2\gamma ).$$
(A4)
Here, $`z`$ is the matrix defined in (12), and $`\beta `$ is given by

$$\mathrm{cos}\beta =\left(\mathrm{cos}^2\varphi _0\mathrm{cos}^2\varphi _1+\mathrm{sin}^2\varphi _0\mathrm{sin}^2\varphi _1+2\mathrm{cos}\varphi _0\mathrm{cos}\varphi _1\mathrm{sin}\varphi _0\mathrm{sin}\varphi _1\mathrm{cos}(\tau _0\tau _1)\right)^{1/2}$$
(A5)
and $`\xi `$ by

$$\xi =\mathrm{arg}\left(\mathrm{cos}\varphi _1\mathrm{sin}\varphi _0e^{i\tau _0}\mathrm{cos}\varphi _0\mathrm{sin}\varphi _1e^{i\tau _1}\right)+\gamma .$$
(A6)
The exceptional point occurs at

$$\lambda _{\mathrm{EP}}=\frac{ϵ_1ϵ_2}{\omega _1\omega _2}e^{\pm 2i\beta }.$$
(A7)