text:
# Stability domains for time-delay feedback control with latency

## I Introduction

Throughout the last decade the use of time delayed signals for controlling unstable periodic orbits (UPOs) has been a field of increasing interest. A method first introduced by Pyragas PYR92 , known as “time-delay autosynchronization” (TDAS), calculates the control force from the difference of the current state to the state one period in the past. Socolar et al.SOC94  have shown that this technique can be improved by using states further in the past. This generalization of TDAS is called “extended time-delay autosynchronisation” (ETDAS). One great advantage of ETDAS over conventional feedback controller schemes is that it can be applied to high frequency oscillators. Since it employs a direct comparison of continuous signals generated by the system itself the only factors limiting the speed of the controller is the bandwidth of the amplifiers and signal propagation times. There is no need to generate a reference signal independently.

Experiments on electronic oscillators have shown that ETDAS can be effective SUK97 . They also show, however, that the latency time associated with signal propagation – the time required to compare the the current signal with its time delayed counterpart and inject the feedback into the system – can have important effects. Just JUS99b  has shown how longer latency times decrease the range of feedback gains over which control is achieved for simple systems controlled by TDAS. Here we extend his analytic formalism, which consists of a first-order perturbation theory in the gain, to ETDAS and note some novel features of the behavior of the individual Floquet multipliers. We then show how the theory provides a qualitative explanation of the shape of the domain of control observed in the experiments of Sukow et al.SUK97 .

## II System equations and perturbation theory

In this section we review the formalism developed by Just, generalizing to the case of ETDAS controllers. We use the same notation as Bleich and Socolar BLE96  for the system equations and feedback signal, and the notation of Just JUS99b  for the Floquet theory parameters.

Consider a $`N`$-dimensional dynamical system defined by

$$\dot{𝐱}(t)=𝐟(𝐱(t),ϵ_0+ϵ(t)),$$
(1)
where $`𝐱(t)`$ denotes a $`N`$-dimensional state vector, $`ϵ_0`$ is a parameter that can be modulated to achieve control, and $`ϵ`$ is a feedback signal. Assume that in the absence of control, $`ϵ(t)=0`$, there exists an unstable periodic orbit $`𝐱_0(t)`$ with period $`\tau `$. Let $`\xi (t)=\widehat{𝐧}𝐱(t)`$ be a component of $`𝐱`$ that can be continuously measured. The ETDAS feedback signal can be written in several equivalent forms:

$`ϵ(t)`$ $`=`$ $`\gamma {\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}R^k\left[\xi (tt_lk\tau )\xi (tt_l(k+1)\tau )\right]`$ (2)
$`=`$ $`\gamma \left(\xi (tt_l)(1R){\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}R^{k1}\xi (tt_lk\tau )\right)`$
$`=`$ $`\gamma \left(\xi (tt_l)\xi (tt_l\tau )\right)+Rϵ(t\tau ),`$
where $`\gamma `$ (the feedback gain) and $`R(1,1)`$ are real parameters, and $`t_l`$ is the latency time. The control force vanishes if the UPO is stabilized, since $`\xi (tk\tau )=\xi (t)`$ for all $`t`$ and any integer $`k`$. TDAS corresponds to the special case $`R=0`$. The last form is the basis for simple implementation of ETDAS in experiments.

In order to determine whether the controlled orbit is stable, Equation (1) is linearized around the UPO. Defining $`𝐲(t)=𝐱(t)𝐱_0(t)`$, we have

$`\dot{𝐲}(t)`$ $`=`$ $`𝐉(t)𝐲(t)+ϵ(t){\displaystyle \frac{𝐟}{ϵ}}|_{𝐱_0(t),ϵ=0}`$
$`=`$ $`𝐉(t)𝐲(t)+\gamma 𝐌(t)\left(𝐲(tt_l)(1R){\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}R^{k1}𝐲(tt_lk\tau )\right),`$
where $`𝐉(t)=𝐉(t)|_{𝐱_0(t),ϵ=0}`$ denotes the Jacobian matrix of the uncontrolled system, and $`𝐌(t)=\frac{𝐟}{ϵ}|_{𝐱_0(t),ϵ=0}\widehat{𝐧}`$ is a $`N\times N`$ matrix containing all information about the control force.

Since $`𝐉(t)`$ and $`𝐌(t)`$ are both periodic with period $`\tau `$, Floquet theory ensures that $`𝐲(t)`$ can be written as

$$𝐲(t)=\underset{m=0}{\overset{\mathrm{}}{}}\underset{n=1}{\overset{N}{}}c_m^{(n)}e^{(\mathrm{\Lambda }_m^{(n)}+i\mathrm{\Omega }_m^{(n)})t}𝐩_m^{(n)}(t),$$
(4)
where $`𝐩_m^{(n)}(t)`$ is a periodic function with period $`\tau `$:

$$𝐩_m^{(n)}(t)=𝐩_m^{(n)}(t+\tau ).$$
(5)
The factorization of the sum into a double sum is done for convenience in the discussions below. In the absence of control, i.e. the absence of time delay terms, there are $`N`$ eigenmodes of the system, indexed by $`n`$. When control is turned on, each of these gives rise to a countably infinite set of eigenmodes indexed by $`m`$. For each set, there is one eigenvalue that begins at the original value $`\lambda ^{(n)}+i\omega ^{(n)}`$ for $`\gamma =0`$ and varies as $`\gamma `$ is increased. The remaining members all have eigenvalues that approach either $`\mathrm{log}|R|/\tau `$ or $`\mathrm{}`$ as $`\gamma `$ approaches $`0`$.

$`\mathrm{\Lambda }_m^{(n)}`$ and $`\mathrm{\Omega }_m^{(n)}`$ are the real and imaginary parts of the Floquet exponent corresponding to the eigenmode $`𝐩_m^{(n)}(t)`$. Inserting Equation (4) into Equation (II) will lead to conditions that must be satisfied by $`\mathrm{\Lambda }_m^{(n)}`$ and $`\mathrm{\Omega }_m^{(n)}`$. The system is linearly stable if and only if all $`\mathrm{\Lambda }_m^{(n)}`$ that satisfy these conditions are negative. Equations (II) and (4) yield the following equation for each of the modes $`𝐩_m^{(n)}(t)`$, where we drop the subscript $`m`$ and superscript $`n`$ for notational convenience:

$$(\mathrm{\Lambda }+i\mathrm{\Omega })𝐩(t)+\dot{𝐩}(t)=𝐉(t)𝐩(t)+\gamma 𝐌(t)e^{(\mathrm{\Lambda }+i\mathrm{\Omega })t_l}\frac{1e^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}{1Re^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}𝐩(tt_l).$$
(6)
This equation is equivalent to

$`(\mathrm{\Lambda }+i\mathrm{\Omega })𝐩(t;\kappa )+\dot{𝐩}(t;\kappa )=\left[𝐉(t)+\kappa 𝐌(t)𝐖(t,t_l)\right]𝐩(t;\kappa ),`$ (7)
where $`𝐖(t,\mathrm{\Delta }t)`$ is the propagator defined by $`𝐩(t;\kappa )=𝐖(t,\mathrm{\Delta }t)𝐩(t\mathrm{\Delta }t;\kappa )`$ and

$$\kappa \gamma e^{(\mathrm{\Lambda }+i\mathrm{\Omega })t_l}\frac{1e^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}{1Re^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}.$$
(8)
Since $`\kappa `$ is proportional to $`\gamma `$, it can be loosely thought of as a measure of the strength of the control gain. One must keep in mind, however, that the value of $`\kappa `$ is ultimately determined by the solutions for the exponent $`\mathrm{\Lambda }+i\mathrm{\Omega }`$. Note also that the Floquet eigenmodes $`𝐩(t;\kappa )`$ themselves depend on $`\mathrm{\Lambda }+i\mathrm{\Omega }`$ through $`\kappa `$, making for a nontrivial modification of the usual eigenvalue problem.

An expression taking the effects of control into account can be derived by perturbation theory. Equation (7) can be written as

$$(\mathrm{\Lambda }+i\mathrm{\Omega })𝐩(t;\kappa )=\left(\frac{d}{dt}+𝐉(t)+\kappa 𝐌(t)𝐖(t,t_l)\right)𝐩(t;\kappa )(t).$$
(9)
We regard $`\frac{d}{dt}+𝐉(t)`$ as an operator with known eigenvalues $`\lambda ^{(n)}=i\omega ^{(n)}`$ and consider $`\kappa 𝐌(t)𝐖(t,t_l)`$ to be a perturbation, a technique familiar from quantum mechanics. The effects of the controller on the Floquet exponents can be expanded in powers of $`\kappa `$ as

$$\mathrm{\Lambda }+i\mathrm{\Omega }=\lambda +i\omega +\chi (t_l)\kappa +o(\kappa ^2),$$
(10)
where the coefficient $`\chi (t_l)`$ is a complex valued function. Any effects of interactions between the Floquet modes enter only at $`\kappa ^2`$. This is why it is convenient to index the modes by $`n`$ and $`m`$ for the purposes of first-order perturbation theory.

As mentioned above, Equation (10) has an infinite number of solutions for $`\mathrm{\Lambda }`$ and $`\mathrm{\Omega }`$ which approach $`\mathrm{log}|R|/\tau `$ PYR01 ; BEC02  or minus infinity as the feedback gain $`\gamma `$ goes to zero. This behavior arises from the essential singularity at $`\mathrm{\Lambda }=\mathrm{log}|R|/\tau `$ and the divergence at $`\mathrm{\Lambda }=\mathrm{}`$, the latter arising only due to the nonzero latency time.

Just JUS99b  has shown how the coefficient $`\chi (t_l)`$ can be calculated. Let $`𝐮(t)=𝐩(t;\kappa =0)`$ and $`𝐯^{}(t)`$ be the right and left Floquet eigenmode in the absence of control, $`\kappa =0`$. Both are periodic with period $`\tau `$, i.e. $`𝐮(t)=𝐮(t+\tau )`$ and $`𝐯^{}(t)=𝐯^{}(t+\tau )`$. Now let $`\widehat{𝐮}(t)=\mathrm{exp}(i\omega t)𝐮(t)`$ and $`\widehat{𝐯}^{}(t)=\mathrm{exp}(i\omega t)𝐯^{}(t)`$. These satisfy the equations

$`\lambda \widehat{𝐮}(t)+\dot{\widehat{𝐮}}(t)`$ $`=`$ $`𝐉(t)\widehat{𝐮}(t)`$
$`\lambda \widehat{𝐯}^{}(t)\dot{\widehat{𝐯}}^{}(t)`$ $`=`$ $`\widehat{𝐯}^{}(t)𝐉(t),`$ (11)
with the boundary conditons:

$`\widehat{𝐮}(t+\tau )`$ $`=`$ $`e^{i\omega \tau }\widehat{𝐮}(t)`$
$`\widehat{𝐯}^{}(t+\tau )`$ $`=`$ $`e^{i\omega \tau }\widehat{𝐯}^{}(t).`$ (12)
The standard first-order perturbation theory result for the coefficient $`\chi (t_l)`$ is

$`\chi (t_l)=e^{i\omega t_l}\rho (t_l),`$ (13)
where

$$\rho (t_l)=\frac{_0^\tau \widehat{𝐯}^{}(t)𝐌(t)𝐖(t,t_l)\widehat{𝐮}(t)𝑑t}{_0^\tau \widehat{𝐯}^{}(t)\widehat{𝐮}(t)𝑑t}.$$
(14)
Now $`\rho (t_l)`$ depends on $`t_l`$ only through $`𝐖`$, and since $`𝐮(t)`$ is $`\tau `$-periodic, $`𝐖(t,t_l+\tau )𝐮(t)=𝐖(t,t_l)𝐮(t)`$. Therefore, $`\rho (t_l)`$ has to satisfy

$$\rho (t_l+\tau )=e^{i\omega \tau }\rho (t_l).$$
(15)
Inserting Equation (13) in Equation (10) and neglecting second order terms yields

$`\mathrm{\Lambda }+i\mathrm{\Omega }=\lambda +i\omega +\rho (t_l)\gamma e^{(\mathrm{\Lambda }+i(\mathrm{\Omega }\omega ))t_l}{\displaystyle \frac{1e^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}{1Re^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}}.`$ (16)
Following the treatment of Just JUS99b  for the $`R=0`$ case, we note that this expression can be simplified in the case of a so-called flip orbit, where $`\omega =\pi /\tau `$. Defining the frequency deviation $`\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Omega }\pi /\tau `$, Equation (16) can be rewritten as

$$\mathrm{\Lambda }+i\mathrm{\Delta }\mathrm{\Omega }=\lambda +\rho (t_l)\gamma e^{(\mathrm{\Lambda }+i\mathrm{\Delta }\mathrm{\Omega })t_l}\frac{1+e^{(\mathrm{\Lambda }+i\mathrm{\Delta }\mathrm{\Omega })\tau }}{1+Re^{(\mathrm{\Lambda }+i\mathrm{\Delta }\mathrm{\Omega })\tau }}.$$
(17)
Moreover, since all coefficients in Equation (II) are real, and for $`\omega =\pi /\tau `$ the boundary conditions of Equation (II) are invariant under complex conjugation, $`\widehat{𝐮}(t)`$ and $`\widehat{𝐯}^{}(t)`$ can always be chosen to be real valued. From Equation (14) it is then clear that $`\rho (t_l)`$ is real and from Equation (15) we have $`\rho (t_l)=\rho (t_l+\tau )`$.

We emphasize that $`\widehat{𝐮}(t)`$, $`\widehat{𝐯}^{}(t)`$, and $`\rho (t_l)`$ are also real valued in the case of zero torsion, i.e. $`\omega =0`$. In this case, $`\chi (t_l)=\rho (t_l)`$ and $`\rho `$ is $`\tau `$-periodic: $`\rho (t_l)=\rho (t_l+\tau )`$. For torsion-free perturbations, Equation (16) becomes

$`\mathrm{\Lambda }+i\mathrm{\Omega }=\lambda +\chi (t_l)\gamma e^{(\mathrm{\Lambda }+i\mathrm{\Omega })t_l}{\displaystyle \frac{1e^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}{1Re^{(\mathrm{\Lambda }+i\mathrm{\Omega })\tau }}}.`$ (18)
Although a nonzero torsion of all unstable eigenmodes is a necessary condition for possible control, JUS97  the case of $`\lambda <0`$ and $`\omega =0`$ might be interesting because an initially stable eigenmode can become unstable when the control force is applied and thus limit the domain of control.

## III Floquet exponents for fixed latency time

To understand the solutions to Equation (17), it is helpful to first consider their behavior for fixed values of $`g\rho (t_l)\gamma `$. A detailed discussion of the effects of variations of $`\rho `$ with $`t_l`$ will be presented in section IV.

The behavior of the real and imaginary parts of the Floquet exponent in Equation (17) can be seen in Figure 1. In each panel, curves are shown for both $`t_l=0`$ (thick lines) and $`t_l=\tau /2`$ (thin lines). For $`g=0`$ (no control), the real part $`\mathrm{\Lambda }`$ is equal to $`\lambda `$. For increasing $`g`$ the value of $`\mathrm{\Lambda }`$ decreases, reaching $`0`$ at $`g=\lambda (1+R)/2`$, then changing its sign; thus the orbit becomes stable. Further increase of $`g`$ usually leads to a collision with an exponent created by the control scheme, forming a complex conjugate pair (see Figure 1(b), (d), and (f), except for the thick line in (f)). After the collision, $`\mathrm{\Lambda }`$ then begins to increase (see Figure 1(a), (c), and (e), except for the thick line in (e)). For $`g`$ sufficiently large, $`\mathrm{\Lambda }`$ becomes positive again and control is lost.

Note that as $`g`$ increases from zero, an infinite number of solutions to Equation (17) emerge from $`\mathrm{log}|R|/\tau `$ as complex conjugate pairs. PYR01 ; BEC02  In order to collide with the single exponent coming from $`(\lambda +i\omega )\tau `$, one pair has to become real and separate (see Figure 1(d)). If this does not happen, a crossing of branches can occur. After the crossing, the complex conjugate pair becomes the branch with largest $`\mathrm{\Lambda }`$ and thus responsible for the stability of the system (see thick lines in Figure 1(e) and 1(f)).

Figure 1 illustrates that increasing the latency time $`t_l`$ and/or decreasing $`R`$ leads to a smaller range of $`g`$ for successful control. In fact, one can compute a maximum latency time $`t_{max}`$ for which control can be achieved, which corresponds to the case where the collision of the two branches occurs at $`\mathrm{\Lambda }=0`$ as in the thin line in Figure 1(a). Using the imaginary part of Equation (17) to eliminate the factor of $`g`$ in the real part of that same equation and setting $`\mathrm{\Lambda }`$ equal to zero, one has to search for nontrivial solutions for $`\mathrm{\Delta }\mathrm{\Omega }`$. A condition for the existence of such a solution is

$$t_lt_{max}=\frac{1}{\lambda }+\frac{\tau }{2}\frac{R1}{R+1}.$$
(19)
This agrees with the result of JUS99b  for $`R=0`$. For larger latency time, control is not possible to first-order in $`\kappa `$. Equation (19) shows that negative $`R`$ decreases the maximum latency time. We therefore focus on $`R0`$ from here on.

Another way to visualize the effects of the control scheme is to consider how the Floquet multipliers (as opposed to exponents) evolve with varying $`g`$. We consider again the case $`\omega =\pi /\tau `$. The Floquet multipliers are defined by $`\mu _m^{(n)}=e^{(\mathrm{\Lambda }_m^{(n)}+i\mathrm{\Delta }\mathrm{\Omega }_m^{(n)})\tau }`$. Dropping the subscript $`m`$ and the superscript $`n`$, Equation (17) can be rewritten:

$$\mu =\mathrm{exp}\left[\lambda \tau g\tau \left(\frac{\mu +1}{\mu +R}\right)\mu ^{t_l/\tau }\right].$$
(20)
Stability is achieved if all multipliers satisfy $`|\mu |<1`$; i.e. the multipliers $`\mu `$ are located inside the unit circle in the complex plane.

Considering the real and imaginary part of Equation (20) and numerically following the roots of each equation in the complex plane (using Mathematica), we observe topologically different cases depending on the signs of $`\lambda `$ and $`g`$ in the absence of latency, $`t_l=0`$, as illustrated in Figure 2.

* $`g>0;\lambda >0`$: (Figure 2(a) and (b)) As $`g`$ increases, the largest Floquet multiplier starts outside the unit circle at $`\mathrm{exp}(\lambda \tau )`$ on the real axis indicated by point A, moves towards the unit circle, and eventually crosses it. Meanwhile an infinite number of complex conjugate pairs spread out from $`\mu =R`$ indicated by point B one of which will determine the stability range. Two scenarios are possible. One possibility is that one pair recombines on the real axis and one of the multipliers collides with the system’s while the other approaches zero. After the collision they form a complex conjugate pair and cross the unit circle again (see Figure 2(a)). The other possibilitiy is that the created pair becomes the largest multiplier, turns around and crosses the unit circle (see Figure 2(b)). The first case corresponds to the thick lines in Figure 1(c) and (d), the second to Figure 1(e) and (f).
* $`g>0;\lambda <0`$: (Figure 2(c) and (d)) Similar to the previous case the largest Floquet multiplier starts at $`\mathrm{exp}(\lambda \tau )`$ on the real axis, this time inside the unit circle. It moves towards the origin for increasing $`g`$ and may either collide with a multiplier created by the control force as in Figure 2(c), or continue towards zero while a complex conjugate pair becomes the largest multiplier as in Figure 2(d). In both cases all multipliers stay inside the unit circle for increasing $`g`$; for $`\lambda <0`$ the system is stable for all $`g>0`$.
* $`g<0;\lambda >0`$: The largest multiplier starts at $`\mathrm{exp}(\lambda \tau )`$ on the real axis and goes to $`\mathrm{}`$, thus control is never successful. All other multipliers created by the control scheme stay inside the unit circle for decreasing $`g`$.
* $`g<0;\lambda <0`$: The largest multiplier starts inside the unit circle, crosses it, and goes to $`\mathrm{}`$. More mulitpliers created by the control scheme cross the unit circle for further decrease of $`g`$.

## IV Shapes of the domain of control in the $`t_l`$-$`\gamma `$ plane

For a discussion of the domain of control in the $`t_l`$-$`\gamma `$ plane let us consider first $`\rho (t_l)=1`$. We will show later how the coefficient $`\rho (t_l)`$ scales the domain at every value of the latency time.

For each $`\lambda `$, the lower branch in Figure 3 is the horizontal line $`\gamma =\lambda (1+R)/2`$, where there is a flip instability associated with the real exponent originating from $`\lambda `$ at $`\gamma =0`$. The upper branch corresponds to a Hopf bifurcation that can arise in two different ways: (1) the relevant complex conjugate pair of exponents originates in a collision between the branch associated with $`\lambda `$ and a real eigenvalue created by the feedback scheme, as in Figure 1(a) or (c); or (2) the relevant complex conjugate pair originates at $`\mathrm{log}|R|/\tau `$, as in the thick lines in Figure 1(d). It can be seen that increasing $`R`$ and decreasing $`\lambda `$ increases the domain of control at fixed $`t_l`$ and increases $`t_{max}`$ (see Equation (19)). The upper and lower curves do not intersect if the system is already stable, i.e. $`\lambda <0`$.

The actual domain of control is a distortion of Figure 3 owing to the variation of $`\rho `$ with $`t_l`$. The distortion is simple to compute, however, since changing the value of $`\rho `$ is entirely equivalent to changing $`\gamma `$. Thus the values of $`\gamma `$ on the upper and lower curves at a particular value of $`t_l`$ are simply multiplied by $`1/\rho (t_l)`$. Note that these variations in $`\rho `$ cannot change the ratio of the upper and lower values of $`\gamma `$.

For eigenmodes with $`\omega =\pi /\tau `$, the antiperiodicity of $`\rho (t_l)`$ and the fact that $`\rho (t_l)`$ is real require that $`\rho (t_l)`$ has at least one root in the interval $`[0,\tau ]`$. We assume for convenience a sinusoidal form $`\rho (t_l)=A\mathrm{sin}(\frac{\pi }{\tau }t_l\varphi )`$. Figure 4 shows the domain of control for different values of the phase $`\varphi `$. In each panel, the shading indicates the stable domain for one choice of $`\varphi `$.

If a root of $`\rho (t_l)`$ appears before the upper and lower branches intersect, the domain of control will include a region with negative feedback gain $`\gamma `$, as shown in Figure 4(a). Since the domain of control at fixed $`t_l`$ scales like $`1/\rho `$, divergences appear at the values of the latency time for which $`\rho (t_l)`$ vanishes. For $`\lambda <0`$, shown in Figure 4(b), the upper and lower branches still do not intersect.

The case $`\lambda <0`$ can be important because it can reduce the domain of control, as shown in Figure 5(c) and (d), when the uncontrolled system has exponents $`\lambda _1>0`$ and $`\lambda _2<0`$. Since second order terms are neglected in the perturbation theory, the different Floquet modes of the system do not interact with each other, so the effective domain of control is just the intersection of the single domains for each exponent.

Domains similar to those shown in Figure 5(c) and (d) have been observed in experiments on high speed diode resonator circuits. The analogous figure obtained from experiments is reproduced here as Figure 5(a) and (b) to facilitate comparison. To construct the theoretical figure, parameter values are adjusted to reproduce several features of the experimental results. From the experiments, three parameters are known: the weighting parameter $`R`$, the largest Lyapunov exponent $`\lambda _1`$, and, since the instability is a flip, the frequency $`\omega _1=\pi /\tau `$. Equation (19) for the $`\lambda _1`$ mode then gives an immediate prediction for $`t_{max}`$, the largest latency time for which control can be achieved. The agreement with the experiment is reasonable, especially given that the very narrow tails of the domains may be hard to detect in experiments.

We make the plausible assumption that the second largest Floquet mode is a stable flip, so $`\omega _2=\pi /\tau `$. In our simple model, there are then five parameters that determine the shapes of the domains of control: $`\lambda _2\tau `$ (the real part of the subleading exponent), $`A_1`$ and $`A_2`$ (the amplitudes of the variation in $`\rho _1`$ and $`\rho _2`$), and the phases $`\varphi _1`$ and $`\varphi _2`$ (which determine where $`\rho _1`$ and $`\rho _2`$ vanish).

To fix these five parameters, we consider the $`R=0.68`$ domain. The phase $`\varphi _1`$ must lie somewhere between $`t_1\pi `$ and $`t_2\pi `$ in order for the divergence in $`1/\rho `$ and associated sign change in the domain of control to be right. From the fact that the onset of divergence is not evident yet at $`t_1`$, we estimate that $`\varphi _1`$ is close to $`t_2\pi `$ and fix it at $`0.6\pi `$. From the fact that the boundary of the subleading mode does not appear to cut off the domain near $`t_l=0`$, we take the divergence of $`\rho _2`$ to occur there, requiring $`\varphi _2=0`$ or $`\pi `$. The remaining parameters, $`A_1`$, $`A_2`$, and $`\lambda _2\tau `$ are adjusted to fit $`t_1`$, $`t_2`$, the latency times (in units of $`\tau `$) corresponding to the limits of the domains, and $`\gamma (t_2)`$, the gain at which the lower domain of control is cut off at $`t_2`$. The parameters determined from the $`R=0.68`$ data are used for the $`R=0`$ plot as well since $`\rho (t_l)`$ is determined purely from the uncontrolled system.

The primary conclusion we draw is that the theory does give qualitative insight into the structure of the stability domains. Even with our crude constraints on the form of $`\rho (t_l)`$, the general shapes of the domains are reproduced surprisingly well. As expected, larger $`R`$ increases the size of the domain of control, especially the part with negative $`\gamma `$.

It may appear that by adjusting the full functions $`\rho _1`$ and $`\rho _2`$ one could fit arbitrary shapes of the stability domains, and it is true that many types of undulations in the domain boundaries could be fit. Moreover, one can always appeal to lack of robustness to noise to explain why very narrow regions of a predicted stability domain would not show up in the experimental data. There are, however, some features of the experimental data that cannot be reproduced by the theory presented here. In particular, consider the width of the region at negative $`\gamma `$ for $`R=0.68`$. The upper and lower boundaries of the right half of this region both represent instabilities in the same mode, the mode associated with $`\lambda _1`$. For a fixed value of $`t_l`$, then, the ratio of the $`\gamma `$’s at these boundaries is independent of $`\rho `$. The substantially larger width of the experimental stability domain cannot be obtained by adjusting any of the parameters in our theory. Including additional modes or making a different assumption about $`\omega _2`$ would not help, since the narrowness of the theoretical domain is determined by $`\lambda _1`$. We therefore conclude that second-order effects are significant in the experimental system. These effects may involve interactions between modes associated with different $`\lambda `$’s or just interactions of modes within the set generated by $`\lambda _1`$ alone.

Our analysis can also be used to explore the possibilities for qualitatively different domain shapes. An interesting example is obtained when $`\varphi _2`$, the phase associated with the subleading mode of the uncontrolled system is taken to be $`0.1\pi `$ (and $`A_2<0`$). As shown in Figure 6, this can lead to a situation in which a nonzero latency time is required for effective control.

## V Conclusion

We have discussed the effects of latency time on a feedback control scheme known as ETDAS. Using Floquet theory and carrying out a first-order perturbation theory in the feedback gain, we have shown that nontrivial domain shapes can arise in the plane parameterized by feeback gain and latency time. Within the first order theory, we find that no control is possible above a maximum latency time determined solely by the Floquet exponent of the most unstable mode in the uncontrolled system. We also find that Floquet modes that are stable in the uncontrolled system contribute significantly to the overall stability picture, reducing the domain of control substantially.

The theory accounts well for qualitative features of the stability domains observed in experiments. As expected, larger values of the ETDAS parameter $`R`$ give larger stability domains. Detailed comparison indicates, however, that second order effects are experimentally observable.

###### Acknowledgements.

We would like to thank E. Schöll, W. Just, A. Amann, I. Harrington, and D. Gauthier for useful conversations. This work was supported by NSF Grant PHY-98-70028 and within the framework of the exchange program between TU Berlin and Duke University. P.H. acknowledges a Fulbright Scholarship.


format: text

text:
# Quantum corrections to microscopic diffusion constants

## 1 Introduction

Basic physical considerations suggest that, in addition to convective mixing – routinely included in stellar evolution computations – additional transport processes are efficient within the stellar interior; they are driven by pressure, temperature and chemical abundance gradients, and by the effect of radiative pressure on the individual ions. These processes are collectively called ’diffusion processes’, and their inclusion in stellar evolution computations is necessary in order to satisfy the helioseismological constraint for the solar models (Bahcall et al. 1995).

In general, individual ions are forced to move under the influence of pressure as well as temperature gradients, which both tend to move the heavier elements toward the centre of the star, and of concentration gradients that oppose the above processes. Radiation, which causes negligible diffusion in the Sun, pushes the ions toward the surface, whenever the radiative acceleration of an individual ion species is larger than the gravitational acceleration. The speed of the diffusive flow depends on the collisions with the surrounding particles, as they share the acquired momentum in a random way. It is the extent of these ’collision’ effects that dictates the timescale of element diffusion within the stellar structure, once the physical and chemical profiles are specified. The most general treatment for the element transport in a multicomponent fluid associated with diffusion is provided by Burgers’ (1969, B69) equations. In these equations the effect of collisions between ions is represented by the so-called resistance coefficients, i.e. the matrices $`K`$, $`z`$, $`z^{}`$, $`z^{\prime \prime }`$, whose precise evaluation is fundamental in order to estimate correctly the diffusion timescales for the various elements111In the original B69 book the elements of the matrix $`K`$ are called resistance coefficients, whereas $`z`$, $`z^{}`$, $`z^{\prime \prime }`$, which are related to the efficiency of the thermal diffusion, have no specific names. In this paper we follow the notation of Paquette et al. (1986), who termed both $`K`$ and $`z`$, $`z^{}`$ and $`z^{\prime \prime }`$ as resistance coefficients..

Recent spectroscopic determinations of Fe and Li abundances in Galactic globular cluster turn-off stars discussed in, e.g., Gratton et al. (2001) or Bonifacio et al. (2002), have shown that the present standard treatment of diffusion is in disagreement with the observed surface abundance of these two elements in stars of globular clusters (for the Li problem see, e.g., Michaud et al. 1984 or Vauclair & Charbonnel 1995). In the light of these results, it is very important to conclusively assess how much of this discrepancy is due to an incorrect treatment of the diffusion process, and how much is due to competing rotationally induced or other non-standard macroscopic mixing phenomena, which inhibit the efficiency of diffusion.

In this paper the state of the art regarding the computation of the resistance coefficients in conditions typical of the stellar plasma is reviewed and a detailed comparison of the results from different authors is performed. The effect of different resistance coefficients on solar models is examined, too. Moreover, we introduce and discuss for the first time in astrophysical computations the effect of quantum corrections to the evaluation of the resistance coefficients. We also provide simple fitting formulae for accurate calculations of resistance coefficients including the appropriate quantum corrections. In § 2 we compare existing determinations of the resistance coefficients, and discuss the differences on solar models; in § 3 we determine the appropriate quantum corrections to these coefficients, and conclusions will follow in § 4. Analytical formulae for the computations of updated resistance coefficients, including quantum corrections, are given in the appendices222FORTRAN77 routines to compute the improved resistance coefficients (including quantum corrections) and their implementation into Schlattl’s (2002) diffusion-constant routine are publicly accessible under http://www.astro.livjm.ac.uk/$`\stackrel{~}{}`$hs..

## 2 Classical treatment of element diffusion

The diffusion of elements in a multicomponent fluid can be treated either according to the Chapman & Cowling (1970) or the B69 formalism. We are following the latter description, which is equivalent to the so-called “second approximation” in Chapman & Cowling (1970). It would be desirable to use a higher order approximation, because an uncertainty of the order of 10% is introduced by using only the second one (Roussel-Dupré 1982), independent of the accuracy in the resistance coefficients. But this can only be done in the scheme of Chapman & Cowling (1970), which becomes very cumbersome for gases with more than two components.

### 2.1 Existing calculations of resistance coefficients

Burgers’ equations are obtained assuming the gas particles have approximate Maxwellian velocity distributions; the temperatures are the same for all particle species; the mean thermal velocities are much larger than the diffusion velocities; magnetic fields are unimportant. Burgers’ scheme involves the resistance coefficients $`K`$, $`z`$, $`z^{}`$ and $`z^{\prime \prime }`$; following Mason et al. (1967, MMS67), they can be expressed in terms of the so-called reduced collision integrals $`\mathrm{\Omega }_{ij}^{(l,s)}{}_{}{}^{}`$ (Hirschfelder et al. 1954)333Note the error in the table at page 44 of B69: Hirschfelder et al.’s (1954) $`8\mu \mathrm{\Omega }_{st}^{(ij)}`$ is equivalent to $`\mathrm{\Sigma }_{st}^{(ij)}`$ in B69 and Chapman & Cowling’s (1970) $`8\mu \mathrm{\Omega }_{st}^{(i)}(j)`$., according to the following relationships:

$`{\displaystyle \frac{K_{ij}}{K_{ij}^0}}`$ $`=`$ $`4{\displaystyle \frac{T_{ij}^{}{}_{}{}^{2}\mathrm{\Omega }_{ij}^{(1,1)}{}_{}{}^{}}{\mathrm{ln}\left(\mathrm{\Lambda }_{ij}^2+1\right)}},`$ (1)
$`z_{ij}`$ $`=`$ $`11.2{\displaystyle \frac{\mathrm{\Omega }_{ij}^{(1,2)}{}_{}{}^{}}{\mathrm{\Omega }_{ij}^{(1,1)}{}_{}{}^{}}},`$ (2)
$`z_{ij}^{}`$ $`=`$ $`2.5{\displaystyle \frac{6\mathrm{\Omega }_{ij}^{(1,2)}{}_{}{}^{}4.8\mathrm{\Omega }_{ij}^{(1,3)}{}_{}{}^{}}{\mathrm{\Omega }_{ij}^{(1,1)}{}_{}{}^{}}}\text{ and}`$ (3)
$`z_{ij}^{\prime \prime }`$ $`=`$ $`2{\displaystyle \frac{\mathrm{\Omega }_{ij}^{(2,2)}{}_{}{}^{}}{\mathrm{\Omega }_{ij}^{(1,1)}{}_{}{}^{}}},`$ (4)
where

$`T_{ij}^{}`$ $`=`$ $`k_\mathrm{B}T{\displaystyle \frac{\lambda _\mathrm{D}}{\left|Z_iZ_je^2\right|}}\text{and}`$ (5)
$`K_{ij}^0`$ $`=`$ $`{\displaystyle \frac{2}{3}}\sqrt{{\displaystyle \frac{2\mu _{ij}\pi }{(k_\mathrm{B}T)^3}}}\left(Z_iZ_je^2\right)^2n_in_j\mathrm{ln}\left(\mathrm{\Lambda }_{ij}^2+1\right)`$ (6)
with $`\lambda _\mathrm{D}=\sqrt{k_\mathrm{B}T/(4\pi n_\mathrm{e}e^2)}`$ being the Debye-length, $`\mathrm{\Lambda }_{ij}=4T_{ij}^{}`$ the plasma parameter, $`\mu _{ij}=m_im_j/(m_i+m_j)`$ the reduced mass, and $`Z_i`$, $`m_i`$ and $`n_i`$ the charge number, mass and particle number density of species $`i`$, respectively. The collisions between the particles of the stellar plasma determine the values of the $`\mathrm{\Omega }_{ij}^{(l,s)}{}_{}{}^{}`$ integrals; the physics of the collisions is specified by some form of the Coulomb interaction which, as a first approximation, can be described by a pure Coulomb potential with a long-range cut-off distance, usually assumed to be $`\lambda _\mathrm{D}`$. According to B69, using this truncated pure Coulomb potential the resistance coefficients become

$`{\displaystyle \frac{K_{ij}}{K_{ij}^0}}`$ $`=`$ $`2{\displaystyle \frac{\mathrm{ln}\mathrm{\Lambda }_{ij}C_\mathrm{E}\pm \frac{\pi ^2}{4}}{\mathrm{ln}\left(\mathrm{\Lambda }_{ij}^2+1\right)}},`$ (7)
$`z_{ij}`$ $`=`$ $`0.6,`$ (8)
$`z_{ij}^{}`$ $`=`$ $`1.3\text{ and}`$ (9)
$`z_{ij}^{\prime \prime }`$ $``$ $`2,`$ (10)
where $`C_\mathrm{E}`$ is Euler’s constant and the alternative signs in Eq. (7) represent repulsive ($`+`$) and attractive ($``$) particles, respectively. For high densities and thus small values of the plasma parameter, the resistance coefficients become negative (for an attractive potential already for $`\mathrm{\Lambda }20`$). Since in stellar plasmas, in particular in case of white dwarfs, $`\mathrm{\Lambda }`$ is often much smaller than this value, more elaborate computations have been performed.

More accurate collision integrals have been obtained by considering a Debye-Hückel type of potential

$$V_{ij}(r)=\frac{Z_iZ_je^2}{r}\mathrm{e}^{r/\lambda _\mathrm{D}}$$
(11)
where $`r`$ is the particle distance. The results are provided either in tabulated form by MMS67 or as fitting formulae by Muchmore (1984, M84), Iben & MacDonald (1985, IM85) and Paquette et al. (1986, P86). In the latter cases the integrals are given as functions of

$$\mathrm{\Phi }_{ij}=\mathrm{ln}\left(\mathrm{ln}(1+\mathrm{\Lambda }_{ij}^2)\right).$$
(12)
Different notations for the same physical quantities have been used by these groups; in Table 1 their relations to the present work are summarized. M84, IM85 and P86 argue that while $`\lambda _\mathrm{D}`$ is an appropriate screening distance at low densities, the Debye sphere loses its significance at high densities, in which case a more appropriate screening distance is the mean interionic distance; thus they suggest to use in the actual stellar model computations the larger of $`\lambda _\mathrm{D}`$ or the mean interionic distance. Whatever one chooses for the actual screening distance, its value has to be employed as $`\lambda _\mathrm{D}`$ in Eq. (5) to obtain the appropriate value for $`\mathrm{\Lambda }_{ij}`$ for determining the collision integrals. We also remark that, at least for the entire evolution of the Sun up to now, $`\lambda _\mathrm{D}`$ has always been larger than the mean interionic distance.

The results of the different groups are compared in Fig. 1. Noticeably, given that they use all the same physical assumptions, they all got very similar values, but they all differ significantly from Burgers’ results at higher densities (i.e., lower values of $`\mathrm{\Phi }`$), where his approximations are not adequate anymore. M84’s formulae for the case of a repulsive potential have been obtained by fitting solely the values represented by dots in Fig. 1. Although the formulae work pretty well overall, $`z^{\prime \prime }`$ deviates considerably at $`\mathrm{\Phi }=2`$ from the mean value of $`2.6`$, and the almost linear behaviour of $`z^{}`$ for $`3<\mathrm{\Phi }<1`$ is only poorly reproduced.

The supposedly most accurate calculation of the collision integrals has been performed in P86, where cubic splines for 50 equally spaced intervals in $`\mathrm{\Phi }`$ are provided; the results agree very well with the tabulation by MMS67. To obtain more manageable but still accurate formulae, we fitted polynomials of at most 5th order to their values for $`2.5<\mathrm{\Phi }<3`$; this enables to compute diffusion of elements up to Fe in stars on the main sequence and to follow with sufficient accuracy the early white dwarf cooling phase. The numerical values for the polynomial coefficients are given in Appendix A.

In Fig. 2 the relative differences between our fitting formulae and the values of P86 are shown. For the case of a repulsive potential ($`e`$-$`e`$ or ion-ion collisions) the results of P86 for all quantities can be reproduced with an accuracy better than 2%; for an attractive potential ($`e`$-ion collisions) the accuracy is still well within 5%, and for a large range of $`\mathrm{\Phi }`$ values it is certainly better than 3%. Taking into account that the results of different groups for the same physical assumptions, e.g., between IM85 and P86 (long-dashed line in Fig. 2a), differ by up to 15%, we consider our polynomial fits to be sufficiently accurate.

In main sequence stars the value of $`\mathrm{\Phi }`$ does not deviate considerably from $`2`$ for H and He (see Fig. 3b), thus constant values for $`z`$, $`z^{}`$ and $`z^{\prime \prime }`$ are usually assumed. Moreover, in stellar plasmas the amount of element diffusion is basically determined by ion-ion collisions, and as the differences for attractive and repulsive potential are small when $`\mathrm{\Phi }2`$, the resistance coefficients for a repulsive potential are often employed for all cases. For instance, the widely used diffusion routine by (Thoul et al. 1994, TBL) uses the values of B69 for $`z`$, $`z^{}`$ and $`z^{\prime \prime }`$ (Eq. 810) while for $`K`$ the fitting formula of IM85 has been implemented. This assumption overestimates $`z`$ and $`z^{}`$, but underestimates $`z^{\prime \prime }`$, which results in a somewhat too high efficiency of the thermal diffusion (see M84).

### 2.2 Effect on solar models

In order to estimate the influence of different choices of the resistance coefficients on the solar structure, we computed various solar models utilizing the stellar evolution code and element diffusion routine described in Schlattl (2002) and references therein. Element diffusion is treated following the scheme by TBL, considering in addition the effect of electron degeneracy and partial ionization of elements; the latter is implemented by defining a mean charged ion per element, instead of computing diffusion for each ion separately. This treatment of partial ionization is sufficient for a large range of stellar masses, including the Sun. We neglect radiative levitation, which leads to a tiny improvement of the theoretical sound-speed profile compared to the Sun of at most 0.025% at $`r0.6R_{}`$ (see Turcotte et al. 1998).

The sound speeds of various solar models computed using different choices of the collision integrals are summarized in Fig. 3a, while the variation of He and metal abundances at the surface and at the centre of the Sun are shown in Fig. 4. The solid line represents a solar model (S1, see Table 2) computed with $`K`$ of IM85 and constant values for the $`z`$’s of B69 corresponding to the values used by TBL. This description has been used, for instance, in solar models of Bahcall et al. (1998) and Schlattl (2002). Dropping the assumption of constant $`z`$’s, and using instead the functions of M84, the sound-speed difference to a seismic model becomes about 4 times higher for $`0.3<r/R_{}<0.65`$ (dotted line in Fig. 3a; S2), caused by the reduced thermal diffusion efficiency. The latter leads also to a diminished surface depletion of He and metals (see Fig. 4). When one employs, in addition, M84’s values for $`K`$ instead of the IM85 ones, the difference to model S1 (dash-dot-dot-dotted line) lowers, because M84 computed an about 7% smaller value for $`K`$ than IM85 (Fig. 1a) for the relevant $`\mathrm{\Phi }`$ range, which results in slightly higher diffusion velocities.

In solar conditions, M84’s value of 2.62 for $`z^{\prime \prime }`$ is however too high with respect to the exact result; in fact $`\mathrm{\Phi }`$ is about 2.2 in solar radiative regions (see Fig. 3b), and $`z^{\prime \prime }`$ has to be between 2.3 and 2.4 (Fig. 1d). Applying the correct $`z^{\prime \prime }(\mathrm{\Phi })`$ would further reduce the sound-speed difference compared to models computed with the value of M84. Indeed, using our formulae for P86’s collision integrals, which gives similar values of $`K`$ as M84 (Fig. 2a), but includes the variation of $`z^{\prime \prime }`$, lowers the sound-speed difference (compare S3 and S4 in Fig. 3a). A small additional increase of the diffusion velocities is obtained by considering different collision integrals for $`e`$-ion (attractive) and ion-ion or $`e`$-$`e`$ (repulsive) interactions (see also Fig. 4), which leads to model S5 (long-dashed line in Fig. 3a).

Typical additional quantities inferred by helioseismological methods are the depth of the convective zone and the surface helium abundance, which are determined to be $`0.713\pm 0.001`$ (Basu & Antia 1997) and $`0.246\pm 0.002`$ (Basu & Antia 1995), respectively. In all models of Table 2 apart from S2 these quantities are well within the observational limits. Nonetheless, a general trend following the overall sound speed behaviour can be observed: The larger the sound-speed difference of the respective model the shallower the convective envelope and the higher the surface helium content.

## 3 Quantum corrections

All computations of the resistance coefficient discussed in the previous section have assumed that the collisions are dominated by the classical interaction between two point-charge particles. For long-range potentials like the Debye-Hückel one, quantum effects become important at high energies, whereas the behaviour is classical at low energies. Hahn et al. (1971, H71) have computed corrections to the classical collision integrals, such that

$`\mathrm{\Omega }_{}^{(p,q)}{}_{}{}^{}`$ $`=`$ $`\mathrm{\Omega }_{\mathrm{class}}^{(p,q)}{}_{}{}^{}(T^{})`$ (13)
$`\left(\mathrm{\Omega }_{\mathrm{qu}}^{(p,q)}{}_{}{}^{}(T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2})\pm {\displaystyle \frac{1}{2s+1}}\mathrm{\Omega }_{\mathrm{ex}}^{(p,q)}{}_{}{}^{}(T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2})\right),`$
where $`\mathrm{\Omega }_{\mathrm{class}}^{(p,q)}{}_{}{}^{}`$ is the classical value of §2, $`\mathrm{\Omega }_{\mathrm{qu}}^{(p,q)}{}_{}{}^{}`$ contains the quantum mechanical correction, and $`\mathrm{\Omega }_{\mathrm{ex}}^{(p,q)}{}_{}{}^{}`$ the exchange contribution, which has to be dropped for distinguishable particles (the indices $`ij`$ denoting different interacting particles have been omitted for the sake of clarity.) The upper sign in the term in parenthesis of Eq. (13) refers to particles obeying Fermi-Dirac statistic, the lower one for particles following Bose-Einstein statistics. The spin of the particle is denoted, as usual, by $`s`$, and $`\mathrm{\Lambda }^{}`$ is the de Boer parameter given by

$$\mathrm{\Lambda }_{ij}^{}{}_{}{}^{2}=\frac{h^2}{2\mu _{ij}Z_iZ_je^2\lambda _\mathrm{D}}.$$

Tabulated values for $`T_{}^{}{}_{}{}^{2}\mathrm{\Omega }_{\mathrm{qu},\mathrm{ex}}^{(p,q)}{}_{}{}^{}`$ are given in H71, which we again fitted by polynomials of 5th order (see Appendix A).

The excellent accuracy of our fitting formulae can be checked in Fig. 5, where the effect of the quantum corrections on the resistance coefficient and the thermal diffusion coefficients in the Sun is shown. The corrections increase toward the centre, and especially for the $`\mathrm{H}^+`$-$`e`$ collisions they considerably alter $`K`$. However, the diffusion velocity of H is basically determined by $`\mathrm{H}^+`$-$`\mathrm{H}^+`$ interactions and in that case the corrections for all quantities are below 1%. The expected weak influence on the solar structure is demonstrated by model S6, where our fitting formulae to the quantum corrections of H71 have been used. Overall, the diffusion velocities for all elements have been increased, more prominently in the centre, where the difference to models S5 (without quantum corrections) is the biggest (Fig. 5). This caused a further reduction of the sound-speed difference of model S6 compared to S5 (Fig. 3a).

It should be mentioned that the spin-dependent term in Eq. (13) would demand to know exactly the number of atoms in each ionization and in each excited state. However, the diffusion velocity of all elements is basically determined by their collision frequency with the most abundant element (H), and thus only for $`\mathrm{H}^+`$-$`\mathrm{H}^+`$ collision this spin-dependent term has to be accounted for (the spin-dependent term appears only in case of indistinguishable particles). Moreover, when quantum corrections become important inside stars most of the elements are already fully ionized. Thus only very tiny errors are introduced when, as in this work, only one mean-charged ion per element is considered.

It is interesting to notice that model S6, which includes the most accurate collision integrals of P86 accounts for the difference between repulsive and attractive interactions, and includes quantum effects, has almost the same sound-speed profile as model S1 (cf. Fig. 3a), which contains the most crude treatment of diffusion among all models considered in this work. Also the depth of the convective envelope and the surface helium content are nearly identical (Table 2). Therefore, just by chance, a very accurate solar model with respect to diffusion is obtained when using the resistance coefficients chosen by TBL.

## 4 Conclusions

In this paper, we have compared resistance coefficients computed by various groups for the calculation of atomic diffusion constants, and we have discussed the inclusion of quantum effects on the otherwise classical determination of these quantities. We provide simple analytical formulae to implement easily the presently most accurate determination of resistance coefficients, i.e., the classical result by P86 plus the effect of quantum corrections determined by H71. (FORTRAN routines for computing updated resistance coefficient and the consequent diffusion constants are publicly available.22footnotemark: 2)

The different classical computations based on a a Debye-Hückel type of potential produce solar sound-speed profiles with significant differences when compared to the current accuracy of helioseismological determinations. Just by chance, our accurate treatment including quantum corrections produces a solar sound-speed profile comparable to the one obtained using the less accurate $`z`$, $`z^{}`$ and $`z^{\prime \prime }`$ coefficients selected by TBL. In this context, we would like to add that with none of the descriptions for the collision integrals – neither with the inclusion of radiative levitation – the bump in the sound-speed difference at $`r0.65R_{}`$ disappears. So, an additional mixing process beyond the formal boundary of the convective zones is still a probable candidate to resolve this discrepancy (see, e.g., Richard et al. 1996).

The introduction of quantum corrections increases the efficiency of diffusion with respect to the classical case, and their effect is more pronounced for higher densities. It will therefore be important to test the effect of our accurate resistance coefficients on models of objects denser than the Sun, like white dwarfs, and main sequence turn-off stars of galactic globular clusters, where also radiative levitation is altering the surface abundance patterns (Richard et al. 2002).

Although the accuracy of the diffusion constants could be improved considerably when using the correct functions for the resistance coefficients, there still remains the limitation that Burgers’ formalism is equivalent only to Chapman & Cowling’s (1970) second approximation. This causes an intrinsic uncertainty when using Burgers’ equation of the order of 10% (Roussel-Dupré 1982), which is difficult to reduce further.

###### Acknowledgements.

H.S. has been supported by a Marie Curie Fellowship of the European Community programme “Human Potential” under contract number HPMF-CT-2000-00951.

## Appendix A

In this appendix we provide analytic formulae for $`T_{}^{}{}_{}{}^{2}\mathrm{\Omega }_{}^{(s,t)}{}_{}{}^{}`$ and those combinations of the indices $`s`$ and $`t`$ which are needed to compute the resistance coefficients $`K`$, $`z`$, $`z^{}`$ and $`z^{\prime \prime }`$ using Eq. (1)(4). The polynomial functions for the classical case have been obtained by a least square-difference fitting of the results of P86.

The best fits for the classical collision integrals could be obtained with

$$T_{}^{}{}_{}{}^{2}\mathrm{\Omega }_{}^{(s,t)}{}_{}{}^{}=\underset{k=0}{\overset{5}{}}c_{st}^\mathrm{A}(k)\mathrm{\Phi }^k$$
(14)
for attractive and

$`T_{}^{}{}_{}{}^{2}\mathrm{\Omega }_{}^{(s,t)}{}_{}{}^{}`$ $`=`$ $`{\displaystyle \underset{k=0}{\overset{3}{}}}c_{st}^\mathrm{R}(k)(\stackrel{~}{\mathrm{\Phi }}_{st})^k\text{with}`$ (15)
$`\stackrel{~}{\mathrm{\Phi }}_{st}`$ $`=`$ $`\mathrm{ln}\left(\mathrm{ln}(1+\mathrm{\Lambda }^{b_{st}})\right)`$ (16)
for repulsive potentials. We dropped the indices $`i,j`$ for the sake of clarity. The numerical values for the coefficients $`c^\mathrm{A}`$ and $`c^\mathrm{R}`$, and the exponents $`b`$ are summarized in Table 3. Notice that it is sufficient to determine $`T_{}^{}{}_{}{}^{2}\mathrm{\Omega }_{}^{(s,t)}{}_{}{}^{}`$, as in the definitions of the resistance coefficients the factor $`T_{}^{}{}_{}{}^{2}`$ always cancels out.

The quantum-mechanical (“qu”) and exchange contribution (“ex”) can both be computed according to the following relationship

$`T_{}^{}{}_{}{}^{2}\mathrm{\Omega }_{\mathrm{qu},\mathrm{ex}}^{(p,q)}{}_{}{}^{}(T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2})`$ $`=`$ $`{\displaystyle \frac{pN_p}{q(q+1)}}I_{\mathrm{qu},\mathrm{ex}}^{(q)}(T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2})\text{with}`$ (17)
$`N_p^1`$ $`=`$ $`1{\displaystyle \frac{1+(1)^p}{2(1+p)}}.`$ (18)
In these formulae the quantum and exchange parts of the collision integrals are expressed as functions of the quantity $`T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2}`$. The best fits to the values tabulated in H71 have been obtained with

$`I_\mathrm{C}^{(q)}(x)`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{4}{}}}d_\mathrm{C}^{(q)}(i)\xi _{\mathrm{C}}^{(q)}{}_{}{}^{i}(x)\text{,\hspace{1em}where}`$ (19)
$`\xi _\mathrm{C}^{(q)}(x)`$ $`=`$ $`\{\begin{array}{ccc}\mathrm{ln}(1+g_\mathrm{C}^{(q)}\times x)\hfill & \text{for}& \mathrm{C}=\text{“qu”},\hfill \\ \mathrm{ln}(1+g_\mathrm{C}^{(q)}/x)\hfill & \text{for}& \mathrm{C}=\text{“ex”}\hfill \end{array}`$ (22)
with $`x=T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2}`$ and the additional constraint that $`\xi _{\mathrm{ex}}^{(q)}(x)=0`$ for $`x0.1`$. The index “C” denotes the appropriate expression for either the quantum or the exchange contribution to the collision integrals. The respective coefficients are provided in Table 4. Note that $`I_{\mathrm{ex}}^{(1)}=0`$ and that Eqs. (19) and (22) are only valid for $`x<\mathrm{1\hspace{0.17em}000}`$, as no values were provided by H71 for larger $`x`$. We suggest not to extrapolate $`T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2}`$ beyond this upper boundary, but to use the values at $`T^{}\mathrm{\Lambda }_{}^{}{}_{}{}^{2}=\mathrm{1\hspace{0.17em}000}`$, which are $`I_{\mathrm{ex}}^{(2)}=0.4913`$, $`I_{\mathrm{ex}}^{(3)}=0.4960`$, $`I_{\mathrm{qu}}^{(1)}=1.5053`$, $`I_{\mathrm{qu}}^{(2)}=1.9581`$ and $`I_{\mathrm{qu}}^{(3)}=2.2010`$.

## Appendix B

Here we describe the modifications to be applied to the routine by TBL, in order to use any preferred formulation of the resistance coefficients.

The first change involves Eqs. (9)–(10) of TBL, where the IM85 definition of $`K`$ is implemented. They have to be replaced with the chosen representation of $`K`$.

The matrix of elements $`Y_{st}`$, introduced at page 830 of TBL, has to be changed to

$`Y_{st}=3y_{st}+z_{st}^{}x_{st}m_t/m_s.`$ (23)
Finally, the following elements of the matrix $`\mathrm{\Delta }_{ij}`$ as defined in Eqs. (33) and (34) of TBL become

$$\mathrm{\Delta }_{ij}=\{\begin{array}{ccc}_{ki}z_{ik}k_{ik}x_{ik}\hfill & \text{for}& j=i+S,\hfill \\ z_{iq}k_{iq}y_{iq}\hfill & \text{for}& j=S+1,\mathrm{},2S\hfill \\ & & \text{and }ji+S,\hfill \end{array}$$
(24)
for $`i=1,\mathrm{},S`$ and

$`\mathrm{\Delta }_{ij}=\{\begin{array}{ccc}2.5_{kj}z_{pk}k_{pk}x_{pk}\hfill & \text{for}& j=p,\hfill \\ 2.5z_{pj}k_{pj}x_{pj}\hfill & \text{for}& j=1,\mathrm{},S\text{ and }jp,\hfill \\ \multicolumn{3}{c}{_{ki}k_{pk}y_{pk}(0.8z_{pk}^{\prime \prime }x_{pk}+Y_{pk})0.4z_{pp}^{\prime \prime }k_{pp}}\\ & \text{for}& j=i,\hfill \\ \multicolumn{3}{c}{k_{pq}y_{pq}x_{pq}(3+z_{pq}^{}0.8z_{pq}^{\prime \prime })}\\ & \text{for}& j=S+1,\mathrm{},2S\hfill \\ & & \text{and }ji,\hfill \end{array}`$ (32)
for $`i=S+1,\mathrm{},2S`$, where $`p=iS`$ and $`q=jS`$. Here, the same symbols have been used as in TBL, to which the reader is referred for their definitions.


text:
# Semi–linear gravitational lens inversion

## 1 Introduction

This paper is concerned with the problem of inversion of a gravitationally–lensed image of an extended source, i.e. a galaxy rather than a star or quasar. This problem is interesting because lensed images of extended sources provide more information than images of point sources, and so potentially are more useful for determining the mass profiles in galaxies and clusters of galaxies. Also, because of the magnification, one can measure structure in the light profile of the source at enhanced resolution. In this paper we show how this problem can be separated in a natural way into linear and non–linear dimensions, and how this provides a number of advantages.

In this introduction we review solutions to the inversion problem and introduce some of the terminology used in the remainder of the paper. In all the solutions described here the mass in the lens is parameterized. Nevertheless the analysis applies equally to a pixelized mass distribution.

When presented with the lensed image of an extended source, the unknowns to solve for are the source light profile and the lens mass profile, and the uncertainties in these quantities. One approach to this problem, suggested by Kayser and Schramm (1988), uses the fact that regions of the source that are multiply–imaged have the same surface–brightness. For a trial mass distribution, the method traces image pixels to the source plane where the counts in different image pixels mapping to the same source pixel are compared. The solution for the mass is obtained by minimizing the dispersion in the image pixel counts for such multiply–imaged source pixels. Kochanek et al. (1989) successfully applied this approach to the inversion of the radio Einstein ring MG1131+0456. The algorithm was refined by Wallington, Kochanek, and Koo (1995) who applied it to the triply–imaged giant arc in the galaxy cluster Cl 0024+1654.

The main shortcoming of this approach is that it does not deal with the image point–spread–function (psf). If psf smearing of the image (either instrumental or atmospheric) is significant, the light profile of the source is not correctly recovered by backward tracing the image, even if the mass distribution is exactly known. To deal with the psf, a forward approach is needed i.e. one chooses a model for the source light profile (parameterized or pixelized), and a model for the mass (parameterized or pixelized), forms the image, convolves it with the psf, and compares it to the actual image, adjusting the source and lens models to minimize a merit function e.g. $`\chi ^2`$.

An argument for choosing to parameterize rather than pixelize the source light profile is that it forces the solution to be smooth. Nevertheless, the source light profile may be complex. This is true, for example, in the cases of MG1131+0456 and Cl 0024+1654, cited above. A large number of parameters might be required to provide a satisfactory description. Without clues to the character of the source it is extremely difficult to select the best parameterization i.e. the one which provides a satisfactory fit with the smallest number of parameters. In the most extreme example Tyson, Kochanski, & dell’Antonio (1998) used 232 parameters to model the source light distribution of the galaxy lensed by the cluster Cl0024+1654.

If the source light profile is complex it is natural to consider pixelizing the source, i.e. the counts in each pixel is a free parameter. This removes the difficulty in finding a good parameterization for the source, and thereby avoids any bias in the fitted mass profile resulting from a poor choice. On the other hand, due to the deconvolution, and because the pixels are independent, the solution can be noisy. It is possible to achieve a smooth pixelized solution by adding a suitable ‘regularizing’ term to the merit function. This term, if minimized on its own, would produce a smooth source light profile. By adding this term to $`\chi ^2`$ the final solution involves a balance between obtaining the best fit to the image (minimizing $`\chi ^2`$), and obtaining a smooth source solution (minimizing the regularizing term). Wallington, Kochanek, and Narayan (1996) apply this approach to the case of the radio Einstein ring MG 1654+134. They use a maximum entropy approach i.e. the regularizing term to be minimized is the negative of the entropy, the negentropy. Labeling the counts in source pixel $`i`$ by $`s_i`$, the source plane negentropy is $`_is_i\mathrm{ln}(s_i)`$, and the merit function they minimize is

$$G=\chi _{im}^2+\lambda \underset{i}{}s_i\mathrm{ln}(s_i)$$
(1)
(here we have followed the notation in Press et al., 2001, §18.7). The purpose of the multiplier $`\lambda `$ is to give more or less weight to the negentropy term.

The inversion proceeds as follows: For a fixed value of $`\lambda `$, the solution is determined by searching through the multi–parameter space for the minimum of the merit function. The number of dimensions of the parameter space to search is the sum of the number of source pixels and the number of parameters used for the mass. This search is most efficiently achieved with a pair of nested cycles. The inner (source) cycle searches for the best source light profile for a fixed mass profile. The outer (mass) cycle adjusts the mass profile. Outside this cycle is a third ($`\lambda `$) cycle where the multiplier is adjusted. Because the negentropy term acts to smooth the source, as $`\lambda `$ increases, the $`\chi _{im}^2`$ term also increases, i.e. the fit becomes worse. The principle for reaching the final solution (e.g. Press et al., 2001, §18.4) is to start with $`\lambda `$ large, then progressively reduce the weight of the regularizing term until the $`\chi _{im}^2`$ becomes satisfactory. In other words the solution has the smoothest source that provides a satisfactory fit to the image. ‘Satisfactory’ is usually interpreted as reaching the criterion for the $`\chi ^2`$ for the image $`\chi _{im}^2=min(\chi _{im}^2)+\sigma (\chi _{im}^2)`$. With three nested cycles, the $`\lambda `$, mass, and source cycles, the routine can be slow.

In this paper we describe a new technique which we suggest simplifies and clarifies the problem in a number of ways. In purely formal terms, the method is very similar to the maximum entropy method of Wallington et al.: Algebraically, we simply replace the negentropy term in the merit function (1) with a linear regularization term. However, the insight we bring is to show that for a fixed mass distribution, the minimization of the merit function is now a linear problem i.e. can be solved by matrix inversion. In other words the source cycle – the innermost of the three minimization cycles – is eliminated. This has major benefits. In the first place the inversion is much quicker, thereby allowing a more thorough search for the best fit mass model. At the same time, the uncertainty of identifying the true minimum has been removed. The method also greatly simplifies calculation of the uncertainties, as we show below. More generally, the formalism clarifies the essence of the problem: The source parameters are linear dimensions and the mass parameters are non–linear dimensions. For this reason we call the method ‘semi–linear’.

At this point it is worth noting that, because of magnification and multiple imaging, the number of constraints to the solution can be much greater than the number of parameters to solve for. In this respect the lens inversion problem differs from many inversion problems encountered in astronomy (for example image deconvolution). We find, as a consequence, that in many circumstances the regularization term can be removed altogether. So the $`\lambda `$ cycle is also eliminated. The merit function is then just $`\chi _{im}^2`$, and this is our starting point for the presentation of the theory. For a fixed mass profile, the pixelized source light distribution that produces the min.$`\chi _{im}^2`$ fit is obtained by linear inversion. The mass profile is then adjusted to find the minimum of these min.$`\chi _{im}^2`$ fits. The advantage, besides speed (only the mass cycle remains), is that the solution is unbiased, since there is no smoothing of the source.

The outline of the remainder of the paper is as follows. In §2 we explain the basic theory, demonstrating that for a fixed mass profile, the problem of obtaining the source light profile by $`\chi ^2`$ minimization in the image plane is a linear one, and obtaining the covariance matrix for the counts in the source pixels. We then extend the basic theory to include a linear regularization term. In §3 we apply the method to a realistic problem, assessing the performance for different psf widths, and different source pixel sizes, with and without linear regularization. In §4 we provide a summary of the main points, together with recommendations for applying the method.

## 2 Theory

In this section we present the theory of semi–linear inversion, firstly without regularization (§2.1), and then with regularization (§2.2). In each sub–section we begin with the case where the mass is fixed, and then treat the general case, minimizing also on the mass parameters.

### 2.1 Semi–linear inversion without regularization

#### 2.1.1 Fixed mass: Eliminating the source cycle

Without any regularizing term, the merit function is simply $`G=\chi _{im}^2`$. The basic problem is to find the counts in the source pixels that, for a given mass distribution, minimize the merit function $`G`$, i.e. give the best fit to the observed image. Pixels in the source plane are labeled $`i=1,I`$. There is no restriction on how the source plane is tessellated. In principle, the pixels could change in both size and shape across the source region, which itself could be of any shape. Pixels in the image plane are labeled $`j=1,J`$. It is assumed in the following that the image pixels include counts from the image of the lensed source only i.e. the images of any lensing galaxies, and the mean sky count, have been subtracted. Also we suppose that the data in each image pixel are independent i.e. are characterized by the counts $`d_j`$, and dispersion $`\sigma _j`$, with no covariance between pixels (appropriate for CCD data).

The inversion proceeds as follows: Choose the mass model parameters, then, for each source pixel $`i`$, in turn, form the image for unit counts (surface brightness) by appropriate ray tracing and convolution with the known point spread function i.e. compute the counts in the $`i`$th image $`f_{ij},j=1,J`$. The problem now is to combine these $`I`$ images with scalings $`s_i,i=1,I`$, to minimize $`G`$. These scalings are the deconvolved intrinsic source surface–brightness distribution.

The problem is of a standard type. The merit function is written

$$G=\chi _{im}^2=\underset{j=1}{\overset{J}{}}\left[\frac{_{i=1}^Is_if_{ij}d_j}{\sigma _j}\right]^2.$$
(2)
Minimizing with respect to each of the source terms we have a set of $`I`$ simultaneous equations of the form

$$\frac{1}{2}\frac{G}{s_i}=0=\underset{j=1}{\overset{J}{}}\left[\frac{f_{ij}_{k=1}^Is_kf_{kj}f_{ij}d_j}{\sigma _j^2}\right]$$
(3)
where the reason for the factor $`\frac{1}{2}`$ will soon become clear. These equations may be written in matrix form

$$\text{F}\text{S}=\text{D}.$$
(4)
Here S is a column matrix of length $`I`$ containing the elements $`s_i`$, to be solved for. F is a symmetric $`I\times I`$ matrix, with elements $`\text{F}_{ik}=_{j=1}^Jf_{ij}f_{kj}/\sigma _j^2`$. Finally D is a column matrix of length $`I`$ containing the elements $`\text{D}_i=_{j=1}^Jf_{ij}d_j/\sigma _j^2`$.

The solution for the counts in the source pixels is then simply obtained by matrix inversion

$$\text{S}=\text{F}^1\text{D}$$
(5)
thus eliminating the source cycle.

The solution for the errors has a particularly simple form. We seek the covariance matrix for the source pixels. Noting that

$$\text{F}_{ik}=\frac{1}{2}\frac{{}_{}{}^{2}G}{s_is_k},$$
(6)
we see that the matrix F is one half times the Hessian matrix of $`\chi _{im}^2`$, which is to say that F is the curvature matrix of the problem (Press et al., 2001, §§15.4, 15.5) – this was the reason for using the factor $`\frac{1}{2}`$ in equation (3). We now show that the matrix $`\text{C}=\text{F}^1`$ is the required covariance matrix of S.

For independent image pixels, the covariance between source pixels $`i`$ and $`k`$ is given by

$$\sigma _{ik}^2=\underset{j=1}{\overset{J}{}}\sigma _j^2\frac{s_i}{d_j}\frac{s_k}{d_j}.$$
(7)
Using equation (5) this becomes

$$\sigma _{ik}^2=\underset{j=1}{\overset{J}{}}\sigma _j^2\underset{l=1}{\overset{I}{}}\text{C}_{il}\frac{f_{lj}}{\sigma _j^2}\underset{m=1}{\overset{I}{}}\text{C}_{km}\frac{f_{mj}}{\sigma _j^2}.$$
(8)
Multiplying this out gives

$$\sigma _{ik}^2=\text{C}_{ik}$$
(9)
as required.

We see that for the case of fixed mass, the covariance matrix for the source pixel counts is computed in the inversion step, without the need for further calculation. We shall refer to this $`I\times I`$ matrix as the source covariance matrix hereafter. Even though it is not the complete solution for the source pixel errors (because the mass parameters have been fixed), the source covariance matrix is extremely useful, for example, in exploring different mass models and pixelizations (§3).

It is worth noting that the semi–linear inversion solution, either with or without regularization, differs in character from the maximum entropy solution. With the semi–linear method the counts in any source pixel are unbounded, so the best–fit value could be negative, since some image pixels contain negative counts (i.e. are below mean sky). With the maximum entropy method all source counts must be positive. The semi–linear method provides the best estimate of the counts in a source pixel, and the solution is satisfactory provided the result is consistent with being positive. If the counts in any source pixel are significantly negative (e.g. $`<4\sigma `$) this would indicate a bad mass model. This possibility of testing the suitability of the mass model with a source–plane statistic can be viewed as an extra positive feature of the semi–linear method.

#### 2.1.2 Mass cycle

The full solution proceeds by searching through the mass–distribution parameter space, at each point minimizing $`\chi _{im}^2`$ by linear inversion, to find the smallest of these min.$`\chi _{im}^2`$ values, the global minimum. Because the number of dimensions of the parameter space for the non–linear search has been greatly reduced, it is now a much simpler problem to locate the true minimum securely.

The solution for the errors is more complicated than above, since we have added in the non–linear mass dimensions. If there are $`L`$ parameters that describe the mass, labeled $`m_l`$, we need to form the $`(I+L)\times (I+L)`$ (symmetric) curvature matrix of the problem. But note that the majority of the terms, the $`I\times I`$ terms $`\frac{1}{2}\frac{{}_{}{}^{2}G}{s_is_k}`$, have already been computed and are the elements of the matrix F at the global minimum. The remaining terms, the $`L`$ rows (and columns) of terms such as $`\frac{1}{2}\frac{{}_{}{}^{2}G}{m_lm_n}`$, and $`\frac{1}{2}\frac{{}_{}{}^{2}G}{m_ls_i}`$, can be filled in by simple measurement of the shape of the $`\chi _{im}^2`$ surface about the global minimum. The covariance matrix for the mass and source parameters is the inverse of this curvature matrix. We shall refer to this $`(I+L)\times (I+L)`$ matrix as the full covariance matrix hereafter.

### 2.2 Semi–linear inversion with regularization

#### 2.2.1 Fixed mass: Eliminating the source cycle

The possibility of replacing the negentropy term in equation (1) by a term (a linear regularization term) which preserves the linearity of the min.$`\chi ^2`$ approach is made evident by the linearity of equation (3) with respect to the source parameters. Clearly we can form a merit function by adding to $`\chi _{im}^2`$ any term $`G_L`$ which is a linear combination of terms $`s_is_k`$

$$G_L=\underset{i,k}{}a_{ik}s_is_k$$
(10)
since the partial differentials of these additional terms will also be linear. One example of a linear regularization term is $`G_L=_{i=1}^Is_i^2`$. The choice of $`G_L`$ is discussed below.

Writing the merit function generally as

$$G=\chi _{im}^2+\lambda G_L$$
(11)
then, following through the same analysis as in §2.1.1, the solution for the counts in the source pixels can be written

$$\text{S}=[\text{F}+\lambda \text{H}]^1\text{D}.$$
(12)
We call the matrix H the regularization matrix. The elements of H are

$$\text{H}_{ik}=\frac{1}{2}\frac{{}_{}{}^{2}G_{L}^{}}{s_is_k}.$$
(13)
For example, if the regularization term is $`G_L=_{i=1}^Is_i^2`$, then we have $`\text{H}=\text{I}`$, the identity matrix.

The form of $`G_L`$ should be chosen to penalize noisy solutions. The choice $`G_L=_{i=1}^Is_i^2`$, termed “zeroth–order” regularization in the literature, is one attempt to achieve this. Other widely–used linear regularization terms include gradient and curvature forms. These three regularization forms correspond, loosely speaking, to the prejudice that the source light profile is, respectively, approximately zero, constant, or planar (see Press et al., 2001, §18.5, for a detailed account of the theory of linear regularization and its implementation). In practice, if $`\lambda `$ is not too large, all three terms serve to smooth the source in a rather similar way, and there is little to distinguish between the solutions.

The gradient and curvature forms consider the differences between counts in neighboring pixels. Until now we have used a one–dimensional numbering scheme for the source pixels. In this case, since we need to take account of the relative spatial locations of pixels in the source plane we use coordinates $`x,y`$. The simplest gradient term is

$$G_L=\underset{x,y}{}[s_{x,y}s_{x+1,y}]^2+\underset{x,y}{}[s_{x,y}s_{x,y+1}]^2.$$
(14)
Another form uses $`[s_{x,y}0.5(s_{x+1,y}+s_{x,y+1})]^2`$. In forming the sum it is necessary to translate the indices $`x,y`$ to the index $`i`$, and then equation (13) is used to form the matrix H. Note that zeroth–order regularization is computationally by far the simplest method, since it does not involve this step of translation of indices.

The simplest curvature form is

$`G_L={\displaystyle \underset{x,y}{}}[s_{x,y}0.5(s_{x1,y}+s_{x+1,y})]^2`$
$`+{\displaystyle \underset{x,y}{}}[s_{x,y}0.5(s_{x,y1}+s_{x,y+1})]^2.`$ (15)
Another form uses $`[s_{x,y}0.25(s_{x1,y}+s_{x+1,y}+s_{x,y1}+s_{x,y+1})]^2`$.

The source covariance matrix for the regularized case is fortunately only slightly more difficult to compute than for the unregularized case. Writing $`\text{R}=[\text{F}+\lambda \text{H}]^1`$, and following the same line of reasoning as in §2.1.1, we obtain the analogous equation to equation (8)

$$\sigma _{ik}^2=\underset{j=1}{\overset{J}{}}\sigma _j^2\underset{l=1}{\overset{I}{}}\text{R}_{il}\frac{f_{lj}}{\sigma _j^2}\underset{m=1}{\overset{I}{}}\text{R}_{km}\frac{f_{mj}}{\sigma _j^2}.$$
(16)
Multiplying out we obtain

$$\sigma _{ik}^2=\text{R}_{ik}\lambda \underset{l=1}{\overset{I}{}}\text{R}_{il}[\text{RH}]_{kl}.$$
(17)
#### 2.2.2 Mass cycle

The procedure for the full solution is the same as for the unregularized case. One searches through the mass–distribution parameter space, at each point minimizing $`G`$ by linear inversion, to find the smallest of these min.$`G`$ values, the global minimum. In the regularized case there is no simple solution for the full covariance matrix however. In the unregularized case, we were able to use the fact that the inverse of the full curvature matrix is the full covariance matrix. But in the regularized case this is no longer true since we are minimizing $`G=\chi _{im}^2+G_L`$. Instead, an alternative approach must be used, for example, a Monte Carlo method which inverts an ensemble of realisations of the image by adding noise to the original image.

## 3 Simulations

In this section we apply the semi–linear inversion method to a realistic test problem. To validate the linear inversion step, we begin with the case of fixed mass. We quantify the effectiveness of the method under variations of the image $`S/N`$, psf width, and source pixel size, for both the unregularized and regularized cases. We then present an analysis of the full problem, allowing the mass parameters to vary. Finally we debate the advantages of the unregularized and regularized approaches, for different practical applications.

### 3.1 Test problem

To make the computations more useful we have based our investigation on a realistic simulation of a deep image of an Einstein ring gravitational lens system observed with the Advanced Camera for Surveys (ACS) aboard HST. The camera has a pixel size of $`0.05\mathrm{}`$. We have used cosmological parameters $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$. The lens is placed at $`z=0.3`$ and the source lies at $`z=3.0`$.

Figure 1 illustrates the test problem. The lens (not shown) is modelled as a singular isothermal ellipsoid with one–dimensional velocity dispersion $`260`$ km s<sup>-1</sup> and ellipticity $`e=1b/a=0.25`$. The semi–major axis of the lens is aligned at $`40^{}`$ counterclockwise from the vertical. The Einstein angle is $`\theta _E=4\pi \sigma ^2D_{ds}/(c^2D_s)=1.58\mathrm{}`$. The source, shown in the top left panel, is contained within a square of size $`0.75\mathrm{}`$ and is modelled as two circular sources of Gaussian profile, binned in $`0.05\mathrm{}`$ pixels, the same as the image pixel size. The peak surface brightness of each source is 1.0, in arbitrary units. One source lies inside the inner caustic, while the second source straddles the inner caustic. This source configuration resembles that inferred for the gravitational lens $`00472808`$ (Wayth et al., 2003). To create a realistic ACS simulation the image was formed by ray tracing, then convolved with the point spread function, and noise added (in the manner described in the following paragraph). For simplicity we modelled the psf as a Gaussian, and chose FWHM $`0.08\mathrm{}`$ which is the resolution of a diffraction–limited telescope of the same diameter as HST, at a wavelength $`\lambda =800`$nm. Because of the slight undersampling, the convolution is made on a sub–pixelized grid and then binned up to the full pixel size.

The data pixels used for the inversion were the 3626 pixels within the annulus in the image plane marked in the figure. This annulus is defined by the region covered by imaging the entire source plane.<sup>1</sup><sup>1</sup>1The region of the central image should also be included for non–singular mass models. An important point to note is that the analysis region must at least cover this annulus, otherwise the counts in some source pixels will be unconstrained and the inversion will fail. A larger region may be used, but if it becomes too large the usefulness of the $`\chi ^2`$ statistic is diminished, as then most of the pixels are in the background. The final step in the simulation is to apply uniform Gaussian random noise over the image plane. The noise level is defined in terms of the total $`S/N_{im}`$ integrated over the annulus. The same noise realisation was used for all the simulations, but scaled in order to vary $`S/N_{im}`$.

The upper middle panel shows the final simulated ACS image. This image, with source pixel size $`0.05\mathrm{}`$, $`S/N_{im}=60`$, and psf FWHM$`=0.08\mathrm{}`$, is the reference test problem to invert. We later vary these three parameters. The parameters of the different models we have run are listed in Table 1. Col. (1) provides the simulation number. The reference problem is numbered 1. Col. (2) gives the source pixel size in arcsec, col. (3) the summed $`S/N`$ in the image, and col. (4) the psf FWHM in arcsec. Col. (5) is a label U or R depending on whether the inversion was unregularized or regularized. Then, in the case of regularized inversion, col. (6) provides the degree of regularization, quantified by the parameter $`N_\lambda `$, which is the increase of $`\chi _{im}^2(\nu )`$ from the minimum in units of the standard deviation $`\sigma (\chi _{im}^2(\nu ))`$. Recalling the discussion in §1, a value $`N_\lambda =1`$ in this column corresponds to the usual criterion for the maximum allowed degree of regularization. The other columns are explained in the next section.

### 3.2 Fixed mass

The main purpose of the simulations is to illustrate the linear inversion step, the ‘inner cycle’, i.e. that part of the semi–linear inversion method that differs from previous methods. For this reason in this sub–section we fix the mass parameters at the input values. The image inversion, therefore, is achieved in a single step using equations (5) and (12), for the unregularized and regularized cases respectively, and using equations (9) and (17) for the source covariance matrix. We consider the full problem, solving also for the mass parameters in §3.3.

#### 3.2.1 Unregularized inversion

The unregularized inversion of the reference problem is provided in the remaining panels of Figure 1. The lower left panel shows the reconstructed source and the upper right panel shows the image of the reconstructed source convolved with the psf, i.e. the min.$`\chi _{im}^2`$ model fit to the simulated image. The bottom right panel shows the source $`\sigma `$ image i.e. the standard deviation in each pixel. This provides a visual impression of the uncertainties – note how the region of lowest $`\sigma `$ is bounded by the inner caustic. However the whole covariance matrix is required for a proper interpretation of the results. The lower middle panel is the source $`S/N`$ image. In all the Figures $`13`$, for the source $`\sigma `$ and $`S/N`$ images the grayscale covers the full range of numbers in the panel. For the other panels the same grayscale range is used in each figure, to allow comparison of the relative noise levels.

We measure several quantities to assess the quality of the inversion, listed in the remaining columns of Table 1. The reduced $`\chi ^2`$ in the image plane, $`\chi _{im}^2(\nu )`$ is provided in col. (7). The quoted uncertainty is given by $`\sqrt{2/\nu }`$ where $`\nu `$ is the number of degrees of freedom i.e. the number of image pixels (3626) minus the number of source pixels (225 or 900). The reduced $`\chi ^2`$ in the source plane, $`\chi _{so}^2(\nu )`$, and its uncertainty, is provided in col. (8). To account for the covariance terms this is computed using

$$\chi _{so}^2=\underset{i,k}{}\mathrm{\Delta }s_i\text{C}_{ik}^1\mathrm{\Delta }s_k=\underset{i,k}{}\mathrm{\Delta }s_i\text{F}_{ik}\mathrm{\Delta }s_k$$
(18)
where $`\mathrm{\Delta }s_i`$ is the residual in the $`i`$th pixel. Here the number of degrees of freedom is the number of source pixels. Col. (9) provides the $`S/N`$ summed over the small box in the source plane shown in Figure 1. The noise is computed as the square root of the sum of the elements in the covariance matrix, formed by stripping out from the source covariance matrix C the rows and columns corresponding to the pixels in the box. Col. (10) provides the absolute value of the significance $`\mathrm{\Delta }s/\sigma `$ of the worst–fit source pixel, and col. (11) lists the $`r.m.s.`$ of the residuals in the source plane.

The results for the reference problem, line (1) in Table 1, are all satisfactory: The reduced $`\chi ^2`$ values in the image and source planes are both consistent with $`1.0`$, and the significance of the worst pixel $`2.80\sigma `$ (col. 10) is not unexpected given that there are 225 source pixels. The summed $`S/N`$ in the source box is an improvement on $`S/N_{im}`$. This might be expected since the box is restricted to the small region of the source plane containing nearly all the signal. At the same time it shows that the $`S/N`$ is not greatly degraded by amplification of noise in the deconvolution step. We return to this issue below. We interpret these results as meaning that the inversion has succeeded and produced the correct solution to the well–posed problem of finding the source–pixel counts that give the best fit to the image.

In simulation 2 we doubled the noise in the image plane. Comparing lines (1) and (2) in the table, the effect of this is to double the noise in the source plane (col. 11), and so halve the S/N of the detected source (col. 9), as expected.

Variation of psf FWHM. We have investigated the effect of varying the width of the psf. Line (3) provides the results for no psf, and line (4) provides the results for a psf FWHM of $`0.16\mathrm{}`$, double the reference value. Comparison of lines (1), (3), and (4) shows that as the psf FWHM increases the noise in the source plane, col. (11), increases and the source detection $`S/N`$, col. (9), decreases. This is as expected: In Fourier space the effect of the psf is to suppress the amplitude of the power spectrum of the source for large wave numbers. Therefore in the deconvolution process the noise on these scales is amplified. As the psf is broadened the power suppression is greater, and so the noise amplification in the deconvolution step is greater. The reduction in $`S/N_{so}`$ from 111.2 (line 3) to 79.9 (line 1), in going from no psf to psf FWHM of $`0.08\mathrm{}`$, is quite modest. This demonstrates that satisfactory inversion of ACS images using $`0.05\mathrm{}`$ source pixels is possible without regularization.

Regardless of the degree of amplification of noise the various statistical quantities in cols (7)–(11) of Table 1, lines (1)–(4), are all reasonable. This shows that in these cases the inversion is well behaved, and in none of the cases is the matrix F singular. This contrasts with the usual inversion problem, for example image deconvolution. With image deconvolution the number of parameters to solve for (the counts in the deconvolved image pixels) is typically the same as or greater than the number of constraints (the number of image pixels). In lensing, because of magnification, the number of image pixels may be much greater than the number of source pixels. This suggests, further, that in regions where the magnification is greatest it would be possible to use source pixels smaller than the image pixels. We consider this issue below.

The results of these first four simulations indicate that, in some circumstances, provided the psf is not too broad, unregularized semi–linear inversion provides a useful solution.

Variation of source pixel size. Figure 2 shows the same problem as Figure 1 but with $`0.025\mathrm{}`$ source pixels rather than $`0.05\mathrm{}`$ pixels. The results are summarized in line (5) of Table 1. The quality of the reconstruction, lower left panel, is now dramatically worse, and outside the central region is clearly unsatisfactory. (The grayscale range of this panel is the same as in Figure 1.) Compared to line (1) the noise in the source (col. 11) has risen by a factor 8, whereas intuitively one would expect only a factor 2 increase (4 times as many pixels). This is indicative of large amplification of noise, because the psf has suppressed the signal on these scales. This is a consequence of the fact that the separation in the image plane of the images of two adjacent source pixels is smaller than the psf size. Put another way, a resolution element in the image plane, traced back to the source plane, is oversampled by the source pixel size. The source covariance matrix now contains large, predominantly negative, covariance terms which correspond to the odd/even appearance in the outer regions of the source plane.

The high noise level in the outer parts of the source–plane belies the usefulness of this image. In fact the source is strongly detected, albeit at reduced $`S/N_{so}=20`$, even though not readily apparent to the eye. The source is clearly visible in the $`S/N`$ image, however. At the same time the $`\chi ^2`$ values in both the image and source planes remain satisfactory. Because of the larger magnification, the reconstruction is much better within the caustic line. This suggests it would be advantageous to use a variable pixel size across the source plane. For example, with reference to Figure 2, a scheme where the pixel size is $`0.05\mathrm{}`$ outside the caustic and $`0.025\mathrm{}`$ inside might be appropriate. We need to identify a criterion for choosing the pixel size that avoids the excessive amplification of noise evident in Figure 2. There are clearly three variables which determine the minimum source pixel size: The image pixel size, $`a`$, the psf FWHM, $`b`$, and the magnification, $`c`$. We have had some success with a scheme which relates the source pixel size to the variable $`max(a,b/2)/c^{1/2}`$. The results will be reported elsewhere (Dye and Warren, in prep.).

#### 3.2.2 Regularized inversion

We now include linear regularization in the inversion. All the results reported here used the gradient form, equation (14). The results are quite similar for the different linear regularizing schemes described in §2.1.1, however.

Lines (6) and (7) in Table 1 are the results for the reference problem with different degrees of regularization. As the regularizing term increases, the source becomes smoother and $`\chi _{im}^2`$ increases. Line (6) is for $`N_\lambda =1`$. Comparing line (6) to line (1) we see that the effect of regularization is to suppress the noise in the reconstructed source and to increase substantially the source $`S/N`$, col. (9). This is at the expense of a poorer match to the true source light profile, as measured by cols (8) and (10). For $`N_\lambda =2`$ the agreement with the input source is no longer acceptable.

In lines (8) to (10) we provide solutions for source pixel size $`0.025\mathrm{}`$, and psf FWHM of $`0.08\mathrm{}`$, and different degrees of regularization, $`N_\lambda =1,3,5`$. These results compare directly to the unregularized solution to the same problem, line (5). The solutions for simulation no. 8, $`N_\lambda =1`$, and no. 9, $`N_\lambda =3`$, are shown in Figure 3. The visual improvement, comparing the sequence of Figure 2 (unregularized), Figure 3 middle row (regularzsed, $`N_\lambda =1`$), and Figure 3 bottom row (regularized, $`N_\lambda =3`$), is dramatic.

Comparing lines (8) to (10) against line (5) we see that, again, regularization successfully suppresses noise, increasing the $`S/N`$ of the detection of the source. As $`N_\lambda `$ increases, in this case $`\chi _{so}^2`$ increases only very slowly, much more slowly than in the case for larger pixels. This is partly due to the fact that we chose a smooth source, and the results would be different for a source with more small–scale structure. Nevertheless, it indicates that the standard criterion for the degree of regularization to apply, $`N_\lambda =1`$, is somewhat arbitrary.

To summarise this sub–section, using a realistic problem, we have validated the theory of the linear inversion step set out in §2. This is the step that differs from the maximum–entropy method of Wallington et al. (1996), and therefore is the main point of the paper.

### 3.3 Mass cycle

In the present sub-section we report the results of solving the complete problem i.e. determining both the mass profile and the source light distribution.

We first consider the unregularized case. Referring back to the example of Figure 1, the problem is to invert the image at upper middle. The free parameters are the five parameters describing the mass: $`x`$, $`y`$, ellipticity, position angle, and velocity dispersion. We searched through the parameter space to find the min.$`\chi ^2`$ fit. At the minimum the matrix F supplies most of the terms of the curvature matrix. Following the precepts of §2.1.2, the remaining terms were filled in by measuring the relevant second partial derivatives of the $`\chi ^2`$ surface. We found the surface to be completely smooth and parabolic near the minimum. The full curvature matrix was inverted to obtain the full covariance matrix for all the parameters, the mass terms as well as the counts in the source pixels. We found the input mass parameters were correctly recovered to within the uncertainties. We checked the full covariance matrix against the results of Monte–Carlo simulations and found excellent agreement. This confirms that, provided the chosen source pixel size is not too small, the unregularized semi–linear inversion method is a practical solution to the problem of inversion of a gravitational lens image with a resolved source.

We also compared the terms in the source covariance matrix $`\text{C}=\text{F}^1`$ against the corresponding terms in the full covariance matrix. The differences are relatively small. Therefore the matrix C, at the global minimum, provides an approximation to the true source–pixel errors that may be very useful in the exploration stage, when considering different mass models and different pixelizations.

In the regularized case we found, generally, that the procedure converged more rapidly than in the unregularized case. Regularized inversion can produce solutions which are not true representations of the source (§3.4, Table 1). Nevertheless, we found, in contrast, that the solution for the mass parameters is very insensitive to the degree of regularization. In the regularized case, the curvature of the merit function cannot be used to obtain the full covariance matrix so that an alternative approach such as a Monte Carlo method must be adopted (§2.2.2). So the source covariance matrix, equation (17), at the global minimum, is particularly useful here as an approximation to the true source–pixel errors.

### 3.4 Regularized vs unregularized

We now include a debate on the relative advantages of the unregularized and regularized approaches. This may seem surprising given the excellent results achieved with regularization (comparing Figures 2 and 3). The weakness of the unregularized inversion is that in deconvolving the psf, the noise at large wavenumbers is boosted. The regularization term in the merit function imposes smoothness on the solution. In effect, the deconvolution (division in Fourier space) is limited to the smaller wave numbers. However, this means that any real structure in the source at large wavenumbers is also suppressed. We are imposing a prejudice that the source is smooth and this might not be justified (see comments in §1). Regularization introduces $`+`$ve covariance between adjacent pixels, forcing the counts to be similar. The regularized solution, then, is not so different to the unregularized solution with larger pixels. In this respect, it is interesting to compare the unregularized solution with $`0.05\mathrm{}`$ source pixels (Figure 1), with the $`N_\lambda =1`$ regularized solution with $`0.025\mathrm{}`$ source pixels (Figure 3). Noting that the two values of the source $`S/N`$ are quite similar (79.9, 64.5, respectively, Table 1), it is a debatable point whether there is more information in the latter figure.

A further point to note is that the regularized inversion can produce solutions which are satisfactory in terms of the fit to the image but which are not true representations of the source in the sense that the $`\chi _{so}^2`$ is unsatisfactory. For example in line (7), Table 1, both $`\chi _{so}^2`$ and $`|\mathrm{\Delta }s/\sigma |`$ are unsatisfactory. Measured by the same statistics, none of the unregularized inversions in the Table is unsatisfactory. The unregularized inversion gives a noisier but unbiased solution for the source light distribution, while the regularized inversion gives a smoother but biased solution.

Despite regularization biasing the source, in the full mass cycle we find that the minimized mass parameters show little sensitivity to the degree of regularization. Furthermore, the regularized solution has the advantage that it converges more quickly. We discuss the practical significance of these two points in §4. Associated with this is the fact that regularization allows source pixel sizes of almost any size, unlike the unregularized case when pixel size must be chosen carefully. This can yield further speed advantages in the initial stages of an analysis, before the solution is refined.

Overall we consider there are important advantages to using both regularized and unregularized inversion in exploring the solution to a particular problem, and the choice will depend on the question being posed and any a priori knowledge concerning the source. Perhaps equally importantly, however, it makes sense to match the source pixel size to the data information content in terms of the $`S/N`$ at different wavenumbers, or, in other words, to vary the pixel size depending on the magnification, as suggested in §3.2.1.

## 4 Summary and recommendations

We have developed a new method for the inversion of gravitationally–lensed images of extended sources for the case where the source light profile is pixelized. The method separates the linear dimensions of the problem (the counts in the source pixels) from the non–linear dimensions (the mass parameters). The method has been extended in a natural way to allow linear regularization of the inversion. The core of the routine is the procedure for inverting an image given a fixed mass profile. We have shown that this step, including deconvolution of the psf, with or without regularization, is a linear one. Since this step is usually achieved by searching the source parameter space for the merit–function minimum, the solution is reached much more quickly. The non–linear part of the problem has been reduced to the search for a minimum in the space of the mass parameters only. In the case of unregularized inversion, the full covariance matrix for all the (source$`+`$mass) parameters can be obtained very quickly. In the case of regularized inversion, a useful approximation to the covariance matrix for the source counts is obtained very simply, but Monte Carlo methods are needed to obtain the full covariance matrix.

How the semi–linear method should be applied in practice depends on the problem posed. If one is interested in the quantitative details of the source light profile, for example, whether some apparent feature is real, then we recommend the unregularized solution. This is because regularization produces source profiles which are too smooth. Without regularization, an optimal source pixel size should be chosen. Having too large a source pixel may cause interesting detail to be lost. However, if the source pixel size is too small, the inverted image may have low $`S/N`$ because of amplification of noise in the deconvolution step. If, on the other hand, one is interested only in the mass parameters, a regularized solution would be the appropriate choice: The mass parameters are rather insensitive to the degree of regularization and one benefits from an increase in inversion speed.

Another consideration is that pixelizing the source uses a large number of parameters. As a rule, one is interested in finding the model with the smallest number of parameters that provides a satisfactory fit to an image. Therefore, in many cases, one might simply use the semi–linear method of inversion to provide an image of the source to guide the choice of parameterization. Here, again, the regularized solution might be the preferred option.

In general, because it is so much easier to implement (§2.2.1), we recommend using zeroth–order regularization. Nevertheless, other considerations may override simplicity. The zeroth–order regularization term, in common with the maximum–entropy regularization term, is a local measure, independent of the counts in adjacent pixels. This can be an advantage or a disadvantage, depending on the actual light profile in the source.

In the simulations presented here, we have used square source pixels which form a regular grid. However, since the resolution across the image plane is fixed while the magnification varies, the resolution across the source varies. Therefore, to maximize the information content in the reconstruction of the source it is necessary to use a variable source pixel size. We will present an analysis of semi–linear inversion with variable source pixel size in a future paper (Dye and Warren, in prep.).

We have benefited from discussions with Paul Hewett, Geraint Lewis, Leon Lucy, and Randall Wayth.