text:
# Fighting cancer with virus

## Abstract

One of the most promising strategies to treat cancer is attacking it with viruses. Virus can kill tumor cells specifically or act as carriers that deliver normal genes into cancer cells. A model for virus therapy of cancer is investigated and some of its predictions are in agreement with results obtained from experimental tumors. Furthermore, the model reveals an oscillatory (periodic or aperiodic) response of tumor cells and virus populations which may difficult clinical prognosis. These results suggest the need for new in vivo and in vitro experiments aiming to detect this oscillatory response.

In the last few years, in which most of the research in molecular biology of cancer focused in new target genes and potential biomarkers, promising roads for tumor treatment such as gene therapy, virus and antiangiogenic drugs have been openned Evan ; Liotta . However, most of the clinically-used cancer therapies have been developed empirically Evan , and therefore mathematical models might be complementary (maybe necessary) tools for the understanding of the dynamics of drug response in the organism. Noticeable examples of such modeling are the recent works of Scalerandi et al. Scalerandi , in which the action of different antiangiogenic drugs was simulated, of Wu et al. Wu  and Wein et al. Wein , concerning the spreading of a virus that replicates selectively in tumor cells, and of Kozusko et al., addressing cancer cell growth and response to treatment with cyclic-specific chemotherapeutic drugs Kozusko .

Current cancer therapy is based on damaging DNA by irradiation or chemicals. However, cells lacking p53 tumor supressor gene are unable to respond to DNA damage by inducing cell cycle arrest or apoptosis, frequently becoming refractary to chemotherapy or radiation. These cells occur in more than $`50\%`$ of all human cancers, and hence it is imperative to develope new forms of treatment. Viruses can be used to trigger a DNA damage response without damaging cellular DNA and to selectively replicate in and lyse p53-deficient human tumor cells Raj ; Bischoff . Other potential oncolytic agents were recently tested. Among them are the reovirus, which can be used in more than $`50\%`$ of all human tumors having an activated Ras signaling pathway Coffey , and the avian Newcastle disease virus Lorence , a genetically altered herpes simplex virus, which kill cancer cells exhibiting an altered p16/pRB tumor suppressor pathway Chase . In particular, clinical trials in humans using Delta-24 RGD adenovirus, an agent with unprecedented success in the treatment of gliomas Suzuki , are currently being designed. Also, since genetically engineered virus can be selected to kill only cancer cells, virus therapy might be more specific than standard DNA-damaging drugs associated to many toxic effects on the normal cells Vogelstein .

In this letter, a model for virus therapy of cancer is investigated and a few of their predictions are compared with results from experimental tumors. The model, based on a previous one recently proposed by us for the avascular growth of tumors Silvio  and applied to chemotherapy modelling Silvio2 , combines macroscopic reaction-diffusion equations for the nutrients and virus with microscopic stochastic rules governing the actions of individual cancer cells. It seems to be appropriated to study the response of cancers accessible to direct intratumoral virus injection (primary brain tumors and cancers of the head and neck, for example), but mainly their distant metastatic outgrowths. Its prime advantages are the use of both, an individual-based framework in which the fate of each cancer cell is traced, and a set of stochastic rules for cell actions as a simple effective kinetic cellular model describing the population dynamics of interacting cells Bellomo ; Lachowicz .

The tissue, represented by a square lattice, is fed by a single capillary vessel at $`i=0`$, the top of the lattice. The model considers four different cell populations: normal $`\sigma _n`$, uninfected $`\sigma _c`$, infected $`\sigma _v`$ and dead $`\sigma _d`$ tumor cells. In contrast to the normal and dead tumor cells, one or more cancer cells can pile up in a given site. Thus, the cell populations in any site, with coordinates $`\stackrel{}{x}=(i,j)`$, $`i,j=0,1,\mathrm{},L`$, assume one of the possible values $`\sigma _n(\stackrel{}{x}),\sigma _d(\stackrel{}{x})=0,1`$ and $`\sigma _c(\stackrel{}{x}),\sigma _v(\stackrel{}{x})=0,1,\mathrm{}`$. Initially, a single cancer cell at $`\stackrel{}{x}=(d,L/2)`$ is introduced in the normal tissue. The row $`i=L+1`$ constitutes the external border of the tissue. The evolution rules implemented in our model are the following:

Nutrient evolution.—The nutrients are divided into two classes: directly used, $`N(\stackrel{}{x},t)`$, or not, $`M(\stackrel{}{x},t)`$, for DNA synthesis. Both, $`N`$ and $`M`$, are the stationary solutions of the simplest (linear with constant coefficients) dimensionless reaction-diffusion equation

$$\frac{H}{t}=^2H\alpha ^2H\sigma _n\lambda _H\alpha ^2H(\sigma _c+\sigma _v)$$
(1)
in which $`H`$ represents the nutrient fields. Distinct nutrient consumption rates for normal and cancer cells, characterized by factors $`\lambda _N`$ and $`\lambda _M`$ ($`\lambda _N>\lambda _M`$), are assumed. (see ref. Silvio  for the complete variable transformations). Also, it is assumed that infected cancer cells sustain their metabolism until lysis.

These equations satisfy periodic and Neumann boundary conditions along the horizontal axis and at the border of the tissue ($`i=L`$), respectively. At the capillary vessel ($`i=0`$) $`N=M=1`$, representing the continuous and fixed supply of nutrients provided by the blood flow.

Cell dynamics.—Each cancer cell, randomly selected with equal probability, can carry out one of three actions.

1.Division. Uninfected cancer cells divide by mitosis with probability

$$P_{div}(\stackrel{}{x})=1\mathrm{exp}\left[\left(\frac{N(\stackrel{}{x})}{\sigma _c(\stackrel{}{x})\theta _{div}}\right)^2\right]$$
(2)
determined by the concentration per cancer cell of the nutrients $`N`$ necessary for DNA synthesis. The parameter $`\theta _{div}`$ controls the shape of $`P_{div}`$. If the chosen cell is inside the tumor, its daughter will pile up at that site, and $`\sigma _c(\stackrel{}{x})\sigma _c(\stackrel{}{x})+1`$. Otherwise, if the selected cell is on the tumor border, its daughter will occupy at random one of their normal or necrotic nearest neighbor sites $`\stackrel{}{x^{}}`$. Thus, $`\sigma _c(\stackrel{}{x^{}})=1`$ and $`\sigma _{n,d}(\stackrel{}{x^{}})=0`$. In contrast, an infected cancer cell does not divide since its slaved cellular machinery is focused on virus replication.

2.Migration. Only uninfected cancer cells migrate with probability

$$P_{mov}(\stackrel{}{x})=1\mathrm{exp}\left[\sigma _c(\stackrel{}{x})\left(\frac{N(\stackrel{}{x})}{\theta _{mov}}\right)^2\right].$$
(3)
A selected cell will move to a randomly chosen nearest neighbor site $`\stackrel{}{x^{}}`$. In terms of cell populations $`\sigma _c(\stackrel{}{x^{}})\sigma _c(\stackrel{}{x^{}})+1`$, $`\sigma _c(\stackrel{}{x})\sigma _c(\stackrel{}{x})1`$, $`\sigma _{n,d}(\stackrel{}{x})=0`$, or $`\sigma _{n,d}(\stackrel{}{x})=1`$ if $`\sigma _c(\stackrel{}{x})=1`$.

3.Death. Uninfected cancer cells die with probability

$$P_{del}(\stackrel{}{x})=\mathrm{exp}\left[\left(\frac{M(\stackrel{}{x})}{\sigma _c(\stackrel{}{x})\theta _{del}}\right)^2\right]$$
(4)
determined by the concentration per cancer cell of the nutrients $`M`$ present on the cell’s microenvironment. This Gaussian distribution has a variance fixed by the parameter $`\theta _{del}`$. Clearly, if the chosen cell dies $`\sigma _c(\stackrel{}{x})\sigma _c(\stackrel{}{x})1`$ and $`\sigma _d(\stackrel{}{x^{}})=1`$ when $`(\sigma _c+\sigma _v)`$ vanishes.

Infected cancer cells mainly die by lysis, and the death due to insufficient nutrient supply is neglected. So, an infected cell die with probability

$$P_{lysis}(\stackrel{}{x})=1\mathrm{exp}\left(\frac{T_{inf}}{T_l}\right)$$
(5)
where $`T_{inf}`$ is the time since the infection and $`T_l`$ is the characteristic period for cell lysis. Again, after lysis $`\sigma _v(\stackrel{}{x})\sigma _v(\stackrel{}{x})1`$ and $`\sigma _d(\stackrel{}{x^{}})=1`$ when $`(\sigma _c+\sigma _v)`$ vanishes.

Virus spreading.—The lysis of each infected cancer cell releases $`v_0`$ virus to the extra-cellular medium. A random fraction of these virus remains on the site of the lysed cell, and the rest is equally distributed among its nearest neighbors. Taken these sources into account, the time evolution of virus concentration is given by the quasi-stationary solutions of the discrete reaction-diffusion equation

$$v(\stackrel{}{x},t+1)=v(\stackrel{}{x},t)+\frac{D_v}{4}\underset{<\stackrel{}{x},\stackrel{}{x^{}}>}{}[v(\stackrel{}{x^{}},t)v(\stackrel{}{x},t)]\gamma _vv(\stackrel{}{x},t).$$
(6)
Here, $`D_v`$ is the virus diffusion constant, $`\gamma _v`$ is the viral clearance rate and the sum, representing the discrete Laplacian on the square lattice, extends over all the nearest neighbors of the site $`\stackrel{}{x}`$. In all simulations, Eq.( 6) was iterated $`100`$ times at each time step, and $`D_v`$ and $`\gamma _v`$ were varied on the interval $`[0,1]`$. Periodic boundary conditions along the y-axis and Neumann boundary conditions at the border of the tissue ($`i=L`$) and at the capillary vessel ($`i=0`$) were used for the virus.

Free virus can infect a tumor cell with probability

$$P_{inf}(\stackrel{}{x})=1\mathrm{exp}\left[\left(\frac{v(\stackrel{}{x})}{\sigma _c(\stackrel{}{x})\theta _{inf}}\right)^2\right],$$
(7)
a sigmoidal function of the local virus concentration per cancer cell controlled by the parameter $`\theta _{inf}`$, the efficacy of virus infection.

The virus therapy begins when the tumor attains a growth $`N_0`$ and consists in a single virus injection. In the direct intratumoral (intravenous) administration an uniforme virus concentration $`V=v_0=1`$ over the entire tumor (tissue) is supplied. This approach corresponds to the experimental protocols used in severe combined immune deficient (SCID) mice Coffey  and in vitro essays Raj ; Bischoff .

Results and discussion.—Simulations in lattices with linear size $`L=500`$ have been done for each morphology (compact, papillary and ramified) reproduced by our previous model Silvio  and commonly observed in real tumors. The corresponding fixed parameters are reported in table I, whereas those controlling virus therapy, namely, $`D_v`$ (virus diffusion coefficient), $`\gamma _v`$ (viral clearance rate), $`\theta _{inf}`$ (efficacy of virus infection) and $`T_l`$ (characteristic time for cell lysis) were varied. Since the main results are qualitatively the same for all morphology types, only those obtained for compact tumors will be reported here for brevity.

Figure 1 shows the total number of uninfected cancer cells as a function of time for different virus diffusion and clearance rates. Without treatment, the progress in time of cancer cell populations for all the simulated patterns grows accordingly Gompertz curves Silvio . In the case of a succesfull therapy the tumor mass is completely eradicated a few time steps after virus injection. However, it is worthwhile to notice that in addition to the virus oncolytic activity ($`\theta _{inf}`$ and $`T_l`$), its spreading properties ($`D_v`$ and $`\gamma _v`$) are key features which determine success of the treatment. Indeed, a smaller virus diffusion constant and a greater clearence rate lead to an unsuccesfull therapy. Moreover, after a significant inicial decrease in uninfected cancer cells following virus injection, the tumor grows faster than without treatment. Such result agrees with the claim that cancer cells facing severe threats can trigger adaptation mechanisms as an integrated defense program similar to those observed in bacterial colonies Israel . A tumor submitted to therapeutic approaches which do not lead to its complete eradication may pogressively become more resistant, aggressive and malignant (more invasive), as it is well known in the clinical practice. These results address the central issue concerning the effect of the host immune response on virus therapy, which for human adenoviruses can not be investigated through any appropriate animal model Bischoff . An immune response directed at free virus and late viral antigens expressed on the surface of infected cancer cells can be modeled by an increase in $`\gamma _v`$ and the introduction of an additional probability to the death of infected cells, respectively. For very small $`\gamma _v`$ values the tumor is eradicated with an unique dose, in agreement with the essays performed for immune deficient (SCDI) Coffey  and athymic Bischoff  mice. In turn, for the same dose but higher values of $`\gamma _v`$, more frequent virus injections are necessary to eradicate the tumor, again as it was observed for syngeneic immune-competent C3H mice in Ref. Coffey .

A striking and novel result shown in Figure 1 is the possibility of an oscillatory behavior for both virus and cancer cell populations. It contrasts to the long-term, spatially uniform profiles for cells and free virus reached through damped oscillations obtained by Wu et al. Wu . Instead, the total number of virus, uninfected and infected cancer cells undergo growing (unstable) oscillations, as shown in Figure 2. Similar patterns exhibiting necrotic regions behind expanding fronts of infected cancer cells have been found by Heise et al. Heise  (Fig. 2(e)). An oscillatory dynamics may turn innappropriate to monitor the therapy at predefined time intervals, difficulting clinical prognosis. A central feature of such oscillations is that they are not generated by whatsoever time delay between virus infection and cell lysis, as required in Ref. Wein . In fact, our simulations demonstrate that for a vanishing $`T_l`$ this oscillatory behavior is sustained and, mainly, the effectiveness of virus therapy is largely increased. Furthermore, in this region of parameter space, a single intratumoral virus injection can result in a complete tumor regression with a continously decreasing probability as the clearance rate increases. Thus, the experimentally observed outcomes for reovirus Coffey  and adenovirus Bischoff  therapies, in which the complete regression of tumors occurs in 60 to 80$`\%`$ of the mice, seems to correspond to that region of the model parameter space. The complete parameter subspace $`\theta _{inf},\gamma _v`$ will be reported elsewhere.

Additionally, our results suggest that virus treatment may be more effective against highly invasive papillary tumors (e. g. glioblastoma multiforme or trichoblastoma), since the cells in such tumors grow slowly due to nutrient limitation. This finding contrast with the conjecture that very low cell densities at the periphery of the ramified tumors may prevent further spreading of the wave of virus infection Wein . The differences between the predictions of our model and those in Refs. Wu ; Wein , namely, oscillations without time delay and sustained wave infection at low cell densities, probably arise from the combined use of macroscopic equations and stochastic microscopic dynamics for discrete tumor cells, instead of partial differential equations (PDEs). The interplay between interactions and fluctuations inherent to the discrete microscopic components of the system may lead to collective behaviors unrevealed by a PDE based model Louzoun .

Figure 3 provides further comparison among model predictions and data for virus therapy of tumors in mice. It shows the evolution in time of the simulated number of uninfected tumor cells after virus infection for distinct virus oncolytic activity. The insets show the tumor population and volume size ratios. These ratios are defined as the uninfected cancer cell numbers and tumor volumes normalized by their values at the time of virus injection, respectively. Also, in order to build this figure a time for cancer cells to pass through the DNA synthesis phase of the cell cycle of about $`10`$Rew  was assumed, meaning that the time steps in the simulations correspond to about $`45`$ h. In the inset (b) of Fig. 3, the experimental data correspond to the effects of three different virus (a wild-type and UV-inactivated adenovirus, and a mutant adenovirus dl1520) on C33A human cervical tumors xenografts grown subcutaneously in SCID mice Bischoff . These experiments may be simulated in our approach by associating distinct values for $`\theta _{inf}`$ to different virus, whereas the virus diffusion constants and clearance rates are kept fixed since only one unique experimental model (mouse and implanted tumor) is considered. Remarkably, the predictions of the present model are in good quantitative agreement with the experimental data.

In summary, a model combining macroscopic reaction-diffusion equations for the nutrients and virus with microscopic stochastic rules governing the actions of the cancer cells was proposed in order to describe the effects of virus on tumor growth. The simulational results of the model are consistent with the experimental findings for virus therapy. Moreover, they reveal an oscillatory (periodic or aperiodic) response of tumor cells and virus populations on an intermediate region of the model parameter space ($`D_v`$, $`\gamma _v`$). In this region, a single intratumoral injection of virus can result in a complete tumor regression with a continously decreasing probability as the clearance rate increases. Such oscillations, which may difficulty clinical prognosis, suggest the need for new in vivo and in vitro experiments aiming to detect this oscillatory response.

This work was partially supported by the Brazilian Agencies CAPES, CNPq, and FAPEMIG.


format: text

text:
# BI-TP 2003/25 SWAT/03/379 Where is the chiral critical point in 3-flavor QCD? This work is supported by BMBF under grant No.06BI102, DFG under grant FOR 339/2-1 and PPARC grant PPA/a/s/1999/00026.

## 1 Introduction

Sufficiently hot and dense hadronic matter undergoes a transition to a deconfined, chirally symmetric medium. It is expected that in QCD with the correct, physically realized spectrum of quark masses this transition is a true phase transition, related to non-analyticities of the QCD partition function, only for specific choices of temperature ($`T`$) and baryon chemical potential ($`\mu _B3\mu _q`$). The transition then is first order for large values of $`\mu _B`$, ends in a second order endpoint – the chiral critical point – at a certain value $`\mu _B^{crit}`$ and is a crossover transition, not related to singular behavior in thermodynamic observables, for smaller values of the chemical potential (Fig. 1).

The chiral critical point is part of a critical surface in the 3-d parameter space of degenerate light ($`m_{u,d}`$) and a strange ($`m_s`$) quark mass as well as $`\mu _B`$. At $`\mu _B=0`$ this surface contains the line of second order phase transitions, which separates the region of first order phase transitions for light quarks from a broad crossover region.

At $`\mu _B=0`$ we determine the location of the 3-flavor critical point on this transition line. Using reweighting in $`\mu _B`$ we determine the curvature of the critical surface and estimate the location of the chiral critical point, $`\mu _B^{crit}`$. All calculations have been performed with improved staggered fermions (p4-action) and a Symanzik-improved gauge action .

## 2 Critical point at $`\mu =0`$

In previous studies with unimproved staggered fermions (see, for instance ) the chiral critical point in 3-flavor QCD has been located at a pion mass value of about 300 MeV. Using Binder cumulants constructed from the chiral condensate, $`B_4(0)=(\delta \overline{\psi }\psi )^4/(\delta \overline{\psi }\psi )^2^2`$, it has been verified that this point belongs to the Ising universality class . Using the knowledge of the universality class as an input, first calculations with improved staggered fermions (p4-action) on a $`12^3\times 4`$ lattice suggested that the critical point shifts to substantially smaller values. In  a preliminary value of $`m_\pi ^{crit}190`$ MeV has been reported.

Our present analysis with the p4-action is based on simulations on $`12^3\times 4`$ and $`16^3\times 4`$ with a bare quark mass of $`ma=0.005`$ which corresponds to a pion mass of about 170 MeV. A finite size scaling analysis should allow to draw conclusions on the nature of the transition. At this light pion mass we find no indication for a significant finite size dependence of the peak height of chiral ($`\chi _{\overline{\psi }\psi }`$) and Polyakov loop ($`\chi _L`$) susceptibilities,

$$\chi _{\overline{\psi }\psi }=(\overline{\psi }\psi )^2\overline{\psi }\psi ^2,\chi _L=L^2L^2.$$
(1)
Furthermore, the Binder cumulant $`B_4(0)`$ is significantly larger than the 3d-Ising value, $`B_4^{Ising}=1.604`$. This indicates that the chiral critical point has to lie at even smaller values of $`m_\pi `$. Reweighting of both observables in the quark mass around the simulation point at $`ma=0.005`$ suggests a critical value of $`ma=0.0007(4)`$. At this value of the quark mass the Binder cumulant is consistent with the Ising value and the peak heights in the chiral susceptibility show a significant volume dependence which is consistent with $`\chi _{\overline{\psi }\psi }^{\mathrm{peak}}V`$. Some results for the reweighted chiral condensate and $`\chi _{\overline{\psi }\psi }`$ obtained on the $`16^3\times 4`$ lattice are shown in Fig. 2. The rather small value found for the critical quark mass translates into a small value of the pion mass, which is only about half its physical value. Such a small value is consistent with the small values typically found in effective model calculations and follows the trend towards a smoother transition found in other calculations with improved staggered fermion actions . It differs significantly from the result obtained with standard, unimproved staggered fermions. The actual location of the chiral critical point at $`\mu _B=0`$ thus seems to be strongly cut-off dependent. In summary we find,

$$m_\pi ^{\mathrm{crit}}=\{\begin{array}{cc}\hfill 290(20)\text{MeV},& \text{standard action}\hfill \\ \hfill 67(18)\text{MeV},& \text{p4 improved action}\hfill \\ \hfill 47\text{MeV},& \text{linear sigma model}\hfill \end{array}$$
(2)
## 3 Reweighting in $`\mu _𝐪`$ at $`\mathrm{𝐦𝐚}=\mathbf{0.005}`$

At fixed value of the quark mass we use  Taylor expanded reweighting to $`𝒪(\mu _q^2)`$ of $`\chi _{\overline{\psi }\psi }`$ and $`B_4(0)`$ to locate the chiral critical point as that point at which $`B_4`$ attains the 3d-Ising value. The reweighting is done only in the 2-flavor sector for (u,d)-quarks, i.e. we keep $`\mu _{strange}=0`$.

Reweighting was possible only for the data on the $`12^3\times 4`$ lattice as reweighting to large enough values of $`\mu _q`$ was prohibited on the larger lattice due to too large fluctuations in the sign of the fermion determinant. Even on the smaller lattice a tendency of $`B_4`$ to decrease to the Ising value only sets in for $`\mu _qa>0.05`$. A linear extrapolation for $`\mu a0.05`$ yields $`\mu _q^{crit}a=0.074(13)`$ or $`\mu _q^{crit}/T=0.296(52)`$. An estimate for the transition temperature at the corresponding pion mass is obtained from  $`T_c/\sqrt{\sigma }=0.40(1)+0.039(4)(m_\pi /\sqrt{\sigma })`$. This sets the scale for our estimate of the critical value of the chemical potential in 3-flavor QCD, $`\mu _q^{crit}=52(10)\mathrm{MeV}`$ at $`m_\pi =170\mathrm{MeV}`$. A quadratic interpolation to the physical pion mass value fixes the location of the chiral critical point in 3-flavor QCD,

$$\mu _q^{crit}=40(9)\mathrm{MeV},m_\pi =140\mathrm{MeV}.$$
(3)
As expected the value, $`\mu _B^{crit}=120\mathrm{MeV}`$, found here in 3-flavor QCD is substantially smaller than the value of about 700 MeV estimated for (2+1)-flavor QCD .

## 4 The $`\mu `$-dependence of the transition line

Following the prescription outlined in  the $`\mu `$-dependence of the transition temperature has been obtained using Taylor expanded reweighting of susceptibilities to $`𝒪(\mu _q^2)`$. From the shift of the peak positions one finds the critical couplings, $`\beta _c(\mu _q)`$, which to leading order depend linearly on $`\mu _q^2`$ and yield an estimate for the slope of $`T_c(\mu _q)`$:

$$\frac{\mathrm{d}^2T_c}{\mathrm{d}\mu _q^2}=\frac{1}{N_\tau ^2T_c(0)}\frac{\mathrm{d}^2\beta _c}{\mathrm{d}\mu _q^2}/\left(a\frac{\mathrm{d}\beta _c}{\mathrm{d}a}\right)$$
(4)
Comparing calculations in 3-flavor QCD for $`ma=0.005`$ and 0.1 we find a significant increase in the slope of $`T_c(\mu _q)`$ with decreasing $`ma`$. Using the asymptotic $`\beta `$-function in (4) gives

$`{\displaystyle \frac{T_c(\mu _q)}{T_c(0)}}=\{\begin{array}{cc}\hfill 10.025(6)\left(\frac{\mu _q}{T_c(0)}\right)^2,& ma=0.1\hfill \\ \hfill 10.114(46)\left(\frac{\mu _q}{T_c(0)}\right)^2,& ma=0.005\hfill \end{array}`$

We note that an additional quark mass dependence is hidden here in the transition temperature at $`m_q=0`$, i.e. $`T_c(0)`$, which is used to normalize the transition temperature at $`\mu _q0`$. Furthermore, taking into account violations of asymptotic scaling in the $`\beta `$-function will lead to a further increase of the curvature of $`T_c(\mu _q)`$. At $`ma=0.01`$ deviations from asymptotic scaling have been analyzed in detail and suggest that the slope increases by almost a factor of 2 .

The transition lines for the two different quark masses together with the estimates for the location of the chiral critical points in 3-flavor QCD (left and right boxes) and the ($`\mu _q`$-independent) transition line in the infinite quark mass limit are shown in Fig. 3.

## 5 Conclusions

At vanishing chemical potential the cut-off dependence of the chiral critical point is large. Calculations with an improved action suggest that the chiral transition in 3-flavor QCD becomes first order at much smaller values of the pion mass than estimated previously. These findings are consistent with other studies of QCD thermodynamics performed with improved staggered fermions . A similar influence of cut-off effects is found for the dependence of the chiral critical point on the chemical potential.


text:
# Hardware and software status of QCDOC Based on talks by K. Petrov and T. Wettig at Lattice 2003, Tsukuba.

## 1 Introduction

Continued advances in commodity processing and networking hardware make PC (or workstation) clusters a very attractive alternative for lattice QCD calculations . Indeed, there are quite a few important problems that can be addressed on PC clusters, and many lattice physicists are taking advantage of this opportunity. However, for the most demanding problems in lattice QCD, e.g. dynamical fermion simulations with realistic quark masses, one would like to distribute the global volume over as many nodes as possible, resulting in a very small local volume per node. PC clusters are inadequate to deal with this case because the communications latency inherent in their networking hardware implies that the local volume must not be chosen too small if a reasonable sustained performance is to be achieved. In other words, for typical lattice QCD problems PC clusters do not scale well beyond a few hundred nodes.

In custom-designed supercomputers such as QCDOC  and apeNEXT , the communications hardware is designed to reduce the latencies and to assist critical operations (such as global sums) in hardware. As a result, these machines are significantly more scalable and allow for much smaller local volumes. In addition, they provide low power consumption, a small footprint, and a very low price/performance ratio per sustained MFlops. On the downside, the development effort is considerably higher than for PC clusters, but this effort is amortized by the unique strengths of these machines.

## 2 Hardware

The QCDOC hardware has been described in detail in several previous publications, see Refs. , therefore we only summarize its most important features here.

The QCDOC ASIC, shown schematically in Fig. 1, was developed in collaboration with IBM research and manufactured by IBM. It contains a standard PowerPC 440 core running at 500 MHz, a 64-bit, 1 GFlops FPU, 4 MBytes of embedded memory (eDRAM), and a serial communications interface (SCU) which has been tailored to the particular requirements of lattice QCD. The SCU provides direct memory access, single-bit error detection with automatic resend, and a low-latency pass-through mode for global sums. Also on the chip are several bus systems, controllers for embedded and external (DDR) memory, an Ethernet controller, a bootable Ethernet-JTAG interface, and several auxiliary devices (interrupt controller, I2C interface, etc.). A picture of one of the first ASICs, delivered in June of 2003, is shown in Fig. 2.

The physical design of a large machine is as follows. Two ASICs are mounted on a daughterboard, together with two standard DDR memory modules (one per ASIC) with a capacity of up to 2 GBytes each. The only other nontrivial components on the daughterboard, apart from a few LEDs, are four physical layer chips for the MII interfaces (two per ASIC) and a 4:1 Ethernet repeater which provides a single 100 Mbit/s Ethernet connection off the daughterboard. A picture of the very first two-node daughterboard is shown in Fig. 3.

A motherboards holds 32 such daughterboards, eight motherboards are mounted in a crate, and a large machine is built from the desired number of crates. A picture of a QCDOC motherboard is shown in Fig. 4.

The physics communications network of QCDOC is a 6-dimensional torus with nearest-neighbor connections. The two extra dimensions allow for machine partitioning in software so that recabling is not required. A 64-node motherboard has a $`2^6`$ topology, with three open dimensions and three dimensions closed on the motherboard (one of which is closed on the daughterboard). The SCU links run at 500 Mbit/s and provide separate send and receive interfaces to the forward and backward neighbors in each dimension, resulting in a total bandwidth of 12 Gbit/s per ASIC (of which 8 Gbit/s will be utilized in a 4-dimensional physics calculation).

In addition to the physics network, there is an Ethernet based network for booting, I/O, and debugging, as well as a global tree network for three independent interrupts. The Ethernet traffic from/to each motherboard proceeds at 800 Mbit/s to a commercial Gbit-Ethernet switch tree, a parallel disk system, and the host machine. The latter will be a standard Unix SMP machine with multiple Gbit-Ethernet cards. See Fig. 5.

As of the writing of this article (September 2003), all major subsystems of the QCDOC ASIC have been tested in single-daughterboard configurations (2 ASICs per daughterboard) using a temporary test-jig. This configuration allows non-trivial communications between the two ASICs in one of the six dimensions; for the remaining dimensions the ASIC communicates with itself in a loop-back mode. Extensive memory tests with different sized commodity external DDR SDRAM modules have been done, tests of the 4 MByte on-chip EDRAM have been performed, and all the DMA units have been used. High-performance Dirac kernels have been run for Wilson and ASQTAD fermion formulations, confirming the performance figures given in Table 1 and Ref. . No problems with the ASIC have been found to date.

With QCDOC motherboards now in hand for use (Fig. 4), tests of 64 and 128 node machines are imminent. In our test-jig tests, the ASICs appear to be working well close to the target frequency of 500 MHz. With a fully populated motherboard, and the more stable electrical environment provided by the motherboard as compared to our simple test-jigs, we will soon be able to test a large number of ASICs at the 500 MHz design frequency. From preliminary measurements, the power consumption per node is about 5 W.

## 3 Software

One of the major goals of the development team was to make QCDOC accessible and deployable to a large scientific community, so every attempt has been made to allow users to use standard software tools and techniques. In order to reliably boot, monitor and diagnose the 10,000+ nodes in a large QCDOC machine, and to allow user code maximum access to the high-performance capabilities of the hardware, we are writing a custom operating system for QCDOC. This is modeled on the QCDSP operating system, which has successfully allowed the use of a machine with a similar number of processing nodes. The QCDOC operating system (QOS) has been improved in myriad ways, continuing to focus on ease-of-use for the end-user.

### 3.1 Operating System

We now give a more detailed list of the requirements of the operating system and discuss how these features have been implemented. In addition to booting the machine, the QOS must diagnose hardware errors and monitor the hardware during user program execution. Monitoring is accomplished by regular inspection of on-chip registers which monitor memory and communication link status, along with interrupt handlers which are invoked by changes in hardware state deemed crucial. After booting, partitions of the machine may be requested for interactive use by a user or the queuing system. Each partition will only run a single user program at a time; all multitasking will be handled by the host front-end. During program execution, the operating system provides host-QCDOC input/output access using commodity Ethernet connections, access to special QCDOC features such as the mesh communication network and on-chip memory and access from QCDOC processing nodes to a parallel disk array. The QOS is designed to coordinate the interactions between the host and QCDOC, insulating the user from many QCDOC-specific details while still providing detailed control over their execution environment.

Our solution is QOS, written largely from scratch in C++, C and some assembler. The software runs on the host computer and the QCDOC nodes and consists of several parts, namely:

* The Qdaemon on the host, which manages and monitors the entire QCDOC machine. All interactions go through the Qdaemon.
* The Qclient layer on the host, which provides access to the Qdaemon for a variety of planned user interfaces, such as the batch system and Web applications, as well as the existing command line interface, the Qcsh.
* The Qcsh on the host, a modified version of tcsh, which is currently in use and allows complete control and use of QCDOC.
* The boot/run kernels on each node.

The Qdaemon is the heart of the operating system, and it is designed and implemented using modern C++ techniques, e.g. it is heavily templated, it employs POSIX threads to allow efficient use of a multi-processor front end, its queue calls and object lists are thread-safe, and it can drive multiple physical network interfaces. The single Qdaemon on the host uses reliable UDP to transport packets between QCDOC and the host, starting with basic UDP packets in the boot procedure and adding an RPC protocol to this when the run kernels are loaded. The Qdaemon controls the IP addressing scheme to the nodes, does the partitioning that users request, monitors machine status, directs I/O to the host and evokes hardware testing when needed. The Qdaemon can be accessed by a “root” user for system management and monitoring, but most users will only have access to a limited set of Qdaemon’s capabilities, those appropriate for the partition where they are running a job. The Qdaemon also provides flexible software partitioning of QCDOC.

The Qclient is a library of interfaces which communicates with the Qdaemon and can accept input from a variety of planned user tools. To access the Qclient, and hence the Qdaemon, access to the host computer must be gained by conventional authentication tools such as ssh or a secure Web connection. Currently, a command-line interface, the Qcsh, communicates with the Qdaemon via the Qclient library.

Following the successful QCDSP model, Qcsh is a command-line interface. Starting from a standard unix shell, tcsh, extra built-in commands to control QCDOC have been added. Users can use normal shell redirection to control I/O from QCDOC, programs can be executed on QCDOC and backgrounded in Qcsh allowing further interaction between the Qcsh and Qdaemon (but not QCDOC directly since it is executing a program), and standard shell scripts can be used. The special QCDOC commands in the Qcsh all begin with the letter q. Some examples are: qinit \- establish a connection between the Qcsh and the Qdaemon; qpartition\_create \- create a partition of QCDOC; qboot \- boot a QCDOC partition; qrun \- run a user program; qhelp \- display available Qcommands.

The Qdaemon will only permit connections from the local host, which allows a Unix socket to be passed to the Qdaemon from a user interface like the Qcsh. The Qdaemon then does user I/O from/to this Unix socket and can determine Unix user and group identities directly from the host OS kernel. This is a major simplification. While the Qdaemon will be handling many users’ access to QCDOC, it does not have to open/close files and concern itself with protections. This is handled by the initial user interface, thereby ensuring correct ownership and permissions.

The first software loaded to a QCDOC node is the boot kernel, which is loaded directly to the data and instruction caches of the 440 core via the Ethernet-JTAG hardware. Only the 440 core must be functional for the boot kernel to execute. When execution begins, tests of the on-chip and DDR memory can be done. Once the memory is known to be working, the boot kernel enables the standard 100 Mbit Ethernet port on the QCDOC ASIC and the run kernel is loaded down. The run kernel handles the activation of machine-global features such as the nearest-neighbor communications network. When these steps complete, which should be on the scale of ten minutes even for a large machine, it is available for general use. All of the steps outlined above are implemented and have been extensively used on our first ASICs.

The run kernel provides support for users, including access to QCDOC features via system calls. Part of the strategy for using RPC for the host–run-kernel communication is to avoid having a separate communications protocol for disk access. We have written an NFS client for the run kernel that uses the standard RPC-based NFS protocol. The client supports two mount points and open/read/write/close functionality. This NFS support has already been tested on QCDOC. While providing a standard user environment for high-performance computing, the run-kernel will not support multi-tasking. We have chosen to keep the run-kernel lean and compact to ensure reliability and to keep the software task bounded.

### 3.2 Communications Software

A major issue for a massively parallel computer is the effectiveness of its message passing. In hardware, QCDOC is a nearest-neighbor mesh, the majority of QCD based communications are also nearest-neighbor and consequently QOS natively supports nearest-neighbor communications calls. The QOS calls are implemented so as not to interfere with the low latency of the QCDOC communications hardware. The QCDOC hardware supports efficient global sums (done via our pass-through hardware), which are also accessible via QOS calls. We have also implemented the SciDAC Lattice QCD Software Committee’s QCD message passing (QMP) protocol on QCDOC. This protocol supports nearest-neighbor communications, which we efficiently map to the native QOS communications calls, as well as arbitrary communications. For the latter, we will implement a Manhattan-style routing.

### 3.3 Existing Software and Benchmarks

A simple (but not trivial) consequence of QCDOC using a standard processor from IBM’s PowerPC line is the availability of a whole arsenal of Open Source and commercial software tools, most notably the GNU tools and IBM’s tools. On QCDOC, users will be able to use the GNU toolset, the closest thing to a standard across computer platforms. Additionally, we have access to the high-performance commercial tools from IBM, such as the xlc/xlC compilers which we have seen outperform the GNU compilers on PowerPC platforms by as much as a factor of two. We expect that users may do initial code development with the GNU toolset and only compile with the more restrictively-licensed IBM compilers when final performance is an issue.

An issue for the larger potential user group for QCDOC is performance of existing physics codes, not written for a QCDOC or QCDSP type of computer architecture. The MILC code was an obvious choice to test, since it is one of the major lattice simulation codes and is written in C. The MILC code has been run on QCDOC, both the simulator and now the actual ASIC, concentrating on the ASQTAD action. Unmodified MILC code gives performance in the few percent range for small lattice volumes, which was easily improved by a few standard C-code modifications as described in . A summary of the performance for single precision MILC code and our double precision assembly code is given in Table 1.

## 4 Schedule

A 128-node prototype machine is currently being assembled at Columbia. Assuming that no major problems are found, two machines of 10 TFlops (peak) each will be available to UKQCD and the RIKEN-BNL Research Center by the summer of 2004. The operating system and user software needed to utilize these machines is progressing in hand with the hardware developments. Approval is pending for a 3 TFlops (peak) machine at Columbia and a $``$ 20 TFlops (peak) installation at BNL for the U.S. lattice community.

Acknowledgments. This work was supported in part by the U.S. Department of Energy, the Institute of Physical and Chemical Research (RIKEN) of Japan, and the U.K. Particle Physics and Astronomy Research Council.


text:
# HISTORY OF STAR FORMATION RATE AND LUMINOSITY DENSITY OF GALAXIES

## INTRODUCTION

In recent years there have been several studies of the cosmological history of star formation rate (SFR) in galaxies. Usually the SFR is derived from the luminosity density in a selected wavelength band which is an indicator of star formation. Some of the SF indicators that have been used are: H$`\alpha `$ line, far-infrared (FIR), submm, ultraviolet (UV) and nonthermal radio continuum emissions. Each indicator is sensitive to a certain range of stellar masses and the total SFR is obtained by assuming an initial mass function (IMF) of stars at birth, usually a power-law form. Further, the observed luminosity depends on the star formation history and the epoch of observation. Even for a constant rate of star formation, there is a time evolution of luminosity which increases, reaches an equilibrium value and declines after the termination of SF. This evolution and decline also depend on the exponent of the IMF. The changes in the luminosity to SFR ratio as a function of these parameters are not the same for all indicators. Thus, apart from the observational errors, uncertainties arise in the SFR derived from luminosity density using different indicators. In this paper, we study these uncertainties by computing the evolution of the luminosity to SFR ratio for different IMFs and scenarios of star formation. It should be noted that there are other parameters like metallicity that can affect the derivation of SFR from luminosity. For example, Hirashita et al (2001) find that the FIR-SFR conversion factor can vary by a factor of up to 4 on this account. When we have better knowledge of factors like metallicity and cirrus contribution in starburst galaxies, it may be possible to quantify these uncertainties.

## COMPUTATION

In this study we compute the time evolution of L/SFR ratio for the following luminosities: Lbol, the bolometric luminosity, Lfir, the far-infrared luminosity which is a fraction of the bolometric luminosity, Luvi, the intrinsic luminosity in an ultraviolet band, Luve, the luminosity of uv photons escaping the star forming region, LH, the H$`\alpha `$ line luminosity and LR, the nonthermal radio continuum luminosity. For UV we use the 0.16 $`\mu `$m band which has been extensively used in HST observations. The submm luminosity is usually converted into a FIR luminosity to derive the SFR and hence we do not consider it separately. The radio continuum is assumed to be from synchrotron emission from cosmic ray electrons accelerated by supernovae arising from stars of mass $`>`$ 8 M. In starburst(SB) galaxies, which are the subjects of our study, the contribution of the earlier disc population is small and is neglected. Their contribution, if present, will add additional uncertainties. The FIR luminosity is dust-reradiated emission and in many studies, it has been taken to be the same as the bolometric luminosity. Though this is likely to be the case at the very early embedded stage of a star, the efficiency of conversion will be less than unity and decrease as time passes because of the disruption and dispersion of the parent clouds by the massive stars formed. The fact that visible OB stars are observed in star forming complexes in the Milky Way supports this hypothesis. Silva et al. (1998) have introduced the concept of residence time and have shown that their models give good fits to observational spectral energy distributions and spectral line data. In this study, we assume that the FIR conversion efficiency decreases exponentially with a time constant $`\tau `$. It then follows that the flux of UV photons escaping the star forming region depends inversely on this conversion efficiency and can be characterized by an extinction which decreases exponentially with the same time constant. The H$`\alpha `$ luminosity will also have a similar dependence with a proper scaling of the extinction.

The total luminosity in a given band is the integral over the history of stellar formation and evolution and over the stellar masses. The IMF is assumed to have a power law form $`\xi `$(m) $``$ m within mass limits of 120 to 1 M. Since the contribution of stars of mass $`<`$ 1 M to the luminosity, in the bands considered, is negligible, the L/SFR ratio can easily be scaled for an IMF extending to lower masses. The stellar evolutionary tracks are taken from Schaller et al. (1992) and the computational procedure follows that of Mayya (1995). The integrated luminosity at time T, the epoch of observation is given by

$$L_{bol}(T)=_1^{120}\xi (m)𝑑m_T^0p(t)L_{bol}(t,m)𝑑t$$
(1)
$$L_{fir}(T)=_1^{120}\xi (m)𝑑m_T^0p(t)L_{bol}(t,m)e^{t/\tau }𝑑t$$
(2)
$$L_{uvi}(T)=_1^{120}\xi (m)𝑑m_T^0p(t)L_{uv}(t,m)𝑑t$$
(3)
$$L_{uve}(T)=_1^{120}\xi (m)𝑑m_T^0p(t)L_{uv}(t,m)10^{0.4A(uv)exp(t/\tau )}𝑑t$$
(4)
$$L_H(T)=K_H_1^{120}\xi (m)𝑑m_T^0p(t)L_{Lyc}(t,m)10^{0.4A(H)exp(t/\tau )}𝑑t$$
(5)
$$L_R(T)=K_R_8^{120}\xi (m)𝑑m_T^{Min(T,t^{}(m))}p(t)e^{t/\tau _e}𝑑t$$
(6)
Here t is the look back time when a star is born, T is the observational epoch from a reference time when the starburst started, p(t) is the production rate of stars at time t, t(m) is the lifetime of a star of mass m, A(uv) and A(H) are the extinctions in magnitudes at 0.16 $`\mu `$m and H$`\alpha `$ wavelength at the time of birth, Lyc refers to Lyman continuum photons and KH is a constant that relates H$`\alpha `$ luminosity to the Lyman continuum luminosity and KR is the nonthermal radio continuum luminosity per supernova and $`\tau _e`$ is the exponential residence time of synchrotron emitting electrons. The production term p(t) represents the time evolution of the starburst. It is a $`\delta `$ function for an instantaneous burst, constant for a continuous and constant star formation and could have other forms like an exponential decay or a gaussian.

The IMF index is neither known accurately nor is it known to be universal. Observations of stellar clusters in the Milky Way and the Magellanic clouds give a range of 2 to 3 ( Scalo 1998). We, therefore, perform computations for three values of $`\gamma `$ = 2, 2.5 and 3. As for the time dependence of star formation, we will consider, in this paper, only the case of p(t) = constant for t = 0 to TB and zero elsewhere and present the results for $`T_B`$ = 95 Myr. As for the residence time of stars in the parent clouds, Silva et al. (1998) find that for the starburst galaxies they model, a value of 10–20 Myr fits the observations. Here, we take $`\tau `$ = 15 Myr. The extinction at 0.16 $`\mu `$m and H$`\alpha `$ are taken as 10 and 3 magnitudes respectively. The UV extinction adopted leads to a luminosity weighted extinction factor of 5-6 at times $`>`$ 50 Myr similar to the value estimated by Steidel et al. (1999). It may be noted that for moderately obscured regions, a more elaborate attenuation model is needed. However, for the luminous starburst galaxies, especially at moderate and large redshifts, our treatment may be adequate. As remarked earlier the presence of older star-forming regions will add further uncertainty to the SFR. For the residence time of high energy electrons, we take $`\tau _e`$ = 10 Myr.

## 1 RESULTS AND DISCUSSION

In Figure 1, we show for different SF indicators, the time evolution of the ratio of luminosity at time T to that at 95 Myr when the SB terminates. We have not shown the curve for H$`\alpha `$ since it is very similar to that of the FIR. For FIR and UVE, we show the ratio with respect to the bolometric and intrinsic uv luminosity respectively. The general trend of the time evolution is an increase of the luminosity to a plateau value and then a decrease after the SB terminates.The plateau values of the ratio of luminosity to SFR for the bolometric, FIR and uv (0.16 $`\mu `$m)emissions are 1.6$`\times `$1010 LM-1 Yr, 8.8$`\times `$109 LM-1 Yr and 5.5$`\times `$1020 W Hz-1 M-1 Yr respectively. These values, after adjustment for the Salpeter IMF and the 0.1–100 M range adopted by Kennicutt ( 1998),are about the same as those given in his review. Further, the ratio of FIR to radio continuum fluxes is similar to that of Bressan et al (2002) if we adopt their value of nonthermal emission per supernova.

The rise- and fall-times vary from band to band. The time profile for the FIR luminosity is the sharpest while it is the broadest for the UVE. This is due to the fact that luminous massive stars have lifetimes $`<<`$ the residence time and almost the total bolometric luminosity is re-radiated in the FIR while for stars with lifetime $`>`$ 15 Myr, the conversion efficiency decreases. The opposite is the case for the escaping uv photons. For $`\tau `$ = 15 Myr, the maximum conversion efficiency is 0.55. It will be lower if $`\tau `$ is less. The fraction of escaping uv photons reaches a maximum of 0.17. The uv extinction is more sensitive to the values of $`\tau `$ and extinction at zero age than the FIR conversion efficiency. The slight increase in the UVE luminosity after the termination of the SB and the subsequent oscillation is due to a competition between the death of massive stars and increased escape of uv photons from lower mass stars as the stars age. It may be noted that the FIR profile is sharper than that of the bolometric luminosity while the UVE approaches that of UVI in the post starburst period. The radio continuum luminosity has a rise- and fall-time of about 30 Myr, the lifetime of stars of mass 8 M. The plateau FIR and UVE luminosities depend on $`\tau `$, the residence time. Hence, there will be additional uncertainties in the SFR resulting from the variation of $`\tau `$ from galaxy to galaxy and its possible dependence on the intensity of starburst. It is generally believed that the FIR emission is least subject to extinction correction applicable to other wavelengths. However, the concept of residence time introduces corrections up to a factor of 2–3. From this point of view, the radio continuum emission is a better probe of SFR since it needs no such correction.

It is clear that the SFR derived for an individual galaxy depends on the duration of the SB and the time of observation and the luminosity does not give the instantaneous current SFR unless one observes in the plateau region of continuous and constant star formation. Even if we observe a large enough sample of galaxies, say in a given redshift bin, and all galaxies have the same duration of starburst and other parameters remain the same, the uncertainty is not fully eliminated. If we assume that the SB duration is the same for all galaxies, at any given time one will observe galaxies distributed uniformly over the evolutionary stage. Though the plateau luminosity may be higher than the observational threshold, in the rising portion of the evolutionary curve, galaxies will be missed when the luminosity is below the threshold. In the post starburst period, one will observe galaxies above the threshold though there is no ongoing star formation. For example, if the luminosity is n times the observational threshold and t1 and t2 are the times when the fraction of plateau luminosity is 1/n in the rising and falling stages respectively ( see Fig. 1), the fraction of galaxies that will be missed is t1/TB and the additional contribution from the post SB period is (t2-TB)/$`T_B`$. The effect depends on the SF indicator used, the duration of the SB and the luminosity function in the chosen band. The FIR emission is a good tracer of ’instantaneous’ SFR since both the rise and fall are steep. For the radio continuum, there is a net over-estimation since the decline in the post SB period is slower. Using the 1.4 GHz radio luminosity function of Haarsma et al (2000) for convolution and data of Fig.1, we find the net over-estimation of the luminosity density to be 26 %. The over-estimation in the case of uv photons is much higher. The corrections also depend on the duration of the SB. For low values of SB duration, the plateau region is a small fraction of the total and hence the correction is more. Thus, the uncertainties are significant as long as the SB duration is $`<`$ 150-200 Myr.

### 1.1 Effect of IMF Variation

In Figure 2 we have plotted, the luminosity to SFR ratio normalized to the value for the IMF index $`\gamma `$ = 2.5 vs $`\gamma `$ for three epochs, 30 Myr, 95 Myr and 135 Myr. It is seen that the variation with index is not the same for the different bands. The FIR luminosity shows the maximum variation, a factor of more than 5 as the index changes from 2 to 3. The escaping uv photons show the least change. Further, the changes also depend on the observational epoch. Thus, for individual galaxies, one can have large uncertainties in SFR if the IMF varies. If a large sample of galaxies is available, the average over varying index will result in reduced uncertainty.

It is not known whether the IMF extends from high mass to low mass stars in the same region at the same time. In intense starbursts it is believed that the IMF is biased towards high mass stars. In the Galactic star forming region W3, Megeath et al (1996) find that on extrapolation of the observed IMF for M $`<`$ 10 M, they expect to see only one star of mass $`>`$ 10 M, but radio continuum observations show the presence of six high mass stars. If in starbursts, the IMF does not extend below, say 10 M, the luminosity computed would not change much, but the total stellar mass would be lower. If IMF bias is present and varies from galaxy to galaxy, there will be additional uncertainties in the derived SFR.

## 2 SUMMARY

We have computed the time evolution of luminosity in different bands of star formation indicators for the case of a continuous and constant star formation. The luminosity rises from the time the starburst begins, reaches a plateau value and declines after the termination of the starburst. The post-burst decline is generally, slower than the growth at early times. We have adopted the concept of residence time to calculate the FIR and escaping uv luminosities. The FIR luminosity has the sharpest profile while the uv luminosity has the broadest. The star formation rate derived from luminosity density is subject to uncertainties depending on the shape of the IMF, the residence time parameter, the duration and distribution of starburst and the observational epoch.

Email: renga@inaoep.mx (T. N. Rengarajan); ydm@inaoep.mx ( Y. D. Mayya)

Manuscript received: 28 October 2002; revised: 13 May 2003; accepted: 14 May 2003


text:
# Observational Studies of Early-type Overcontact Binaries: TU Muscae

## 1 Introduction

The variability of TU Muscae (catalog ) (HD 100213) was discovered by Oosterhoff (1928) and the system was first studied in detail by Andersen & Grønbech (1975, hereafter AG75) who list references to earlier work on this early-type overcontact (Wilson, 2001) binary. Given the paucity of eclipsing, double-lined O-type binaries, TU Mus has received surprisingly little attention from observers. AG75 published fully covered $`uvby`$ light curves and fourteen radial velocities for each star near quadratures from photographic spectra of moderate resolution. Stickland, et al. (1995, hereafter S95) presented twenty-four radial velocities for the primary star and twenty-three for the secondary from high resolution International Ultraviolet Explorer (IUE) spectra. TU Mus was bright enough ($`V_T=8.4`$) to be observed by Hipparcos and a light curve with about a hundred points was published in the Hipparcos Epoch Photometry Annex (ESA, 1997). Two analyses of these data using modern light curve synthesis tools have been published: a preliminary report on the current work by Terrell (2002) and an analysis of the AG75 photometry by Wilson & Rafert (1981).

The two sets of radial velocities give quite discrepant results for the absolute dimensions. AG75 obtained photographic spectra with a dispersion of 20 Å mm-1 and found $`K_1`$=251.3 $`\pm `$4.8 km sec-1 and $`K_2`$=370.5 $`\pm `$4.3 km sec-1. S95 found $`K_1`$=216.7 $`\pm `$2.7 km sec-1 and $`K_2`$=345.4 $`\pm `$3.1 km sec-1 from their IUE spectra that had signal-to-noise ratios (S/N) of 10-20. Adopting an inclination of $`76^{}`$, these numbers result in masses of $`M_1=17.2M_{\mathrm{}}`$ and $`M_2=10.8M_{\mathrm{}}`$ for the IUE velocities and $`M_1=23.5M_{\mathrm{}}`$ and $`M_2=15.8M_{\mathrm{}}`$ for the optical velocities, a difference much larger than the error limits allow. In order to explore the nature of this discrepancy, we obtained high-resolution, high S/N spectra of TU Mus near the quadratures. Although we have a limited number of spectra due to the difficulty of obtaining observing time on large telescopes for the study of eclipsing binaries, we find very good agreement with the AG75 results.

Given the good agreement between our velocities and those of AG75, we performed a simultaneous solution of those velocities, the AG75 $`uvby`$ photometry, the Hipparcos photometry and new $`BV`$ photometry from the 2002 observing season. We used the 2003 version of the Wilson-Devinney (WD) program (Wilson & Devinney (1971),Wilson (1979),Wilson (1990)) and present new elements for the binary, showing that the period has remained remarkably constant over the roughly thirty year span of all the observations.

## 2 Spectroscopy

Nine ESO 2.2m + FEROS spectra have been secured, in three groups of three spectra each (Table 1). The wavelength range of the spectra extends from 3900 Å to 9200 Å with $`R=48,000`$ and a 400 second exposure time. The S/N is evaluated around H$`\beta `$ and the orbital phase is computed with the ephemeris given in section 4.

### 2.1 Rotational velocity

AG75 estimated the rotational velocities of the primary and secondary stars in TU Mus as 285 and 240 km sec-1 ($`\pm `$10%), while S95 derived significantly slower rotations: 250 and 195 km sec-1 ($`\pm `$10%), respectively. To derive the rotational velocity we used the relation

$$V_{\mathrm{rot}}=42.42\times HIW35\mathrm{km}\mathrm{sec}^1$$
(1)
calibrated by Munari & Tomasella (1999) on the width at half maximum of He I 5876 Å in high resolution spectra of O and B stars. The decovolution of the He I 5876 Å line profile in the two TU Mus quadrature spectra gives a HIW of 7.26 and 7.35 Å for the primary, and 6.37 and 6.36 Å for the secondary, which corresponds to 273 and 277, and 234 and 236 km sec-1 respectively. The Munari and Tomasella relation was calibrated using single, and thus axially symmetric, stars. The components of TU Mus depart quite strongly from axial symmetry so we investigated its influence on the HIW using the line profile feature of WD. We computed the HIW for TU Mus and for a spherical star with radius equal to the $`r_{back}`$ radii of the TU Mus components. We found that the HIW for TU Mus was about 5% smaller than for the spherical star. Applying this correction to our data, we adopt 290 and 248 km sec-1 as the rotational velocities of primary and secondary. Of course, the concept of a rotational velocity is of less obvious usefulness for stars like TU Mus than for axially symmetric ones. The question arises as to where on the equator of TU Mus the above rotational velocities apply. Since our spectra were taken at quadrature phases, $`r_{back}`$ seems appropriate and we find values of 291 and 242 km sec-1 there from our light curve solution (which assumed synchronism), in excellent agreement with the line profile results. Our data support the AG75 values and rule out the smaller values measured by S95.

To confirm the MT method we have applied it to the He I 5876 Å line of the B1.5Vn star HIP 77635 that we observed under identical instrumental conditions soon after TU Mus on Feb. 22. Its rotational velocity is reported as 306 km sec-1 in the Bernacca & Perinotto (1970) catalogue of rotational velocities. The HIW of the line turned out to be 8.16 Å which corresponds to 311 km sec-1 according to Equation 1. This excellent agreement supports our confidence in the rotational velocities we have measured for the components of TU Mus. We have also computed synthetic, broadened profiles for He I 4471, 6678 and 7065 Å and the agreement with the observed profiles is also excellent.

### 2.2 Distance

Assuming a luminosity ratio $`L_2/L_1`$=0.55 from our light curve solution, the $`V_T`$=8.41 magnitude transforms into Johnson V magnitudes of $`V_1`$=8.88 and $`V_2`$=9.52, following Bessell (2000). With $`E_{BV}`$=0.19 (see below) and assuming an O7.5V classification for the primary, the distance to TU Mus amounts to 4.77 kpc using the absolute magnitudes calibrated by Houk, Branch & Winfrey (1994) ($`M_V`$=–5.1, –4.5 for O7.5V, O9.5V stars) and a ratio of total to selective absorption of $`R_V=3.1`$. The same calculation for the O9.5V secondary results in a distance of 4.86 kpc. We therefore adopt 4.8 kpc as the distance to TU Mus.

### 2.3 Reddening

The interstellar lines of Na I offer a means of estimating the reddening affecting TU Mus independent of the knowledge of the intrinsic colors. The profile of Na I D1 and D2 lines on our FEROS spectra of TU Mus are shown in Figure 1. They are clearly composed of three separate components. Their average position, width and equivalent width have been derived from multi-Gaussian fitting on all 9 available spectra, and the results are summarized in Table 2, where the errors are the standard errors. The 9 individual spectra gave extremely consistent results. The multi-Gaussian mean fit of Table 2 is overplotted on the D2 line of spectrum #1629 in Figure 1 to show the accuracy of the fit.

The equivalent widths of the three components have been transformed in Table 2 into the corresponding color excess using the relation calibrated by Munari & Zwitter (1997), with a total color excess amounting to $`E_{BV}`$=0.19$`\pm `$0.01. Similar profiles are observed for K I 7698 Å and H & K Ca II interstellar lines. The Na I widths in Table 2 are not corrected for instrumental width, which amount to 5.0 km sec-1 from unresolved telluric lines in the vicinity of Na I lines (thus intrinsic width does not exceed 3 km sec-1). The line of sight to TU Mus therefore crosses three kinematically very sharp and well separated sources of reddening, of pure interstellar origin. In fact none of the components show radial velocity change with orbital phase and none share the -4 km sec-1 systemic velocity.

Of interest is the heliocentric radial velocity of the three components. The accuracy reported in Table 2 is not fictitious because nearby telluric absorptions show the same wavelength stability from one spectrum to the other, and it matches the expected high FEROS spectrograph stability suitable for extra-solar planet searches. The limited resolution of the AG75 spectra did not allow them to recognize the multi-component structure of the interstellar lines (H and K lines of Ca II), for which they estimated a $`7.3\pm 0.8`$ km sec-1 radial velocity. Averaging the three components in Table 2 weighted according to equivalent width, a radial velocity of $`0.7\pm 0.2`$ km sec-1 is derived for an unresolved interstellar Na I profile. However, a $`7.2\pm 0.4`$ km sec-1 value averaged by equivalent width is obtained on our spectra for both H & K Ca II as well as K I interstellar lines. While showing exactly the same two major components of Na I D1 & D2, they fail to show the third, reddest one.

The radial velocities of the three components are plotted in Figure 2 where they are compared with the expected values of the mean galactic rotation along the line of sight to TU Mus ($`l`$=284.8, $`b`$=–4.1 deg) from the Hron (1987) formalism with $`A`$=–17 km sec-1 and $`\alpha `$=–2.0 km sec-1 pc-2 and Pont, et al. (1994) values for ($`u_{}`$, $`v_{}`$, $`w_{}`$). It is evident that the mean Galaxy rotation curve does not reproduce the expected velocities up to the $``$4.8 kpc distance to TU Mus (note that the line of sight to the target is still mostly within the $`z=\pm 150`$ pc dust layer around the galactic plane at the estimated $``$4.8 kpc distance). The discrete velocity field maps of Brand & Blitz (1993), based on H II regions and reflection nebulae, do not perform much better, which could be anticipated given the scant and widely spaced data used in their analysis in that part of the Galaxy. Only the stronger component can be accounted for and located at $``$2.5 kpc distance by the Brand & Blitz (1993) map. The limited individual reddening implied by the equivalent widths of the three interstellar lines suggests that only a few limited and kinematically scattered bubbles of the interstellar medium are crossed by the line of sight, with a diffuse component playing a marginal role, if any. The kinematic information from the interstellar lines cannot therefore be used to constrain the distance to TU Mus.

### 2.4 Radial velocities

The radial velocity of TU Mus has been measured on the two groups of spectra in Table 1 obtained close to quadrature phases. Those around orbital phase 0.03 are very close to conjunction and thus have insufficient velocity separation between the two components to obtain reliable radial velocities. Line splitting at quadrature is wide and allows an easy and firm determination of the radial velocity, as shown for H$`\beta `$ in Figure 3. The radial velocities have been obtained by Gaussian de-convolution of the profiles into two components, which provide an excellent overall fit to the observed line profile. The measured lines are He I 4016.218, 4471.507, 4921.929, 5015.675, 5875.651, 6678.149, 7065.276 Å; He II 4199.831, 4541.589, 4685.682, 5411.524 Å; and H$`\gamma `$, H$`\beta `$, H$`\alpha `$. The mean of the radial velocities (and associated standard error) derived from these 14 lines is in Table 3.

Our high dispersion, high S/N radial velocities are very close to the photographic ones by AG75 and definitively rule out the smaller ones obtained by S95 on low S/N IUE SWP high resolution spectra. This is clearly demonstrated by Figure 4 that shows the He I 5876 Å line of spectrum # 1629 superimposed with the two component fit using our radial velocities in Table 3 and the corresponding values for AG75 and S95.

### 2.5 Spectral classification and metallicity

Both stars are O-type, showing the whole series of expected He II lines. Their spectral types are somewhat different, as the comparison of H$`\beta `$ and He I 4471 Å profiles in Figures 3 and 4 suggests. The ratio of the two H$`\beta `$ lines is 1.8 while those of the He I is 1.3, favoring a two spectral subtype higher temperature for the primary. A reasonable estimate seems to be O7.5 for the primary and O9.5 for the secondary.

Other than the hydrogen and helium lines, only two C III lines are clearly visible at 4650 Å and 4069 Å(̇Si IV lines are unsuitable in this regard.) Both C III lines are actually blends of three individual lines with complicated profiles caused by wide wavelength separation within the triplets, fast rotation, and the SB2 nature of the spectra. However, a comparison with synthetic spectra does not support a carbon content, and therefore metallicity, far from the solar value.

## 3 Photometry

TU Mus was observed by RHN on 10 nights from 2002 April 13 to June 4 at the Mount John University Observatory while he was a guest of the University of Canterbury in Christchurch, New Zealand. The f/13.5 Cassegrain focus on the 0.61 meter Optical Craftsmen telescope was used together with a Santa Barbara Instrument Group ST-9E CCD camera, utilizing the KAF-0261E chip (512x512) with 20 micron pixels. A telecompressor lens operated at 2.11x compression reduced the focal ratio to 6.4 and gave a field of view of 9’ x 9’. The B and V (Johnson) filters were by Schuler.

Reduction of the CCD images and aperture photometry was done with the MIRA software package.111Available at http://www.axres.com Bias and dark removal plus flat fielding (dome flats) were done in the usual way. The comparison star was GSC 8984:0579 and the check stars were GSC 8984:1473 and GSC 8984:1623. Weather and time constraints did not permit the determination of transformations to the standard system. There were 1973 frames taken in B and 1962 in V. However, wispy clouds often caused scattered measurements. In order to remove errant measurements in an objective way, the raw comparison magnitudes were plotted versus time. Points adjacent to rapidly changing values (greather than 0.2 magnitudes per sampling interval of 30 seconds) were rejected. In this way, the data set was reduced to 1148 points in B and 1493 points in V, and the light curves were cleaned up considerably.

## 4 Data Analysis

We performed a simultaneous analysis of the AG75 photometry and radial velocities, the Hipparcos photometry, and our new photometry and velocities with the 2003 version of the Wilson-Devinney (WD) code. This version of WD has the capability of modeling the radiation of the stars by Kurucz (1993) atmospheres and we used this feature in our solutions.

In early exploratory fits with the light curve program of WD, it became clear that the AG75 photometry, the Hipparcos photometry, and our new photometry were not consistent in terms of the depths of the eclipses. The AG75 and Hipparcos photometry matched well, but our new photometry had slightly deeper eclipses, requiring an inclination about one degree higher. Rather than indicating some real change in the inclination of the system, the differences result from the fact that the older photometry suffers from third light contamination. The Digitized Sky Survey image of the field shows a companion about $`15\mathrm{"}`$ southeast of TU Mus. Our aperture photometry was done with a $`10\mathrm{"}`$ aperture and therefore did not include the companion. The effective aperture for the Hipparcos observations was about $`30\mathrm{"}`$ (ESA, 1997), thus including the companion. The AG75 data also included the companion since their aperture was $`30\mathrm{"}`$. The Guide Star Catalog lists the brightness of the companion as $`13_.^m0\pm 0_.^m4`$. In our light curve solutions below, we allowed for third light in all of the photometry but the values for the $`BV`$ curves were always very tiny, thus we conclude that third light does not affect our $`BV`$ data. The solution did find third light in the AG75 and Hipparcos data and the results are consistent with the GSC brightness of the companion. Figure 5 shows the DSS image of the area around TU Mus.

Previously, WD had been modified to allow for the use of either binary phase or time as the independent variable (Wilson & Terrell, 1998). In the latter case, one can solve for the orbital period ($`P`$), its first time derivative ($`\dot{P}`$), and the reference epoch ($`HJD_0`$, usually the time of primary minimum). In order to investigate the behavior of the period of TU Mus, we used time as the independent variable for our simultaneous solutions, adjusting the period and reference epoch. Attempts to adjust Ṗ always resulted in values statistically indistinguishable from zero, indicating that the period of TU Mus has remained constant over the thirty year span of the observations being analyzed. Our estimate of the ephemeris is

$$HJD_{min}=2448500.3080(1)+1.38728653(2)\times E$$
(2)
for primary eclipse where the quantities in parentheses are the $`1\sigma `$ uncertainties in the last digits of the parameters.

Other parameters adjusted in the simultaneous solution were the semi-major axis of the relative orbit ($`a`$), the binary center of mass radial velocity ($`V_\gamma `$), orbital inclination ($`i`$), secondary mean effective temperature ($`T_2`$), common envelope surface potential ($`\mathrm{\Omega }`$), mass ratio ($`q`$), and the bandpass-specific luminosity of the primary ($`L_1`$). Other parameters, such as the bolometric albedos and gravity brightening exponents, were held fixed at their expected theoretical values. The square root limb darkening law was used with coefficients from Van Hamme (1993). The mean effective temperature of the primary was set to 35,000 K based on the O7.5 spectral type. Data set weights were determined by the scatter of the observations. WD gives output at each step from which the standard deviations of the data sets, $`\sigma `$, can be calculated. It is sufficient to determine reasonably good starting values for the standard deviations by examining plots of the observations and then adjust them as the solution progresses. Within a given light curve or radial velocity curve, WD can also apply individual weights. For the radial velocities, we set the base curve weights for the AG75 velocities and then gave our new velocities higher individual weights based on their estimated errors. In order to eliminate potential local minimum problems, we began the solutions from several starting points and found that the same minimum was always recovered.

Table 4 shows the results of the simultaneous solution. TU Mus is just barely in an overcontact configuration, having a surface potential of 3.137 $`\pm `$ 0.002 while the inner critical potential is 3.156 for the estimated mass ratio of 0.651. Using the degree of contact index defined by Wilson & Rafert (1981), where point contact has a value of 1.0, we find a value of 1.006 for TU Mus. Figure 6 shows the fit to the AG75 radial velocities and the ones we obtained. The fit to the Hipparcos photometry is shown in Figure 7 and since all of the various observations were phased with the ephemeris given in Equation 2, it shows that the period has remained constant to the precision with which we can measure it. Figure 8 shows the fit to the AG75 $`u`$ observations and exhibits a small but noticeable asymmetry in the maxima, a phenomenon frequently seen in the light curves of late-type overcontact (W UMa) systems. As can be seen in Figure 9, the asymmetry is slightly smaller in the $`y`$ data, indicating that the source of the asymmetry is a high temperature phenomenon, perhaps a superluminous area as suggested in CE Leo (catalog ) by Samec, et al. (1993) and YZ Phe (catalog ) by Samec & Terrell (1995). Figure 10 shows the fits to our $`B`$ and $`V`$ data. An asymmetry is present in these newer data but in the opposite sense from the AG75 data, indicating that the source of the asymmetry is variable in time.

## 5 Age and Evolutionary Status

Very little is known about the structure and evolution of overcontact binaries with radiative envelopes. To our knowledge, no models of this type have been published. Although obviously limited, a comparison with single-star models can be instructive in estimating the age and evolutionary status of TU Mus. Using a stellar evolution code (Han, et al. (1994), Eggleton (1973), Eggleton (1972), Eggleton (1971)) kindly supplied to us by P. Eggleton, we have constructed evolutionary models of $`23.5M_{\mathrm{}}`$ and $`15.3M_{\mathrm{}}`$. The model for the primary reaches the observed radius of the TU Mus primary at about 2.6 Myr and the model for the secondary reaches the observed radius at about 6.2 Myr. Because of the luminosity exchange from the primary to the secondary, the radius of the primary is probably less than it would be without the luminosity transfer. Similarly, the secondary is probably larger than it otherwise would be. Thus the true age of TU Mus probably lies between these two values. Our ignorance of the structure and evolution of these stars is so great, however, that one cannot discount the possibility that the current configuration is a result of a much more complicated evolutionary history.

Two results from this work, the overcontact configuration and the lack of any period change, have interesting implications for the structure of early-type overcontact systems. In the late-type overcontact (W UMa) systems, the luminosity transfer drives mass transfer from the less massive to the more massive star, resulting in period changes. But if TU Mus truly is in an overcontact configuration, the same level of mass transfer cannot be taking place because of the constant period. What is clearly needed is more work in this area with modern codes like Djehuty (Eggleton, et al., 2003). More sophisticated modeling by 3D stellar evolution codes will certainly give us a better understanding of the structure and evolution of early-type overcontact systems and the well-determined properties of TU Mus makes it a good target of such modeling.

## 6 Conclusions

Our analysis of new and existing data on TU Mus shows it to be an overcontact binary consisting of a $`23.5M_{\mathrm{}}`$ primary and a $`15.3M_{\mathrm{}}`$ secondary with a small degree of contact. Our results are in very good agreement with the results of AG75 and contradict the smaller masses found by S95. The source of the discrepancy between the UV and optical velocities is not clear. One idea we explored was the effect of the change in the center of light due to the reflection effect between the UV and optical, but our experiments indicated that this effect is much too small to explain the large discrepancy. Given our agreement with the AG75 velocities, one is left with two possible conclusions: that the UV velocities suffer from some unrecognized systematic error or there is some fundamental effect that causes a difference between UV and optical velocities for stars of this type. We note that Stickland and colleagues have seen a similar discrepancy for LY Aur (Stickland, et al., 1994). Our current understanding of the structure and evolution of early-type overcontact binaries is quite limited and we hope that our results for the observed properties of TU Mus may serve as a good test of future 3D models of these stars.

RHN thanks the American Association of Variable Star Observers for the loan of the ST-9E camera. RHN would also like to thank specifically the two technicians at the observatory, Alan Gilmore and Pam Kilmartin, for their superb help, and the members of the Physics and Astronomy Department at the University of Canterbury for their warm hospitality. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.


text:
# Non-monotonic angular magnetoresistance in asymmetric spin valves

## Abstract

The electric resistance of ferromagnet/normal-metal/ferromagnet perpendicular spin valves depends on the relative angle between the magnetization directions. In contrast to common wisdom, this angular magnetoresistance is found to be not necessarily a monotone function of the angle. The parameter dependence of the global resistance minimum at finite angles is studied and the conditions for experimental observation are specified.

The discovery of the giant magnetoresistance (GMR) Baibich et al. (1988) has driven much of the current research to enrich the functionalities of electronic circuits and devices employing ferromagnetic elements. The current perpendicular to plane (CPP) transport techniqueBass and Pratt Jr. (1999); Sun et al. (2002); Yang et al. (1995) turned out to be especially suited to study the physics of spin dependent transport. Nanostructured and perpendicular spin valves are ideal devices to investigate the current-induced magnetization reversal,Slonczewski (1996) which has potential applications for magnetic random access memories. These structures allow the measurement of the angular magnetoresistance (aMR)Dauguet et al. (1996); Giacomo et al. (2003) introducing an analogue degree of freedom between the conventional parallel vs. antiparallel digital configurations. A semiclassical theoretical treatment of the aMR leads to the concept of a spin-mixing conductanceBrataas et al. (2000) that turned out useful for phenomena like the spin torqueXia et al. (2002); Bauer et al. (2003) and interface-enhanced Gilbert damping. Tserkovnyak et al. (2002)

This Rapid Communication addresses the aMR of asymmetric perpendicular spin valves. We show that the parallel magnetization configuration of ferromagnet(F)/normal metal(N)/ferromagnet heterostructures does not necessarily correspond to the minimal resistance. This non-monotonic behavior requires a redefinition of the GMR ratio in terms of the global maximum and minimum resistances instead of those for parallel and antiparallel configurations. We discuss how to optimize the conditions for an experimental observation and demonstrate that the spin torque is strongly affected by the asymmetry as well.

First, we summarize necessary concepts from Ref. Brataas et al., 2001 for resistive elements such as an interface between a monodomain ferromagnet with magnetization parallel to the unit vector $`𝒎`$. The charge and spin current, $`I_\mathrm{c}`$ and $`𝑰_\mathrm{s}`$, driven by a potential and spin accumulation bias, $`\mathrm{\Delta }\mu _\mathrm{c}`$ and $`\mathrm{\Delta }𝝁_\mathrm{s}`$, read

$$I_\mathrm{c}=\frac{e}{h}\left[(g^{}+g^{})\mathrm{\Delta }\mu _\mathrm{c}+(g^{}g^{})𝒎\mathrm{\Delta }𝝁_\mathrm{s}\right],$$
(1)
$$\begin{array}{c}𝑰_\mathrm{s}=\frac{1}{4\pi }𝒎\left[(g^{}g^{})\mathrm{\Delta }\mu _\mathrm{c}+(g^{}+g^{})𝒎\mathrm{\Delta }𝝁_\mathrm{s}\right]\hfill \\ +\frac{1}{4\pi }2\mathrm{R}\mathrm{e}(g^{})𝒎\times (\mathrm{\Delta }𝝁_\mathrm{s}\times 𝒎).\hfill \end{array}$$
(2)
where $`g^{}`$ and $`g^{}`$ are the conductances for electrons with majority and minority spin, respectively, and $`g^{}`$ is the mixing conductance for a spin current polarized transverse to the magnetization. We disregarded $`\mathrm{Im}(g^{})`$, which for metallic interfaces is usually smaller than 10% of $`\mathrm{Re}(g^{})`$.Xia et al. (2002); Stiles and Zangwill (2002) It is convenient to introduce $`g=g^{}+g^{}`$, $`p=(g^{}g^{})/g`$ and $`\eta =2g^{}/g`$, where $`g`$ is the total conductance, $`p`$ the polarization and $`\eta `$ the relative mixing conductance.

Let us examine a two terminal system (F-N-F) as shown in Figure 1. The contacts need not be identical; the conduction parameters are summarized as $`G_\mathrm{L}`$ and $`G_\mathrm{R}`$. The electric resistance as function of the angle between the magnetization directions of the reservoirs, can simply be calculated using Eqs. (1) and (2), assuming charge and spin conservation on the normal metal node. For a symmetric structure ($`G_\mathrm{L}=G_\mathrm{R}=G`$) the resistance $`R(\theta )`$ reads:Brataas et al. (2001)

$$\frac{e^2}{h}R(\theta )=\frac{2}{g}\frac{\mathrm{tan}^2\theta /2+\eta }{(1p^2)\mathrm{tan}^2\theta /2+\eta }.$$
(3)
If necessary, spin flip processes in the normal metal can be included.Brataas et al. (2001) A finite angle between the magnetizations causes a spin accumulation on the normal metal node. Since we disregard the imaginary part of the mixing conductance, it lies in the plane of the magnetization vectors. The resistance increases with increasing spin accumulation, whose creation costs energy, and thus with $`\theta `$. Therefore the resistance is minimal when the magnetizations are parallel and maximal for $`\theta =\pi `$. The mixing conductance can be interpreted as an additional channel for dissipating the spin accumulation on the normal metal node for $`0<\theta <\pi `$; an increasing mixing conductance will therefore reduce the total resistance. This is the mechanism behind deviations of the aMR from a simple $`\mathrm{cos}^2\theta /2`$ behavior, which can be used to determine the mixing conductance from experimental curves.Bauer et al. (2003); Giacomo et al. (2003)

In the following, we focus on an asymmetric configuration with $`G_\mathrm{L}G_\mathrm{R}`$. The asymmetry in conductance ($`g_\mathrm{L}g_\mathrm{R}`$) causes a charge accumulation on the normal metal node. Similarly, when $`p_\mathrm{L}p_\mathrm{R}`$, a spin accumulation is excited on the normal metal node even for $`\theta =0`$. We find here that configurations with $`\theta 0`$ may correspond to a spin accumulation that is smaller than that of the parallel one, and therefore a global resistance minimum at finite angles. The recipe for a significant effect is a large polarization of the current by the source contact (e.g. $`p_\mathrm{L}1`$) and efficient dissipation of the spin accumulation for a finite angle $`\theta `$ by a large mixing conductance ($`\eta _\mathrm{R}>1`$) of the drain. The polarization direction of the spin current differs from the magnetization directions for finite angles $`\theta `$ as in the symmetric case, but the noted asymmetry forces it to be close to the magnetization direction of the source contact. A large mixing conductance $`g^{}`$ favors the transverse over the longitudinal spin current. Spins on the normal metal node therefore escape easily and the reduced spin accumulation is equivalent to a decrease of the total resistance. This interplay between spin accumulation and magnetization angles strongly modifies the total aMR profile.

Unfortunately, the exact equations for $`R(\theta )`$ are not very transparent. A perturbation approach to these equations is not helpful because no small parameters can be identified for the experimentally relevant metallic structures. However, we did find relatively simple analytical expressions for the angle $`\theta _\mathrm{m}`$ of the global resistance minimum as well as a simple expression for the maximal aMR in the limit $`\eta _\mathrm{R}1`$.

We derive that next to $`\theta =0,\pi `$, the resistance may have extrema at two additional angles:

$`\mathrm{cos}\theta _{\mathrm{m1}}=\left({\displaystyle \frac{p_\mathrm{R}}{p_\mathrm{L}}}\right)\left({\displaystyle \frac{1+\left(\frac{g_\mathrm{L}}{g_\mathrm{R}}\right)\frac{1p_\mathrm{L}^2}{\eta _R}}{1\frac{1p_\mathrm{R}^2}{\eta _R}}}\right),`$ (4)
$`\mathrm{cos}\theta _{\mathrm{m2}}=\left({\displaystyle \frac{p_\mathrm{L}}{p_\mathrm{R}}}\right)\left({\displaystyle \frac{1+\left(\frac{g_\mathrm{R}}{g_\mathrm{L}}\right)\frac{1p_\mathrm{R}^2}{\eta _L}}{1\frac{1p_\mathrm{L}^2}{\eta _L}}}\right).`$ (5)
where the absolute value of $`\mathrm{cos}\theta _\mathrm{m}`$ must be smaller than unity, which is clearly not the case for a symmetrical spin valve. The condition for one extra extremum is easily fulfilled. Two  additional extrema are not consistent with the condition $`\eta _\mathrm{L},\eta _\mathrm{R}>1`$, which rigorously holds for high contact resistances,Brataas et al. (2001) but not necessarily for highly transparent interfaces.Bauer et al. (2003) It can be proven that when one extremum exists and $`\eta _\mathrm{L},\eta _\mathrm{R}>1`$, the extremum is the global minimum and located in the interval $`0<\theta <\pi /2`$. When $`\eta _\mathrm{R}<1p_\mathrm{R}^2`$ (that does not seem very likely for metals), an additional extremum may exist in the interval $`\pi /2<\theta <\pi `$. It turns out to be a maximum that can be understood in the same way as the minimum. Additional minima and maxima may even coexist for specific parameter combinations, which do not appear relevant for metallic spin valves, however.

The position of the global minimum does not depend on $`\eta `$ of the source contact. The source polarizes the current through the total structure parallel to its magnetization, therefore the source $`\eta `$ does not play a role at all. The component of the spin current orthogonal to the magnetization is called the spin torqueXia et al. (2002) acting on this magnetization, since it is absorbed by the magnetic order parameter and may excite the magnetization when exceeding a threshold value.Slonczewski (1996); Oezyilmaz et al. (2003); Kiselev et al. (2003) In the global minimum the spin torque on the source magnetization vanishes with the transverse component of the spin current. The spin torque on the drain is large, but not at the maximum as a function of $`\theta `$.

Let us choose the left lead to be the polarizing source $`(p_\mathrm{L}>p_\mathrm{R})`$ and the right lead to be the dissipating drain $`(\eta _\mathrm{L}>\eta _\mathrm{R})`$. The condition for a non-collinear resistance minimum is now:

$$\left|\left(\frac{p_\mathrm{R}}{p_\mathrm{L}}\right)\left(\frac{1+\left(\frac{g_\mathrm{L}}{g_\mathrm{R}}\right)\frac{1p_\mathrm{L}^2}{\eta _R}}{1\frac{1p_\mathrm{R}^2}{\eta _R}}\right)\right|<1.$$
(6)
When the second factor is larger than one (true for $`\eta _\mathrm{R}>1`$), only a polarization ratio $`p_\mathrm{L}/p_\mathrm{R}>1`$ can save this inequality. The condition is never fulfilled when the left hand side of Equation (6) diverges:

$$\frac{1p_\mathrm{R}^2}{\eta _\mathrm{R}}=\frac{2\left(\frac{1}{g_\mathrm{R}^{}}+\frac{1}{g_\mathrm{R}^{}}\right)^1}{g_\mathrm{R}^{}}1.$$
(7)
Therefore $`g_\mathrm{R}^{}`$ should be considerably larger than $`g_\mathrm{R}^{}`$ and $`g_\mathrm{R}^{}`$. When the average conductance of the source is smaller than the mixing conductance of the drain, the numerator

$$\left(\frac{g_\mathrm{L}}{g_\mathrm{R}}\right)\frac{1p_\mathrm{L}^2}{\eta _R}=\frac{2\left(\frac{1}{g_\mathrm{L}^{}}+\frac{1}{g_\mathrm{L}^{}}\right)^1}{g_\mathrm{R}^{}},$$
(8)
reduces $`\mathrm{cos}\theta _\mathrm{m}`$, and hence increases $`\theta _\mathrm{m}`$.

The GMR ratio is usually defined in terms of the resistance in parallel or antiparallel configurations, in terms of the previously introduced parameters:

$`GMR^{}={\displaystyle \frac{R^{\mathrm{ap}}R^\mathrm{p}}{R^{\mathrm{ap}}}}=`$ (9)
$`{\displaystyle \frac{g_\mathrm{L}g_\mathrm{R}\left[(p_\mathrm{L}+p_\mathrm{R})^2(p_\mathrm{L}p_\mathrm{R})^2\right]}{(g_\mathrm{L}+g_\mathrm{R})^2(g_\mathrm{L}p_\mathrm{L}g_\mathrm{R}p_\mathrm{R})^2}},`$
where $`R^{\mathrm{ap}}`$ and $`R^\mathrm{p}`$ are the resistances for antiparallel and parallel configuration, respectively. Hence, the GMR ratio increases when the total polarization increases, as expected. However, when the difference between the polarizations of both sides is large, $`GMR^{}`$ decreases, because of a larger spin accumulation on the normal metal node for $`\theta =0`$, as noted above. Since now $`R^{\mathrm{ap}}R^\mathrm{p}`$ is no longer the maximal resistance difference, a new definition for the magnetoresistance is appropriate in terms of the global maximum (for which we still take the antiparallel configuration) and the newly found global minimum. In the limit of large $`\eta _\mathrm{R}1`$ we arrive at the simple result

$$GMR=\frac{R^{\mathrm{ap}}R^\mathrm{m}}{R^{\mathrm{ap}}}=GMR^{}\frac{(p_\mathrm{L}+p_\mathrm{R})^2}{4p_\mathrm{L}p_\mathrm{R}},$$
(10)
where $`R^\mathrm{m}`$ is the global minimum of the resistance. It can easily be verified that $`GMR`$ is indeed larger than $`GMR^{}`$.

Next we investigate the conditions under which this enhanced magnetoresistance can be measured in magnetic spin valves with a current perpendicular to plane geometry.Bass and Pratt Jr. (1999) Even for identical magnetic layers an asymmetry can be realized by a spin independent resistance or tunnel barrier at the outside of one of the magnetic films, as long as the spin diffusion length is larger than the total bilayer. Such an additional series resistor then effectively decreases $`p`$, $`g^{}`$ and $`g^{}`$ of this magnet. Because the spin current normal to the magnetization is absorbed by the magnetSlonczewski (1996) over a couple of monolayers, $`g^{}`$ is not modified by the extra resistance.Xia et al. (2002) $`g^{}`$ thus can indeed be engineered to be larger than $`g^{}`$ for a given contact such that a non-monotonic aMR can be expected.

Spin dependent bulk resistances contribute to the aMR over thicknesses smaller than the spin diffusion length. Copper and cobalt have relatively large spin flip lengths, respectively $`250\mathrm{nm}`$ and $`50\mathrm{nm}`$, which makes them useful materials to explore this effect. $`\mathrm{Al}_2\mathrm{O}_3`$ tunnel barriers are routinely used for tunnel MR studies and suitable materials for the present purposes.

The full aMR profile can best be calculated numerically. We consider a structure consisting of two identical cobalt layers (thickness is $`3\mathrm{nm}`$) separated by a thin copper layer ($`10\mathrm{nm}`$). Both magnets are sandwiched by tunnel junctions, another copper layer, and finally normal metal reservoirs as sketched in Fig. 2. Bulk resistances of copper and cobalt are disregarded because they are relatively very small. $`G_\mathrm{F}`$ symbolizes all conductance parameters of a copper-cobalt interface, $`G_\mathrm{L}`$ and $`G_\mathrm{R}`$ stand for the outer normal resistances including the tunnel junctions. For interfaces between a ferromagnet and a normal metal, $`g^{}`$ ( $``$ number of modes in the normal metal) usually lies between $`g^{}`$ and $`g^{}`$. For a Co/Cu interface $`g`$ is typically $`1413\times 10^3`$ (for an interface cross section of $`140\times 90\mathrm{nm}^2`$), $`p=0.75`$ (Ref. Yang et al., 1995) and $`\eta =0.38`$ (Ref. Xia et al., 2002); these values include the Boltzmann corrections for transparent interfaces.Bauer et al. (2003) In order to compare configurations with different values for $`G_\mathrm{L}`$ and $`G_\mathrm{R}`$, its series resistance is assumed constant at $`1/G_\mathrm{L}+1/G_\mathrm{R}=0.37\mathrm{\Omega }`$.

The computed aMR is presented in Figure 3 for three different ratios $`G_\mathrm{L}:G_\mathrm{R}`$. When $`G_\mathrm{L}`$ differs sufficiently from $`G_\mathrm{R}`$, the global minimum shifts away from the parallel configuration, as predicted. The position of the global minimum, $`\theta _\mathrm{m}`$ increases with increasing polarization contrast. The GMR ratio increases as well, which is in qualitative agreement with Eq. (10).

Finally we compute the spin torque, $`i.e.`$ the transverse component of the spin current. The spin torque in spin valves is governed by similar expression as the charge currentBauer et al. (2003) and is strongly affected by the asymmetry as well. It is convenient to normalize the spin torque by the charge current:

$$i_\mathrm{s}=\frac{|𝒎\times (𝑰_\mathrm{s}\times 𝒎)|}{|I_\mathrm{c}|}.$$
(11)
In Figure 4, $`i_\mathrm{s}`$ of the left magnetization is plotted as a function of $`\theta `$ and different $`G_\mathrm{L}/G_\mathrm{R}`$ ratios. The zeroes in the intervals $`0<\theta <\pi /2`$ illustrate that when the left side is the polarizing source, the spin torque at the global minimum vanishes, which agrees with the finding that the resistance minimum is not a function of $`\eta _\mathrm{L}`$. We observe that the spin torque is strongly enhanced when $`G_\mathrm{L}/G_\mathrm{R}0`$ because the relative mixing conductance of the left hand side is then highly increased, which physically means that the spin accumulation on the normal metal node easily can be dissipated.

To summarize, we have shown that the angular magnetoresistance (aMR) of perpendicular spin valves can be a non-monotonic function when the contacts between the central normal metal node and the outer ferromagnets differ. An analytical expression is derived for the angle $`\theta _\mathrm{m}`$ at which the magnetoresistance has its global minimum. A new definition for the GMR ratio is proposed to take this effect into account. This GMR ratio is now larger than the conventional definition in terms of the resistance of parallel and antiparallel configurations. The spin torque in asymmetric structures is also importantly modified.

J.M. thanks Martin Grønsleth and Jan Petter Morten for fruitful discussions. This work has been supported by ERASMUS, the FOM and the NEDO joint research program “Nano-Scale Magnetoelectronics”.