text:
# Forcing a unitary transformation by an external field: comparing two approaches based on optimal control theory.

## I Introduction.

Quantum computation is based on implementing selected unitary transformations representing algorithms nielsen00 . In many physical implementations the unitary transformation is generated using an external driving field. This driving field has to perform the quantum gate between, for example, two qubits, without altering the other levels which represent additional qubits. This means that the specific unitary transformation has to address a set of levels in an environment where other energy levels are present. An approach to correct this undesired coupling to other levels has been suggested for specific cases tian00  but a general solution is not known.

The presence of a large number of levels coupled to the external driving field is specially relevant in the implementation of quantum computing in molecular systems zadoyan01 ; amitay02 ; tesch01 . Tesh et al. tesch01  proposed the use of optimal control theory (OCT) as a possible remedy. OCT rice00  is a well developed approach that allows to obtain the driving field which induces a desired transition between preselected initial and final states. However, the dependence on the particular transitions (characterized by the initial and final states) makes the traditional formulation of OCT inappropriate for that purpose. For example, if the unitary transformation relates the initial states $`\phi _{ik}`$ with the final states $`\phi _{fk}`$ (the index $`k`$ denotes all the relevant states involved in the transformation), the traditional OCT approach derives an optimal field $`ϵ_k`$ for each pair $`\{\phi _{ik},\phi _{fk}\}`$. But the fields $`ϵ_k`$ obtained are in general different, and then the transition induced by $`ϵ_k`$ for the initial state $`\phi _{ik^{}}`$ won’t necessarily give the right final state $`\phi _{fk^{}}`$. To implement a given unitary transformation a single field $`ϵ`$ that relates all the relevant pairs $`\{\phi _{ik},\phi _{fk}\}`$ is needed.

An approach to overcome this problem was suggested in palao01poster ; palao02 . The idea is to generalize OCT to deal directly with the evolution operator. Recently, a different approach using simultaneous optimization of several state to state transitions has been suggested by Tesch and Vivie-Riedle tesch02 . The purpose of this work is to compare the two approaches and two point the similarities and differences between them. This comparison has led to new insight into the use of optimal control theory for quantum systems.

Unitary transformation optimization can be formulated in the following way. We consider a quantum system with a Hilbert space of dimension $`M`$, expanded by an orthonormal basis of states $`\{|k\}`$ ($`k=1,\mathrm{},M`$). In the following we will assume that the $`k`$ states correspond to the free Hamiltonian proper states. The selected unitary transformation is imposed on the subspace of the first $`N`$ energy levels of the system ($`NM`$). For example, the $`N`$ levels could correspond to the physical implementation of the qubit(s) embedded in a larger system. The additional levels ($`k=N+1,\mathrm{},M`$) are not directly involved in the target unitary transformation and they are generally considered as “spurious levels” coupled to the field. However, it is not always the case: an example is the proposal of implementing quantum computation using rovibronic molecular levels. In the simplest description two electronic surfaces are considered. Two rovibronic states of one of the electronic surfaces are chosen as the implementation of a qubit (in this case $`N=2`$) and the unitary transformation is implemented using field induced transitions between the electronic surfaces. Out of the relevant subspace ($`k=3,\mathrm{},M`$) there are levels residing on both electronic surfaces, spurious in the sense that any leakage to them would destroy the desired final result, but at intermediate times they are used as intermediate storage which allow to carry out the desired unitary transformation palao02 ; note1 .

The optimization objective is to implement a selected unitary transformation in the relevant subspace at a final time $`T`$. The target unitary transformation is represented by an operator in the system Hilbert space and denoted by $`\widehat{𝐎}`$. For $`N<M`$, the matrix representation of $`\widehat{𝐎}`$ in the basis $`\{|k\}`$ has two blocks of dimension $`N\times N`$ and $`(MN)\times (MN)`$. The elements connecting these blocks are zero. This structure means that population at the target time is not transferred between the two subspaces. Only the $`N\times N`$ block is relevant for the optimization procedure, and the other block remains arbitrary.

The dynamics of the system is generated by the Hamiltonian $`\widehat{𝐇}`$,

$$\widehat{𝐇}(t;ϵ)=\widehat{𝐇}_0\widehat{𝝁}ϵ(t),$$
(1)
where $`\widehat{𝐇}_0`$ is the free Hamiltonian, $`ϵ(t)`$ is the driving field and $`\widehat{𝝁}`$ is a system operator describing the coupling (transition dipole operator in molecular systems). Eq. (1) can be generalized to more than one independent driving field, for example, controlling separately the two components of the polarization of an electro-magnetic field brixner01 . The dynamics of the system at time $`t`$ is fully specified by the evolution operator $`\widehat{𝐔}(t,0;ϵ)`$. An optimal field $`ϵ_{opt}`$ induces the target unitary transformation $`\widehat{𝐎}`$ on the system, at time $`T`$ if

$$\widehat{𝐔}(T,0;ϵ_{opt})=e^{i\varphi (T)}\widehat{𝐎}.$$
(2)
Eq. (2) implies a condition only on the $`N\times N`$ block of the matrix representation of $`\widehat{𝐔}`$. The phase $`\varphi (T)`$ is introduced to point out that in some cases the target unitary transformation $`\widehat{𝐎}`$ can be implemented up to an arbitrary global phase. The phase $`\varphi `$ can be decomposed into two terms, $`\varphi _1(T)+\varphi _2(T)`$. The term $`\varphi _1`$ originates from the arbitrary choice of the origin of the energy levels which formally means that a term proportional to the identity operator can always we added to the Hamiltonian. The phase $`\varphi _1`$ is given by,

$$\varphi _1(T)=\frac{_{k=1}^ME_kT}{M\mathrm{}},$$
(3)
where $`E_k`$ is the energy of the level $`k`$. By $`\varphi _2`$ we denoted other contributions to the global phase due to the structure of the unitary transformation and its arbitrariness for the levels $`k=N+1,\mathrm{},M`$.

## II Evolution equation approach.

The optimization approach proposed in palao01poster ; palao02  is based on Eq. (2) by defining a complex parameter $`\tau `$ as

$$\tau (\widehat{𝐎};T;ϵ)=\underset{k=1}{\overset{N}{}}k|\widehat{𝐎}^{}\widehat{𝐔}(T,0;ϵ)|k.$$
(4)
As $`\widehat{𝐎}`$ is a unitary transformation in the relevant subspace, $`\tau `$ is a complex number inside a circle of radius $`N`$. Its modulus is equal to $`N`$ only when the unitary transformation generated by the field $`\widehat{𝐔}`$ is equal to the target unitary transformation in the relevant subspace, except for a possible global phase. The modulus of $`\tau `$ is a measurement of the fidelity of the target unitary transformation implementation by the field palao02 . For $`N=M`$, the sum in Eq. (4) is the trace of the operator product.

In palao01poster ; palao02  the optimization of the real part of $`\tau `$, or the imaginary part, or a linear combination of both was suggested as a method to find the optimal field. For simplicity we will consider the real case. The maximization of $`\mathrm{Re}[\tau ]`$ can be formulated as a functional optimization. In this work we use the following form

$`\overline{J}(\widehat{𝐔},\widehat{𝐁},\mathrm{\Delta }ϵ)`$ $`=`$ $`\mathrm{Re}\left[{\displaystyle \underset{k=1}{\overset{N}{}}}k|\widehat{𝐎}^{}\widehat{𝐔}(T,0;ϵ_0+\mathrm{\Delta }ϵ)|k\right]\lambda {\displaystyle _0^T}|\mathrm{\Delta }ϵ|^2𝑑t`$ (5)
$``$ $`\mathrm{Re}\left[{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle _0^T}k|\widehat{𝐁}\left({\displaystyle \frac{}{t}}+{\displaystyle \frac{i}{\mathrm{}}}\widehat{𝐇}(t;ϵ_0+\mathrm{\Delta }ϵ)\right)\widehat{𝐔}|k𝑑t\right],`$
where $`\widehat{𝐔}`$, $`\widehat{𝐁}`$, $`ϵ_0`$, and $`\mathrm{\Delta }ϵ`$ depends on $`t`$. The first term in the right-hand-side is the original objective. In our formulation, the other two terms are constrains depending on a reference field $`ϵ_0`$ and a field correction $`\mathrm{\Delta }ϵ`$. The term including $`|\mathrm{\Delta }ϵ|^2`$ minimizes the total energy of the correction. The last term introduces the dynamics of the system under the field $`ϵ_0+\mathrm{\Delta }ϵ`$. $`\widehat{𝐁}`$ is an operator Lagrange multiplier and $`\lambda `$ a scalar Lagrange multiplier palao02 . The common OCT form for the functional rice00  is recuperated setting $`ϵ_0=0`$ and interpreting $`\mathrm{\Delta }ϵ`$ as the driving field. However, Eq. (5) offers some advantages in the interpretation of the equations derived from the functional. More elaborated constrains are possible sundermann99 ; hornung01 , but they are not relevant for the following discussion.

Applying the calculus of variations, $`\delta \overline{J}=0`$, with respect to $`\widehat{𝐁}`$, the Schrödinger equation for the evolution operator of the system is obtained,

$$\frac{\widehat{𝐔}}{t}=\frac{i}{\mathrm{}}\widehat{𝐇}(t;ϵ_0+\mathrm{\Delta }ϵ)\widehat{𝐔},$$
(6)
with the condition $`\widehat{𝐔}(0,0;ϵ_0+\mathrm{\Delta }ϵ)=𝟙`$. The variation of $`\widehat{𝐔}`$ give the Schrödinger evolution equation for the operator $`\widehat{𝐁}`$,

$$\frac{\widehat{𝐁}^{}}{t}=\frac{i}{\mathrm{}}\widehat{𝐇}(t;ϵ_0+\mathrm{\Delta }ϵ)\widehat{𝐁}^{},$$
(7)
with the condition $`\widehat{𝐁}^{}(T,T;ϵ_0+\mathrm{\Delta }ϵ)=\widehat{𝐎}`$. $`\widehat{𝐁}^{}`$ can be interpreted as the backwards propagation in time of the target unitary transformation $`\widehat{𝐎}`$. It is related to $`\widehat{𝐔}`$ by

$$\widehat{𝐁}^{}(t,T;ϵ)=\widehat{𝐔}(t,T;ϵ)\widehat{𝐎}.$$
(8)
Eq. (6) and (7) represent the propagation forward and backwards in time of the boundary conditions of the problems, that is, the identity, $`𝟙`$, at time $`t=0`$ and the target unitary transformation, $`\widehat{𝐎}`$, at time $`t=T`$. The variation of $`\mathrm{\Delta }ϵ`$ leads to an equation for the correction to the field,

$$\mathrm{\Delta }ϵ(t)=\frac{1}{2\lambda \mathrm{}}\mathrm{Im}[\underset{k=1}{\overset{N}{}}k|\widehat{𝐁}(t,T;ϵ_0+\mathrm{\Delta }ϵ)\widehat{𝝁}\widehat{𝐔}(t,0;ϵ_0+\mathrm{\Delta }ϵ)|k].$$
(9)
Eq. (9), more than Eq. (5), is the central result of the optimal control procedure and constitutes the starting point of the iterative algorithms devoted to determine the optimal field palao02 . When an optimal field $`ϵ_{opt}`$ is found $`\mathrm{\Delta }ϵ=0`$, and from Eq. (8) and Eq. (9),

$$\mathrm{Im}[\underset{k=1}{\overset{N}{}}k|\widehat{𝐎}^{}\widehat{𝐔}^{}(t,T;ϵ_{opt})\widehat{𝝁}\widehat{𝐔}(t,0;ϵ_{opt})|k]=0.$$
(10)
Eq. (10) constitutes a condition for the optimal fields for the first approach. It is the base for the following analysis. Let us denoted by $`\stackrel{~}{ϵ}`$ a field that generated the target unitary transformation up to a global phase, $`\widehat{𝐔}(T,0;\stackrel{~}{ϵ})=e^{i\varphi }\widehat{𝐎}`$. Using

$$\widehat{𝐔}(t,0,ϵ)=\widehat{𝐔}(t,T,ϵ)\widehat{𝐔}(T,0,ϵ),$$
(11)
the diagonal block structure of the matrix representation of $`\widehat{𝐎}`$ in the basis $`\{|k\}`$, and the relation $`k|\widehat{𝝁}|k=0`$, it is found that the left-hand-side of Eq. (10) gives

$$\mathrm{Im}[\underset{k=1}{\overset{N}{}}e^{i\varphi }k|\widehat{𝐔}(t,T;\stackrel{~}{ϵ})\widehat{𝐎}\widehat{𝐎}^{}\widehat{𝐔}^{}(t,T;\stackrel{~}{ϵ})\widehat{𝝁}|k]=0.$$
(12)
This result implies that any field inducing the target unitary transformation up to a global phase is a possible optimal solution of the optimization algorithm based on Eq. (9). The convergence to such a solution will depend on the particular numerical implementation. It must be remarked that Eq. (10) is a necessary but not a sufficient condition for the optimal field. For example, let us consider a target unitary transformation diagonal in the basis $`\{|k\}`$ denoted by $`\widehat{𝐎}_D`$. The unitary transformation generated by the free Hamiltonian $`\widehat{𝐔}(t_1,t_2;ϵ=0)`$ is also diagonal in that basis, and then,

$$\mathrm{Im}[\underset{k=1}{\overset{N}{}}k|\widehat{𝐎}_D^{}\widehat{𝐔}^{}(t,T;ϵ=0)\widehat{𝝁}\widehat{𝐔}(t,0;ϵ=0)|k]=0.$$
(13)
However, in general $`\widehat{𝐔}(T,0;ϵ=0)`$ is not equal to the target $`\widehat{𝐎}_D`$. These spurious solutions to the optimization can be avoided with a different choice of the initial guess for the algorithm.

In optimization procedures based on Eq. (9) the full operator propagation in Eq. (6) and (7) is not needed since only the action of $`\widehat{𝐔}`$ and $`\widehat{𝐎}`$ on the states $`|k`$ in the relevant subspace ($`k=1,\mathrm{},N`$) appears. Then only the first $`N`$ rows of the matrix operator representations are propagated. Denoting by $`U^k`$ ($`(B^{})^k`$) the $`k`$ row of the matrix representation of $`\widehat{𝐔}`$ ($`\widehat{𝐁}^{}`$) in the basis $`\{|k\}`$, and being $`U_j^k`$ ($`(B^{})_j^k`$) the $`j`$ element of the row ($`j=1,\mathrm{},M`$), the evolution equations (6) and (7) take the form,

$`{\displaystyle \frac{U^k(t)}{t}}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}H(t)U^k(t),`$
$`{\displaystyle \frac{(B^{})^k(t)}{t}}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}H(t)(B^{})^k(t),`$ (14)
with the conditions $`U_j^k(t=0)=\delta _{jk}`$ and $`(B^{})_j^k(t=T)=O_j^k`$ respectively, being $`O_j^k`$ the matrix elements of the target unitary transformation $`\widehat{𝐎}`$ and $`H(t)`$ the matrix representation of $`\widehat{𝐇}`$. The $`2N`$ evolution equations ($`k=1,\mathrm{},N`$) are equivalent to the propagation of $`2N`$ states of the system. The advantage is that in this case only $`2(N\times M)`$ elements are propagated instead of the $`2(M\times M)`$ in Eq. (6) and Eq. (7).

## III State to state approach.

The state to state approach tesch02  is based on the simultaneous optimization of $`N`$ transitions between pairs of initial and final states. These pairs of states $`\{\phi _{il},\phi _{fl}\}`$, ($`l=1,\mathrm{},N`$) are related by the target unitary transformation, $`|\phi _{fl}=\widehat{𝐎}|\phi _{il}`$. The objective is formulated as

$$\eta (\widehat{𝐎};T;ϵ)=\underset{l=1}{\overset{N}{}}|\psi _{il}(T;ϵ)|\phi _{fl}|^2,$$
(15)
where $`\psi _{il}(T;ϵ)`$ is the state at the target time that evolves with the driving field $`ϵ(t)`$ and the initial condition $`\psi _{il}(t=0)=\phi _{il}`$. $`\eta `$ is a positive real number and its maximum value $`\eta =N`$ is reached when all the initial states $`\phi _{il}`$ are driven by the field to the correct final states $`\phi _{fl}`$. The task of obtaining the optimal field is equivalent to the maximization of $`\eta `$. The set of initial states $`\phi _{il}`$ must be chosen carefully. In order to account for all the possible transitions the states $`\phi _{il}`$ have to represent the relevant subspace. However, the choice of an orthonormal basis could produce undesired results. Let us denoted by $`\{\stackrel{~}{\phi }_{il}\}`$ an orthonormal basis of the relevant subspace and by $`\widehat{𝐃}`$ an arbitrary unitary transformation diagonal in that basis. The product $`\widehat{𝐎}\widehat{𝐃}`$ is also a unitary transformation. If $`ϵ_O`$ and $`ϵ_{OD}`$ are fields that generate $`\widehat{𝐎}`$ and $`\widehat{𝐎}\widehat{𝐃}`$ at time $`T`$ respectively, it is found that the same optimization objective is obtained,

$$\eta _{}(\widehat{𝐎};T;ϵ_O)=\eta _{}(\widehat{𝐎};T;ϵ_{OD}),$$
(16)
where $``$ denotes that $`\eta `$ was evaluated using an orthonormal basis. Then any algorithm based on $`\eta `$ and using an orthonormal basis would find the optimal field corresponding to any of the possible targets $`\widehat{𝐎}\widehat{𝐃}`$ ($`\widehat{𝐎}`$ is a particular case when $`\widehat{𝐃}`$ is the identity operator). The reason is that $`\eta `$ is sensitive to the modulus of the projection of each pair $`\{\phi _{il},\phi _{fl}\}`$ but not to the relative phases between them. The phase problem can be overcome with a careful choice of the states $`\phi _{il}`$. Keeping the first $`N1`$ states of the basis and substituting the $`N`$ state by $`_{l=1}^N\stackrel{~}{\phi }_{il}/\sqrt{N}`$, the maximum condition is achieved only when the field induces the target unitary transformation up to a possible global phase.

The optimization is formulated as the maximization of the functional

$`\overline{K}(\psi _{ik},\psi _{fk},\mathrm{\Delta }ϵ)`$ $`=`$ $`{\displaystyle \underset{l=1}{\overset{N}{}}}|\psi _{il}(T;ϵ_0+\mathrm{\Delta }ϵ)|\phi _{fl}|^2\lambda {\displaystyle _0^T}|\mathrm{\Delta }ϵ(t)|^2𝑑t`$
$``$ $`2\mathrm{Re}\left[{\displaystyle \underset{l=1}{\overset{N}{}}}\psi _{il}(T;ϵ_0+\mathrm{\Delta }ϵ)|\phi _{fl}{\displaystyle _0^T}\psi _{fl}|\left({\displaystyle \frac{}{t}}+{\displaystyle \frac{i}{\mathrm{}}}\widehat{𝐇}(t;ϵ_0+\mathrm{\Delta }ϵ)\right)|\psi _{il}𝑑t\right],`$
where $`\psi _{il}`$, $`\psi _{fl}`$, $`ϵ_0`$ and $`\mathrm{\Delta }ϵ`$ depend on time. The form in Ref. tesch02  is recuperated setting $`ϵ_0=0`$ and interpreting $`\mathrm{\Delta }ϵ`$ as the driving field.

The variations with respect to $`\psi _{li}`$ and $`\psi _{lf}`$ ($`l=1,\mathrm{},N`$) lead to the following $`2N`$ equations

$$\frac{}{t}|\psi _{il}=\frac{i}{\mathrm{}}\widehat{𝐇}(t;ϵ_0+\mathrm{\Delta }ϵ)|\psi _{il},$$
(18)
with the condition $`|\psi _{il}(t=0)=|\phi _{il}`$, and

$$\frac{}{t}|\psi _{fl}=\frac{i}{\mathrm{}}H(t;ϵ_0+\mathrm{\Delta }ϵ)|\psi _{fl},$$
(19)
with the condition $`|\psi _{fk}(t=T)=|\phi _{fk}`$. The variation respect to $`\mathrm{\Delta }ϵ`$ leads to the equation for the correction to the field,

$$\mathrm{\Delta }ϵ_{ss}(t)=\frac{1}{\lambda \mathrm{}}\mathrm{Im}\left[\underset{l=1}{\overset{N}{}}\psi _{fl}(T;ϵ_0+\mathrm{\Delta }ϵ)|\phi _{fl}\psi _{fl}(t;ϵ_0+\mathrm{\Delta }ϵ)|\widehat{𝝁}|\psi _{il}(t;ϵ_0+\mathrm{\Delta }ϵ)\right].$$
(20)
Eq. (18) and (19) are formally equivalent to Eq. (II), propagating forward or backward on time the boundary conditions of the problem. Then the main difference between the two approaches is found in the expression for the correction to the field, Eq. (9) and Eq. (20). The connection between the two expressions is established by rewriting the states $`\phi `$ and $`\psi `$ as

$`|\phi _{il}`$ $``$ $`|l,`$
$`|\phi _{fl}`$ $``$ $`\widehat{𝐎}|l,`$
$`|\psi _{il}(t;ϵ)`$ $``$ $`\widehat{𝐔}(t,0;ϵ)|l,`$
$`|\psi _{fl}(t;ϵ)`$ $``$ $`\widehat{𝐁}^{}(t,T;ϵ)|l`$ (21)
and Eq. (20) as,

$$\mathrm{\Delta }ϵ_{ss}(t)=\frac{1}{\lambda \mathrm{}}\mathrm{Im}\left[\underset{l=1}{\overset{N}{}}l|\widehat{𝐔}^{}(T,0;ϵ_0+\mathrm{\Delta }ϵ)\widehat{𝐎}|ll|\widehat{𝐁}(t,T;ϵ_0+\mathrm{\Delta }ϵ)\widehat{𝝁}\widehat{𝐔}(t,0;ϵ_0+\mathrm{\Delta }ϵ)|l\right].$$
(22)
This expression is formally equivalent to Eq. (9) except for the factor $`l|\widehat{𝐔}^{}(T,0,ϵ_0+\mathrm{\Delta }ϵ)\widehat{𝐎}|l`$. (The different factor $`1/2`$ can be removed with a redefinition of $`\lambda `$). Omitting this factor the two approaches would be completely equivalent: in that case the careful choice of the initial state set wouldn’t be necessary and the two approaches would lead to the same set of equations. In general the equivalence is only formal due to the additional factor and the different set of states $`\{|k\}`$ and $`\{|l\}`$. In the following we will chose $`|l=|k`$ for $`l=1,\mathrm{},N1`$ and $`|l=N=_{k=1}^N|k/\sqrt{N}`$. The condition for the optimal field ($`\mathrm{\Delta }ϵ=0`$) is obtained using Eq. (20) and (8),

$$\mathrm{Im}\left[\underset{l=1}{\overset{N}{}}l|\widehat{𝐔}^{}(T,0;ϵ_{opt})\widehat{𝐎}|ll|\widehat{𝐎}^{}\widehat{𝐔}^{}(t,T;ϵ_{opt})\widehat{𝝁}\widehat{𝐔}(t,0;ϵ_{opt})|l\right]=\mathrm{\hspace{0.17em}0}.$$
(23)
Eq. (23) is a necessary but not sufficient condition for the optimal field, as is also Eq. (10). Denoting as before by $`\stackrel{~}{ϵ}`$ a driving field generating the target unitary transformation up to a global phase, the left-hand-side of (23) gives,

$$\mathrm{Im}\left[\underset{l=1}{\overset{N}{}}l|\widehat{𝐎}^{}\widehat{𝐔}^{}(t,T;\stackrel{~}{ϵ})\widehat{𝝁}\widehat{𝐔}(t,T;\stackrel{~}{ϵ})\widehat{𝐎}|l\right]=\mathrm{\hspace{0.17em}0}.$$
(24)
(The property $`\mathrm{Im}[\mathrm{\Psi }|\widehat{𝝁}|\mathrm{\Psi }]=0`$ for Hermitian operator $`\widehat{𝝁}`$ was used). Contrary to Eq. (12), the phase $`\varphi `$ doesn’t appears in Eq. (24) due to the factor $`l|\widehat{𝐔}^{}(T,0,\stackrel{~}{ϵ})\widehat{𝐎}|l`$ correcting each term in the sum. Using the same arguments leading to Eq. (13), it can be shown that the left-hand-side of Eq. (23) is null for $`ϵ=0`$ when the target unitary transformation is diagonal. The previous results imply that any field that generates the target unitary transformation up to a global phase fulfills the condition in Eq. (23) and in addition spurious solutions could be found, as in the first approach.

## IV Conclusions

In this study we have employed a modified formulation of OCT which uses the correction to the driving field $`\mathrm{\Delta }ϵ`$ as central element. $`\mathrm{\Delta }ϵ=0`$ is a necessary but not sufficient condition for obtaining the objective. The relations derived from this condition allow a better analysis of the control equations, in particular, the similarities and differences between the two approaches. Optimal solutions of the evolution equation approach will be also optimal solutions of the state to state approach. In this sense both approaches are equivalents. The difference between them is a term in the state to state approach which modified the phases. This term can cause a phase ambiguity in the target unitary transformation. A careful choice of the initial set of states can solve this problem.

## V Acknowledgments

We want to thank Lajos Diosi and Zohar Amitay for many useful discussions. This work was supported by the Spanish MCT BFM2001-3349 and the Israel Science Foundation. The Fritz Haber Center is supported by the Minerva Gesellschaft für die Forschung, GmbH München, Germany.


format: text

text:
# References

Reply to “Comments on Kullback-Leibler and renormalized entropies: Applications to electroencephalograms of epilepsy patients”.

R. Quian Quiroga1,2, J. Arnhold2,3, K. Lehnertz3 and P. Grassberger2

1Sloan-Swartz Center for Theoretical Neurobiology.
California Institute of Technology, MC 216-76. 91125 Pasadena, CA. USA
2John von Neumann Institute for Computing,
Forschungszentrum Jülich GmbH, D - 52425 Jülich, Germany
3Department of Epileptology, University of Bonn,
Sigmund-Freud Str. 25, D - 53105 Bonn, Germany

## Abstract

Kopitzki et al (preceeding comment) claim that the relationship between Renormalized and Kullback-Leibler entropies has already been given in their previous papers. Moreover, they argue that the first can give more useful information for e.g. localizing the seizure-generating area in epilepsy patients.

In our reply we stress that if the relationship between both entropies would have been known by them, they should have noticed that the condition on the effective temperature is unnecessary. Indeed, this condition led them to choose different reference segments for different channels, even if this was physiologically unplausible. Therefore, we still argue that it is very unlikely that renormalized entropy will give more information than the conventional Kullback-Leibler entropy.

We thank the authors of the preceding comment for pointing out a misprint in  (in the line following Eq. (8) it should read $`p\stackrel{~}{q}`$ instead of $`pq`$), and a numerical inconsistency in Figs. 2-4 of ref. . The latter resulted from an error in our code. The corrected data is shown in Fig. 1 below, where we added the renormalized entropy values calculated with a pre-seizure reference for completeness. Note that Eq. (10) of , i.e. $`|\mathrm{\Delta }H|K(p|q)`$, is now verified in all cases . Despite this correction, the data are qualitatively similar to the ones presented in , and we still conclude that renormalized entropy does not give more information than standard Kullback-Leibler (KL) entropy.

Apart from this, we do not agree with any of the other claims raised in the preceding comment (and we still have some discrepancy in details with the numerical results shown in  whose origin is unclear to us).

The first main point of ref.  was to show that the “renormalized entropy” (RE) proposed in  and applied to EEG data in  was indeed a KL entropy, but taking an unusual “renormalized” reference. We maintain, in contrast to claims made in the comment, that this relation (Eq. (9) in ) was not mentioned in , and not in  either. Indeed, due to it, the condition $`T_{eff}1`$ postulated in  is not needed to obtain the inequality $`\mathrm{\Delta }H0`$. The fact that the latter was claimed in  to hold only for $`T_{eff}1`$ indicates that the authors were not aware of the relation to KL (or “relative”) entropy. Apart from this, we also wanted to give a simple treatment free of all allusions to statistical thermodynamics, the latter making the treatments in  hard to understand.

Our second point was that RE is very unlikely to be more useful than the usual (un-renormalized) KL entropy for the analysis of EEGs from epileptic patients, as claimed in . On the one hand this was based on the numerical similarity between RE and standard KL entropies, which is enforced by several inequalities and which makes it unlikely a priori that either is superior. On the other hand, we verified this explicitly by detailed numerical calculations which indeed showed that both behaved very similar. It is clear from Fig. 1 that major differences are due to the choice of the reference window. In contrast to what is suggested in the comment, we did not base this conclusion entirely on theoretical arguments.

Finally, we also stressed that the condition $`T_{eff}1`$ – which is not needed at all – has led the authors of  to choose reference points which are physiologically very unfortunate. Again, we remark that it is very unreliable to compare a relative entropy measure obtained from EEG recordings at different electrodes by using different references (from the pre-seizure stage, during the seizure, or from the post-seizure stage) for each electrode in order to localize an epileptic focus. Thus the failure of realizing that RE is a sort of KL entropy – or at least of drawing the obvious consequences from this observation – has hampered its application to EEGs.


text:
# New methods to measure phase transition strength WJ and DJ were partially supported by the EC IHP network HPRN-CT-1999-000161: Discrete Random Geometries: From Solid State Physics to Quantum Gravity, and DJ and RK were partially supported by an Enterprise Ireland/British Council Research Visits Scheme. RK also wishes to thank the TrinLat collaboration for its hospitality during an extended stay at Trinity College Dublin.

## 1 INTRODUCTION

Lattice regularization of a quantum field theory renders the system a statistical mechanical one. The issues of phase transitions and their properties become, therefore, of central importance. While a true phase transition can only occur for a system of infinite extent, the non-perturbative computational approach to lattice field theory and statistical physics accesses only systems of limited size. The partition function of such a system can be written as a polynomial in an appropriate temperature-like or field-like variable and the complex zeroes of such a polynomial encode all of the information on the behaviour of standard thermodynamic quantities.

Traditional statistical mechanical techniques involving partition function zeroes are mainly confined to analysing the zeroes closest to the real axis, as these are the strongest contributors to the critical or pseudocritical behaviour. However, a full understanding of the critical properties of the infinite-size system requires knowledge of the density of zeroes too. It has long been known that the density of zeroes for a *finite*-size system would provide a lucrative source of information but a reliable technique for the extraction of this quantity from numerical data proved elusive.

Recently, however, some of us have succeeded in providing just such a technique . The general idea is to focus on the integrated density of zeroes rather than directly on the density itself, which, for a finite system, is a series of delta-functions. This new method has been seen to be quite reliable and robust and is applicable to phase transitions of the temperature-driven  and field-driven  types and to transitions of first and higher order. For a comparison between this and other approaches, see .

For systems hitherto analysed, the zeroes of the partition function had two special properties, which seem to be common to the bulk of standard models encountered in statistical physics. These are (i) the zeroes are all simple zeroes (zeroes of order one) and (ii) they lie on a curve called the *singular line*, which impacts on to the real axis at the phase transition point.

The question now arises as to the generality of the techniques developed in . Here, we show how the methods can indeed be extended to systems for which some or all of the zeroes occur in degenerate sets and/or they are not confined to a singular line, but instead form a two-dimensional pattern in the complex plane. Such two-dimensional patterns of zeroes have been observed in various lattice field theory and statistical physics models . Here, we demonstrate the extended technique by application to the case of free Wilson fermions. The zeroes of this lattice field theory display the two new features we wish to address.

The partition function for a lattice of finite extent, $`L`$, is $`Z_L(z)_j\left(zz_j(L)\right)`$, where $`z`$ is an appropriate coupling parameter. In the case where the zeroes, $`z_j`$, are on a singular line impacting on to the real axis at the critical point, $`z_c`$, they can be parameterised by $`z_j=z_c+r_j\mathrm{exp}(i\phi )`$. The density is then defined as $`g_L(r)=L^d_j\delta (rr_j(L))`$. The cumulative distribution function of zeroes is $`G_L(r)=_0^rg_L(s)𝑑s`$, which is $`j/L^d`$ if $`r(r_j,r_{j+1})`$. At a zero one assumes the cumulative density is given by the average

$$G_L(r_j)=(2j1)/2L^d.$$
(1)
In the thermodynamic limit, for a second-order transition, the integrated density is, in fact ,

$$G_{\mathrm{}}(r)r^{2\alpha },$$
(2)
where $`\alpha `$ is the usual critical exponent associated with specific heat. Standard finite-size scaling emerges quite naturally from this approach .

## 2 GENERAL DISTRIBUTIONS OF ZEROES

A departure from smooth linear sets of zeroes was found in 1984 when it was shown that for anisotropic two-dimensional lattices there can exist a two-dimensional distribution (area) of zeroes . Since then, a host of systems have been discovered with this feature. A common characteristic of all two-dimensional distributions of zeroes is that they cross the physically relevant real axis at only one point, in the thermodynamic limit, corresponding to the phase transition.

For such two-dimensional distributions, the density of zeroes in the infinite-volume limit has been shown to be

$$g(x,y)=y^{1\alpha n}f\left(\frac{x}{y^n}\right),$$
(3)
where $`(x,y)`$ give the location of zeroes with the critical point as the origin. Integrating out the $`x`$-dependence yields

$$g(y)=_{x_1}^{x_2}g(x,y)𝑑xy^{1\alpha },$$
(4)
where $`x_1`$ and $`x_2`$ mark the extremities of the distribution of zeroes at a distance $`y`$ from the $`x`$ axis. Integrating again, to determine the cumulative density of zeroes at the point $`r`$ in the $`y`$-direction, yields an expression identical to (2). Thus, the strength of the transition, as measured by $`\alpha `$, can be determined by similar methods to those previously used. Rather than counting the zeroes along the singular line, one now counts them up to a line $`y=r`$ within the two-dimensional domain they inhabit.

The second new feature we wish to accommodate is the existence of degeneracies in the set of zeroes. If a number of zeroes coincide, $`G_L`$, as defined in (1), is multivalued and is no longer a proper function. A more appropriate density function is determined as follows. Suppose, in general, that $`z_j=z_{j+1}=\mathrm{}=z_{j+n1}`$ are $`n`$-fold degenerate. It is easy to convince oneself that the densities to the left and right of an actual zero are

$$G_L(r)=\{\begin{array}{cc}j+n1\hfill & \text{for }r(r_{j+n1},r_{j+n})\hfill \\ j1\hfill & \text{for }r(r_{j1},r_j).\hfill \end{array}$$
(5)
The density at the $`n`$-fold degenerate zero, $`r_j`$, is again sensibly defined as an average:

$$G_L(r_j)=\frac{1}{L^d}\left(j+\frac{n}{2}1\right).$$
(6)
This is the most general formula for extracting the density of any distribution of zeroes and deals with two-dimensional spreads and degeneracies.

## 3 APPLICATION

We wish to test the above technique in a situation with the two new features of interest – namely a two-dimensional set of degenerate zeroes. The free Wilson fermion model is an ideal testing ground and its zeroes in two dimensions are easily generated exactly . The zeroes for a system of size $`L=40`$ are depicted in Fig. 1 in the complex $`z=1/2\kappa `$ plane. Here $`\kappa `$ is the hopping parameter.

We generated zeroes for a wide range of lattices, and their distributions, as given by (6), are plotted in Fig. 2. The data collapses onto a universal curve which goes through the origin, indicating that (6) is indeed a suitable form for the density of zeroes and indicating the occurence of a phase transition. Fits close to the origin yield an exponent compatible with the expected value, $`\alpha =0`$. The error estimates appropriate to such a fit are non-trivial and we leave their discussion to a separate publication .

## 4 CONCLUSIONS

A new method to extract the (continuous) density of zeroes from (discrete) finite-size data has been extended to deal with the case of two-dimensional distributions of zeroes and systems in which the zeroes occur in degenerate sets. The method has been demonstrated in an application to the free Wilson fermion model and seen to be capable of direct determination of the strength of the phase transition as measured by the critical exponent $`\alpha `$.


text:
# Structure and Dynamics of Candidate O Star Bubbles in N44

## 1 Introduction

Mechanical feedback from massive stars is a fundamental driver of galaxy evolution. The supernovae (SNe) and supersonic stellar winds from these stars generate bubbles and shells in the interstellar medium (ISM), which may dominate the ISM structure formation in many star-forming galaxies. The standard, adiabatic model for these shells and superbubbles (e.g., Weaver et al., 1977) predicts the production of hot, $`10^610^7`$ K, low-density gas within these shells, which bear the products of SN and massive star nucleosynthesis. Therefore, these bubble structures are thought to be the source of the diffuse, hot ionized component of the ISM. Vigorous star-forming regions may generate supergiant shells that blow out of galaxy disks, thereby dispersing metals and mass far from the parent star formation event. In addition, this mechanical energy contributes to ISM kinematics and is thought to be a primary source of turbulence (e.g., Norman & Ferrara, 1996; Goldman, 2000).

The effects of mechanical feedback have been clearly demonstrated, for example, in starburst galaxies (Johnson et al., 2000; Strickland et al., 2000) and in superbubbles around OB associations (e.g., Saken et al., 1992; Oey, 1996). However, several troubling discrepancies remain unsolved. For example, almost all clean examples of bubbles blown by Wolf-Rayet stars and of superbubbles present shells that are too small for their observed parent stellar population (e.g., Treffers & Chu, 1982; Drissen et al., 1995). Many superbubbles also show anomalous kinematics and X-ray emission (e.g., Rosado et al., 1982; Chu & Mac Low, 1990; Wang & Helfand, 1991). To resolve these puzzles and more clearly understand the mechanical feedback process, it is therefore necessary to examine simpler systems.

In active star-forming regions, high concentrations of OB stars collectively produce superbubbles. But in their peripheries, where stars are loosely distributed, discrete bubbles may be produced by individual stars. H$`\alpha `$ images of H ii regions ionized by OB associations in the Large Magellanic Cloud (LMC) have revealed large (50–150 pc) gas shells around the OB associations and numerous small ($`<`$ 15 pc), ring-like nebulae outside the large shell (see, e.g., Davies, Elliott, & Meaburn, 1976). These small nebulae are candidate wind-blown bubbles of single O or B stars (Weis & Duschl, 1999), and offer perhaps the simplest test of the standard, adiabatic bubble model.

We have investigated the kinematics and physical nature of three ring-like nebulae in N44 (Henize, 1956) in the LMC. N44 is a bright H ii complex with a 44 pc$`\times `$67 pc superbubble at the center and several compact H ii regions along the shell rim. In the surrounding area, small dense H ii regions as well as extended diffuse nebulae are also present at distances up to 160 pc. Figure 1 shows these features and their identifications given by Henize (1956). Three OB associations exist in N44 (Lucke & Hodge, 1970): LH 47 in the central shell, LH 48 in the compact H ii region N44 I, and LH 49 partially embedded in the bright H ii region N44 D. The three nebulae that we have studied are N44 F, N44 J, and N44 M. N44 F is a bright, circular H ii region at the northwest rim of the superbubble and is adjacent to a bright filament connected to the main nebula. In contrast, N44 J and N44 M, appearing respectively on the northern and eastern parts of N44, are fainter and relatively isolated. None of these three nebulae have been studied since their original identification by Henize (1956).

We have obtained ground-based CCD images and high-dispersion echelle spectra of N44 F, N44 J, and N44 M to study the morphology and kinematics of these nebulae. Additionally, we have used Hubble Space Telescope (HST) WFPC2 images of N44 F in the H$`\alpha `$ and \[S ii\] lines to analyze its ionization structure. In this paper, we first describe these observations in §2, then present the analysis and results for N44 F in §3 and N44 J and N44 M in §4. We present our conclusions in §5.

## 2 Observations

The datasets used in this study include: (1) CCD images taken with the Curtis Schmidt and 0.9 m telescopes at Cerro Tololo Inter-American Observatory (CTIO), (2) CCD images taken with the HST WFPC2, and (3) high-dispersion echelle spectra taken with the CTIO 4 m telescope. The high-dispersion spectra are useful to diagnose expanding shells that are not morphologically identifiable.

### 2.1 Ground-based CCD Images

CCD images of N44 were taken with the Curtis Schmidt telescope at CTIO in 1993 December. The detector was a front-illuminated Thomson 1024$`\times `$1028 CCD camera with 19 $`\mu `$m pixels, or 1$`\stackrel{}{\mathrm{.}}`$835 pixel-1. Images were obtained with H$`\alpha `$ ($`\lambda _c`$=6561Å, $`\mathrm{\Delta }\lambda `$=26Å), \[O iii\] ($`\lambda _c`$=5010Å, $`\mathrm{\Delta }\lambda `$=50Å), and \[S ii\] ($`\lambda _c`$=6720Å, $`\mathrm{\Delta }\lambda `$=50Å) filters for exposure times of 1800 s, 2700 s, and 2700 s, respectively. The H$`\alpha `$ filter includes not only the H$`\alpha `$ line but also the neighboring \[N ii$`\lambda \lambda `$6548,6583 lines. At the LMC’s radial velocity, $``$300 km s-1, the \[N ii$`\lambda `$6583 line is red-shifted further away from the red edge of filter, but the \[N ii$`\lambda `$6548 line is red-shifted further into the H$`\alpha `$ bandpass. Using the \[N ii\]/H$`\alpha `$ intensity ratio measured from our echelle spectra (see Section 2.3), we estimate that the \[N ii$`\lambda `$6548 line emission contributes 2.4%, 4.3% and 3.3% of the total flux of N44 F, N44 J and N44 M, respectively. This \[N ii\] contamination has been corrected in our H$`\alpha `$ flux measurements.

Higher-resolution H$`\alpha `$ images were taken with the 0.9 m telescope at CTIO in 2000 December. The detector was a thinned, AR coated SITe 2048$`\times `$2048 CCD with 24 $`\mu `$m pixels, or 0$`\stackrel{}{\mathrm{.}}`$4 pixel-1. The H$`\alpha `$ filter ($`\lambda _c`$=6559Å, $`\mathrm{\Delta }\lambda `$=64Å) used transmits H$`\alpha `$ and both \[N ii\] lines, therefore the “H$`\alpha `$” images are really H$`\alpha `$+\[N ii\] images, and thus are only used for morphology analysis. The angular resolution, determined by the FWHM of the stellar images in the field of view, was approximately 1$`\stackrel{}{\mathrm{.}}`$3. To map the entire N44 complex, five fields were observed with exposure times of 2$`\times `$600 s for each field.

All CCD images were bias-subtracted and flat-fielded in a standard way using IRAF software, and multiple frames were combined to remove cosmic rays and improve the S/N ratio. Figure 1 shows the H$`\alpha `$ mosaic of the N44 complex, and Fig. 2 presents close-ups of N44 F, N44 J and N44 M in the H$`\alpha `$, \[O iii\], and \[S ii\] lines.

### 2.2  HST WFPC2 Images

HST WFPC2 images of N44 F were acquired via the Hubble Observation Problem Report 111The original observations of N44 were made with an incorrect orientation, resulting in a loss of 1/3 of the field of view. To compensate the loss, HOPR observations were granted. To facilitate the scheduling of these HOPR observations and to maximize the scientific yield, a different pointing center in N44 suitable for arbitrary orientations was chosen. (HOPR) observations of the Cycle 6 program 6698. In these observations, made with the $`F656N`$ filter for 2$`\times `$500 s and the $`F673N`$ filter for 2$`\times `$600 s, N44 F was included in the Wide-Field Camera 3 (WF3). The $`F656N`$ filter, centered at 6563.7 Å with a FWHM of 21.4 Å, includes the H$`\alpha `$ line and the neighboring \[N ii\] lines. Using the \[N ii\]/H$`\alpha `$ intensity ratio measured from our echelle spectra (see Section 2.3), we estimate that the \[N ii$`\lambda `$6548 line emission contributes 2.4% of the flux in N44 F. This \[N ii\] contamination has been corrected in our H$`\alpha `$ flux measurements. The $`F673N`$ filter, centered at 6732.2 Å with a FWHM of 47.2 Å includes the nebular lines of \[S ii$`\lambda \lambda `$ 6716,6731.

The calibrated WFPC2 images were produced by the standard HST pipeline. We processed them further with IRAF and STSDAS routines. The images taken with the same filter were combined to remove cosmic rays and to produce a total-exposure map. The combined \[S ii\] and H$`\alpha `$ images were then each corrected for the intensity- and position-dependent charge transfer efficiency as explained by Holtzman et al. (1995). Following the recommended procedures for narrowband WFPC2 photometry222available at http://www.stsci.edu/instruments/wfpc2/Wfpc2\_faq/wfpc2\_nrw\_phot\_faq.html, we find the rectangular widths of the $`F656N`$ and $`F673N`$ filters to be 28.3 Å and 63.3 Å, respectively, and use these widths to convert the flux densities into fluxes.

Figure 3 shows the WFPC2 images of N44 F. To show the excitation variations, the \[S ii\]/H$`\alpha `$ ratio map of N44 F is also presented, in addition to the H$`\alpha `$ and \[S ii\] images.

### 2.3 CTIO 4 m Echelle Spectra

Echelle spectra were obtained with the CTIO 4 m telescope in 2000 December. The spectrograph was used with a 79 line mm-1 echelle grating and the long-focus red camera. We observed a single order centered on the H$`\alpha `$ line by inserting a post-slit H$`\alpha `$ interference filter and replacing the cross-dispersing grating with a flat mirror. A Tek 2048$`\times `$2048 CCD with 24 $`\mu `$m pixel-1 was used to record the images. This provided a spectral sampling of $``$0.08 Å pixel-1 and a spatial sampling of 0$`\stackrel{}{\mathrm{.}}`$26 pixel-1. The wavelength coverage, limited by the H$`\alpha `$ filter and the echelle order, was 125 Å; the spatial coverage, limited by the optics, was $``$200′′. The angular resolution, determined by the FWHM of the seeing, was approximately 1″. A slitwidth of 1$`\stackrel{}{\mathrm{.}}`$65 was used. The resultant spectral resolution, measured from the Th-Ar lamp lines, was about 13 km s-1 FWHM.

Each of the three nebulae, N44 F, N44 J, and N44 M, was observed with an east-west oriented slit. The exposure times were 600 s for N44 F, but only 300 s for N44 J and N44 M. The spectral lines detected include nebular H$`\alpha `$ and \[N ii\] lines, and telluric H$`\alpha `$ and OH lines (Osterbrock et al., 1996). We have adopted rest wavelengths of 6562.7885 Å, and 6548.080 Å and 6583.454 Å for the H$`\alpha `$ and \[N ii$`\lambda \lambda `$6548,6583 nebular lines, respectively (Spyromilio, 1995).

## 3 Physical Structure of N44 F

We will first describe the overall morphology of N44 F, examine its excitation and ionization structure, and derive its rms electron density. Special emphasis is given to the prominent dust “pillars” where star formation might be taking place. Next, we describe the kinematic features detected in the echelle spectra, fit Gaussian components to the line profiles, and use the resultant velocities to identify expanding structures and to determine expansion velocities. Finally, we combine the morphological and kinematic information and present an integrated view of this H ii region. Throughout this paper, we will adopt a distance of 50 kpc to the LMC (Feast, 1999).

### 3.1 Morphology and Density

N44 F is a bright circular nebula nestled on a bright filament at the northwest rim of the superbubble in N44. In low-resolution H$`\alpha `$ images, N44 F appears to have a double-ring structure with a very bright region in the northwest. The details of these features are resolved in the HST WFPC2 images. The outer ring, of radius $``$ 22″, is surrounded by dusty features. Finger- or pillar-like dust columns, similar to those in the Eagle Nebula (=M16; Hester et al., 1996), are also seen at the periphery, with two especially prominent ones at the west, close to the brightest region of N44 F. These dust columns are revealed in detail by the HST images, and will be discussed more extensively in Section 3.2. The inner ring, of radius $``$ 6″, is resolved into diffuse emission without sharp features. The western part of this ring merges into the brightest region of N44 F. A group of bright stars are located within the central ring, but they are substantially offset towards the bright emission region at the northwest. Unfortunately, none of the stars in the vicinity of N44 F have been studied photometrically or spectroscopically.

The \[S ii\] images of N44 F show a very different morphology: only the outer ring is visible and it is mostly composed of sharp filaments. The pillar-like features are also present and stars at the tips of two of them are clearly seen (see Fig. 4). The \[S ii\]/H$`\alpha `$ ratio in the central region is very low, about 0.1. The outer ring possesses a higher \[S ii\]/H$`\alpha `$ ratio, with values ranging from 0.2 to 0.4. Enhanced \[S ii\]/H$`\alpha `$ ratios are often associated with shell structures or high velocity shocks. As we do not find the latter in the nebula (see Section 3.3), the somewhat elevated \[S ii\]/H$`\alpha `$ ratio is likely to result from the low ionization parameter333The ionizing flux derived further can be used to estimate the ionization parameter $`U(Q(H^0)nf^2)^{1/3}`$, where $`Q(H^0)`$ is the ionizing flux in photons s-1, $`n`$ the density in cm-3, and $`f`$ the filling factor. Using $`Q(H^0)`$= 1.0$`\times `$1049 photons s-1, $`n`$=40 cm-3 and $`f`$=$`1\left(1\frac{\mathrm{\Delta }R}{R}\right)^3`$0.6, we find a value of $`10^3`$ for the ionization parameter. caused by the shell morphology. The values of \[S ii\]/H$`\alpha `$ are in fact similar to those seen in other shell structures (e.g., Lasker, 1977; Hunter, 1994).

Using the Curtis Schmidt images, we also examine the \[O iii\] morphology of N44 F. Because of the poorer resolution, only the overall morphology can be described. The \[O iii\] image of N44 F looks very similar to that in H$`\alpha `$, except for the bright H$`\alpha `$ outer ring, which disappears almost completely in \[O iii\]. We have used the HST WFPC2 observations of N44 C (Garnett, Galarza, & Chu, 2000) to calibrate the observed \[O iii\]/H$`\alpha `$ ratio derived with the Curtis Schmidt images. The \[O iii\]/H$`\alpha `$ ratio of N44 F is relatively high, about 0.5–1.2 inside the ring, with the highest value near the star(s). The outer ring has the lowest \[O iii\]/H$`\alpha `$ ratio, about 0.25. The integrated \[O iii\]/H$`\alpha `$ ratio for the whole nebula is about 0.55. The high ratio near the central region indicates that the ionizing source must be an early-type massive star.

The WFPC2 H$`\alpha `$ image can be further analyzed to estimate the rms density of N44 F. A color excess of E(B$``$V)=0.11 mag for N44 (Oey & Massey, 1995) is adopted in the extinction correction. The emission measure (EM) derived from the extinction-corrected surface brightness (Peimbert, Rayo, & Torres-Peimbert, 1975) ranges from 7$`\times `$103 cm-6 pc in the faintest interior regions of N44 F to 4$`\times `$104 cm-6 pc at the tips of the “pillars”. To derive the density from the EM, we can assume two simple geometries. First, the density can be estimated by assuming the depth to be equal to the diameter of the nebula: for the central region, we find a value ranging from 25 cm-3 to 35 cm-3. Second, N44 F is limb-brightened and we can therefore estimate the density by assuming a constant density spherical shell with radius $`R`$ and thickness $`\mathrm{\Delta }R`$. This geometry predicts a unique radial surface brightness profile for any $`\mathrm{\Delta }R/R`$ based on the varying path length through the shell material. Comparing these simple modeled profiles to the observed H$`\alpha `$ surface brightness profiles near the shell rim, we find a shell thickness $`\mathrm{\Delta }R/R`$ of $``$0.25. Using the peak surface brightness and the length of the longest line of sight, $`2\sqrt{2R\mathrm{\Delta }R(\mathrm{\Delta }R)^2}`$, we then derived a density in the shell of 40 cm-3. On average, the diffuse gas of N44 F probably possesses a density of a few tens cm-3.

Using the Curtis Schmidt H$`\alpha `$ image and the flux calibration from Kennicutt & Hodge (1986), we measured over the entire N44 F nebula an H$`\alpha `$ flux of 3.6$`\times `$10-11 erg cm-2 s-1. For a distance of 50 kpc, the H$`\alpha `$ luminosity of N44 F is then 1.4$`\times `$1037 erg s-1, after correcting for an extinction of E(B$``$V)=0.11 mag. This H$`\alpha `$ luminosity requires an ionizing flux of $`Q(H^0)`$=1.0$`\times `$1049 photons s-1. A second evaluation of the extinction correction has been made by comparing radio data and the Curtis Schmidt H$`\alpha `$ image: a mean value E(B$``$V) of 0.28 mag for N44 F, N44 J, and N44 M is found (J. R. Dickel, private communication), and this results in a $`Q(H^0)`$ of 1.4$`\times `$1049 photons s-1 for N44 F. This value suggests that the total ionizing power of the star(s) within N44 F is equivalent to a single O6.5 V–O7.5 V star (Panagia, 1973; Vacca, Garmany, & Shull, 1996; Schaerer & de Koter, 1997). This is compatible with the rather high \[O iii\]/H$`\alpha `$ ratio: the photoionization models of Stasińska & Schaerer (1997) for CoStar stellar atmosphere models suggest that the observed integrated \[O iii\]/H$`\alpha `$ ratio of N44 F could be produced by the equivalent of an O8 V–O9 V star444We have compared N44 F, N44 J, & N44 M to the closest model available in Stasińska & Schaerer (1997), i.e. a spherical nebula of density 10 cm-3 and metallicity 0.25 Z, ionized by one single star.. The comparison to the photoionization models should however be interpreted with caution because the models do not assume a shell structure, which affects the excitation and ionization structure. The dominant spectral type could therefore be somewhat earlier than inferred from the modeled excitation.

### 3.2 Pillar-like Dust Columns

The dark clouds of N44 F exhibit new examples of the intriguing finger- or pillar-like structures, similar to those reported in the Eagle Nebula (Hester et al., 1996). The most prominent two are embedded in the bright northwest region (see Figs. 3 and 4). Several additional features of this type are present in the south and possibly in the southwest, but they are fainter. In Fig. 4, we present the most striking of these structures in the H$`\alpha `$ and \[S ii\] images as well as in the ratio map. The dust columns in N44 F are 5″ to 10″ long and 1″ to 2 ″ wide, corresponding to 1 to 2 pc long and 0.25 to 0.5 pc wide at the LMC’s distance. These sizes and morphologies are similar to those seen in 30 Doradus (Scowen et al., 1998) and M16 (Hester et al., 1996). Two of these columns harbor stars at their tips, as indicated by arrows in Fig. 4. Hester et al. (1996) suggested that stars at the tips of pillars could be young stellar objects which were formed in the molecular cloud and finally emerge after the photoevaporation of most of their parent cloud. However, the Hester et al. (1996) model of pre-existing condensations in the molecular cloud does not provide the unique conditions for producing photoionized pillars. Hydrodynamical simulations by Williams, Ward-Thompson, & Whitworth (2001) illustrate that these structures can be formed by a variety of initial conditions in the structure with inhomogeneities in density and/or radiation field, but current observations do not yet permit discrimination among these different models. It is therefore of great interest to identify additional examples of such objects.

Other important characteristics of these dust columns are their ionization structure and their density. All of the dust columns point towards the central ionizing source, i.e., a bright star or a small cluster between the center and the northwest edge of N44 F. As can be seen in Fig. 5, the \[S ii\] peak is separated from the H$`\alpha `$ peak by approximately 0$`\stackrel{}{\mathrm{.}}`$1, i.e., 0.024 pc, with H$`\alpha `$ peaking slightly closer to the source of UV radiation which is photoevaporating the cloud. The shape of these profiles is very similar to those observed in the pillars of 30 Doradus, or in ground-based images of M16 (Scowen et al., 1998). Unfortunately, no WFPC2 \[O iii\] images of N44 F are available, thus we can not undertake a complete analysis using photoionization models or constrain further the properties of the ionizing source. Following the method used in section 3.1, we estimate an rms density larger than 300 cm-3 at the tip of one of the pillars. This density is compatible with that determined by Scowen et al. (1998) in 30 Doradus, 700 cm-3. These pillars are thus photoevaporating clouds, whose evaporated matter enhances the density and brightness of the diffuse material seen in the northwest part of the nebula.

### 3.3 Kinematic Properties

The kinematic properties of N44 F are derived from the high-dispersion echelle observations of the H$`\alpha `$ and \[N ii$`\lambda `$6583 lines. The \[N ii\] line has a smaller thermal width than the H$`\alpha `$ line, thus it shows the line splitting better although it is noisier. We have extracted H$`\alpha `$ and \[N ii\] line profiles and fitted them with Gaussian components. The resultant velocity-position plot (Fig. 6) shows that N44 F is an expanding shell with a radius of 22″, coincident with the outer ring described in §3.1. Both lines indicate a shell expansion velocity of $``$12 km s-1. An example of the \[N ii\] and H$`\alpha `$ profiles are presented in Fig. 7. The approaching (blueshifted) component is fainter, suggesting a lower density. Further east, outside N44 F, an additionnal blue component is present, but it is probably associated with the complex filamentary structure at the outskirts of N44 and is not directly related to N44 F.

### 3.4 Test of the Standard Bubble Model

N44 F has a classical ring-nebula morphology, and its kinematics confirms an expanding shell, most likely created by the wind of the early-type massive star(s) within this bubble. We can now test whether the observed expansion dynamics of this bubble are consistent with those expected from the standard, adiabatic bubble model (e.g., Weaver et al., 1977). With the observed radius $`R`$ in units of pc and expansion velocity $`V`$ in units of km s-1, we can compute the dynamical timescale $`t_6`$ in units of 106 yr using $`t_6=0.6R/V`$, and further determine the ratio of stellar wind luminosity $`L_{36}`$ in units of 1036 ergs s-1 to ambient density $`n_0`$ in units of H-atom cm-3, using $`R=27(\frac{L_{36}}{n_0})^{1/5}t_6^{3/5}`$.

For a radius of 5.3 pc and an expansion velocity of $``$12 km s-1, N44 F has a dynamical timescale of 2.7$`\times `$105 yr and a $`L_{36}/n_0`$ ratio of 0.016 ergs s-1 cm3. If we assume that the shell consists of mass swept up from the volume, the ambient density, $`n_0`$, can be estimated by $`n_0=3n\frac{\mathrm{\Delta }R}{R}`$, where $`n`$ is the rms density of the shell, $`\mathrm{\Delta }R`$ the shell thickness, and $`R`$ the shell radius. For $`n`$40 cm-3 and $`\frac{\mathrm{\Delta }R}{R}=0.25`$ (see Section 3.1), we derive an ambient density of $``$30 cm-3 and a wind luminosity of $``$5$`\times `$1035 ergs s-1. This wind luminosity is compatible with that expected from the spectral type determined in Section 3.1, i.e., an O7V star (de Jager et al., 1988; Prinja, Barlow, & Howarth, 1990). One uncertainty in these estimates is the shell thickness; however, even if the shell thickness is off by a factor of two, the predicted wind luminosity will still be compatible, within the errors, with an O7V star.

Thus, the dynamics of N44 F appear to be consistent with the standard adiabatic model for bubble evolution, within uncertainties generated by the complex interstellar environment. This contrasts with the results from superbubbles, which show evidence of retarded growth (e.g., Oey, 1996). To investigate further this discrepancy, it is worth considering an alternative bubble model.

The presence of dust pillars around N44 F clearly indicates a highly inhomogeneous ISM. The photoevaporation of the dense clumps may serve as an energy sink and govern the dynamic evolution of this H ii region. Since photoevaporation processes take place in the early evolutionary stages of a bubble/superbubble, the energy loss by photoevaporation may explain the smaller-than-expected sizes of previously studied superbubbles and Wolf-Rayet ring nebulae (see §1). Therefore, we will now consider models of bubbles formed in a cloudy medium by McKee, Van Buren, & Lazareff (1984). In their models, only bubbles blown by stars with weak winds will follow the Weaver et al. model, while stars with stronger winds will produce bubbles whose expansion is governed by the dynamics of photoevaporating clouds. They define a dimensionless wind luminosity $`L_W^{}L_W/(1.26\times 10^{36}(S_{49}^2/n_m)^{1/3})`$ and a homogenization radius $`R_h(\mathrm{pc})0.009t^{4/7}\left(\frac{S_{49}}{n_m^2}\right)^{1/7}`$, where $`L_W`$ is the wind luminosity in erg s-1, $`S_{49}`$ the ionizing flux in photons s-1, $`n_m`$ the mean cloud density in cm-3 and $`t`$ the dynamical timescale in years. Strong winds possess $`L_W^{}1`$ and blow bubbles of radius $`R_h`$.

Using the observed parameters of N44 F (with $`n_m`$=$`n_0`$) and assuming a stellar wind luminosity for an O7V star, we estimate $`L_W^{}`$1 and $`R_h`$4.5 pc, close to the observed bubble’s radius. This seems to suggest that the dynamics of photoevaporation actually governs the bubble’s expansion. The agreement of the data with both models could be a simple coincidence, or maybe the bubble has just encountered the cloudy medium and begun the photoevaporation processes, in which case there may not be a large deviation from the Weaver et al. model. Note, however, that all these conclusions hinge upon the validity of the stellar wind luminosity, which is based on an assumed spectral type of the central star. Future study of the actual stellar content of N44 F is needed to place stronger constraints on the applicable bubble models.

## 4 Physical Structure of N44 J and N44 M

### 4.1 N44 J

Situated in the northernmost part of the N44 complex, N44 J is a slightly elongated, ring-like nebula with a size of 37″$`\times `$32″(or 9.0 pc$`\times `$7.8 pc). The 0.9 m H$`\alpha `$ image shows a bright rim around dusty features at the northwest boundary. A bright central star is visible in the H$`\alpha `$ and \[S ii\] images (see Fig. 2). In H$`\alpha `$ and \[S ii\], the nebula exhibits modest limb-brightening, with the edge-to-center brightness ratio $`<`$2, while its central part displays a rather uniform surface brightness. The \[S ii\]/H$`\alpha `$ ratio is very similar to that of N44 F, ranging from 0.25 to 0.4. In contrast, the \[O iii\] image of N44 J shows emission only near the central star. The \[O iii\]/H$`\alpha `$ ratio is 0.4 near the star, and only 0.05–0.08 elsewhere, significantly lower than the corresponding value in N44 F. The integrated \[O iii\]/H$`\alpha `$ ratio of N44 J is about 0.12.

The observed H$`\alpha `$ flux of N44 J is 5.7$`\times `$10-12 erg cm-2 s-1. This requires an ionizing flux $`Q(H^0)`$ of 1.6$`\times `$1048 photons s-1 when adopting E(B$``$V)=0.11 mag, or 2.3$`\times `$1048 photons s-1 for E(B$``$V)=0.28 mag. These values are expected for an O9V–B0V type ionizing star (Panagia, 1973; Vacca, Garmany, & Shull, 1996; Schaerer & de Koter, 1997), which is fully consistent with that suggested by the \[O iii\]/H$`\alpha `$ ratio of the nebula (Stasińska & Schaerer, 1997). Assuming a spherical geometry, we derive an rms electron density of $``$25 cm-3, slightly less dense than N44 F. The ionized gas mass of N44 J is then $``$200 M.

The echelle spectrum of N44 J shows a kinematically quiescent H ii region. The centroids of the H$`\alpha `$ and \[N ii\] lines are centered at 292 km s-1, with an rms variation of 3 km s-1 (see Figs. 6 and 7). No line splitting is seen. The observed FWHMs are 27$`\pm `$1 km s-1 for H$`\alpha `$ and 16$`\pm `$2 km s-1 for \[N ii\]. These observed line widths imply a turbulent FWHM of about 10 km s-1 for H$`\alpha `$ and 8 km s-1 for \[N ii\] if we take into account the instrumental profile and the thermal broadening at 104 K. The turbulence can thus explain the small variations of the centroid velocity. The difference between the turbulence widths is probably caused by a difference in the emission regions of the two lines. If N44 J is actually expanding, the expansion velocity should be very small, $`<`$8 km s-1. For such small expansion velocity, no strong compression is expected, hence no sharp filaments should be seen.

Finally, we may compare the observations with the standard, adiabatic bubble model. As the density in the ring-like nebula is a few tens cm-3, we will assume in our calculations a density of 10 cm-3 in the ISM. For an expected wind luminosity of $``$3$`\times `$1034 erg s-1 for an O9.5V star (de Jager et al., 1988; Prinja, Barlow, & Howarth, 1990) and with the observed 17$`\stackrel{}{\mathrm{.}}`$5 radius of N44 J, the Weaver et al. (1977) model predicts an expansion velocity of about 8 km s-1 for N44 J. This expected expansion velocity is just at the limit implied by the velocity widths of the H$`\alpha `$ and \[N ii\] lines. Raising the ambient density by a factor of 2 can lower the expected expansion velocity to 6.5 km s-1, which is more consistent with the observed limit. We conclude that the observed bubble dynamics of N44 J does not greatly contradict the model predictions.

### 4.2 N44 M

On the southeast part of N44, another small ring-like nebula is present: N44 M. It is the faintest nebula of our sample. It shows a ring with a radius of $``$22″ (or 5.3 pc), but its north boundary is not well defined (see Fig. 2). The sharp southwest boundary suggests the presence of dust clouds. Several bright stars appear in its interior. N44 M possesses an excitation structure very similar to N44 J. The \[S ii\]/H$`\alpha `$ ratio of N44 M is comparable to that in N44 F, but the \[O iii\]/H$`\alpha `$ ratio of N44 M is much lower, ranging from 0.07 to 0.16 in the nebula and rising to 0.4 near the stars. The integrated \[O iii\]/H$`\alpha `$ ratio of N44 M is about 0.18.

With an observed H$`\alpha `$ flux of 5.3$`\times `$10-12 erg cm-2 s-1 and the two evaluations of the extinction, the required ionizing flux $`Q(H^0)`$ is 1.5–2.1$`\times `$1048 photons s-1. This suggests an O9-B0V type ionizing star (Panagia, 1973; Vacca, Garmany, & Shull, 1996; Schaerer & de Koter, 1997), which is again compatible with the \[O iii\]/H$`\alpha `$ ratio (Stasińska & Schaerer, 1997). If we assume a spherical geometry, we find that the rms density is about $``$17 cm-3 and the estimated mass of the nebula $``$275 M.

The echelle spectrum of N44 M shows only one component, centered at 298 km s-1 with an rms variation of 1 km s-1, and a FWHM of 28$`\pm `$1 km s-1 in the H$`\alpha `$ line and 20$`\pm `$4 km s-1 in the \[N ii\] line (see Figs. 6 and 7). The turbulent FWHM are then 12 km s-1 for H$`\alpha `$ and 14 km s-1 for \[N ii\]. With an observed radius of $``$22″, a density of 10 cm-3 for the ISM and a O9.5V ionizing star (see Section 4.1.), the standard bubble model predicts an expansion velocity of $``$7 km s-1 for N44 M, which appears consistent with the upper limit of the echelle observations.

The kinematic characteristics of N44 J and N44 M are very similar, despite the somewhat different morphology. Neither H ii region contains shells expanding at velocities greater than $``$8 km s-1. In the future, additional observations may enable a better comparison of these nebulae with the standard model. Precise photometry and/or stellar spectroscopy would give accurate informations about the stellar content of these nebulae, while nebular spectra with higher resolution should more strongly constrain the kinematics of N44 J and N44 M. Thus, improved data is needed to more stringently test the standard model for N44 J and N44 M.

## 5 Conclusion

Massive stars inject a large amount of energy into the ISM through their ionizing radiation, fast stellar winds, and SNe. The stellar winds, which dominate the early mechanical feedback, are able to sweep the ambient medium into an expanding shell and create a cavity around the parent stars. H ii regions ionized by massive stars are therefore thought to evolve from compact H ii regions to ring-like nebulae.

We have studied three candidate wind-blown bubbles in the N44 complex: N44 F, N44 J, and N44 M. They all have diameters of $``$10 pc. Their ionizing fluxes and \[O iii\]/H$`\alpha `$ ratios suggest that their ionizing stars have spectral types of roughly O7V, O9.5V and O9.5V, respectively. N44 F has an expansion velocity of 12 km s-1, while the other two have expansion velocities $`<`$ 8 km s-1. The stellar wind luminosity implied by the ionizing stars and the observed bubble dynamics appear consistent with that expected in the standard, adiabatic bubble model within the limits of observational uncertainties. Similar results were found by Oey & Massey (1994) for two bubbles around single O star in M33. This is in sharp contrast to the anomalous kinematics seen in superbubbles and Wolf-Rayet bubbles. We caution, however, that our observed expansion velocities may not provide the most stringent constraints when comparing with models, because the expected expansion velocities of N44 J and N44 M are comparable to or smaller than the convolution of the isothermal sound velocity ($``$10 km s-1) and the instrumental FWHM (13 km s-1). Deeper, higher-resolution spectroscopic observations of the \[N ii\] or \[O iii\] lines are needed to pin down the expansion velocity for a more definitive test of models.

We also report the discovery of several dust pillars in N44 F, which are similar to those found in the Eagle Nebula (Hester et al., 1996) and 30 Doradus (Scowen et al., 1998). These structures are indicative of the evaporation of molecular clouds by the ionizing source. The examples in N44 F are especially well-suited to investigating the process and role of photoevaporation in star formation and shell evolution, because of their location within such a well-defined, wind-blown bubble. We found that photoevaporation, as evidenced by the pillar structure, may be dynamically important for the shell evolution in N44 F, slowing down the expanding shell and inhibiting its growth. These photoevaporation processes may thus undermine the seeming agreement between the observed bubble dynamics and the simple adiabatic bubble model.

We thank Charles Danforth for assisting in the echelle observations and Rosie Chen for help in the HST WFPC2 image reduction. This work has been supported by the NASA grant STScIGO-6698.01-95A. Y.N. acknowledges support from the F.N.R.S., contracts P4/05 and P5/36 ‘Pôle d’attraction Interuniversitaire’ (SSTC-Belgium) and the PRODEX XMM-OM and Integral Projects. R.C.S. acknowledges the support of the Dean B. McLaughlin fellowship at the University of Michigan and NSF grant AST 95-30747.


text:
# Structural Determinant of Protein Designability

## Abstract

Here we present an approximate analytical theory for the relationship between a protein structure’s contact matrix and the shape of its energy spectrum in amino acid sequence space. We demonstrate a dependence of the number of sequences of low energy in a structure on the eigenvalues of the structure’s contact matrix, and then use a Monte Carlo simulation to test the applicability of this analytical result to cubic lattice proteins. We find that the lattice structures with the most low-energy sequences are the same as those predicted by the theory. We argue that, given sufficiently strict requirements for foldability, these structures are the most designable, and we propose a simple means to test whether the results in this paper hold true for real proteins.

Successful protein design relies in part on knowledge of how a polypeptide chain’s native structure varies as a function of its amino acid sequence. Yet even if this “protein folding problem” were solved in its entirety, the would-be protein designer would still face the formidable task of finding those sequences which folded to the target structure he or she wished to engineer B.Dahiyat and S.Mayo (1997); Shakhnovich (1998). It is therefore vital to know how many solutions to this search of sequence space exist for a given target structure, i.e. how designable the target structure is A.V.Finkelstein et al. (1993); Li et al. (1996); Shakhnovich (1998); Kussell and Shakhnovich (1999); Buchler and Goldstein (2000); Koehl and Levitt (2002). The question of what makes a particular protein structure more designable also bears fundamentally on our understanding of how proteins first evolved.

Past studies of designability have been limited in large part by their lack of generality. The previous contributions in Li et al. (1996); Kussell and Shakhnovich (1999), for example, rely heavily on their study of two-letter monomer alphabets and Cartesian lattice-polymers, and can therefore offer no clear implications for the 20-letter, off-lattice world of real proteins Shakhnovich (1998). In contrast, studies like that of Koehl and Levitt Koehl and Levitt (2002) come closest to probing real protein designabilities, but have no theoretical foundations from which to extrapolate beyond their numerical results. Finally, a number of investigations have assumed that the distribution of amino acid sequence energies is either nearly Tiana et al. (2001) or totally Buchler and Goldstein (2000) independent of the structural topology of the target fold, a premise which is flatly contradicted by the findings in Kussell and Shakhnovich (1999). The intention of this communication is therefore to identify a theoretically motivated, generally applicable quantitative measure of structural topology which we expect to be a good predictor of designability.

In this Letter, we develop an approximate analytical theory of the spectrum of possible monomer sequence energies for a heteropolymer in any given conformation. Our theory leads us to identify a novel topological property of a conformation, its contact trace, as a determinant of the shape of the conformation’s sequence spectrum. We then confirm the predicted usefulness of the contact trace by running Monte Carlo simulations under conditions which are more realistic than those for which the contact trace was originally identified. Finally, we use the results of these simulations to suggest a connection between the contact trace and designability, and we propose a way to test this hypothesis on real protein structures.

As a preliminary, we must draw a distinction between what we term “strong designability” and “weak designability”. We define a structure’s “strong designability” in the same way that previous studies have defined its “designability”: as the total number of amino acid sequences that fold to that structure. Thereotical arguments based on the random energy model suggest that in order for a sequence to be foldable, its native state energy must lie below a certain low-energy threshold $`E_c`$ Shakhnovich and Gutin (1990a); Bryngelson and Wolynes (1987); A.V.Finkelstein et al. (1995); Kussell and Shakhnovich (1999); Sali et al. (1994), which represents the lowest energy range accessible to misfolded conformations. We therefore define a structure’s “weak designability” to be the fraction of its sequences which lie below a low energy cutoff $`E_c`$. It is important to note that though the strong and weak versions of designability are set apart a priori, it is possible that they will turn out to be the same thing in practice. Indeed, we will argue later on that if natural thermodynamic and kinetic requirements for folding are sufficiently stringent, then a structure’s strong and weak designability become essentially indistinguishable.

Our first aim is to determine whether contact topology effects a structure’s distribution of possible sequence energies, and to this end we derive a closed-form sequence partition function for a special class of amino acid alphabets. We begin by considering a polymer of $`N`$ monomers, where each monomer can be one of $`2M`$ different possible kinds. We may construct for any polymer configuration an $`N\times N`$ traceless contact matrix $`C`$ whereby lattice nearest neighbors that are not sequence neigbors are considered to be in contact Shakhnovich and Gutin (1990b). Next, we represent the amino acid type of the $`i`$th monomer in the chain as a $`2M`$ dimensional unit-vector $`\stackrel{}{s}^{\text{ }(i)}=[0,\mathrm{},0,1,0,\mathrm{},0]`$, where the non-zero vector element is the $`k`$th element of the vector if the $`i`$th monomer is of type $`k`$.

We define the Hamiltonian to be a standard nearest-neighbor contact potential, i.e.

$$=\frac{1}{2}\underset{i,j}{\overset{N,N}{}}C_{i,j}\stackrel{}{s}^{\text{ }(i)}(B\stackrel{}{s}^{\text{ }(j)})$$
(1)
where $`B`$ is the $`2M\times 2M`$ matrix of interaction energies for the different pairs of monomer types. A closed form expression for the partition function of this Hamiltonian may be obtained if we let $`B`$ take on the special form

$$B=\left[\begin{array}{ccccc}V_{1,1}& V_{1,1}& V_{1,2}& V_{1,2}& \mathrm{}\\ V_{1,1}& V_{1,1}& V_{1,2}& V_{1,2}& \mathrm{}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \end{array}\right]$$
(2)
In other words, $`B`$ is the direct product of a matrix $`V`$ with an Ising-like potential. Though we have now constrained the nature of our contact potential, we still retain unlimited freedom to choose the types of interactions represented in $`V`$.

Now let us define the $`M`$-vector $`\stackrel{}{\sigma }`$ through $`\sigma _{k+1}=s_{2k+1}s_{2k+2}`$. Our Hamiltonian becomes

$$=\frac{1}{2}\underset{i,j}{\overset{N,N}{}}C_{i,j}\stackrel{}{\sigma }^{\text{ }(i)}(V\stackrel{}{\sigma }^{\text{ }(j)})$$
(3)
where $`\stackrel{}{\sigma }^{(i)}`$ is a vector of unit length whose single non-zero element may be either $`1`$ or $`1`$.

At this point we define the $`MN\times MN`$ block matrix $`𝐔`$ through $`𝐔_{M(i1)+k,M(j1)+l}=C_{i,j}V_{k,l}`$. In other words, we turn every element of $`C`$ into an $`M\times M`$ block which couples the vectors $`\{\stackrel{}{\sigma }^{\text{ }(i)}\}`$ through $`V`$ wherever there is a contact. Finally, if we write $`𝝈[\stackrel{}{\sigma }^{\text{ }(1)},\mathrm{},\stackrel{}{\sigma }^{\text{ }(N)}]`$ then the Hamiltonian can be expressed as

$$=\frac{1}{2}𝝈(𝐔𝝈)$$
(4)
This form of the Hamiltonian will allow us to calculate the sequence space partition function, for any contact map and any potential matrix $`V`$. To do so, we employ a continuous-spin approximation, allowing each vector $`\stackrel{}{\sigma }^{\text{ }(i)}`$ to swing anywhere on the $`M`$-dimensional unit sphere instead of restricting it to one of the $`2M`$ available unit-integer lattice points. This approximation not only smears the discrete sequence spectrum into a continuous one, but also distorts the spectral width by altering the relative sizes of $`M`$ and $`M1`$ letter sequence spaces, thereby preventing the theory from making quantitatively accurate predictions. Our assumption is that the model nevertheless retains important information about the *effect* that variations in contact topology have on the shape of the spectrum. The partition function now becomes

$$Z(\beta )=d^{MN}𝝈\left(\underset{i=1}{\overset{N}{}}\delta \left(|\stackrel{}{\sigma }^{\text{ }(i)}|1\right)\right)\mathrm{exp}\left[\frac{\beta }{2}𝝈(𝐔𝝈)\right]$$
(5)
Assuming hereafter that $`M1`$, we have

$$d^M\stackrel{}{\sigma }\mathrm{exp}\left[\frac{(M1)}{2}|\stackrel{}{\sigma }|^2\right]d|\stackrel{}{\sigma }|\text{ }\delta (|\stackrel{}{\sigma }|1)$$
(6)
Defining $`𝐌=(M1)𝐈`$, where $`𝐈`$ is the $`MN\times MN`$ identity matrix, and $`𝐮=𝐌^1𝐔`$ we may write

$$Z(\beta )=d^{MN}𝝈\mathrm{exp}\left[\frac{1}{2}𝝈(𝐌+\beta 𝐔)𝝈\right]$$
(7)
and, since $`𝐔`$ is a real symmetric matrix, we finally obtain

$$\frac{Z(\beta )}{Z(0)}z(\beta )=\sqrt{\frac{det[𝐌]}{det[𝐌+\beta 𝐔]}}=\frac{1}{\sqrt{det[𝐈+\beta 𝐮]}}$$
(8)
If we use $`det𝐀=\mathrm{exp}\text{Tr}\mathrm{ln}𝐀`$ to expand $`\mathrm{ln}z`$ to $`𝒪(\beta ^2)`$ and inverse-Laplace transform the resulting approximation of the partition function to obtain a distribution of sequence energies $`n(E)`$, we get

$`n(E){\displaystyle \frac{1}{2\pi }}{\displaystyle e^{i\beta E}z(i\beta )𝑑\beta }\mathrm{exp}\left[{\displaystyle \frac{E^2}{\text{Tr }𝐮^2}}\right]`$ (9)
$`\mathrm{exp}\left[{\displaystyle \frac{E^2}{2N_C\sigma _B^2}}\right]`$
where $`N_C`$ is the number of contacts and $`\sigma _B^2`$ is the variance of the interaction energies in $`B`$. A similar result was obtained in Shakhnovich (1998) by a straightforward expansion of the sequence partition function. Equation (9) says that the Gaussian approximation of the density of states leads to a naïve sequence space random energy model (REM) for $`n(E)`$. It is this REM picture, in which all structures with the same number of contacts have identical sequence spectra, which has been used implicitly in Buchler and Goldstein (2000). We demonstrate below, however, that the consideration of higher order terms in the expansion of $`\mathrm{ln}z`$ can have a profound impact on how we understand designability.

Let us now consider the free energy $`F=\frac{1}{\beta }\mathrm{ln}z`$ of our sequence space partition sum. Defining the matrix $`v=\frac{1}{M1}V`$ and expanding about high “design temperature” (i.e. requiring that$`|\beta \lambda _i|<1`$, for all $`\lambda _i`$ which are eigenvalues of $`𝐮`$), we get

$`F={\displaystyle \frac{\beta }{4}}(\text{Tr }v^2)(\text{Tr }C^2)+{\displaystyle \frac{\beta ^2}{6}}(\text{Tr }v^3)(\text{Tr }C^3)`$ (10)
$`{\displaystyle \frac{\beta ^3}{8}}(\text{Tr }v^4)(\text{Tr }C^4)+𝒪(\beta ^4)`$
Those structures which minimize $`F`$ will be the ones with the greatest number of amino acid sequences that have low energy when threaded onto that structure. Minimization of $`F`$ is therefore a means to maximize weak designability, which we recall is determined by the fraction of low-energy sequences in a structure’s sequence spectrum. Thus, we now consider which choices of the contact matrix $`C`$ would serve to make the free energy $`F`$ as negative as possible. Because $`\text{Tr }C^n`$ is equal to the number of $`n`$-step paths along the contact map which return to their starting place, we know that all such contact traces must be positive. Thus, the exact behavior of the series in (10) will hinge on whether the largest eigenvalues of $`v`$ are positive or negative.

For either type of potential matrix $`v`$, however, we expect that there will be some positive correlation between the trace of an even power of a structure’s contact matrix and the number of low energy monomer sequences in that structure. Furthermore, the dependence of the free energy expansion in (10) on such coarse quantities as the traces of powers of $`v`$ suggests that the impact of the contact matrix on the spectrum of sequence energies should be relatively insensitive to the detailed features of the potential. We therefore determined to empirically test whether the above results remained valid for a discrete monomer alphabet which violated the special form of the potential assumed in (2). We first calculated the contact matrices for all $`103346`$ different compact conformations of 27-mers on a cubic lattice Shakhnovich and Gutin (1990b). Next, we calculated $`E`$ vs. $`T`$ annealing curves for random starting sequences on different structures for a standard Monte Carlo search of sequence space with a move set containing composition-preserving two-monomer and three-monomer permutations. The energy of each sequence was determined using a potential set given by Table 6 of Miyazawa and Jernigan (1985). This set of interactions, where average interactions are subtracted out, is one of the most diverse potentials possible for a 20-letter alphabet, and therefore provides the most general empirical test of the predicted relationship between sequence energies and contact topology. From the annealing curves, we then calculated the entropy in sequence space $`S(E)`$ according to the prescription given by eq. (11) of Shakhnovich (1998).

Figure 1 plots the sequence space entropy difference between low and near-modal energy versus the largest eigenvalue of the structure’s contact matrix for $`86`$ randomly-selected lattice structures. As predicted, the entropy difference between the peak and the left tail decreased as the largest contact matrix eigenvalue increased (corr. $`=0.92`$), indicating that more sequences have low energy in high trace structures. Figure 2 illustrates that the effect observed in Fig. 1 results from global differences in the shapes of the sequence spectra of high trace and low trace structures. The higher the contact trace, the more gradually the number of sequences falls off as energy decreases, and therefore the greater the relative number of sequences of low energy. Clearly, the contact trace of the target structure controls how low in energy a Monte Carlo sequence optimization algorithm running at fixed temperature $`T_{des}`$ will be able to search. The greater the contact trace, the larger the $`S(E)`$ at low energies, i.e. the greater the weak designability.

Interestingly, the most designable 27-mer structures identified using our maximum eigenvalue determinant are similar to the one identified in Li et al. (2002) using random sampling of sequences and a different, ”solvation-like” Miyazawa- Jernigan potential. This attests to a generality of our proposed structural determinant of designability with respect to potentials.

Structures of high contact trace are weakly designable, but are they strongly designable? In order to address this question, we examined the stability of sequences designed on two structures of maximal and minimal contact trace. For each target structure, we determined how many of its designed sequences were “on-target”, that is, had the target structure as their unique energetic ground state determined over all compact conformations, and calculated the gap in energy between that of the ground state and those of the lowest-energy structures with low structural similarity to the ground state conformation. We found that while only $`12`$% of the sequences designed for the low trace structure were on-target, $`24`$% were in the case of the high trace structure. Furthermore, as Figure 3 illustrates, the sequences designed successfully on the high trace structure tended on average to have larger energy gaps than their low trace cousins, consistent with earlier observations that low energy in the native conformation also leads on-average to a larger gap and greater stability Shakhnovich and Gutin (1993); Abkevich et al. (1994)

Fig. 3 therefore gives us a means to unify strong and weak designability. Past studies have suggested that in order for a protein sequence to fold rapidly to its native state, it must have a larger-than-average energy gap Goldstein et al. (1992); Sali et al. (1994); Sato et al. (2001). Fig. 3 demonstrates that if the conditions for folding stably and quickly in nature demand sufficiently high energy gaps, weak designability will be one and the same with strong designability, since only by having very low energy in its target structure will a sequence have an appreciable chance of exhibiting the gaps which are thermodynamically and kinetically required by the natural environment. We speculate that protein evolution under such conditions would lead to the emergence of natural protein folds with near-optimal contact traces.

We have presented an analytical theory which identifies the first instance of a generally applicable, well-defined, numerical measure of a protein structure’s topology which is expected to correlate with the structure’s designability. Using a Monte-Carlo search in the sequence space of a lattice model with $`20`$-letter energetics, we have shown that the implications of the theory extend beyond the special assumptions under which they were first developed. The finding that higher contact trace may lead to greater potential for thermal stability leads us to hypothesize that thermophilic organisms, whose natural environment makes uncommonly strict demands for protein stability, exhibit a genomic bias towards folds of higher contact trace. We have recently found that this bias exists (unpublished results), providing further encouragement that the contact trace may yet offer new insight into the laws governing structural diversity in natural proteins.

We thank Edo Kussel, William Chen, Nikolay Dokholyan, and Lewyn Li for valuable discussions, and NIH and Pfizer Inc. for support.


text:
# NTUTH-02-505E August 2002 Quenched chiral logarithms in lattice QCD with overlap Dirac quarksTalk presented by T.W. Chiu at Lattice 2002. This work was supported in part by National Science Council, R.O.C. under the grant numbers NSC90-2112-M002-021

## 1 Introduction

In quenched chiral perturbation theory (q$`\chi `$PT) , the pion mass to one-loop order reads as

$`m_\pi ^2=Cm_q\{1\delta [\text{ln}(Cm_q/\mathrm{\Lambda }_\chi ^2)+1]\}+Bm_q^2`$ (1)
where $`m_q`$ denotes the bare ($`u`$ and $`d`$) quark mass, $`\mathrm{\Lambda }_\chi `$ is the chiral cutoff which can be taken to be $`2\sqrt{2}\pi f_\pi `$ ( $`f_\pi 132`$ MeV ), $`C`$ and $`B`$ are parameters, and $`\delta `$ is the coefficient of the quenched chiral logarithm.

Theoretically, $`\delta `$ can be estimated to be

$`\delta ={\displaystyle \frac{m_\eta ^{}^2}{8\pi ^2f_\pi ^2N_f}}`$ (2)
where $`m_\eta ^{}`$ denotes the $`\eta ^{}`$ mass in q$`\chi `$PT, and $`N_f`$ is the number of light quark flavors. For $`f_\pi =132`$ MeV, $`N_f=3`$, and $`m_\eta ^{}=\sqrt{𝐦_\eta ^{}^2+m_\eta ^22m_K^2}=853`$ MeV, (2) gives

$`\delta 0.176`$ (3)
Evidently, if one can extract $`\delta `$ from the data of $`m_\pi ^2`$, then $`m_\eta ^{}`$ in q$`\chi `$PT can be determined by (2).

Besides from the data of $`m_\pi ^2`$, one can also obtain $`\delta `$ via (2) by extracting $`m_\eta ^{}`$ from the propagator of the disconnected hairpin diagram. However, to compute the propagator of the hairpin is very tedious.

Fortunately, with the realization of exact chiral symmetry on the lattice , the quark propagator coupling to $`\eta ^{}`$ is $`(D_c+m_q)^1`$ , thus only the zero modes of $`D_c`$ can contribute to the hairpin diagram, regardless of the bare quark mass $`m_q`$. Therefore one can derive an exact relation  between the $`\eta ^{}`$ mass in q$`\chi `$PT and the index susceptibility of any Ginsparg-Wilson lattice Dirac operator, without computing the hairpin diagram at all. Explicitly, this exact relation reads as

$`(m_\eta ^{}a)^2={\displaystyle \frac{4N_f}{(f_\pi a)^2}}{\displaystyle \frac{(n_+n_{})^2}{N_s}}`$ (4)
where $`N_s`$ is the total number of sites, and $`\chi (n_+n_{})^2/N_s`$ is the index susceptibility of any Ginsparg-Wilson lattice Dirac operator in the quenched approximation. Then (2) and (4) together gives

$`\delta ={\displaystyle \frac{1}{2\pi ^2(f_\pi a)^4}}{\displaystyle \frac{(n_+n_{})^2}{N_s}}`$ (5)
A salient feature of (5) is that $`\delta `$ can be determined at finite lattice spacing $`a`$, by measuring the index (susceptibility) of the overlap Dirac operator, and with $`f_\pi a`$ extracted from the pion propagator.

Now it is clear that, in order to confirm the presence of quenched chiral logarithm in lattice QCD, one needs to check whether the coefficient $`\delta `$ obtained by fitting (1) to the data of $`m_\pi ^2`$, agrees with that (5) from the index susceptibility. This is a requirement for the consistency of the theory, since the quenched chiral logarithm in $`m_\pi ^2`$ (1) is due to the $`\eta ^{}`$ loop coupling to the pion propagator through the mass term (in the chiral lagrangian), thus $`\delta `$ ($`m_\eta ^{}`$) must be the same in both cases. We regard this consistency requirement as a basic criterion for lattice QCD (with any fermion scheme) to realize QCD chiral dynamics in continuum.

Recently, we have determined $`\delta =0.203(14)`$, from the data of $`m_\pi ^2`$, as well as $`\delta =0.197(27)`$ from the index susceptibility, for the $`8^3\times 24`$ lattice . Their excellent agreement suggests that lattice QCD with overlap Dirac quarks indeed realizes the QCD chiral dynamics in the continuum.

In this paper, we extend our studies to other lattices of larger volumes and/or smaller lattice spacings. Our computations are performed with a Linux PC cluster of 30 nodes at the speed $`30`$ Gflops .

## 2 Results

With the Wilson $`SU(3)`$ gauge action and Creutz-Cabbibo-Marinari heat bath algorithm, we generate ensembles of gauge configurations as follows111The 56 configurations on the $`16^3\times 32`$ lattice are retrieved from the gauge connection (http://qcd.nersc.gov/)..

| $`L^3\times T`$ | $`\beta `$ | \# configs. |
| --- | --- | --- |
| $`8^3\times 24`$ | 5.8 | 100 |
| $`10^3\times 24`$ | 5.8 | 100 |
| $`12^3\times 24`$ | 5.8 | 100 |
| $`16^3\times 32`$ | 6.0 | 56 |

For each configuration, quark propagators are computed for 10-12 bare quark masses. Then the pion propagator and its time correlation function $`G(t)`$ are obtained, and $`G(t)`$ is fitted by the usual formula

$`G_\pi (t)={\displaystyle \frac{Z}{2m_\pi a}}[e^{m_\pi at}+e^{m_\pi a(Tt)}]`$
to extract the pion mass $`m_\pi a`$ and the pion decay constant

$`f_\pi a=2m_qa{\displaystyle \frac{\sqrt{Z}}{m_\pi ^2a^2}}`$
The data of $`f_\pi a`$ can be fitted by straight line. Thus taking $`f_\pi a`$ at $`m_qa=0`$ equal to $`132`$ MeV times the lattice spacing $`a`$, then $`a`$ can be determined. Fixing $`\mathrm{\Lambda }_\chi a=2\sqrt{2}\pi f_\pi a`$, we fit (1) to the data of $`(m_\pi a)^2`$, and obtain $`Ca`$, $`\delta `$ and $`B`$. Then $`m_\eta ^{}`$ is computed from (2). Our results are summarized in Table 1. Evidently, all $`\delta `$ values are in good agreement with the theoretical estimate $`\delta 0.176`$ from q$`\chi `$PT. Furthermore, the $`\eta ^{}`$ masses are also in good agreement with the theoretical estimate $`m_\eta ^{}853`$ MeV.

Even though the quenched chiral logarithm may not be easily detected in the graph of $`(m_\pi a)^2`$ vs. $`m_qa`$, it can be unveiled by plotting $`(m_\pi a)^2/(m_qa)`$ vs. $`m_qa`$, as shown in Fig. 1. Further, plotting $`(m_\pi a)^2/(m_qa)B(m_qa)`$ vs. $`\mathrm{log}(m_qa)`$, the presence of quenched chiral logarithm is evident, as shown in Fig. 2.

Note that, in Fig. 1, for the $`16^3\times 32`$ lattice at $`\beta =6.0`$, the minimum of $`(m_\pi a)^2/(m_qa)`$ occurs at $`m_qa=Ca\delta /B0.1`$, which corresponds to $`m_q200`$ MeV for $`a0.095`$ fm (Table 1). Similarly, for $`8^3\times 24`$, $`10^3\times 24`$, and $`12^3\times 24`$ lattices at $`\beta =5.8`$, the minima of $`(m_\pi a)^2/(m_qa)`$ occur at $`m_qa`$ 0.21, 0.15 and 0.16 (i.e., $`m_q`$ 280, 190, and 217 MeV) respectively. Now, except for the $`8^3\times 24`$ lattice which has the smallest volume, the minima of $`(m_\pi a)^2/(m_qa)`$ for the other three lattice sizes almost occur at the same $`m_q200\pm 20`$ MeV. This provides rather substantial evidence that the quenched chiral logarithm parameters extracted from our data of $`m_\pi ^2`$ is indeed genuine, not due to any finite size effects or discretization errors.

Next we measure the index susceptibility of the overlap Dirac operator by the spectral flow method , and then determine $`\delta `$ via (5). The $`\eta ^{}`$ mass in q$`\chi `$PT is computed from (4). Our results are summarized in Table 2, including our earlier results for the $`8^3\times 24`$ lattice. Evidently, each $`\delta `$ is in good agreement with the corresponding $`\delta `$ extracted from $`m_\pi ^2`$ (Table 1), as well as with the theoretic estimate $`0.176`$.

## 3 Conclusion

Our results in this paper provide strong evidences that quenched QCD with overlap Dirac quarks realizes quenched QCD chiral dynamics, as depicted by quenched chiral perturbation theory.