text:
# Nonlinear Wave Equations

## 1. Introduction

Consider the constant and variable coefficient wave operators in $`\times ^n`$,

$$\mathrm{}=_t^2\mathrm{\Delta }_x,\mathrm{}_g=g^{ij}(t,x)_i_j.$$

In the variable coefficient case the summation occurs from $`0`$ to $`n`$ where the index $`0`$ stands for the time variable. To insure that the equation is hyperbolic in time we assume that the matrix $`g^{ij}`$ has signature $`(1,n)`$ and that the time level sets $`t=const`$ are space-like, i.e. $`g^{00}>0`$. We consider semilinear wave equations,

$$\mathrm{}u=N(u)(SLW),\mathrm{}u=N(u,u)(GSLW)$$

and quasilinear wave equations,

$$\mathrm{}_{g(u)}u=N(u)(u)^2(NLW),\mathrm{}_{g(u,u)}u=N(u,u)(GNLW).$$

To each of these equations we associate initial data in Sobolev spaces

$$u(0)=u_0H^s(^n),_tu(0)=u_1H^{s1}(^n).$$

There are two natural questions to ask: (i) Are the equations locally well-posed in $`H^s\times H^{s1}`$? (ii) Are the solutions global, or is there blow-up in finite time?

### Local well-posedness.

In a first approximation we define it as follows:

###### Definition 1.

A nonlinear wave equation is well-posed in $`H^s\times H^{s1}`$ if for each $`(v_0,v_1)H^s\times H^{s1}`$ there is $`T>0`$ and a neighborhood $`V`$ of $`(v_0,v_1)`$ in $`H^s\times H^{s1}`$ so that for each initial data $`(u_0,u_1)V`$ there is an unique solution $`uC(T,T;H^s)`$, $`_tuC(T,T;H^{s1})`$ which depends continuously on the initial data.

In practice in order to prove uniqueness one often has to further restrict the class of admissible solutions. In most problems, the bound $`T`$ from below for the life-span of the solutions can be chosen to depend only on the size of the data.

It is not very difficult to prove that all of the above problems are locally well-posed in $`H^s\times H^{s1}`$ for large $`s`$. The interesting question is what happens when $`s`$ is small. One indication in this regard is given by scaling. At least in the case when the nonlinear term has some homogeneity, for instance $`N(u)=u^p`$ or $`N(u)=u^p(u)^q`$, one looks for an index $`\alpha `$ so that all transformations of the form $`u(x,t)\lambda ^\alpha u(\lambda x,\lambda t)`$, $`\lambda >0`$ leave the equation unchanged. Correspondingly one finds an index $`s_0=\frac{n}{2}\alpha `$ so that the norm of the initial data $`(u_0,u_1)`$ in the homogeneous Sobolev spaces $`\dot{H}^s\times \dot{H^{s1}}`$ is preserved by the above transformations.

Below scaling ($`s<s_0`$) a small data small time result rescales into a large data large time result. Heuristically one concludes that local well-posedness should not hold. Still, to the author’s knowledge there is no proof of this yet.

###### Conjecture 2.

Semilinear wave equations are ill-posed below scaling.

This becomes much easier to prove if one strengthens the definition of well-posedness, e.g. by asking for uniformly continuous or $`C^1`$ dependence of the solution on the initial data.

If $`s=s_0`$ then for small initial data local well-posedness is equivalent to global well-posedness. The same would happen for large data if we were to strengthen the definition of well-posedness and ask for a lifespan bound which depends only on the size of the data. This is the only case where this distinction makes a difference.

If $`s>s_0`$ then a local well-posedness result gives bounds for life-span $`T_{max}`$ of the solutions in terms of the size of the data,

$$(u_0,u_1)_{H^s\times H^{s1}}MT_{max}M^{s_0s}.$$

The better localization in time makes the problems somewhat easier to study. However, besides scaling there are also other obstructions to well-posedness. These are related to various concentration phenomena which can occur depending on the precise structure of the equation.

### Global well-posedness.

We briefly mention that there is a special case in which the global well-posedness is well understood, namely when the initial data is small, smooth, and decays at infinity. This is not discussed at all in what follows.

Consider first the case when $`s`$ is above scaling, $`s>s_0`$, and local well-posedness holds in $`H^s\times H^{s1}`$. Then any solution can be continued as long as its size does not blow-up. Hence the goal of any global argument should be to establish a-priori bounds on the $`H^s\times H^{s1}`$ norm of the solution. All known results of this type are for problems for which there are either conserved or quasi-conserved positive definite quantities. Such conserved quantities can often be found for equations which are physically motivated or which have some variational structure. For simplicity suppose that there is some index $`s_c`$ and an energy functional $`E`$ in $`H^{s_c}\times H^{s_c1}`$ which is preserved along the flow. The index $`s_c`$ needs not be equal to the scaling index $`s_0`$. There are three cases to consider:

(i) The subcritical case $`s_c>s_0`$. Then a local well-posedness result at $`s=s_c`$ implies the global result for $`ss_c`$. Furthermore, in recent years there has been considerable interest in establishing global well-posedness also for $`s_0<s<s_c`$. This is based on an idea first introduced by Bourgain  in a related problem for the Schröedinger equation, and followed up by a number of authors.

(ii) The critical case $`s_c=s_0`$. Here the energy is not needed for small data, when local and global well-posedness are equivalent. For large data, however, the energy conservation is not sufficient in order to establish the existence of global solutions. In addition, one needs a non-concentration argument, which should say that the energy cannot concentrate inside a characteristic cone.

(iii) The supercritical case, $`s_c<s_0`$. No global results are known:

###### Open Problem 3.

Are supercritical problems globally well-posed for $`ss_0`$?

A simple example is the equation (NLW) with $`N(u)=|u|^{p1}u`$. The energy is

$$E(u)=|u_t|^2+|_xu|^2+\frac{1}{p+1}|u|^{p+1}.$$

Then $`s_c=1`$, while $`s_0=\frac{n}{2}\frac{2}{p1}`$. In $`3+1`$ dimensions, for instance, $`p=3`$ is subcritical, therefore one has global well-posedness in $`H^1\times L^2`$. The exponent $`p=5`$ is critical and in this case the problem is known to be globally well-posed in $`H^1\times L^2`$; the non-concentration argument is due to Grillakis . The exponent $`p=7`$ is supercritical.

### Blow-up.

Not all nonlinear wave equations are expected to have global solutions. Quite the contrary, generic equations are expected to blow up in finite time; only for problems with some special structure it seems plausible that global well-posedness may hold. A simple way to produce blow-up is to look for self-similar solutions, $`u(x,t)=t^\gamma u(\frac{x}{t})`$. If they exist, self-similar solutions disprove global well-posedness. Because they must respect the scaling of the problem, they are not so useful when trying to disprove local well-posedness.

Another way to produce blow-up solutions is the so-called ode blow-up. In the simplest setting this means looking at one dimensional solutions (say $`u(x,t)=u(t)`$) which solve an ode and blow up in finite time. Then one can truncate the initial data spatially and still retain the blow-up because of the finite speed of propagation. This is still not very useful for the local problem.

A better idea is to constructs blow-up solutions which are concentrated essentially along a light ray, see Lindblad , and Alihnac . In this setup the actual blow-up occurs either because of the increase in the amplitude, in the semilinear case, or because of the focusing of the light rays, in the quasilinear case. As it turns out, the counterexamples of this type are often sharp for the local well-posedness problem.

## 2. Semilinear wave equations

Usually, a fixed point argument is used to obtain local results for semilinear equations. We first explain this for the case when $`s=s_0`$. We define the homogeneous and inhomogeneous solution operators, $`S`$ and $`\mathrm{}^1`$ by

$$S(u_0,u_1)=u\{\mathrm{}u=0,u(0)=u_0,_tu(0)=u_1\},$$

$$\mathrm{}^1f=u\{\mathrm{}u=f,u(0)=0,_tu(0)=0\}.$$

Then the equation (NLW) for instance can be recast as

$$u=S(u_0,u_1)+\mathrm{}^1N(u).$$

To solve this using a fixed point argument one needs two Banach spaces $`X`$ and $`Y`$ with the correct scaling and the following mapping properties:

$$S:H^s\times H^{s1}X,\mathrm{}^1:YX,N:XY.$$

The first two are linear, but the last one is nonlinear. The small Lipschitz constant is always easy to obtain provided the initial data is small and that $`N`$ decays faster than linear at $`0`$. The solutions given by the fixed point argument are global.

In the case $`s>s_0`$ the scaling is lost, and with this method one can only hope to get results which are local in time. To localize in time one chooses a smooth compactly supported cutoff function $`\chi `$ which equals $`1`$ near the origin. The fixed point argument is now used for the equation

$$u=\chi S(u_0,u_1)+\chi \mathrm{}^1N(u).$$

A solution to this solves the original equation only in an interval near the origin where $`\chi =1`$. The modified mapping properties are

$$\chi S:H^s\times H^{s1}X,\chi \mathrm{}^1:YX,N:XY.$$

How does one choose the spaces $`X`$, $`Y`$? One approach is to use the energy estimates for the wave equation and set

$$X=\{uL^{\mathrm{}}(H^s),uL^{\mathrm{}}(H^{s1})\},Y=L^1(H^{s1}).$$

The first two mapping properties are trivial. However, if the third holds then we must also have $`N:XL^{\mathrm{}}(H^{s1})`$. The one unit difference in scaling between $`L^1`$ and $`L^{\mathrm{}}`$ implies that this can only work for $`ss_0+1`$.

What is neglected in the above setup is the dispersive properties of the wave equation. Solutions to the linear wave equation cannot stay concentrated for long time intervals. Instead, they will disperse and decay in time (even though the energy is preserved). In harmonic analysis terms, this is related to the restriction theorem (see ) and is a consequence of the nonvanishing curvature of the characteristic set for the wave operator, namely the cone $`\xi _0^2=\xi _1^2+\mathrm{}+\xi _n^2`$. Here $`\xi `$ stands for the Fourier variable. One way of quantifying the dispersive effects is through the Strichartz estimates. They apply both to the homogeneous and the inhomogeneous equation (see  and references therein):

$$S:H^\rho \times H^{\rho 1}L^pL^q,|D|^{1\rho _1\rho }\mathrm{}^1:L^{p_1^{}}L^{q_1^{}}L^pL^q$$

where $`(\rho ,p,q)`$ and $`(\rho _1,p_1,q_1)`$ are subject to

$$\frac{1}{p}+\frac{n}{q}=\frac{n}{2}\rho ,\frac{2}{p}+\frac{n1}{q}\frac{n1}{2},2p,q\mathrm{},(\rho ,p,q)(1,2,\mathrm{}).$$

The worst case in these estimates occurs for certain highly localized approximate solutions to the wave equation, which are called wave packets. A frequency $`\lambda `$ wave packet on the unit time scale is essentially a bump function in a parallelepiped of size $`1\times \lambda ^1\times (\lambda ^{\frac{1}{2}})^{n1}`$ which is obtained from a $`\lambda ^1\times (\lambda ^{\frac{1}{2}})^{n1}`$ parallelepiped at time zero which travels with speed $`1`$ in the normal direction. Because of the uncertainty principle, this is the best possible spatial localization which remains coherent up to time $`1`$. Of course one can rescale and produce wave packets on all time scales.

In low dimension $`n=2,3`$ the Strichartz estimates provide a complete set of results for generic equations of both (NLW) and (GNLW) type. Consider the following two examples, of which the second is wrong but almost right:

$$\mathrm{}u=u^3,n=3,s=s_0=\frac{1}{2},X=L^4,Y=L^{\frac{4}{3}},$$

$$\mathrm{}u=uu,n=3,s_0=\frac{1}{2},s=1X=|D|^1L^{\mathrm{}}L^2L^2L^{\mathrm{}}Y=L^2.$$

For $`n4`$, however, the Strichartz estimates no longer provide all the results. The reason is as follows. The worst nonlinear interaction in both (NLW) and (GNLW) occurs for wave packets which travel in the same direction. One can use the Strichartz estimates to accurately describe the interaction of same frequency wave packets. But in the interaction of two wave packets at different frequencies, the low frequency packet is more spread, and only a small portion of it will interact with the high frequency packet. However, unlike in low dimension, the Strichartz estimates do not provide sharp bounds for this smaller part of a wave packet.

A more robust idea due to Bourgain  and Klainerman-Machedon  is to use the $`X^{s,b}`$ spaces associated to the wave equation very much in the same way the Sobolev spaces are associated to the Laplacian:

$$u_{X^{s,b}}=(1+|\xi |)^s(1+||\xi _0||\xi ^{}||)^b\widehat{u}_{L^2}.$$

Then one chooses $`X=X^{s,\frac{1}{2}}`$ and $`Y=X^{s1,\frac{1}{2}}`$. The Strichartz information is not lost since for $`\rho ,p,q`$ as above we have the dual embeddings

$$X^{\rho ,\frac{1}{2}+}L^pL^q,L^p^{}L^q^{}X^{\rho ,\frac{1}{2}}.$$

Within the framework of the $`X^{s,b}`$ spaces one can prove bilinear estimates which provide a better description of the interaction of high and low frequencies, see  and references therein. The bilinear estimates are obtained as weighted convolution estimates in the Fourier space, by using the above embeddings, or by combining the two methods. Sometimes even this setup does not suffice and has to be modified further, see .

###### Conjecture 4.

The equation $`\mathrm{}u=u^p`$ is locally well-posed in $`H^s\times H^{s1}`$ for $`n4`$, $`0s\frac{1}{2}`$, $`p(\frac{n+1}{4}s)(\frac{n+5}{4}s)`$. (see  for more details)

### The null condition.

A natural question to ask is whether there are equations which behave better than generic ones. This may happen if the worst interaction (between parallel wave packets) does not occur in the nonlinearity. A good example is (GNLW) with a quadratic nonlinearity $`Q(u,u)=q^{ij}_iu_ju`$. The cancellation condition, called null condition, asserts that

$$q^{ij}\xi _i\xi _j=0\text{in the characteristic set}g^{ij}\xi _i\xi _j=0.$$

All such null forms are linear combinations of

$$Q_{ij}(u,v)=_iu_jv_iv_ju,Q_0(u,v)=g^{ij}_iu_jv.$$

###### Open Problem 5.

Study semilinear wave equations corresponding to variable coefficient wave operators for $`n4`$ (generic case) or $`n2`$ (with null condition).

In the constant coefficient case one can easily use the null condition in the context of the $`X^{s,b}`$ spaces. This is done using inequalities of the following form:

$$|q_0(\xi ,\eta )|c(|p(\xi )|+|p(\eta )|+|p(\xi +\eta )|)$$

respectively

$$|q_{ij}(\xi ,\eta )|c|\xi |^{\frac{1}{2}}|\eta |^{\frac{1}{2}}|\xi +\eta |^{\frac{1}{2}}(|p(\xi )|^{\frac{1}{2}}+|p(\eta )|^{\frac{1}{2}}+|p(\xi +\eta )|^{\frac{1}{2}})$$

where by $`p(\xi )`$ we denote the symbol of the constant coefficient wave operator, given by $`p(\xi )=\xi _0^2\xi _1^2\mathrm{}\xi _n^2`$. Combining this with the embeddings above one can lower the $`s`$ in the local theory whenever the null condition is satisfied. Unfortunately, this does not always give optimal results. The problem of obtaining improved $`L^pL^q`$ estimates for null forms has also been explored, see  , but without immediate applications to semilinear wave equations. We limit the following discussion to two of the more interesting models.

### Wave maps.

These are functions from $`^n\times `$ into a complete Riemannian manifold $`(M,g)`$ which are critical points for

$$I(\varphi )=_{^n\times }|_t\varphi |_g^2|_x\varphi |_g^2dxdt.$$

In local coordinates the equation for wave maps has the form

$$\mathrm{}\varphi ^k=\mathrm{\Gamma }_{ij}^k(\varphi )Q_0(\varphi ^i\varphi ^j)$$

where $`\mathrm{\Gamma }_{ij}^k`$ are the Riemann-Christoffel symbols. The energy functional is

$$E(u)=_^n|_t\varphi |_g^2+|_x\varphi |_g^2dx.$$

The scaling index is $`s_0=\frac{n}{2}`$ and $`s_c=1`$. Local well-posedness for $`s>s_c`$ can be obtained using the $`X^{s,b}`$ spaces. For $`s=s_c`$, using some modified $`X^{s,b}`$ spaces, local (and therefore small data global) well-posedness was established first in homogeneous Besov spaces $`B_{2,1}^{\frac{n}{2}}\times B_{2,1}^{\frac{n}{2}1}`$ in Tataru  and then in Sobolev spaces by Tao  (for the sphere, $`n2`$) and other authors (general target manifold, $`n3`$). Large data global well-posedness is false in the supercritical case $`n3`$, where self-similar blowup can occur. This leaves open problems in the critical case $`n=2`$:

###### Conjecture 6.

(i) The two dimensional wave maps equation is globally well-posed for small data in $`H^1\times L^2`$ for any complete target manifold.

(ii) The two dimensional wave maps equation is globally well-posed for large data in $`H^1\times L^2`$ for “good” target manifolds.

### The Yang Mills equations.

Given a compact Lie group $`𝒢`$ whose Lie algebra $`𝐠`$ admits an invariant inner product $`,`$ one considers $`𝐠`$ valued connection $`1`$-forms $`A_jdx^j`$ in $`^n\times `$. The covariant derivatives of $`𝐠`$ valued functions are defined by

$$D_jB=_jB+[A_j,B].$$

The ($`𝐠`$ valued) curvature of the connection $`A`$ is

$$F_{ij}=_iA_j_jA_i+[A_i,A_j].$$

This is invariant with respect to gauge transformations

$$A_jOA_jO^1_jOO^1,O𝒢.$$

A Yang-Mills connection is a critical point for the Yang-Mills functional

$$I(A)=_{^n\times }F_{ij},F^{ij}𝑑x𝑑t$$

where indices are lifted with respect to the Minkovski metric. Then the Yang-Mills equations have the form

$$D^jF_{ij}=0$$

and the energy functional is

$$E(A)=_{^n\times }F_{ij},F_{ij}𝑑x𝑑t.$$

A Yang-Mills connection is not a single connection, but instead it is a class of equivalence with respect to the above gauge transformation. In order to view the Yang-Mills equations as semilinear wave equations and solve them one has to fix the gauge, i.e. select a single representative out of each equivalence class. Common gauge choices include: (i) the temporal gauge $`A_0=0`$, (ii) the wave gauge $`_jA^j=0`$ and (iii) the Coulomb gauge $`_{j=1}^n_jA_j=0`$. To understand the equation better it may help to look first at an oversimplified version, namely

$$\mathrm{}u=(u_x)u+_xp,_xu=0.$$

This exhibits a $`Q_{ij}`$ type null condition. The scaling index is $`s_0=\frac{n2}{2}`$ and $`s_c=1`$. Using the $`X^{s,b}`$ spaces one can improve the local theory somewhat, but certain more subtle modifications of this are needed in order to handle high-low frequency interactions. In  such an approach is used to prove that local well-posedness holds for $`s>s_0`$, $`n4`$.

###### Open Problem 7.

Is the Yang-Mills equation well-posed for $`s>s_0`$, $`n=2,3`$? (Likely not for $`n=2`$. For $`n=3`$ one can obtain $`s>\frac{3}{4}`$ using the $`X^{s,b}`$ spaces.)

###### Conjecture 8.

(i) The Yang-Mills equation is globally well-posed for small data in $`H^{s_0}\times H^{s_01}`$ for $`n4`$.

(ii) The Yang-Mills equation is globally well-posed for large data in $`H^{s_0}\times H^{s_01}`$ for $`n=4`$.

## 3. Nonlinear wave equations

We consider (NLW), since (GNLW) reduces to it by differentiation. The fixed point argument in the semilinear case cannot be applied in the nonlinear case, because the wave equation parametrix is not strongly stable with respect to small changes in the coefficients. Instead, one must adopt a different strategy: (i) show that local solutions exist for smooth data, (ii) obtain a-priori bounds for smooth solutions uniformly with respect to initial data in a bounded set in $`H^s\times H^{s1}`$ and (iii) prove continuous dependence on the data in a weaker topology, and obtain solutions for $`H^s\times H^{s1}`$ data as weak limits of smooth solutions. Steps (i) and (iii) are more or less routine, it is (ii) which causes most difficulties. A good starting point is Klainerman’s energy estimate

$$u(t)_{H^{s1}}u(0)_{H^{s1}}\mathrm{exp}(_0^tu(s)_L^{\mathrm{}}𝑑s).$$

This shows that all Sobolev norms of a solution remain bounded for as long as $`u_{L^1L^{\mathrm{}}}`$ stays bounded. It remains to see how to obtain bounds on $`u_{L^1L^{\mathrm{}}}`$. The classical approach uses energy estimates and Sobolev embeddings, but, as in the semilinear case, it only yields results one unit above scaling, namely for $`s>\frac{n}{2}+1`$.

Better results could be obtained using the Strichartz estimates instead. However, this is very nontrivial as one would have to establish the Strichartz estimates for the operator $`\mathrm{}_{g(u)}`$, which has very rough coefficients. Compounding the difficulty, the argument is necessarily circular,

coefficients

regularity

$``$  Strichartz estimates  $``$  solution regularity  $``$  coefficients regularity

One can get around this with a bootstrap argument of the form

$$\begin{array}{c}(u_0,u_1)_{H^s\times H^{s1}}ϵ\\ |g(u)|2\end{array}\}\{\begin{array}{c}\text{Strichartz estimates for }\mathrm{}_{g(u)}\text{ in }[1,1]\\ |g(u)|1\end{array}$$

where the (possibly nonlinear) triple norm contains the needed information about the metric. Still, apriori there is no clear way to determine exactly how it should be defined. A starting point is to set $`|g(u)|=g_{L^1L^{\mathrm{}}}`$, but this only leads to partial results. Following partial results independently obtained by Bahouri-Chemin , and Tataru  and further work of Klainerman-Rodnianski , the next result represents the current state of the problem:

###### Theorem 9.

(Smith-Tataru ) The equation (NLW) is locally well-posed in $`H^s\times H^{s1}`$ for $`s>\frac{n}{2}+\frac{3}{4}`$ ($`n=2`$) and $`s>\frac{n}{2}+\frac{1}{2}`$ ($`n=3,4,5`$). In addition, the Strichartz estimates with $`q=\mathrm{}`$ hold for the corresponding wave operator $`\mathrm{}_{g(u)}`$.

Lindblad’s counterexamples correspond to $`s=\frac{n+3}{4}`$ and show that this result is sharp for $`n=2,3`$. The restriction to $`n5`$ is not central to the problem, it can likely be removed with some extra work.

###### Open Problem 10.

Improve the above result in dimension $`n4`$.

### Wave equation parametrices.

In most approaches, the key element in the proof of the Strichartz estimates is the construction of a parametrix for the wave equation. There are many ways to do this for smooth coefficients, however, as the regularity of the coefficients decreases, they start to break down. Let us begin with the classical Fourier integral operator parametrix, used in the work of Bahouri-Chemin:

$$K(x,y)=a(x,y,\xi )e^{i\varphi (x,y,\xi )}𝑑\xi .$$

The phase $`\varphi `$ is initialized by $`\varphi (x,y,\xi )=\xi (xy)`$ when $`x_0=y_0`$ and must solve an eikonal equation, while for the amplitude $`a`$ one obtains a transport equation along the Hamilton flow. The disadvantage is that all spatial localization comes from stationary phase, which seems to require too much regularity for the coefficients.

One way to address the issue of spatial localization is to begin with wave packets, which have the best possible spatial localization on the unit time scale. In the variable coefficient case the frequency $`\lambda `$ wave packets are bump functions on curved parallelepipeds of size $`1\times \lambda ^1\times (\lambda ^{\frac{1}{2}})^{n1}`$. These parallelepipeds are images of $`\lambda ^1\times (\lambda ^{\frac{1}{2}})^{n1}`$ parallelepipeds at the initial time, transported along the Hamilton flow for $`\mathrm{}_g`$ corresponding to their conormal direction. Then one can seek approximate solutions for $`\mathrm{}_g`$ as discrete superpositions of wave packets, $`u=_Tu_T`$. It is not too difficult to construct individual wave packets, the more delicate point is to show that the wave packets are almost orthogonal. This approach, which is used in , was originally introduced by Smith  and used to prove the Strichartz estimates in $`2`$ and $`3`$ dimensions for operators with $`C^2`$ coefficients.

Another parametrix with a better built in spatial localization can be obtained by doing a smooth phase space analysis:

$$K(\stackrel{~}{y},y)=_Ca(x,\xi )e^{i(\varphi (y,x,\xi )\overline{\varphi (\stackrel{~}{y},x_t,\xi _t)}}𝑑x𝑑\xi 𝑑t\varphi (y,x,\xi )=\xi (xy)+i|\xi |(xy)^2.$$

Here $`(x,\xi )(x_t,\xi _t)`$ is the Hamilton flow for $`\mathrm{}_g`$ on the characteristic cone $`C=\{g^{ij}(x)\xi _i\xi _j=0\}`$. One can factor this into a product of three operators, namely an FBI transform, a phase space transport along the Hamilton flow and then an inverse FBI transform. Neglecting the first one, i.e. setting $`x=y`$ above, produces an operator which is similar to the Fourier integral operators with complex phase. However, it seems to be more useful to keep the Gaussian localizations at both ends. Parametrices of this type were introduced in Tataru  and used to prove Strichartz estimates for operators with $`C^2`$ coefficients in all dimensions. The $`C^2`$ condition was later relaxed in  to $`^2gL^1L^{\mathrm{}}`$. Localization and scaling arguments lead also to weaker estimates for operators whose coefficients have less regularity. Such estimates are known to be sharp, see the counterexamples in Smith-Tataru .

### The null condition.

As in the semilinear case, one may ask whether better results can be obtained for equations with special structure. However, unlike the semilinear case, little is known so far. We propose the following

###### Definition 11.

We say that the equation (GNLW) satisfies the null condition if

$$\frac{g^{ij}(u,p)}{p_k}\xi _i\xi _j\xi _k=0\text{in}g^{ij}(u,p)\xi _i\xi _j=0.$$

###### Conjecture 12.

If the null condition holds then the equation (GNLW) is well-posed in $`H^s\times H^{s1}`$ for some $`s<\frac{n}{2}+\frac{3}{4}`$ ($`n=2`$) respectively for some $`s<\frac{n}{2}+\frac{1}{2}`$ ($`n=3`$).

In $`3+1`$ dimensions a problem which does not quite fit into the above setup but still satisfies some sort of null condition is the Einstein’s equations in general relativity. It is similar to the Yang Mills equations in that it has a gauge invariance, and the null condition is only apparent after fixing the gauge. Klainerman-Rodnianski have obtained a different proof of Theorem 9 for this special case of (NLW).


format: text

text:
# Netons: Vibrations of Complex Networks

## References


text:
# Measurement induced entanglement and quantum computation with atoms in optical cavities

## Abstract

We propose a method to prepare entangled states and implement quantum computation with atoms in optical cavities. The internal state of the atoms are entangled by a measurement of the phase of light transmitted through the cavity. By repeated measurements an entangled state is created with certainty, and this entanglement can be used to implement gates on qubits which are stored in different internal degrees of freedom of the atoms. This method, based on measurement induced dynamics, has a higher fidelity than schemes making use of controlled unitary dynamics.

An essential ingredient in the construction of a quantum computer is the ability to entangle the qubits in the computer. In most proposals for quantum computation this entanglement is created by a controlled interaction between the quantum systems which store the quantum information fortschritte . An alternative strategy to create entanglement is to perform a measurement which projects the system into an entangled state cabrillo ; brian , and based on this principle a quantum computer using linear optics, single photon sources, and single photon detectors has recently been proposed linear . Here we present a similar proposal for measurement induced entanglement and quantum computation on atoms in optical cavities, which only uses coherent light sources and homodyne detection, and we show that high fidelity operation can be achieved with much weaker requirements for the cavity and atomic parameters than in schemes which rely on a controlled interaction.

Several schemes for measurement induced entanglement have been proposed for atoms in optical cavities bose ; plenio ; hong ; duan ; plenioto , but compared to these schemes our proposal has the advantage that entanglement can be produced with certainty even in situations with finite detector efficiency. Furthermore, our procedures are similar to methods already in use to monitor atoms in cavities hood , and we therefore believe that the current proposal should be simpler to implement. Indeed the atom counting procedure in the experiment in Ref. mckeever  could be sufficient to implement this scheme.

We consider atoms with two stable ground states $`|g`$ and $`|f`$ and an excited state $`|e`$ as shown in Fig. 1 (a). Both atoms are initially prepared in an equal superposition of the two ground states, $`(|g+|f)/\sqrt{2}`$. By performing a Quantum Non-Demolition (QND) measurement which measures $`N`$, the number of atoms in state $`|f`$, the state vector is projected into the maximally entangled state $`|\mathrm{\Psi }_{\mathrm{EPR}}=(|gf+|fg)/\sqrt{2}`$ provided we get the outcome $`N=1`$. If the outcomes $`N=0`$ or $`N=2`$ are achieved, the resulting state is $`|gg`$ or $`|ff`$, so that the initial state can be prepared once again, and we can repeat the QND-measurement until the desired entangled state is produced. On average two measurements suffice to produce the state. With $`M`$ atoms in the cavity, the scheme can be extended to the generation of multiparticle entangled states where the population of the $`|f`$ state is distributed symmetrically on all atoms. If, e.g., $`M=3`$ we may produce with certainty the W-states w-state  $`(|ffg+|fgf+|gff)/\sqrt{3}`$ by performing the detection $`4/3`$ times on average. Also, one could entangle a small subset of the atoms by leaving all other atoms in $`|g`$ so that they do not contribute to $`N`$.

The implementation of the QND detection is a modification of a scheme presented in Ref. reflection  where QND detection was achieved by measuring single photons reflected from a cavity. To generate entangled states with a high fidelity this method is undesirable because it is very sensitive to imperfections in the cavity and to imperfect mode matching. Instead we propose to measure the light transmitted through the cavity with homodyne detection. The proposed experimental setup is shown in Fig. 1 (b). The atoms are trapped inside the cavity, and we consider a single field mode described by the annihilation operator $`\widehat{c}`$. Photons in the cavity can decay through two leaky mirrors with decay rates $`\kappa _a`$ and $`\kappa _b`$, and the incoming and outgoing fields at the mirror with decay rate $`\kappa _a`$ ($`\kappa _b`$) are described by $`\widehat{a}_{\mathrm{in}}`$ and $`\widehat{a}_{\mathrm{out}}`$ ($`\widehat{b}_{\mathrm{in}}`$ and $`\widehat{b}_{\mathrm{out}}`$). Light is shined into the cavity in the $`\widehat{a}_{\mathrm{in}}`$ mode, and the transmitted light in the $`\widehat{b}_{\mathrm{out}}`$ mode is measured by balanced homodyne detection. In the cavity the light couples the state $`|f`$ of the $`k`$th atom to the excited state $`|e`$ with a coupling constant $`g_k`$. We mainly consider light far from resonance with the atoms, where the interaction changes the phase of the field by an angle proportional to $`N`$. By measuring the phase of the transmitted light with homodyne detection, we thus obtain the desired QND detection of $`N`$. Spontaneous emission from the excited state limits the amount of light we can send through the cavity, and below we evaluate the fidelity $`F`$, obtainable for a given set of cavity parameters ($`F=\mathrm{\Psi }_{\mathrm{EPR}}|\rho |\mathrm{\Psi }_{\mathrm{EPR}}`$, where $`\rho `$ is the atomic density matrix). The excited state $`|e`$ decays with a rate $`\gamma `$, and for simplicity we assume that the atoms end up in some other state $`|o`$ after the decay. Because these states have a vanishing overlap with the desired state our results will lead to a lower fidelity than if we had assumed that the atoms decay back to the state $`|f`$, but it does not affect the scaling with cavity parameters, which is our main interest here.

The interaction of the cavity light with the atoms is described by the Hamiltonian

$$H=\underset{k}{}g_k|ef|_k\widehat{c}+g_k^{}\widehat{c}^{}|fe|_k+\delta |ee|_k,$$
(1)
where $`\delta `$ is the detuning of the atoms from the cavity resonance, and where the sum is over all atoms in the cavity. The decay of the $`k`$th atom is described by a Lindblad relaxation operator $`d_k=\sqrt{\gamma }|oe|_k`$, and the input/output relations are given by $`d\widehat{c}/dt=i[\widehat{c},H](\kappa /2)\widehat{c}+\sqrt{\kappa _a}\widehat{a}_{\mathrm{in}}+\sqrt{\kappa _b}\widehat{b}_{\mathrm{in}}+\sqrt{\kappa ^{}}\widehat{F}`$, $`\widehat{a}_{\mathrm{out}}=\widehat{a}_{\mathrm{in}}\sqrt{\kappa _a}\widehat{c}`$, and $`\widehat{b}_{\mathrm{out}}=\widehat{b}_{\mathrm{in}}\sqrt{\kappa _b}\widehat{c}`$, where $`\kappa =\kappa _a+\kappa _b+\kappa ^{}`$ is the total decay rate, and $`\kappa ^{}`$ and $`\widehat{F}`$ is the decay rate and noise operator associated with cavity loss to other modes than the two modes considered. For convenience $`\widehat{c}`$ has the standard normalization for a single mode $`[\widehat{c},\widehat{c}^{}]=1`$, whereas the normalization for the free fields are such that, e.g., $`[\widehat{a}_{\mathrm{in}}(t),\widehat{a}_{\mathrm{in}}(t^{})^{}]=\delta (tt^{})`$. We assume that the atoms are separated by more than the resonant optical wavelength, and we hence ignore the dipole-dipole interaction between the atoms. We solve the equations of motion by assuming that the fields are sufficiently weak that only the lowest order terms in $`\widehat{c}`$ are important. Note that this limits the total number of photon in the cavity at any time, but the total number of photons in a pulse may still be large. By taking the Fourier transform we find the cavity field to be given by reflection

$$\widehat{c}(\omega )=\frac{\sqrt{\kappa _a}\widehat{a}_{\mathrm{in}}(\omega )}{\frac{\kappa }{2}i\omega +\frac{g^2\widehat{N}}{\gamma /2+i(\delta \omega )}}+\mathrm{noise},$$
(2)
where we have introduced the operator $`\widehat{N}=_k|ff|_k(t=0)`$ which counts the number of atoms in state $`|f`$ before the pulse, and we have assumed that the magnitude of all coupling constants are identical $`|g_k|^2=g^2`$. Using Eq. (2) we find an expression for the transmitted light

$$\widehat{b}_{\mathrm{out}}(\omega )=\frac{\sqrt{\kappa _a\kappa _b}}{\frac{\kappa }{2}i\omega +\frac{g^2\widehat{N}}{\gamma /2+i(\delta \omega )}}\widehat{a}_{\mathrm{in}}(\omega )+\mathrm{noise}.$$
(3)
We determine the Heisenberg equations of motion for the atoms from Eq. (1) and the relaxation operators and solve them by assuming that the atomic ground state operators vary only little on the time scale $`1/\gamma `$. Here we shall only need the projection operator onto the $`|f`$ state, which after the interaction with the light pulse $`(t=\mathrm{})`$ is given by $`|ff|(t=\mathrm{})=:|ff|(t=0)\mathrm{exp}(\widehat{\theta }):`$, where $`::`$ denotes normal ordering, and where we have introduced the scattering probability operator

$$\widehat{\theta }=\gamma 𝑑\omega \frac{g^2}{\frac{\gamma ^2}{4}+(\delta \omega )^2}\widehat{c}(\omega )^{}\widehat{c}(\omega ).$$
(4)
The incoming light is assumed to be in a coherent state with frequency $`\omega `$ and since Eq. (3) is linear the outgoing field will also be in a coherent state $`|\beta _N`$. Using Eqs. (2)-(4) we see that the difference between the coherent state output of an empty cavity and that of a cavity with $`N`$ atoms can be expressed as

$$|\beta _N\beta _0|^2=N^2\frac{g^2\kappa _b}{\gamma (\frac{\kappa ^2}{4}+\omega ^2)}\theta _N,$$
(5)
The scattering probability per atom $`\theta _N`$ \[obtained by inserting Eq. (2) in Eq. (4)\] depends on the number of atoms, because the field amplitude inside the cavity depends on $`N`$. From Eq. (5) we see that the best signal to noise ratio ($`|\beta _N\beta _0|^2/\theta _N`$) is obtained by use of light which is resonant with the cavity, $`\omega =0`$. For resonant light, the transmittance of the output mirror should match all other losses $`\kappa _b=\kappa _a+\kappa ^{}=\kappa /2`$, and the QND detection then only depends on the cavity parameters through the combination $`g^2/\kappa \gamma `$. In particular, $`|\beta _N\beta _0|^2=N^2\theta _Ng^2/\kappa \gamma `$. If $`\kappa _b`$ deviates from $`\kappa /2`$ the effective $`g^2/\kappa \gamma `$ is changed by $`2\kappa _b/\kappa `$. If transmitted photons are lost, or if they are only detected with a probability $`\eta `$, the effective $`g^2/\kappa \gamma `$ is reduced by this factor.

If we need to distinguish an empty cavity from a cavity with a certain number of atoms, the signal to noise ratio is independent of the detuning of the cavity from the atomic resonance, but here we concentrate on the far detuned case $`\delta g^2/\kappa `$, $`\gamma `$, where the $`N=1`$ and $`N=2`$ components are easier to distinguish. In this limit the intra cavity field and therefore the scattering probability $`\theta _N`$ is independent of $`N`$ and we replace $`\theta _N`$ by the same $`\theta `$ for all $`N`$. If the light is on resonance with the cavity $`\omega =0`$, and far detuned from the atomic resonance, we see from Eq. (3) that the interaction with the atoms changes the phase of the field by an angle proportional to $`N`$. By adjusting the phase of the local oscillator we can configure the measurement so that the photocurrent difference, integrated over the duration of the pulse, can be normalized to a dimensionless field quadrature operator $`x`$ which, in turn, is proportional to the phase change imposed by the atoms and, hence, it provides a QND detection of $`N`$. If the atoms are in the $`N=0`$ state, light is transmitted through the cavity without being affected by the atoms, and the homodyne measurement provides a value for $`x`$ in accordance with a Gaussian probability distribution $`p_0(x)=\mathrm{exp}(x^2)/\sqrt{\pi }`$ with a variance $`(\mathrm{\Delta }x)^2=1/2`$ due to vacuum noise in the coherent state. If $`N=1`$ or 2, Eq. (5) predicts the Gaussian to be displaced so that it is centered around $`x_N=2N\sqrt{\theta g^2/\kappa \gamma }`$.

So far we have ignored that the atoms are slowly pumped out of the state $`|f`$ due to decay of the upper state, and this changes the distribution. If we start with $`N`$ atoms in state $`|f`$, the probability distribution $`p_N(x)`$ can be written as a sum of a part $`p_{\mathrm{ND},N}(x)`$ where there has been no atomic decay and a part $`p_{\mathrm{D},N}(x)`$ where at least one of the atoms have decayed: $`p_N(x)=p_{\mathrm{ND},N}(x)+p_{\mathrm{D},N}(x)`$. With no decay Eq. (5) is valid and we have $`p_{\mathrm{ND},N}(x)=\mathrm{exp}(N\theta )p_0(xx_N)`$. If there is a decay, Eq. (3) is valid until the time of the decay and we find $`p_{\mathrm{D},1}=_0^\theta 𝑑\theta ^{}\mathrm{exp}(\theta ^{})p_0(xx_1\theta ^{}/\theta )`$ and $`p_{\mathrm{D},2}=2_0^\theta 𝑑\theta ^{}\mathrm{exp}(\theta \theta ^{})p_0(xx_1(\theta ^{}+\theta )/\theta )+_0^\theta 𝑑\theta ^{}_0^\theta 𝑑\theta ^{\prime \prime }\mathrm{exp}(\theta ^{}\theta ^{\prime \prime })p_0(xx_1(\theta ^{}+\theta ^{\prime \prime })/\theta )`$. If the homodyne measurement produces an outcome $`\stackrel{~}{x}`$ in a region $`x_a<\stackrel{~}{x}<x_b`$ around $`x_1`$, the atoms are most likely in the desired entangled state and the measurement is successful. The probability for this measurement outcome is $`P_s=_{x_a}^{x_b}𝑑\stackrel{~}{x}P_{\mathrm{tot}}(\stackrel{~}{x})`$, where the total probability distribution for $`\stackrel{~}{x}`$ is given by $`P_{\mathrm{tot}}(\stackrel{~}{x})=p_0(\stackrel{~}{x})/4+p_1(\stackrel{~}{x})/2+p_2(\stackrel{~}{x})/4`$.

The fidelity of the entangled state vanishes if one of the atoms has decayed to a different subspace, and the fidelity of the state produced is therefore given by the conditional probability of starting in the state with $`N=1`$ and not having decayed. By averaging over the interval of accepted values we get the fidelity, $`F=_{x_a}^{x_b}𝑑\stackrel{~}{x}p_{\mathrm{ND},1}(\stackrel{~}{x})/2/P_s`$.

For a given value of $`g^2/\kappa \gamma `$ the success probability and fidelity depend on $`\theta `$, $`x_a`$ and $`x_b`$. To get a high fidelity one must choose a small interval around $`x_1`$ because it minimizes the contribution from $`p_0`$ and $`p_2`$, but a small interval also reduces the success probability. In Fig. 2 we show the smallest error probability $`1F`$ achievable as a function of $`g^2/\kappa \gamma `$, when we vary $`x_a`$ and $`\theta `$ and adjust $`x_b`$ to retain a fixed value of $`P_s`$. In the figure the full curves show $`1F`$ for $`P_s=0.1\%`$, $`P_s=30\%`$ and $`P_s=50\%`$.

In Fig. 2 we observe that the fidelity scales as

$$1F\frac{\kappa \gamma }{g^2}\mathrm{log}\left(\frac{g^2}{\kappa \gamma }\right).$$
(6)
This behavior can be understood from Eq. (5): with a fixed distance $`|\beta _1\beta _0|=\sqrt{2}x_1`$, the error due to scattering decreases as $`\kappa \gamma /g^2`$, but due to the overlap of the Gaussians there is a probability of accepting a wrong value for $`N`$, which will dominate over the error due to scattering if we keep the distance fixed and increase $`g^2/\kappa \gamma `$. Because the overlap of the Gaussians decreases exponentially with the square of the distance, optimum fidelity is obtained if the distance grows logarithmically with $`g^2/\kappa \gamma `$ and this gives the scaling in Eq. (6).

The scaling (6) is valid up to a success rate $`P_s=50\%`$, where $`P_{\mathrm{tot}}(\stackrel{~}{x})`$ consists of three well separated Gaussians, so that even when the scheme is not successful we know in which state ($`N=0`$ or $`N=2`$) the atoms are left. From these states we can deterministically produce the initial product state, and we can repeat the measurement procedure until it is successful. A lower limit on the fidelity with this procedure can be shown to be $`F>_{x_a}^{x_b}𝑑\stackrel{~}{x}p_{\mathrm{ND},1}(\stackrel{~}{x})/2/(1_{\mathrm{}}^{x_a}𝑑\stackrel{~}{x}p_0(\stackrel{~}{x})/4_{x_b}^{\mathrm{}}𝑑\stackrel{~}{x}p_{\mathrm{ND},2}(\stackrel{~}{x})/4)`$. This lower limit is shown with a dashed line in Fig. 2, and the scaling still follows (6).

The deterministic entanglement scheme can be used to implement gates between qubits. Suppose that the levels $`|f`$ and $`|g`$ are both degenerate with two substates. For convenience the four states $`|g_0,|g_1,|f_0,|f_1`$ can be described as tensor products of the ’level’ degree of freedom ($`|g`$ and $`|f`$) and ’qubit’ degree of freedom ($`|0`$ and $`|1`$). (This would be the natural representation if $`|0`$ and $`|1`$ represent, e.g., different motional state, but they may also be different internal states in the atoms, e.g., states with different magnetic quantum numbers). The qubits are initially stored in the states $`|g_0`$ and $`|g_1`$ and, we assume that the interactions during the entanglement preparation are symmetric with respect to the qubit degrees of freedom, so that the qubits are decoupled from the preparation of the entangled state of the levels, i.e, we prepare the state $`|\mathrm{\Phi }_q(|gf+|fg)/\sqrt{2}`$, where $`|\mathrm{\Phi }_q`$ is the initial state of the qubits. The entanglement in the levels $`|f`$ and $`|g`$ can now be used to implement a control-not operation on the qubits as discussed in Ref. cnot : we interchange states $`|f_1`$ and $`|g_1`$ in the control atom, and then measure if this atom is in the $`|f`$ or the $`|g`$ level. If the measurement outcome is $`|g`$ ($`|f`$) the target atom is known to be in $`|g`$ ($`|f`$) if the control atom is in $`|1`$, and the control-not operation amounts to interchanging the states $`|g_0`$ and $`|g_1`$ ($`|f_0`$ and $`|f_1`$) in the target atom. To complete the gate a $`\pi /2`$ pulse is applied between the levels $`|f`$ and $`|g`$, and a QND-measurement of the level of the target atom, followed by single atom transitions, conditioned on the outcome of the measurement leaves the system in $`|\stackrel{~}{\mathrm{\Phi }}_q|gg`$, where $`|\stackrel{~}{\mathrm{\Phi }}_q`$ is the desired outcome of the control-not operation on $`|\mathrm{\Phi }_q`$. The single atom measurements required for the gate operation can be done by homodyne detection of light transmitted through the cavity, and the procedure is accomplished with a fidelity scaling as in Eq. (6).

A different strategy for quantum computation in optical cavities is to construct a controlled interaction between the atoms and the cavity field pellizzari ; domokos ; zheng : If an atom emits a photon into the cavity mode this photon can be absorbed by another atom, and this constitutes an interaction which can be used to engineer quantum gates. With controlled interactions, a first estimate of the error probability $`ϵ`$ in a single gate is given by $`ϵ\mathrm{max}\{\gamma ,\kappa \}/g`$, but this can be improved by use of Raman transitions or a large detuning between the atoms and the cavity mode. If the optical transition is replaced by a Raman transition where a classical beam with Rabi frequency $`\mathrm{\Omega }`$ and detuning $`\mathrm{\Delta }`$ induces a transition between two ground states and the simultaneous creation/annihilation of a photon in the cavity, the effective coupling constant is $`g_{\mathrm{eff}}g\mathrm{\Omega }/\mathrm{\Delta }`$ and the effective decay rate is $`\gamma _{\mathrm{eff}}\gamma \mathrm{\Omega }^2/\mathrm{\Delta }^2`$. By decreasing $`\mathrm{\Omega }/\mathrm{\Delta }`$ we can thus reduce the error probability $`\gamma _{\mathrm{eff}}/g_{\mathrm{eff}}`$ due to spontaneous emission. Similarly it has also been proposed to reduce the effect of cavity decay by using transitions which are detuned from the cavity zheng . If the Raman transitions are detuned by an amount $`\delta `$ from the cavity mode the effective coupling constant between the atoms is $`\chi g_{\mathrm{eff}}^2/\delta `$ and the leakage rate is $`\kappa _{\mathrm{eff}}\kappa g_{\mathrm{eff}}^2/\delta ^2`$, so that the leakage can be reduced by making $`\delta \kappa `$. But, detuning from the cavity mode increases the effect of spontaneous emission because the gate takes longer time. Since the gate duration is $`t1/\chi `$ the probability for spontaneous emission is $`ϵ_\gamma \gamma \delta /g^2`$ and the probability for a cavity decay is $`ϵ_\kappa \kappa /\delta `$. The minimum error probability $`1F=ϵ_\gamma +ϵ_\kappa `$ is then found to scale as

$$1F\sqrt{\frac{\kappa \gamma }{g^2}}.$$
(7)
Although Eq. (7) has been derived for a particular setup, we believe the scaling in Eq. (7) is characteristic for all existing proposals which use a controlled interaction between the atoms. (Also the scheme in the Ref. beige  has the scaling (7) if the photodetectors have finite efficiency, or if the scheme is required to implement gates deterministically). The fundamental problem is that the atomic decay and the cavity loss provide two sources of errors which combine to give the scaling in Eq. (7). For a given cavity, i.e., for given parameters $`g`$, $`\kappa `$, and $`\gamma `$, we obtain the better scaling presented in Eq. (6) by using a scheme where classical light is injected into the cavity and where leakage of a single photon is not a critical event.

In conclusion, we have presented a scheme to entangle atoms and implement gates between qubits which are stored in the atoms. The error probability scales more favorably than it does in schemes which use controlled interactions and unitary dynamics, and the present scheme is therefore advantageous if a high fidelity is required. The best optical cavities currently reach $`g^2/\kappa \gamma 100`$ hood  so that Eq. (7) predicts $`F90\%`$ which is less than the value in Fig. 2. More importantly, to decrease the error rate in the measurement induced scheme, e.g., by a given factor we only need to improve the cavity finesse, which is proportional to $`g^2/\kappa \gamma `$, by the same factor. With the scaling in Eq. (7) the cavity finesse should be increased by the same factor squared.

This work was supported by the Danish Natural Science Research Council and by the Natural Science Foundation through its grant to ITAMP.


text:
# A Code for Stellar Binary Evolution and its Application to the Formation of Helium White Dwarfs

## 1 Introduction

Stellar evolution in binary systems is a classical topic of stellar astrophysics. The study of stellar evolution in this kind of systems has great importance, since an important fraction of know stars belong to binary systems. Inasmuch as, if we want to build a stellar evolution theory able to account for the observations, it is essential to study stellar evolution in binary systems.

According to the standard model, the initially most massive star, usually named the primary evolves faster and, as consequence of its nuclear evolution, begin to inflate. If stars are sufficiently close each other (close binary systems, hereafter CBS) the primary will overflow of the “Roche lobe”. This lobe is found by considering the problem of a binary system in which the orbits are circular, and rotation of both objects is synchronized with orbital motion. Therefore we have the Largange solution to the restricted three body problem, considering that each gas particle is of infinitesimal mass. The Roche lobe corresponds to an equipotential surface that surrounds the stars and represents the limiting volume that can be attained by the stars of the pair before the beginning of mass transfer either, onto the companion, and/or outside the system. From this moment on, the evolution of each star is markedly different from the one they would have in isolation.

According to Kippenhahn & Weigert (1967) mass transfer may begin during core hydrogen burning (case A) and also during central contraction after central hydrogen exhaustion (case B). If it begins after core helium burning it is commonly referred to as case C (Lauterborn 1970). Mass transfer can occur in conservative conditions (total mass and orbital angular momentum conserved) or with mass and/or angular momentum losses from the system. These losses can drastically affect the evolution of the star in a CBS. For reviews on stellar evolution in binary systems, see, e.g., Paczyński (1971), Iben (1985), Iben & Livio (1993), and Taam & Sandquist(2000).

In view of the large number of situations that can be encountered, it is not surprising that CBSs are considered as favorable places for the occurrence of a number of phenomenologies with high astrophysical interest: X-ray sources (van der Klis 2002), low mass helium white dwarfs (Kippenhahn, Kohl & Weigert 1967), planetary nebulae with binary nucleus (Iben & Livio 1993), binary radio pulsars (Bhattacharya & van den Heuvel 1991), supernova progenitors (Podsiadlowski, Joss & Hsu 1992). In this last case, binary evolution naturally accounts for the evolution of the progenitor of SN 1987A. This object was a blue supergiant, a fact difficult to explain in the frame of single stellar evolution, but not in binary evolution (Podsiadlowski, et al. 1992).

In view of the importance of this line of research, in our Observatory we have begun the study of CBSs. To this end, we have developed a code to compute stellar evolution in CBSs, including the self consistent calculation of the mass transfer rate $`\dot{M}`$ (hereafter MTR). The main purpose of the present paper is to describe the numerical techniques employed in our program for the computation the evolution of stars in a CBS. This code is based on a modification of the Henyey technique presented in Kippenhahn, Weigert & Hofmeister (1967) to solve the set of difference equations of stellar evolution with the inclusion of hydrodynamical effects and the MTR. As a first application of this new code to the case of CBSs that lead to the formation of helium white dwarfs.

As stated above, the main characteristic of our code is that, during mass transfer episodes, it finds the MTR in a fully implicit way by means of a single iterative procedure. This is in contrast with the usual procedure employed by other researchers (see, e.g., Podsiadlowski, Rappaport, & Pfahl 2001; Wellstein, Langer & Braun 2001). In these works, in the construction of a stellar model a double iterative procedure is applied: a MTR is estimated, and the stellar structure is relaxed for this MTR value. Then, with this new structure, a new MTR is computed and the stellar structure is relaxed again. This is repeated until consistency is achieved. In calculating the MTR it is usual to employ the physical treatment presented by Ritter (1988) which gives

$$\dot{M}=\dot{M_0}\mathrm{exp}\left(\frac{R_LR}{H_P}\right);$$
(1)
where $`\dot{M_0}`$ is the MTR for a star that exactly fills the Roche lobe, $`R_L`$ is the equivalent radius of the Roche lobe (see below), $`R`$ is the stellar radius and $`H_P`$ is the photospheric pressure scale height. Notice that the double iterative loop is necessary to get numerical stability because of the steep dependence of the MTR upon the sizes of the star and the Roche lobe.

As the main goal of this paper is to present the numerical scheme we have devised and not to produce start - of - the - art evolutionary models here, for simplicity, we shall simply impose $`R_L=R`$. In doing so, we shall get a good global computation (but not the details of the beginning and end) of each mass transfer episode.

The rest of the paper is organized as follows: In Section 2 we briefly describe the differential equations of single stellar evolution we have to solve. In the following Section 3, we describe in detail the Henyey technique and also the code employed to compute binary stellar evolution. We describe our treatment for the orbit of the binary, including the modification we have considered for total mass and angular momentum losses from the system in §4. In Section 5 we discuss the critical problem of how to manage the outer boundary conditions that, in turn, are essential in determining the MTR during mass transfer episodes. Then, in Section 6 we discuss some numerical techniques we have applied, especially the rezoning of the stellar models. In Section 7, we present the results for three different binary systems that lead to the formation of helium white dwarfs. Finally, in Section 8 we comment on our results, further improvements we plan to incorporate to the code and the astrophysical problems we plan to study with.

## 2 The Equations of Stellar Evolution

Here, we briefly summarize the equations of stellar evolution to be solved. As usual, we consider spherically symmetric objects, neglecting rotation and magnetic fields. In these conditions, the equations of stellar structure are (for derivation of these equations see, e.g., Clayton 1968, Kippenhahn & Weigert 1990. For a detailed treatment of hydrodynamic stellar codes, see, e.g., Kutter & Sparks 1972):

i) the Euler equation of fluid motion

$$\frac{v}{t}+v\frac{v}{r}=\frac{1}{\rho }\frac{P}{r}\frac{Gm_r}{r^2},$$
(2)
ii) the definition of velocity

$$\frac{r}{t}=v,$$
(3)
iii) the equation of mass conservation

$$\frac{m_r}{r}=4\pi r^2\rho ,$$
(4)
iv) the equation of energy balance

$$\frac{l_r}{r}=4\pi r^2\rho \left(\epsilon _{nuc}\epsilon _\nu T\frac{S}{t}\right),$$
(5)
v) the equation of energy transport for the radiative case

$$\frac{T}{r}=\frac{3\pi }{4ac}\frac{\kappa \rho }{T^3}\frac{l_r}{4\pi ^2r^2},$$
(6)
and

vi) the equation of energy transport for the convective case

$$\frac{T}{r}=_{conv}\frac{T}{P}\frac{P}{r}.$$
(7)
We employ the Schwarzschild criterium for the onset of convection. The symbols have their standard meaning.

## 3 The Hydrodynamic Code for Binary Stellar Evolution

In order to incorporate the specific phenomena occurring in binary evolution we have to make supplementary assumptions apart from those quoted in the previous section. We shall handle the members of our system as spherical objects, neglecting the departure from spherical symmetry of the equipotentials (e.g., the Roche lobe) and its evolutionary consequences. Moreover, we shall consider that the objects move along circular orbits, and that they influence each other only through gravitational attraction (we neglect irradiation).

As usual we shall consider the problem in Lagrangian coordinates, considering $`\xi `$ the independent variable, defined as

$$\xi =\mathrm{ln}\left(1\frac{m_r}{M}\right).$$
(8)
Radius, pressure, and temperature are handled by means of logarithmic transformations

$`p=\mathrm{ln}P`$,
$`\theta =\mathrm{ln}T`$,
$`x=\mathrm{ln}r`$,

whereas $`l_r`$, $`v`$ are considered linearly.

For simplicity, we have written the difference equations in a centered fashion. It means that we have chosen to write a generic differential equation

$$\frac{dy_i}{dx}=F(x,y_1,\mathrm{},y_5);i=1,\mathrm{},5$$
(9)
as a difference equation

$$\frac{y_{i,j+1}y_{i,j}}{x_{j+1}x_j}F(x_{j+1/2},y_{1,j+1/2},\mathrm{},y_{5,j+1/2})$$
(10)
where $`\eta _{j+1/2}=(\eta _{j+1}+\eta _j)/2`$, being $`\eta `$ any quantity. The second subindex $`j`$ indicates the shell of the star for which the difference equation is written. The results obtained in this way are completely satisfactory for our present purposes. However, this would have not been the case if we were studying stellar objects suffering shock waves somewhere. For calculations of more violent phases we shall have to include non - centered equations together with some dissipative effect (e.g., artificial viscosity, see Kutter & Sparks 1972). Temporal derivatives have been written in the standard backwards differenced form.

In problems with variable stellar mass we need to be careful at computing temporal derivatives at constant mass. We have found it very convenient to rewrite the derivative operator as

$$\frac{}{t}|_{m_r}=\frac{}{t}|_\xi +\frac{\xi }{t}|_{m_r}\frac{}{\xi }|_t.$$
(11)
In this formulation we get the dependence of these derivatives with the MTR. This is important because we shall treat the MTR as a new variable to be relaxed in our code. Then we have to design a Henyey scheme for relaxing the internal structure of the star, the total luminosity, the effective temperature and (in the case of the occurrence of mass transfer) the MTR.

As usual, the calculation of the corrections is based in the solution of the lineal system

$$H\delta =h$$
(12)
where $`H`$ is a matrix with a structure similar to those corresponding to the cases of single star evolution or codes in which the MTR is not treated self consistently. It has blocks on the diagonal, but now, an important difference appears: now there is a non vanishing column due to the derivatives of the structure equations with respect to the MTR. The vector $`\delta `$ represents the corrections vector (see below) and $`h`$ is the inhomogeneity vector that must vanish for a correct stellar model.

Here, some words are in order about the necessity of such a deep modification of the standard numerical treatments of binary evolution. Savonije (1978) has treated MTR as another variable to be relaxed by iterations but only considered it to be relevant in the outer layers integration where he included most of the change in luminosity due to mass transfer (See Section 5 for details). This is so only in the case of excluding a very thick portion of the star from the grid. In this way we are neglecting the temporal derivatives, and so, the inertia of the outermost layers of the star. Other scheme intended for the computation of single and binary hydrostatic stellar evolution has been presented by Ziolkowski (1970), who incorporated an iterative procedure for getting the MTR. However, the outer envelope integrations are performed in a very thick portion of the star, and because of the reasons quoted above, this is not entirely adequate for our purposes

In order to perform a realistic simulation of the evolution of the stars, it is quite desirable to handle most of the mass of the star inside finite - difference part of the Henyey scheme. Notice that there are many relevant situations in which the outer layers structure is a key ingredient in determining the evolution of the stars. Regarding binary evolution, we may cite for example that in order to compute the mass transfer episodes occurring in low mass pre - white dwarf stars we do need to consider very thin outer layers integrations (see below for more details).

Consequently, we have to consider an outer layers integration in a portion of the star so thin that a fraction, or even most, of the luminosity change (see Section 5 for details) due to mass loss during mass transfer episodes occurs inside the finite - differenced portion of the model. Thus, the form of the equation of energy conservation we have to consider is

$$\frac{l_r}{m_r}=\epsilon _{nuc}\epsilon _\nu T\left(\frac{S}{t}|_\xi +\frac{d\xi }{dt}|_{m_r}\frac{S}{\xi }|_t\right),$$
(13)
where, from Eq. (8), we have

$$\frac{d\xi }{dt}|_{m_r}=\frac{\dot{M}}{M}\left(e^\xi 1\right).$$
(14)
Thus, we are forced to consider $`\dot{M}`$ as a extra Henyey unknown.

Notice that, for consistency, we should consider the total mass value of the star as another unknown to be found during the relaxation process by means of the equation

$$M=M^{prev}+\dot{M}\mathrm{\Delta }t,$$
(15)
where $`M^{prev}`$ is the mass of the previous model. Because of the transformation of variables we have considered, the total mass value occurs in most of the equations of structure. Thus, the matrix we shall have to handle in finding the corrections will have the usual blocks on the diagonal but also a non - vanishing column. We shall describe in the rest of this section the technique we employed to solve the numerical problem.

The first block of the matrix includes the four boundary conditions equations that link the values of the structural quantities at the first meshpoint with the effective temperature, luminosity and MTR value. These are constructed by meand of linear interpolations among the outer layers integrations. The resulting block has the form

$$\left[\begin{array}{ccccccc}\frac{h_1}{\dot{M}}& \frac{h_1}{L}& \frac{h_1}{T_{eff}}& \frac{h_1}{x_1}& \frac{h_1}{l_1}& \frac{h_1}{p_1}& \frac{h_1}{\theta _1}\\ \frac{h_2}{\dot{M}}& \frac{h_2}{L}& \frac{h_2}{T_{eff}}& \frac{h_2}{x_1}& \frac{h_2}{l_1}& \frac{h_2}{p_1}& \frac{h_2}{\theta _1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{h_5}{\dot{M}}& \frac{h_5}{L}& \frac{h_5}{T_{eff}}& \frac{h_5}{x_1}& \frac{h_5}{l_1}& \frac{h_5}{p_1}& \frac{h_5}{\theta _1}\end{array}\right]\left[\begin{array}{c}\delta \dot{M}\\ \delta L\\ \delta T_{eff}\\ \delta x_1\\ \delta l_1\\ \delta p_1\\ \delta \theta _1\end{array}\right]$$

$$=\left[\begin{array}{ccc}h_1& & \\ h_2& & \\ \mathrm{}& & \\ h_5& & \end{array}\right],$$
(16)
where $`h_i,i=1,\mathrm{},5`$ are functions known from linealizacion of stellar structure equations. Hereafter in this section, $`L`$ is $`log(L/L)`$ and $`T_{eff}`$ the logarithm of effective temperature.

We define vectors $`U_i`$, $`V_i`$, $`W_i`$, $`X_i`$, $`i=1,\mathrm{},5`$, such that

$$\left[\begin{array}{c}\delta \dot{M}\\ \delta L\\ \delta T_{eff}\\ \delta x_1\\ \delta l_1\end{array}\right]=\left[\begin{array}{cccc}U_1& V_1& W_1& X_1\\ U_2& V_2& W_2& X_2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ U_5& V_5& W_5& X_5\end{array}\right]\left[\begin{array}{c}\delta v_1\\ \delta p_1\\ \delta \theta _1\\ 1\end{array}\right].$$
(17)
For the first block are $`U_i=0`$, $`i=1,\mathrm{},5`$ and $`\delta v_1=0`$.
From Eq. (16) and Eq. (17) we can write

$$\left[\begin{array}{ccccc}\frac{h_1}{\dot{M}}& \frac{h_1}{L}& \frac{h_1}{T_{eff}}& \frac{h_1}{x_1}& \frac{h_1}{l_1}\\ \frac{h_2}{\dot{M}}& \frac{h_2}{L}& \frac{h_2}{T_{eff}}& \frac{h_2}{x_1}& \frac{h_2}{l_1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{h_5}{\dot{M}}& \frac{h_5}{L}& \frac{h_5}{T_{eff}}& \frac{h_5}{x_1}& \frac{h_5}{l_1}\end{array}\right]\left[\begin{array}{cccc}U_1& V_1& W_1& X_1\\ U_2& V_2& W_2& X_2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ U_5& V_5& W_5& X_5\end{array}\right]$$

$$=\left[\begin{array}{ccc}\frac{h_1}{p_1}& \frac{h_1}{\theta _1}& h_1\\ \frac{h_2}{p_1}& \frac{h_2}{\theta _1}& h_2\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{h_5}{p_1}& \frac{h_5}{\theta _1}& h_5\end{array}\right],$$
(18)
and from this we can obtain vector components $`V_i`$, $`W_i`$, $`X_i`$, $`i=1,\mathrm{},5`$. In a similar way we can generalize this mechanism for any interior block. If we propose

$$\delta \dot{M}=A_n\delta v_n+B_n\delta p_n+C_n\delta \theta _n+D_n,$$
(19)
we find, generalizing the equations of the interiors blocks,

| $`A_n=A_{n1}U_{5n4}+B_{n1}U_{5n3}+C_{n1}U_{5n2}`$ |
| --- |
| $`B_n=A_{n1}V_{5n4}+B_{n1}V_{5n3}+C_{n1}V_{5n2}`$ |
| $`C_n=A_{n1}W_{5n4}+B_{n1}W_{5n3}+C_{n1}W_{5n2}`$ |
| $`D_n=A_{n1}X_{5n4}+B_{n1}X_{5n3}+C_{n1}X_{5n2}+D_{n1}`$. |

(20)
The equation that leads to intermediates blocks are in the form

$$\left(\delta v_n\alpha _i+\delta p_n\beta _i+\delta \theta _n\gamma _i\right)+\delta x_{n+1}\frac{g_i}{x_{n+1}}+\delta l_{n+1}\frac{g_i}{l_{n+1}}=$$

$$\delta v_{n+1}\frac{g_i}{v_{n+1}}\delta p_{n+1}\frac{g_i}{p_{n+1}}\delta \theta _{n+1}\frac{g_i}{\theta _{n+1}}\xi _i,$$
(21)
$`i=1,2,\mathrm{},5n=2,3,\mathrm{},N2`$.

The $`g_i`$ functions are know from linealizacion of stellar structure equations in the intermediate shells, and $`\alpha _i,\beta _i,\gamma _i,\xi _i`$ are defined as

| $`\alpha _i`$ = $`A_n\frac{g_i}{\dot{M}}`$ \+ $`U_{5n1}\frac{g_i}{x_n}`$ \+ $`U_{5n}\frac{g_i}{l_n}`$ \+ $`\frac{g_i}{v_n},`$ |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
|  |  |  |  |  |  |  |  |  |
| $`\beta _i`$ = $`B_n\frac{g_i}{\dot{M}}`$ \+ $`V_{5n1}\frac{g_i}{x_n}`$ \+ $`V_{5n}\frac{g_i}{l_n}`$ \+ $`\frac{g_i}{p_n},`$ |  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |  |
| $`\gamma _i`$ = $`C_n\frac{g_i}{\dot{M}}`$ \+ $`W_{5n1}\frac{g_i}{x_n}`$ \+ $`W_{5n}\frac{g_i}{l_n}`$ \+ $`\frac{g_i}{\theta _n},`$ |  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |  |
| $`\xi _i`$ = $`D_n\frac{g_i}{\dot{M}}`$ \+ $`X_{5n1}\frac{g_i}{x_n}`$ \+ $`X_{5n}\frac{g_i}{l_n}`$ \+ $`1.`$ |  |  |  |  |  |  |  |  |
|  |  |  |  |  |  |  |  |  |

Written in matrix form Eq. (21), with

$$\left[\begin{array}{c}\delta v_n\\ \delta p_n\\ \delta \theta _n\\ \delta x_{n+1}\\ \delta l_{n+1}\end{array}\right]=\left[\begin{array}{cccc}U_{5n+1}& V_{5n+1}& W_{5n+1}& X_{5n+1}\\ U_{5n+2}& V_{5n+2}& W_{5n+2}& X_{5n+2}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ U_{5n+5}& V_{5n+5}& W_{5n+5}& X_{5n+5}\end{array}\right]$$

$$\left[\begin{array}{c}\delta v_{n+1}\\ \delta p_{n+1}\\ \delta \theta _{n+1}\\ 1\end{array}\right]$$
(22)
and substituting Eq. (22) in matrix form of Eq. (21), we obtain an equation that allow us to find the components for vectors $`U,V,W,X`$:

$$\left[\begin{array}{ccccc}\alpha _1& \beta _1& \gamma _1& \frac{g_1}{x_{n+1}}& \frac{g_1}{l_{n+1}}\\ \alpha _2& \beta _2& \gamma _2& \frac{g_2}{x_{n+1}}& \frac{g_2}{l_{n+1}}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \alpha _5& \beta _5& \gamma _5& \frac{g_5}{x_{n+1}}& \frac{g_5}{l_{n+1}}\end{array}\right]$$

$$\left[\begin{array}{cccc}U_{5n+1}& V_{5n+1}& W_{5n+1}& X_{5n+1}\\ U_{5n+2}& V_{5n+2}& W_{5n+2}& X_{5n+2}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ U_{5n+5}& V_{5n+5}& W_{5n+5}& X_{5n+5}\end{array}\right]=$$

$$\left[\begin{array}{cccc}\frac{g_1}{v_{n+1}}& \frac{g_1}{p_{n+1}}& \frac{g_1}{\theta _{n+1}}& \xi _1\\ \frac{g_2}{v_{n+1}}& \frac{g_2}{p_{n+1}}& \frac{g_2}{\theta _{n+1}}& \xi _2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ \frac{g_5}{v_{n+1}}& \frac{g_5}{p_{n+1}}& \frac{g_5}{\theta _{n+1}}& \xi _2\end{array}\right].$$
(23)
In the last block we have $`\delta l_n=\delta x_n=\delta v_n=0`$, thus, rewriting the previous matrix equation in the usual way and proposing, like we have done in this derivation

$$\left[\begin{array}{c}\delta v_{N1}\\ \delta p_{N1}\\ \delta \theta _{N1}\end{array}\right]=\left[\begin{array}{ccc}V_{5N4}& W_{5N4}& X_{5N4}\\ V_{5N3}& W_{5N3}& X_{5N3}\\ V_{5N2}& W_{5N2}& X_{5N2}\end{array}\right]\left[\begin{array}{c}\delta p_N\\ \delta \theta _N\\ 1\end{array}\right]$$
(24)
we can compute all components of vectors $`U,V,W,X`$, and then all the corrections. Finally, at the surface we have total luminosity, effective temperature and the MTR for each model.

As stated above, in this work we assume that mass transfer occurs if stellar radius $`R`$ overflows the effective radius $`R_L`$ of the Roche lobe, calculated with the formula of Eggleton (1983),

$$\frac{R_L}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+\mathrm{ln}(1+q^{1/3})}$$
(25)
where $`a`$ is the orbital separation and $`qM/M_2`$ the mass ratio of the binary components.

When a star losses mass, either by a stellar wind in the case of isolated star or by transfer in CBS evolution, not only the external layers, but all stellar structure is affected by this phenomenon. An outer layer that was internal, after a certain time interval, it will be in the stellar surface. Therefore, we need to consider in the code this important effect. We have found it enough to simulate this outflow by means of a linear interpolation in the abundances per mass fractions. Then we should check the meshpoint distribution for consistency.

It is important to stress here that in doing so, we are introducing noise in the chemical profile of the donor star. Presumably, this is responsible (at least in part) for the noise in the MTR vs. time relations shown below for the specific systems we have computed. This is especially noticeable in the case of donor stars with convective envelopes. Evidently, a better strategy should be to employ a moving grid tailored to follow the outwards motion of the stellar layers. We plan to incorporate this refinement in the near future.

## 4 The Orbital Evolution

As a first approximation, we can perform conservative binary evolution calculations. In this case we consider total mass and orbital angular moment as constants. However, we expect the occurrence of total mass and angular momentum losses to be very important in some astrophysical situations of interest. Thus, we have included in our code, these phenomena.

From the definition of total angular orbital moment, and using Kepler’s third law, we can write

$$\frac{\dot{a}}{a}=2\frac{\dot{J}}{J}+\dot{M}\left(\frac{1}{M+M_2}\frac{2}{M}\right)+\dot{M}_2\left(\frac{1}{M+M_2}\frac{2}{M_2}\right)$$
(26)
where $`\dot{J}/J`$ represents energetic losses due to different processes, $`M`$ and $`M_2`$ the masses, and $`\dot{M},\dot{M}_2`$ the MTRs of the lossing - accreting stars.

If we consider non - conservative mass transfer, this causes an angular loss from the system. We follow the method developed for Rappaport, Joss & Webbink (1982) and Rappaport, Joss & Verbunt (1983). It is specified by two free parameters, the fraction $`\beta `$, of the mass lost by the primary star that is accreted by the secondary star $`2`$ and the specific angular momentum, $`\alpha `$, of matter lost from the system in units of $`2\pi a^2/P`$, so that

$$\delta J=\alpha \delta M(1\beta )\frac{2\pi a^2}{P}$$
(27)
where $`\delta M`$ is an incremental mass lost by the donor star, $`\delta J`$ the incremental angular momentum of the matter lost by the system, and $`a`$ and $`P`$ are the orbital semi major axis and period, respectively. We assume that the orbit is always well approximated by a circle of radius $`a`$ (where $`a`$ is a function of time). If we consider only angular momentum losses due to the ejection of matter from the system, the last equation can be rewritten as the differential equation

$$\delta J=\alpha (1\beta )\sqrt{G\left(M+M_2\right)a}\delta M.$$
(28)
Using Kepler’s third law to eliminate $`P`$, and combined with an expression for the total systemic angular momentum,

$$J=\sqrt{\frac{Ga}{M+M_2}}MM_2.$$
(29)
Eq. (28) constitutes a differential equation for $`J`$ (or $`a`$) as a function of $`M`$ (we have neglected the rotational angular momentum of the components stars). Then, we write the loss of angular momentum for the matter ejection as

$$\frac{dlnJ_{ME}}{dt}=\frac{\alpha (1\beta )\sqrt{\left(M+M_2\right)Ga}}{J}\dot{M}.$$
(30)
Angular momentum loss due gravitational radiation is calculated according to the standard formula (Landau & Lifshitz 1959)

$$\frac{d\mathrm{ln}J_{GR}}{dt}=\frac{32G^3\mu }{5c^5}\frac{M_{tot}^2}{a^4}$$
(31)
where $`M_{tot}=M+M_2`$ and $`\mu =MM_2/(M+M_2)`$, $`G`$ and $`c`$ are the gravitational constant and vacuum speed of light, respectively.

To calculate the angular momentum loss due to magnetic braking, we use the prescription of Rappaport, Joss & Verbunt (1983), which is based on the magnetic - braking law of Verbunt & Zwaan (1981),

$$\frac{dJ_{MB}}{dt}=3.8\times 10^{30}MR^4\omega ^3dyncm.$$
(32)
where $`\omega `$ is the angular rotation frequency of the donor star, assumed to be synchronized whit the orbit. We include full magnetic braking when the star has a sizeable convective envelope embracing a mass fraction $`0.02`$. If we substitute Eqs. (30) - (32) in (26), and consider in this last equation

$$\dot{M}_2=\beta \dot{M},$$
(33)
in view of the definition of $`\beta `$, we obtain a differential equation for the orbital separation, which has no analytical solution.

Let us remark that the semiaxis of the orbit is dependent upon the MTR and so, the equivalent radius of the Roche lobe, which in turn, has been assumed to be the radius of the star. Thus, during mass transfer episodes, to be consistent with the iterated value of MTR, we need to perform an orbital integration for each iteration. As the equations for the orbital evolution are well behaved, we decided to handle them with a standard Runge - Kutta technique (Press, et al. 1992). This is very fast and represents a tiny increment of the total numerical effort.

## 5 Handling Outer Layers Integrations

Here we shall treat the problem of handling the outer layers integrations adequately for the case in which there is mass transfer. If not, outer layers are managed as in Kippenhahn, et al. (1967). We refer the reader to that paper for further details in this simple situation and only describe here its modifications to account for the occurrence of mass transfer.

In the outer layers, we integrate the equations of structure neglecting temporal derivatives and the velocity but taking the acceleration of these layers into account. Then, we have a system of ordinary differential equations. In particular, the equation for the luminosity we consider is (see Eq 13)<sup>1</sup><sup>1</sup>1However, we warn the reader that other authors have employed different procedures. See especially, Podsiadlowski, et al. (2001).

$$\frac{dl_r}{dm_r}=T\frac{d\xi }{dt}|_{m_r}\frac{S}{\xi }|_t.$$
(34)
In the case of constant mass evolution, it simply implies that the luminosity is constant in the outer layers. However, in the case of mass transfer/loss episodes, it modifies the profile of luminosity in a sizeable way. For example, for radiative envelopes, which have an entropy increase outwards, this equation predicts a drop in the luminosity. This is an important effect to be considered (see the previous sections). If we take a little amount of mass $`M_{OLI}`$ in the outer layers integrations (e.g. $`M_{OLI}/M10^6`$), it is present in the outer layers integrations as well as in the finite - differenced portion of the star.

In order to perform an adequate treatment of the outer layers of the mass losing component of the pair we need to be careful. Because of the way we have planned our iterative scheme, we need to compute the boundary condition equations and its derivatives with respect to the values of the dependent variables at the first meshpoint, the luminosity and effective temperature of the star and also with respect to the MTR $`\dot{M}`$.

In view of the need for values of the derivatives with respect to the MTR, we have found it very convenient to generalize the triangles method presented in Kippenhahn, et al. (1967) in the following way. Let us consider a MTR $`\dot{M}^{(1)}`$ and a timestep $`\mathrm{\Delta }t`$; then, the mass of the star will be $`M^{(1)}=M^{prev}+\dot{M}^{(1)}\mathrm{\Delta }t`$. Now, we construct a triangle in the HR diagram as proposed by Kippenhahn, et al. (1967), i.e., we perform three integrations with

i) $`\mathrm{log}𝙻/L_{\mathrm{}},\mathrm{log}𝚃_{eff},M^{(1)},\dot{M}^{(1)}`$,
ii) $`\mathrm{log}𝙻/L_{\mathrm{}}+\mathrm{\Delta }\mathrm{log}L/L_{\mathrm{}},\mathrm{log}𝚃_{eff},M^{(1)},\dot{M}^{(1)}`$,
iii) $`\mathrm{log}𝙻/L_{\mathrm{}},\mathrm{log}𝚃_{eff}+\mathrm{\Delta }\mathrm{log}T_{eff},M^{(1)},\dot{M}^{(1)}`$,

where $`\mathrm{log}𝙻/L_{\mathrm{}}`$, $`\mathrm{log}𝚃_{eff}`$ corresponds to one vertex of the triangle in the HR diagram and $`\mathrm{\Delta }\mathrm{log}L/L_{\mathrm{}}`$, $`\mathrm{\Delta }\mathrm{log}T_{eff}`$ are fixed values. We also need to perform integrations with another MTR value $`\dot{M}^{(2)}=\dot{M}^{(1)}+\mathrm{\Delta }\dot{M}`$ so that, for the same timestep, the mass of the star is now $`M^{(2)}=M^{prev}+\dot{M}^{(2)}\mathrm{\Delta }t`$ and the triangle in the HR diagram is

iv) $`\mathrm{log}𝙻/L_{\mathrm{}},\mathrm{log}𝚃_{eff},M^{(2)},\dot{M}^{(2)}`$,
v) $`\mathrm{log}𝙻/L_{\mathrm{}}+\mathrm{\Delta }\mathrm{log}L/L_{\mathrm{}},\mathrm{log}𝚃_{eff},M^{(2)},\dot{M}^{(2)}`$,
vi) $`\mathrm{log}𝙻/L_{\mathrm{}},\mathrm{log}T_{eff}+\mathrm{\Delta }\mathrm{log}𝚃_{eff},M^{(2)},\dot{M}^{(2)}`$.

As the model is iterated, $`\dot{M}`$ changes. We require the MTR to be in between the values at which the envelopes were computed, i.e., $`\dot{M}^{(1)}\dot{M}\dot{M}^{(2)}`$. If not, we change the values of $`\dot{M}^{(1)}`$, $`\dot{M}^{(2)}`$. Also, if $`(\mathrm{log}𝙻/L_{\mathrm{}},\mathrm{log}𝚃_{eff})`$, were outside the assumed HR triangle we change it with the algorithm presented in Kippenhahn, et al. (1967) but now with the same mass and MTR values.

Now, we perform a linear interpolation in $`\dot{M}`$ for the values of the dependent variables at the bottom of the envelope. Notice that in this way we are automatically interpolating for the correct value of the mass of the star. Also, we compute the derivatives of the dependent variables with respect to the MTR at each vertex of the triangle. Then, all the relevant quantities are found by means of a linear bidimensional interpolation inside the triangle in the HR as in Kippenhahn, et al. (1967).

The experience we have gained with our code indicates that a very convenient transformation for the MTR variable is

$$\sigma \mathrm{ln}\left(\frac{\dot{M}}{1M_{\mathrm{}}/y}\right).$$
(35)
Notice that, in this way, the MTR has its correct sign every time. The envelopes are computed for given values $`\sigma _1`$ and $`\sigma _2=\sigma _1+\mathrm{\Delta }\sigma `$. For most purposes we have found it enough to set $`\mathrm{\Delta }\sigma 0.5`$, $`\mathrm{\Delta }\mathrm{log}L/L_{\mathrm{}}=0.02`$, and $`\mathrm{\Delta }\mathrm{log}T_{eff}=0.005`$. This represents a good compromise between the amount of required envelope integrations and the precision in the interpolations and allows the whole iteration to converge in few steps.

We compute a set of six outer envelope integration in the case that the chemical composition of the outermost portion of the star has changed. Otherwise, if the values of $`\sigma `$, $`\mathrm{log}L/L_{\mathrm{}}`$ $`\mathrm{log}T_{eff}`$ fall outside the prism corresponding to the outer boundary integrations employed for the previous model, we change the necessary vertex in order to minimize the number of outer integrations. This is important because they are time consuming.

Let us describe how to compute the moment of the onset of mass transfer. We assume that we have begun our sequence of models in conditions at which the radius of the star is smaller than the corresponding to the sphere equivalent to the Roche lobe, given by Eq. (25). At each moment we compute the structure of the star together with the orbital evolution. Notice that the computation of the orbit is necessary in view that there may be angular momentum dissipation mechanisms operating even without mass transfer (e.g. gravitational radiation). If the components of the pair are close enough each other, at some moment the radius of the most massive component of the pair<sup>2</sup><sup>2</sup>2Notice, however, that there are situations in which, e.g., a canonical 1.4 $`\mathrm{M}_{}`$ neutron star has, e.g., a lower mass main sequence companion. In this case, obviously, the fastest evolving object will be the low mass one. will overflow the Roche radius. At this stage we discard the evolutionary model. Then, we guess a MTR low enough and try to converge the full algorithm presented above. If the MTR value guessed is too high and/or the time step is too long we shall exchange too much mass. This usually makes the iteration divergent, but even if it converges it should be discarded. Then, we halve the timestep, and perform a constant mass integration sequence up to Roche lobe overflow. Usually this procedure is enough to get a plausible initial mass lossing model. From then on, the code works in a very similar fashion compared to the case of constant mass sequences. This procedure for getting the structure of the first model in the mass transfer stage is certainly more straightforward than the one described by Ziolkowski (1970).

A more difficult point is to find a self consistent end of the mass transfer stage. Usually, we set a low MTR threshold, much smaller than the value we guessed for the beginning of the mass transfer stage. As the star gets near the end of the mass exchange epoch, the MTR goes down, sometimes very suddenly. Meanwhile the MTR is above the referred threshold, we assume that the next model will be still loosing mass, otherwise we perform a constant mass evolutionary step and then we compare the radius of the star with that of the Roche lobe. If the star is smaller, we assume that the mass exchange stage has finished and from there on, we continue the evolutionary calculation with constant mass. Otherwise, we discard the mass constant integration and perform a mass exchanging integration with a halved timestep. We have found it to work in most situations.

## 6 Other Relevant Numerical Details

Another critical, and certainly non trivial problem is the rezoning. This is unavoidably necessary during the course of the evolution. Here, rezoning is more problematic compared to the case of the evolution of isolated stars. Notice that, the stellar radius represents a global characteristic of a star that, in the case of isolated stars it appears only in the boundary condition relating effective temperature with luminosity. However, in the case of binary evolution during mass transfer episodes, it also appears in the condition that equates the radius of star with the radius of the Roche lobe. This equation, in turn, largely determines the MTR, a quantity that deeply affects the structure of the whole star.

At rezoning we unavoidably introduce numerical noise that affects the radius of the star. Thus, it is not surprising that, in some critical stages of mass transfer, rezoning spoils the convergence of the code. In most of these cases we recovered convergence simply performing no rezoning. Even in well behaved cases, the noise in the radius of the star produces large oscillations in the value of the MTR. This is especially true in the case of donor stars with convective envelopes.

In the present version of the code we have simply considered the rezoning criteria described in Kippenhahn, et al. (1967): We simply test if the variation of the independent and dependent variables in a zone is larger than a given threshold, if so, we add a new meshpoint dividing the old zone in two. At the new meshpoint we compute the value of the functions by means of linear interpolations. Conversely, if the variation of all variables are below another prefixed value, we remove a meshpoint merging two zones in one. As it will be shown below, this procedure reveals as adequate for the purpose of computing the evolution of low mass CBS that give rise to the formation of helium white dwarfs and also for the case of high mass CBSs. However, in some cases, we have been not able to follow the evolution in very fast stages, like a born again phenomenon of a $`0.8`$ $`\mathrm{M}_{}`$ star formed from a 5 $`\mathrm{M}_{}`$ object in a 3 days period system. We presume such problems to be due to the rezoning.

Notice that the new meshpoints introduced by linear interpolations are not solution of the structure equations and participate in the calculation in the temporal derivatives of pressure, temperature, radius and velocity. We expect to improve significantly the present rezoning scheme by introducing a Hermite interpolation for the values of the dependent variables at inserting new meshpoints as suggested by Wagenhuber & Weiss (1994). Hermite interpolation in a zone can be performed employing a cubic expression whose coefficients are computed using the values of the function and their derivatives at the edges of the zone. (see Wagenhuber & Weiss 1994 for more details). We plan to perform these changes in the near future.

## 7 Application of the Code to Some Low Mass Close Binary Systems

Before describing the calculations we have performed as a test of our new binary evolutionary code, we have to describe the physical ingredients we have included in the present version. Because of convenience in the process of preparing the code, we have chosen to include some very simplified ingredients, mainly in the nuclear description of the problem. We have included the equation of state, OPAL opacities, conductive and molecular opacities, and neutrino emission processes as in Althaus, Serenelli & Benvenuto (2001). We have employed a very simple nuclear reaction network in which we have considered equilibrium expressions for proton - proton, CNO cycles and helium burning cycles given in Clayton (1968). We have followed explicitly the abundances of $`{}_{}{}^{1}H`$, $`{}_{}{}^{4}He`$, $`{}_{}{}^{12}C`$, and $`{}_{}{}^{16}O`$ employing the standard Arnett & Truran (1969) implicit matrix algorithm. We have neglected diffusion processes.

In order to test our numerical scheme we have chosen to compute the evolution of three low mass, CBSs that undergo non - conservative, Class A and B mass transfer episodes from the main sequence to a highly evolved white dwarf configuration.

In all cases we have assumed $`\alpha =0.5`$ and $`\beta =0`$, i.e. that the mass lost by the primary star leaves the system (no accretion onto the secondary) and half of the specific angular momentum is gone away carried by such material. In particular we have considered

1) $`M=2.0`$ $`\mathrm{M}_{}`$ , $`q_0=1.5`$, $`P_0=1.5`$ days
2) $`M=2.0`$ $`\mathrm{M}_{}`$ , $`q_0=1.5`$, $`P_0=0.7`$ days
3) $`M=1.4`$ $`\mathrm{M}_{}`$ , $`q_0=1.0`$, $`P_0=1.0`$ days.

Here, $`q_0M/M_2`$ is the initial the mass ratio of the binary components and $`P_0`$ the initial orbital period. The evolutionary track for the primary component of system 1 is shown in Fig. 1 and the conditions at the beginning (points labelled with odd numbers) and end (points labelled with even numbers) of mass transfer episodes are included in Table 1. The star fills its Roche lobe soon after the end of core hydrogen burning. From then on, the star begins to transfer mass on a timescale of $`10^7`$ y (see Fig. 2). The star ends the first mass transfer episode with a much lower mass value (see Table 1) and begins to evolve bluewards to the white dwarf stage. This evolution is stopped, because the hydrogen envelope is still thick enough to ignite in a flash fashion, forcing the star to inflate on a very short timescale. Then, the star swells very suddenly and fills the corresponding Roche lobe again. Then, the star suffers from a second mass transfer episode depicted in the second panel of Fig. 2. In this second episode the total amount of mass transferred is very tiny and occurs on a very short timescale. As a consequence the Roche lobe size is almost constant and the star is forced to evolve along a constant - radius line. Now, after the end of the second mass transfer episode, the star has a thinner outer hydrogen layer and evolves again bluewards to the white dwarf stage but again it suffers from a second hydrogen thermonuclear flash. Then, the history is repeated: the star suddenly swells up and fills the Roche lobe by third time. Mass transfer begins and again a very little amount of mass is transferred but now in a much shorter timescale compared to the case of the immediately previous mass transfer episode. After ending this third mass transfer episode, the star evolved again bluewards to the white dwarf stage and cools down without suffering any other flash.

As it is well known, hydrogen thermonuclear flashes occur very near the stellar surface. As the hydrogen rich outer layer is very thin, if we want to compute all these stages properly, we do need to consider a very fine zoning for these critical layers. As previously mentioned, at flash - induced mass transfer episodes, a very tiny fraction of the star is lost before the star detaches and starts to evolve bluewards again. Then, we do also need a very thin zoning in order to account properly for the exact amount of hydrogen remaining in the star, as well as the large gravitational energy release when the star undergo large radius excursions. This is the kind of physical situation that suggested us the convenience of designing our binary evolution code in the way we have presented in the previous sections, especially regarding the treatment of the outer layers and the amount of matter we allow to be above the outermost finite - differenced zone.

We remark that, assuming this fully implicit scheme including the MTR we have found that the most difficult quantity to compute is just the MTR. This is especially true when the outer layers become in convective equilibrium. This is reflected in some strong oscillations shown in the first panel of Fig. 2. In any case we feel this not to affect the main course of the evolution. Remarkably, such oscillations are not present in the flash - induced mass transfer episodes.

System 2 is a Class A system (see Fig. 4) which transfers most of its mass when it is still burning hydrogen in the core. Then, it suffers three thermonuclear flashes with its corresponding mass transfer episodes in a way very similar to the corresponding to System 1 (see Table and Figs. 5 and 6).

System 3 (Class B) suffers from four mass transfer episodes, three of them induced by each thermonuclear flash. The main characteristics of these evolutionary calculations are included in Figs. 7 - 9 and in Tables 2 and 3 we include the conditions at the onset and end of mass transfer episodes.

System 3 evolution can be compared with the results presented by Podsiadlowski et al. (2001) in their Fig. 14. They evolved a system composed initially by a 1.4 $`M_{}`$ normal star and a neutron star of the same mass. Mass transfer begins near the end of core hydrogen burning. In doing so, they employed a full nuclear reaction network, convective overshooting, and different assumptions about mass transfer from the ones we have done here. In particular, they used $`\beta =0.5`$, i.e., that half of the mass lost by the normal star is accreted by the neutron star as long as the rate is smaller than the Eddington limit ($`10^8M_{}/y`$). Also, they considered $`\alpha =1`$ for the angular momentum losses. Clearly, these are large differences compared with our System 3 and also with the physical ingredients we considered in the present paper. Notice that Podsiadlowski et al. (2001) found three hydrogen thermonuclear flashes, two of which leading to Roche lobe overflow. In our system, we also have three hydrogen flashes, all of them leading to lobe overflow. This is an important difference that we suspect it to be related with the different assumptions in the orbital evolution of the system. We also consider that the differences in the treatment of nuclear energy release are of key importance. Clearly, our equilibrium cycles predict stronger flashes. However, notice that the mass of the final helium white dwarf is very similar in both calculation (0.199 $`M_{}`$ for Podsiadlowski et al. 2001).

Finally, in Fig. 10 we show the evolution of the central part of the three computed objects in the temperature - density plane. As expected, the evolution is rather similar and is limited to temperatures below helium ignition. The flash episodes suffered by these objects are responsible for the oscillation in the temperature of the star before it gets near the final density value. Then, the objects cool down almost at constant central density, as it should be expected.

These calculations indicate that the above described scheme is capable of handling the evolution of close binary systems in a very adequate way. In our opinion, the here described code should be a valuable tool in handling binary evolution. Let us remark that, we have tested the code with low mass systems because are the ones we are most interested in. Also, we have found very good convergence of the method for the case of massive CBSs.

In closing, we should remark that the above described results should not be considered a state - of - the - art evolutionary calculations but a serious test of the numerical scheme as a previous step for future works. This is so mainly because of the lack of a detailed nuclear reactions network and also for the absence of a treatment of diffusion. Notice that diffusion is a key process in the evolution of objects in the range of masses we have here considered and even largely determine the total number of thermonuclear flashes the object suffers from (see Althaus, et al. 2001).

## 8 Discussion and Future Applications

In this paper we have presented a Henyey-type code tailored for computing stellar evolution in close binary systems (CBSs). This code is a modification of the scheme presented long ago by Kippenhahn, Weigert & Hofmeister (1967) for the case of the evolution of isolated objects.

Here we considered that mass transfer occurs when the star overflows its Roche lobe. The main characteristic of the present code is that it automatically computes the beginning and end of mass transfer stages, computing the MTR of the object in a fully implicit, self consistent fashion. Perhaps, main shortcoming of the scheme we have presented is that it is not able to compute common envelope evolution, simply because at these stages it is not possible to equate the radius of the star to that of the Roche lobe.

In the present version of the code, we have considered radiative OPAL and conductive opacities, neutrino emissivities, convection following the mixing length prescription, and a detailed equation of state. The nuclear reactions have been considered only by means of the equilibrium cycles. Regarding the physics of mass transfer, we have considered the possibility that a fraction of the mass transferred from the primary is accreted by the secondary. Also we have included angular momentum dissipation due to gravitational radiation and magnetic braking.

We have found that the scheme detailed above is able to handle most of the evolutionary stages suffered by a donor star belonging to a CBS. Apart from the results we have presented as a first application, we have also tested it for the case of massive CBSs finding it to be very adequate also for these purposes. Generally speaking, the code is able to perform single stellar evolution and also to find the onset and end of each mass transfer episode. In all the test calculations we have performed we found it much more easy to compute a mass transfer rates (MTRs) when outer layers are in radiative equilibrium. In the case of convective layers convergence is attained but sometimes it is necessary a larger number of iterations.

As a first application of this code to astrophysically relevant conditions, we have applied it to the case of the formation of low mass helium white dwarfs. These computations should not be considered as a state - of - the - art ones but as a test of the numerical scheme as a previous step for future works. This is so mainly because of the lack of a detailed nuclear reactions network and also for the absence of a treatment of diffusion.

We expect that the physical and numerical performance of the code should be largely improved if we replace the outer boundary condition for the stellar radius by the formulation of Ritter (1988) for mass loss transfer (see the Introduction). This should be especially important at the beginning and end of each mass transfer episode and/or in the case of very low mass transfer rates. Also, the implementation of a moving grid should help in producing smooth curves of mass loss vs. time.

As a future application for the present code we plan to study the problem of the formation of helium white dwarfs including the complete physical ingredients we have considered in the previous works of our group on this topic. Also, we plan to study the relevance of the irradiation (Tout, et al. 1989; D’Antona & Ergma 1993) in the evolution of the members of CBSs. This should be especialy relevant for the case of LMXBs which in some evolutionary stages are expected to undergo mass transfer driven by irradiation (Podsiadlowski 1991).

Fig. 4 - 10 are available upon request to the authors at their e-mail addresses.

## 9 Acknowledgements

The authors want to warmly acknowledge our referee Prof. Phillip Podsiadlowski for his very constructive criticism and advice. This has enabled us not only to improve the original version of this work significantly but also has helped us on how to upgrade the performance of the code we presented here from a physical and numerical point of view.

Also, OGB wants to acknowledge Prof. Francesca D’Antona for discussion about the problem of computing binary evolution.