text:
# Saddles and dynamics in a solvable mean-field model

## I Introduction

The landscape paradigm landscape  is based on the idea that the relevant behavior of a given interacting system can be understood by studying the general characteristics of the potential energy surface (PES), i.e. the hypersurface defined by the potential energy as a function of the total configurational degrees of freedom. This way of looking at the dynamics in liquids has received in the last few years new attention due to computational improvement and new theoretical approaches. Within the large variety of studies based on the PES characteristics pes  it has recently pointed out the relevant role played by stationary points (saddles) of the PES in the analysis of the supercooled liquid regime of simple model liquids noi\_sad ; cav\_sad ; doye\_wales ; noi\_jcp . The main result of the saddle-based approach has been the finding of a close relationship between the parameters describing the structural arrest and the saddles visited at a given temperature. More specifically, the saddle order (fraction of negative curvatures -corresponding to negative eigenvalues of the Hessian- at the saddle points) turns out to have a well defined temperature dependence and extrapolates to zero at the mode-coupling temperature $`T_{MCT}`$ mct , where the inverse diffusion coefficient and the structural relaxation time appear to diverge as a power law. This result has been very useful to give a landscape-based interpretation of the dynamic critical temperature for the model liquids investigated: $`T_{MCT}`$ marks a crossover between a saddle-to-saddle dominated dynamics at “high” temperatures and a minimum-to-minimum dominated dynamics below $`T_{MCT}`$. Moreover, to further emphasize the role of the saddles in determining the slow dynamics in liquids, the temperature dependence of saddle order has been shown to be quantitatively related to that of diffusivity noi\_jcp , allowing to give a simple interpretation of the relevant diffusive processes in terms of “open” directions in the $`3N`$ dimensional configuration space. The analysis of saddles has also permitted to obtain important information about the shape of the landscape itself: it has been observed an high regularity in the distances among saddles and in their energy locations. This suggested a surprising simple organization of those points of the PES, i.e. the saddles, that are relevant for the description of the slow dynamics. The observed regularity allows us to introduce simplified models of the landscape reproducing some of the observed PES characteristics. In a previous work noi\_jcp  we explored the features of one of these simple models, the trigonometric model (TM) introduced by Madan and Keyes madan\_keyes , made of $`N`$ independent sinusoidal degrees of freedom. Despite its non-cooperative nature, the trigonometric model seems to capture many important features of the PES of simple liquids.

In this work we use the saddle-approach to analyze another model system that generalizes the TM introducing a cooperative interaction among the different degrees of freedom. More specifically, in the model studied here -that we call $`k`$-trigonometric model ($`k`$TM)- there is a $`k`$-body mean field interaction. Our aim is to extend the validity of the relationship between landscape (saddle properties in our case) and dynamic behavior, besides the model liquids investigated before, to systems where a cooperative behavior exists. The choice of the $`k`$TM has been motivated for two main reasons: it is a generalization of a simple liquid PES model, and it is a new model for testing in an analytical way the ideas underlying the landscape-dynamics relationship. Indeed the thermodynamics of $`k`$TM, as the dynamics (diffusion coefficient), and the saddles visited at a given temperature are analytically computable. This allows us to obtain directly a clear interpretation of the relationship between characteristics of the visited landscape and dynamics in a truly cooperative system, avoiding the computational problems usually encountered in the numerical analysis of PES in model liquids. The main results are: (i) we confirm the existence of a close relationship between saddle order and diffusivity: both have a similar temperature behavior, and approach an asymptotic Arrhenius law at low temperature; (ii) the energy barrier parameter observed in the Arrhenius law of low temperature expansion of diffusivity is the mean elevation $`\delta e`$ of saddles of order one from the minima, and this value grows linearly with the parameter $`k`$ governing the number of interactions in the Hamiltonian. Therefore, our findings confirm the generality of the relationship between landscape (saddle properties) and dynamics, first observed numerically in simple liquids noi\_sad ; cav\_sad ; noi\_jcp  and analytically in spin glasses cav2 .

The paper is organized as follow: in Section II we introduce the model and its characteristics; in Section III the thermodynamics is exactly solved, using a mean-field effective Hamiltonian; in Section IV we explore the properties of saddles sampled by the system at different temperatures; in Section V the diffusion coefficient is determined through a Langevin dynamics, and in Section VI we report the conclusions.

## II The model

The model we study is the mean-field $`k`$-trigonometric model ($`k`$TM), defined by the configurational Hamiltonian:

$`H`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{N^{k1}}}{\displaystyle \underset{i_1,\mathrm{},i_k}{}}[1\mathrm{cos}(\phi _{i_1}+\mathrm{}+\phi _{i_k})]`$ (1)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{N^{k1}}}{\displaystyle \underset{i_1,\mathrm{},i_k}{}}[1\mathrm{}(e^{i\phi _{i_1}}\mathrm{}e^{i\phi _{i_k}})],`$
and by the Langevin dynamics:

$$\gamma \dot{\phi }_i=\frac{H}{\phi _i}+\eta _i,$$
(2)
where $`\{\phi _i\}`$ are angular variables, $`\phi _i[0,2\pi )`$, $`i=1,\mathrm{},N`$, $`\gamma `$ is the friction constant and $`\eta _i`$ is a Gaussian white noise with $`\eta _i(t)=0`$ and $`\eta _i(t)\eta _j(t^{})=2T\gamma \delta _{ij}\delta (tt^{})`$. Defining the variable:

$$z(\{\phi \})=\frac{1}{N}\underset{i}{}e^{i\phi _i}=\xi (\{\phi \})e^{i\psi (\{\phi \})},$$
(3)
the Hamiltonian of the $`k`$TM can be written in a more compact way:

$$H=N\mathrm{\Delta }[1\mathrm{}(z^k)]=N\mathrm{\Delta }(1\xi ^k\mathrm{cos}k\psi ).$$
(4)
The model is invariant under the discrete group $`C_{kv}`$, generated by the transformations: $`\phi _i\phi _i+2\pi /k`$ and $`\phi _i\phi _i`$. It can be considered as a model for a system of bidimensional Heisenberg spins with an interaction that is not invariant under the rotation group but only under the $`C_{kv}`$ group. Depending on the value of the parameter $`k`$ governing the number of interactions, the system manifests a wide spectrum of thermodynamic behavior acprz :

* For $`k=1`$ the $`k`$TM reduces to the trigonometric model introduced by Madan and Keyes madan\_keyes  (the relation between saddles and diffusivity for $`k=1`$ has been investigated in Ref. noi\_jcp ): the system reduces to independent degrees of freedom and all the thermodynamic quantities are smooth functions of temperature.
* For $`k=2`$ a second order phase transition takes place, with a spontaneous breaking of the $`C_{kv}`$ symmetry and a spontaneous magnetization.
* For $`k3`$ the system undergoes a first order phase transition, again with a magnetization breaking the symmetry group.

While in a previous paper we have analyzed the relationship between phase transition points and topological changes in the landscape for this model acprz , here we are interested in the investigation of a link between characteristic of the landscape explored by the system at a given temperature and its dynamic behavior. As the high temperature paramagnetic phase is quite trivial (the effective potential energy is constant and the system can be viewed as an ensemble of free particles), the interesting region for our purpose will be the low $`T`$ one, in which the magnetization is different from zero.

## III Thermodynamics

The thermodynamics of the $`k`$TM is exactly solvable. In a previous paper acprz  we reported the exact derivation of the partition function using the saddle point technique. Here we report a different and more simple derivation based on the introduction of an effective Hamiltonian for each single degree of freedom. Using the standard approach, in order to obtain a mean-field form of the Hamiltonian, we have to make the substitution in Eq. 1

$`e^{i\phi _{i_1}}\mathrm{}e^{i\phi _{i_k}}`$ $``$ $`ke^{i\phi _{i_1}}e^{i\phi _{i_2}}\mathrm{}e^{i\phi _{i_k}}`$ (5)
$``$ $`(k1)e^{i\phi _{i_1}}\mathrm{}e^{i\phi _{i_k}}.`$
We introduce the parameter $`\zeta =e^{i\phi }`$, that has to be determined self-consistently; the mean-field Hamiltonian is then:

$$=\mathrm{\Delta }[1+(k1)\zeta ^kk\zeta ^{k1}\mathrm{cos}\phi ],$$
(6)
where we have chosen $`\zeta `$ to be real, without loss of generality (this corresponds to choose a particular magnetization of the low temperature state). In the canonical ensemble the self consistency equation for $`\zeta `$ turns out to be:

$$\zeta =\mathrm{cos}\phi _{}=\frac{I_1(\beta \mathrm{\Delta }k\zeta ^{k1})}{I_0(\beta \mathrm{\Delta }k\zeta ^{k1})},$$
(7)
where $`\beta =(K_BT)^1`$, $`I_0(\alpha )=(2\pi )^1_0^{2\pi }𝑑\phi e^{\alpha \mathrm{cos}\phi }`$ and $`I_1(\alpha )=I_0^{^{}}(\alpha )`$ are the modified Bessel functions of order 0 and 1 respectively.

For each $`\beta `$ the Eq. 7 gives the thermodynamic value of the mean magnetization $`\zeta (\beta )`$. The value $`\zeta =0`$ always solve Eq. 7, but is a stable solution only at low $`\beta `$ (high temperature). As $`\beta `$ is increased others solutions appear at a given $`\beta `$, and the one with lower free energy is the stable one. In Fig. 1 we report the function $`\zeta (\beta )`$ for $`k=1,2,3`$. For $`k=1`$ the curve is smooth, as no phase transition occurs. For $`k=2`$ there is a singular point at which a second order phase transition takes place, separating an high temperature paramagnetic phase ($`\zeta =0`$) and a low temperature ordered phase, with two possible magnetizations corresponding to $`\psi =0`$ and $`\psi =\pi `$. For $`k=3`$ (and also for $`k>3`$, not reported in the figure) the system undergoes a first order phase transition, corresponding to the appearance of three solutions of Eq. 7: the paramagnetic one $`\zeta =0`$ (always present), the magnetic one $`\zeta >0`$ (with degeneration $`k`$, i.e. $`k`$ different possible value of $`\psi `$) and that corresponding to the maximum of the free energy separating the two previous solutions. For the $`k=3`$ case we report all the solutions of Eq. 7 at different temperatures: the three magnetization values at a given low temperature then correspond to these three different solutions. We report with an arrow in the figure the transition temperature $`T_0=0.63`$, marking the transition point at which the $`k=3`$ model exhibits a first order phase transition (in the figure, as well as in all the other figures, the canonical thermodynamically inaccessible region is represented by dot-dashed lines).

We can write the canonical probability distribution for the variable $`\phi `$ at a given temperature in the form:

$$𝒫(\phi ;\beta )=\frac{e^\beta }{𝒵}=\frac{e^{\beta \mathrm{\Delta }k\zeta ^{k1}\mathrm{cos}\phi }}{2\pi I_0(\beta \mathrm{\Delta }k\zeta ^{k1})},$$
(8)
where $`𝒵`$ is the canonical partition function for the effective Hamiltonian $``$ and the parameter $`\zeta (\beta )`$ is a function of temperature as obtained from Eq. 7.

Using Eq. 8 it is easy to calculate the different static quantities we are interested in. For example, the mean potential energy $`e(\beta )`$ results to be:

$$e(\beta )=\mathrm{\Delta }(1\zeta ^k),$$
(9)
as one can also see directly substituting the mean magnetization value $`\zeta `$ in Eq. 4. In Fig. 2 the temperature dependence of potential energy $`e`$ (full lines in the figures) is shown for $`k=1,2,3`$. For a discussion of the different behaviors see the comments of Fig. 1 reported above.

## IV Saddles

The saddles are the stationary points of the potential energy, satisfying the zero force equation: $`H=0`$. We are interested in locating the saddles of energy $`e_s`$ and in finding their order $`n`$. The latter quantity is defined as the fractional number of negative eigenvalues of the Hessian matrix calculated at the saddle point, i.e. the fractional number of downward curvature directions. For the $`k`$TM it turns out that the saddles located at a given energy level all have the same order $`n`$, and the relationship between $`e_s`$ and $`n`$ is analytical (see Ref. acprz ):

$$e_s(n)=\mathrm{\Delta }[1(12n)^k],$$
(10)
valid for $`n<1/2`$ (or $`e<1`$), condition satisfied in the whole thermodynamically accessible phase space. Similarly to the case of liquids noi\_sad ; cav\_sad , we introduce the concept of saddles visited during the dynamical evolution. To this aim we define a partition of the configurational space, such that each point in $`^N`$ is associated to a saddle point (i.e. we define a basin of attraction for a saddle). The partition of $`^N`$ in basins of saddles of $`H`$ is realized by introducing the pseudo-potential $`W=|H|^2`$ and identifying the basin of attraction of the minima of $`W`$ with that of a saddle of $`H`$. Indeed, $`W`$ has the property that absolute minima correspond to saddles of $`H`$ nota\_w . As for the thermodynamic calculation, we can use also in this case an “effective” pseudo-potential $`𝒲`$ given by:

$$𝒲\left|\frac{}{\phi }\right|^2=\mathrm{\Delta }^2k^2\zeta ^{2(k1)}\mathrm{sin}^2\phi .$$
(11)
The fact that the “effective” pseudo-potential can be obtained differentiating the effective Hamiltonian is not obvious “a priori”. In ref. zampo  we will show that Eq. 11 is correct at small enough temperature, that is the temperature region we are interested in. The search of the minima of $`W`$ goes through a minimization procedure, i.e. through the solution of the steepest descent equation of the form:

$$\{\begin{array}{cc}\dot{\phi }=𝒲=\mathrm{\Delta }^2k^2\zeta (t)^{2(k1)}\mathrm{sin}2\phi ,\hfill & \\ \zeta (t)=\mathrm{cos}\phi (t),\hfill & \end{array}$$
(12)
where now the parameter $`\zeta `$ is a function of time $`\zeta (t)`$, determined by a self consistency dynamic equation (in the above expression $``$ means average over initial conditions). We are interested in the infinite time solutions of Eq. 12, $`\phi (t\mathrm{}|\phi _0)`$, as a function of initial conditions $`\phi _0\phi (t=0)`$. Without explicitly solving Eq. 12, we observe that the sign of $`𝒲`$ at fixed $`\phi `$ does not change during time, due to the fact that the time dependent factor in Eq. 12 is always positive.

This implies that the specific time dependence of $`\zeta `$ does not affect the final point $`\phi (t\mathrm{}|\phi _0)`$ reached from a given initial condition, rather it controls the rapidity of approaching the final point. It it easy to see that the solutions are the following:

$$\phi (t\mathrm{}|\phi _0)=\{\begin{array}{cc}0\hfill & \text{ if }\phi _0(\frac{\pi }{2},\frac{3\pi }{2})\text{ },\hfill \\ & \\ \pi \hfill & \text{ if }\phi _0(\frac{\pi }{2},\frac{3\pi }{2})\text{ }.\hfill \end{array}$$
(13)
The solution $`\phi =0`$ is a minimum of $`𝒲`$ corresponding to a minimum of the effective potential energy (Eq. 6), while the solution $`\phi =\pi `$ is a minimum of $`𝒲`$ that corresponds to a maximum of the effective potential. The mean saddle order $`n`$ as a function of temperature is then obtained as the probability to be in the basin of attraction (defined through the pseudo-potential $`𝒲`$) of the maximum of the potential energy, or, in other words, weighting the function $`\delta (\phi (t\mathrm{}|\phi _0)\pi )`$ over the equilibrium initial conditions:

$`n(\beta )`$ $`=`$ $`{\displaystyle _0^{2\pi }}𝑑\phi _0𝒫(\phi _0;\beta )\delta (\phi (t\mathrm{}|\phi _0)\pi )`$ (14)
$`=`$ $`{\displaystyle _{\pi /2}^{3\pi /2}}𝑑\phi _0𝒫(\phi _0;\beta ).`$
Using Eq. 8 we finally obtain

$$n(\beta )=\frac{1}{2}\left[1\frac{L_0(\beta \mathrm{\Delta }k\zeta ^{k1})}{I_0(\beta \mathrm{\Delta }k\zeta ^{k1})}\right],$$
(15)
where we have introduced the modified Struve function of order $`0`$: $`L_0(\alpha )=2\pi ^1_0^{\pi /2}𝑑\phi \mathrm{sinh}(\alpha \mathrm{cos}\phi )`$. The energy of saddles sampled by the system at different temperatures is then obtained using Eq. 10, in which the temperature dependence is in the order $`n(\beta )`$, as expressed in Eq. 15: $`e_s(\beta )=e_s[n(\beta )]`$.

In Fig. 2 the saddle energies $`e_s`$ (dashed lines) are reported as a function of temperature for $`k=1,2,3`$: qualitatively $`e_s`$ reproduces the shape of the potential energy $`e`$, and it is always below $`e`$, but coincident in the paramagnetic region for $`k2`$. The map $`e_s`$ vs. $`e`$ is shown in Fig. 3, where one observes that, growing the parameter $`k`$, the energy difference between instantaneous configurations and saddles becomes more and more pronounced.

In Fig. 4 the saddle order $`n`$ is reported as a function of temperature for $`k=1,2,3`$. For $`k=3`$ the negative slope branch corresponds to the similar branch in Fig. 1, that is the solution of Eq. 7 which is a maximum of the free energy (in the figure we do not report the paramagnetic branch). In order to visualize better the low temperature regime, in Fig. 5 we plot the saddle order $`n`$ in a semilogarithmic scale as a function of inverse temperature $`\beta `$. Dashed lines in the figure are the low temperature (high $`\beta `$) approximations to Eq. 15:

$$n(\beta 1)\sqrt{\frac{2e^{k1}}{\pi k\beta \mathrm{\Delta }}}e^{\beta \mathrm{\Delta }k},$$
(16)
which corresponds to an Arrhenius behavior.

An important characteristic of saddles is the relationship between energy and order, which is a feature of the energy landscape itself, independent on temperature. In Fig. 6 we report $`e_s`$ vs. $`n`$, as expressed by Eq. 10. The slope of the curves at $`n=0`$ gives information about the energy elevation $`\delta e`$ of saddles of order 1 from underlying minima (in this model all the minima have zero energy), an interesting quantity in the analysis of activated processes in simple liquids. For the $`k`$TM, using Eq. 10, we have

$$\delta e=\frac{de_s}{dn}|_{n=0}=2\mathrm{\Delta }k,$$
(17)
indicating an energy barrier value growing linearly with the parameter $`k`$.

## V Diffusivity

As we want to establish a connection between landscape and dynamics, we now turn our attention to the dynamics of the $`k`$TM. For simple liquids a very interesting relation has been established between saddles and diffusion coefficient, so, also in our case, we focus on this dynamical quantity. In order to determine the diffusivity $`D`$, we calculate the mobility in linear response. We can use again the effective Hamiltonian Eq. 6, as one can show zampo , using standard techniques lc , that the single particle Langevin dynamics is the same if computed with the full Hamiltonian $`H`$ (Eq. 1) or the effective Hamiltonian $``$ (Eq. 6). The final result is then risken  :

$$D(\beta )=(\gamma \beta )^1[I_0(\beta \mathrm{\Delta }k\zeta ^{k1})]^2,$$
(18)
where $`\gamma `$ is the friction in the Langevin equation.

In Fig. 7 the diffusion coefficient $`D`$ (Eq. 18) is reported as a function of inverse temperature $`\beta `$ in a semilogarithmic scale, together with the (Arrhenius) low temperature expansion of Eq. 18:

$$D(\beta 1)\frac{2\pi \mathrm{\Delta }ke^{k1}}{\gamma }e^{\beta 2\mathrm{\Delta }k}.$$
(19)
It is interesting to note that the energy barrier value $`\delta e`$ in the Arrhenius low temperature expansion is exactly the energy elevation of saddles of order $`1`$ from minima (see Eq. 17) $`\delta e=2\mathrm{\Delta }k`$, as expected in the presence of activated processes.

For completeness in Fig. 8 we report the diffusivity $`D`$ against the saddle order $`n`$ in a logarithmic scale. From this plot one can conclude that the slow diffusive dynamics is controlled by the order of the saddles visited at each temperature.

## VI Conclusions

In conclusion we have studied an exactly solvable mean-field model, the $`k`$TM, that generalizes the non-interacting trigonometric model introduced in Ref. madan\_keyes . The useful framework of landscape saddle-approach, introduced in the numerical investigation of simple liquids, has been applied in our case, evidencing how the relation between sampled saddles and dynamics works well also in our model: diffusion coefficient and saddle order have a very similar temperature dependence, both approaching an Arrhenius law at low temperature and having a similar shape in the whole temperature range. The analysis of the $`k`$TM seems to evidence a quite general relationship between dynamics and saddles visited in the potential energy surface, beyond the simple liquid models in which it was firstly observed.

###### Acknowledgements.

We thank L. Cugliandolo and J. Kurchan for pointing out the analytical techniques described in lc  and for very useful and interesting discussions. We also thank F. Sciortino for carefully reading the manuscript. We acknowledge financial support from INFM and MURST COFIN2000.


format: text

text:
# 1 Averaged correlation function 𝐶^{𝑥𝑥}(𝐿/2) as a function of the system size 𝐿 in a log-log-plot for W=0, 0.25, 0.5, 0.75, 1.0 (top to bottom). The data are averaged over 50000 for 𝐿≤1024, 3000 for L=2048 and 500 for L=4096, the statistical error is smaller than the symbol sizes. The data for the pure system (𝑊=0) follow 𝐶^{𝑥𝑥}(𝐿/2)∝𝐿^{-1/2}, the full line with slope -2 is the expected asymptotic behavior according to the IRFP scenario [2]. Inset: Scaling plot according to Eq. (3) for the data of the main figure with 𝜉_𝑊=600, 140, 54, 20 for 𝑊=0.25, 0.5, 0.75 and 1.0, respectively

Comment on ”Disorder Induced Quantum Phase Transition in Random-Exchange Spin-1/2 Chains”

In a recent Letter Hamacher, Stolze and Wenzel  studied the disordered spin-1/2 AFM XXZ chain

$$H=J\underset{i=1}{\overset{L1}{}}\left[\lambda _i(S_i^xS_{i+1}^x+S_i^yS_{i+1}^y)+\mathrm{\Delta }S_i^zS_{i+1}^z\right]$$
(1)
with $`\lambda _i`$ i.i.d. random variables uniformly distributed over the interval $`[1W;1+W]`$, the parameter $`W`$ controlling the strength of the disorder. The authors claim that for $`\mathrm{\Delta }<1`$ they found numerical evidence for non-universal behavior for weak disorder ($`W<1`$) manifested in a continously varying exponent $`\eta (W)`$ describing the asymptotic decay of the transverse spin correlations

$$C^{xx}(r)=S_i^xS_{i+r}^xr^{\eta (W)},$$
(2)
where $`\mathrm{}`$ denotes the ground state expectation value averaged over the disorder and the sites $`i`$. They concluded that there is no universal infinite randomness fixed point (IRFP) as predicted by D. Fisher .

In this comment we show that these conclusions are wrong and that the numerical data presented in  are not in the asymptotic regime. This misconception is due to the fact that the authors ignored the existence of a $`W`$-dependent crossover length scale $`\xi _W`$ that describes the crossover from the pure fixed point to the only relevant IRFP: For $`L\xi _W`$ one observes the critical behavior of the pure system ($`\lambda _i=const.`$), and only for $`L\xi _W`$ the true asymptotic critical behavior ($`\eta (W)=2`$, independent of $`W>0`$) of the disordered chain becomes visible. Even for strong disorder ($`W=1.0`$) $`\xi _W`$ is of the same order of magnitude as the system sizes considered in .

In order to be able to reach sufficiently large system sizes we restrict ourselves to $`\mathrm{\Delta }=0`$ in which case (1) reduces to a free fermion model and the ground state computations are done following . In Fig. 1 we show the averaged bulk correlation function $`C^{xx}(L/2)`$ for different strengths of the disorder. We observe that asymptotically (i.e. for $`L\mathrm{}`$) the data follow the behavior $`C^{xx}(L)L^2`$ as predicted by the real space renormalization group . Only for small $`L`$ the data appear to follow a non-universal (i.e. $`W`$-dependent) power law, and this is the region on which  reports.

Due to the presence of a crossover length scale $`\xi (W)`$ the correlation function obeys the scaling form

$$C^{xx}(L/2)=L^{1/2}\stackrel{~}{c}(L/\xi _W)$$
(3)
where $`\stackrel{~}{c}(x)\mathrm{const}`$ for $`x0`$, and $`\stackrel{~}{c}(x)x^{3/2}`$ for $`x\mathrm{}`$. This implies $`C^{xx}(L/2)L^{1/2}`$ for $`L\xi _W`$ (the pure behavior) and $`C^{xx}(L/2)L^2`$ for $`L\xi _W`$ (the IRFP behavior). In the inset of Fig. 1 we show such a scaling plot of the data in the main figure. We have chosen $`\xi _W`$ for $`W=1`$ such that the crossover region is centered around $`L/\xi _W1`$, the other estimates for $`\xi _W`$ are then chosen to give the best data collapse. We see that for all disorder strength the maximum system sizes used in  are still well within the crossover region and not in the asymptotic regime. We do not expect that this situation will improve for $`\mathrm{\Delta }>0`$. In general $`\xi _W`$ diverges when $`W0`$, and we found that our estimates for $`\xi _W`$ obey $`\xi _W\delta _W^\mathrm{\Phi }`$ where $`\delta `$ is the ($`W`$-dependent) variance of the random variable $`\mathrm{ln}\lambda _i`$, which determines according to  the RG-flow of the model (1). We obtain $`\mathrm{\Phi }1.24\pm 0.02`$.

To conclude we have shown that contrary to what has claimed in  there is no non-universal behavior in model (1) for weak disorder and that by taking into account the existence of a crossover length scale the numerical data, ours as well as theirs, for $`\mathrm{\Delta }<1`$ are compatible with the IRFP scenario predicted in .

Nicolas Laflorenciea and Heiko Riegerb
a Laboratoire de Physique Quantique;   IRSAMC; Université Paul Sabatier;   31062 Toulouse; France;  b Theoretische Physik; Universität des Saarlandes;   66041 Saarbrücken; Germany

Date: 2 December 2002
PACS numbers: 75.10.Jm, 05.70.Jk; 75.40.Mg


text:
# Long-term spectroscopic monitoring of the Luminous Blue Variable HD 160529 Based on observations collected at the European Southern Observatory at La Silla, Chile (proposals 69.D-0378, 269.D-5038)

## 1 Introduction

A-type hypergiants and Luminous Blue Variables (LBVs), which often show spectra similar to A-hypergiants, are the visually brightest stars in galaxies. The early A-hypergiant HD 160529, which has long been known to exhibit the spectroscopic signatures of extreme luminosity (Merrill & Burwell 1943), has therefore long been considered as one of the visually brightest stars of the Galaxy, comparable to the brightest stars in external galaxies.

The variability of HD 160529 was already established 30 – 40 years ago through studies of its variable light and its spectrum at different epochs, which indicated changes in its spectral type. The spectral variability of HD 160529 was studied extensively by Wolf et al. (1974), who found variations of the Balmer line profiles, radial velocity variations with an amplitude of about 40 km s<sup>-1</sup> and line splitting of some metal lines. Intensity variations of absorption lines of the order of 20 % and considerable photometric variations were also found. The light variability was studied and described in detail by Sterken (1977) and Sterken et al. (1991).

HD 160529 was classified as a new galactic LBV by Sterken et al. (1991) due to its brightness decrease of 0.5 mag from 1983 to 1991 and an apparent spectral type change from A 9 to B 8. Apart from the long-term changes, the authors found pulsation-like variations in their photometric data with a quasi-period of 57 days and peak-to-peak amplitudes of 0.1 mag in $`b`$ and $`y`$, while previous analyses suggested a possible 101.3 day period (Sterken 1981). More recent photometry has been published by Manfroid et al. (1991, 1995) and Sterken et al. (1993, 1995).

From comparison with HD 269662 (=R 110), a close photometric and spectroscopic counterpart of HD 160529 in the Large Magellanic Cloud, Sterken et al. (1991) derived an absolute visual magnitude of $`M_\mathrm{V}=8.9`$ and a distance of 2.5 kpc. They estimated the stellar parameters characterizing the phase of maximum visual brightness to be $`T_{\mathrm{eff}}\mathrm{8\hspace{0.17em}000}`$ K, $`\mathrm{log}g\mathrm{\hspace{0.17em}0.55}`$, $`R_{}330`$ R and $`M_{}13`$ M. The derived mass is in good agreement with the mass determined from evolutionary tracks (Lamers et al. 1998).

Sterken et al. (1991) concluded that HD 160529 is located near the lower luminosity limit of the LBV instability strip and possibly in an evolutionary phase after the Red Supergiant state, i.e. in a post-RSG evolutionary phase.

HD 160529 is of much lower luminosity than more typical LBVs such as AG Car (Stahl et al. 2001) and shows much smaller variations of about 0.5 mag in the visual. According to van Genderen (2001), it is at the low amplitude limit of the “strong-active” S Dor variables. The empirical relation between luminosity and amplitude of LBVs (Wolf 1989) is in agreement with this. According to the period-luminosity relation of Stothers & Chin (1995) a “period” of well above 20 yr would be expected. This is in rough agreement with the light curve published by Sterken et al. (1991).

HD 160529 can be considered an intermediate case between normal supergiants and LBVs. The closest counterpart may be the LBV HD 6884 (=R 40) in the SMC (Szeifert et al. 1993) or the LBV HD 269662 (=R 110) in the LMC (Stahl et al. 1990). The variability of HD 160529 is therefore of interest both in connection with typical LBVs such as AG Car (Stahl et al. 2001) and with normal late B – early A supergiants (Kaufer et al. 1996a, b).

Detailed spectroscopic studies of LBVs, which cover both the short and and long time scales, are still very rare and only available for AG Car (Stahl et al. 2001). Studies of more objects with different physical parameters are needed in order to distinguish the general behaviour of LBVs from peculiarities of single objects.

## 2 Observations

### 2.1 Spectroscopy

We have observed HD 160529 with several different echelle spectrographs and several telescopes from 1991, January to 2002, June. We have monitored HD 160529 with high time resolution over a time span of 2 – 4 months in four successive years, i.e. 1992–1995. In addition, a number of snapshot observations have been taken to extend the total time base covered. The epochs of observations and instruments used are summarized in Table 1.

Most of the observations were obtained with the fiber-linked echelle spectrographs Flash (Mandel 1988) and its modified version Heros (Stahl et al. 1996). These instruments were mainly used at the ESO 50cm telescope, but for a few shorter runs also at the ESO 1.52m and 2.2m telescopes. In the Flash configuration, a setting of the spectrograph was used that allowed the observation of the H$`\delta `$ and Hei$`\lambda `$6678 lines with one exposure, i.e. the spectral range covered the wavelength range from about 4 050 Å to 6 780 Å. The Flash data have been published on CD-ROM (Stahl et al. 1995). Heros is a modified version of Flash, where a beam splitter is used to divide the beam after the echelle grating into two channels, each with its own cross-disperser, camera and CCD detector. The red channel is identical to the Flash instrument. The two spectral ranges of the Heros configuration cover the range from 3 450 Å to 5 560 Å in the blue channel and from 5 820 Å to 8 620 Å in the red channel. The spectral resolution is about 20 000 for both Flash and Heros. The signal-to-noise ratio (S/N) of our spectra strongly depends on the wavelength and is lowest in the blue spectral range. For HD 160529 we typically used an exposure time of one to two hours for Flash and Heros with the ESO 50cm telescope. In good conditions, a S/N-ratio of at about 100 is reached in the red spectral range.

Further observations with significantly better S/N and spectral resolution were carried out with Feros (Kaufer et al. 2000) at the ESO 1.52m telescope between 1999 and 2002. The Feros spectra obtained at the ESO 1.52m telescope cover the spectral range from about 3 600 Å to 9 200 Å in 39 orders. The spectral resolution of Feros is about 48 000. The S/N depends strongly on wavelength and is highest in the red spectral region. The typical exposure times were about 10 minutes.

In addition, a few spectra have been obtained with the Caspec echelle spectrograph at the 3.6m ESO telescope and the Ucles echelle spectrograph at the 3.9m AAT telescope. The Caspec spectra cover the spectral range from 3 800 Å to 5 400 Å with a spectral resolution of about 20 000. With Ucles, observations with three grating positions were obtained and the resulting spectrum covers the total spectral range from 3 600 to 6 880 Å. The spectral resolution is about 50 000.

For all instruments, a built-in high-temperature incandescent lamp and a Th-Ar lamp were used to obtain flat-field and wavelength-calibration exposures, respectively.

In addition, we obtained on Sep. 24, 2001, low resolution spectra (3 Å/pix in the wavelength range from 4 000 to 8 000 Å) with the Dfosc focal reducer spectrograph at the Danish 1.54m telescope at ESO, La Silla. In this case, a He-Ne lamp was used for the wavelength calibration. The Dfosc data are of much lower resolution than the echelle data and allow to study the equivalent width variations only. Since the spectra were obtained in non-photometric conditions, they can also not be flux-calibrated.

The Ucles spectra were reduced with IRAF. All other spectra were reduced with ESO-Midas. The Caspec spectra were reduced with the standard echelle package of ESO-Midas (Ponz & Brinks 1983). For the Flash and Heros spectra a modified version of this package (Stahl et al. 1993) was used. The Feros package (Stahl et al. 1999), also running within ESO-Midas, was used for the reduction of the Feros spectra. The Dfosc spectra have been reduced with the long-slit package of ESO-Midas.

In the years 1992–1995, the star was observed at least once per week and in some seasons up to once per night. In these runs, the sampling is sufficiently dense to cover variations with time scales of a few days.

### 2.2 Photometry

Differential photometry was carried out in the framework of the LTPV project (Sterken 1983) using the ESO 50cm telescope (PMT photometry), the Dutch 90 cm telescope (CCD) and the Danish 1.54m telescope (CCD) at ESO, La Silla. Zero point corrections were applied but no proper standardization was attempted. The results based on the measurements in the Strömgren $`y`$ band are displayed in Fig. 1.

## 3 Variability

### 3.1 Changes in spectral type

The large photometric variations of LBVs are typically accompanied by changes of the spectral type, which can change from Of/WN- to A-hypergiant (Stahl 2001). For HD 160529, spectral types between B8Ia and A9Ia have been published in the literature. However, the spectrum of HD 160529 is peculiar and, since the spectral type depends on the exact classification criteria used, it is not clear if the differences in the published spectroscopic classifications are due to real variations. Most authors find a spectral type of A2–A3 (Wolf et al. 1974), only two authors find a strongly deviating type: Houk (1982) finds A9Ia and Sterken et al. (1991) infer spectral type variations from B8 to A9.

It should be noted, that the spectral types published by Sterken et al. (1991) are not based on spectroscopic criteria, but they are inferred from two arguments:

* The observed color variations variations indicate temperature changes from 8 000 – 10 000 K. This assumes that the theoretical colors apply to the extreme supergiant HD 160529.
* The visual amplitude of 0.5 mag is due to a change in the bolometric correction at constant bolometric luminosity. This assumes that the bolometric corrections (Schmidt-Kaler 1982) for supergiants apply to HD 160529.

The classification (A9Ia) of Houk (1982) is therefore the only spectroscopic classification strongly deviating from the A2–A3 range. The remark of Houk (1982): “H lines seem even weaker than in standards; partly filled in?” seems to indicate that this classification is partly based on the weakness of the Balmer lines, which are indeed influenced by emission.

We therefore re-evaluated the evidence for strong changes in the spectral type of HD 160529. From our spectra we infer only minor changes in spectral type for the spectra obtained until 1999. For this period, we determined the spectral type of HD 160529 by comparing with the supergiant stars $`\beta `$ Ori (B8Ia), HD 96919 (B9Ia), HD 92207 (A0Ia) and HD 100262 (A2Ia) from the sample of Kaufer et al. (1997). The match of most lines is reasonably good with an approximate spectral type of A2Ia. The spectral type as determined from the strength of the metal lines is certainly always later than A0 and earlier than A5Ia. Since most of the metal lines are forming partly in the expanding atmosphere of HD 160529, it is not clear how this spectral type relates to the star’s effective temperature. This spectral type of A2Ia is in disagreement with the strength of Hei$`\lambda `$5876, which would indicate a spectral type of B9Ia – A2Ia. The given range reflects the strong variations of this line (cf. Fig. 2).

We find that the spectrum before 1999 was significantly different from the maximum phase spectra of Wolf et al. (1974). The spectral type was earlier (hotter) than in maximum phase. In addition, the emission lines are stronger in minimum phase. The temperature change is inferred mainly from the weakening of Tiii and Crii lines in photometric minimum. However, since these lines form also in the stellar wind, the spectral type changes are unreliable indicators for the effective temperature. The Hei$`\lambda `$5876 line may be a better temperature indicator than the metal lines, although the variable profile of Hei$`\lambda `$5876 indicates that this line is also not of purely photospheric origin.

Nevertheless, part of the spectral changes can certainly be attributed to changes in temperature. The strength of most metal lines strongly increased in 2002. In particular the Cai$`\lambda `$4227 line, which had equivalent widths around 50 mÅ in earlier years, increased to about 270 mÅ in 2002 (cf. Fig.3. This line is a good temperature indicator for A stars. The strength of this line is slightly smaller than in the LBV HD 269662 (=R 110) in 1989, which at that time was classified as about F0Ia (Stahl et al. 1990). The evidence is strong that the star in 2002 was considerably cooler than in the years before. The spectral type was earlier than F0, but much later than A2 and probably late A, confirming the spectral type variations reported earlier.

### 3.2 Variations of line strengths and radial velocity

Since HD 160529 developed a new S Dor phase during our observing campaign, strong spectroscopic changes are to be expected.

The general appearance of the spectrum of HD 160529 changes very little within one observing season. Much stronger variations were observed during the whole campaign from 1992 to 2002. The lines of Crii and Tiii are very temperature-sensitive in the temperature range of HD 160529 (i.e. about 8 000 to 12 000 K) and temperature changes should be reflected in the strength of these lines. Indeed, these lines exhibit particularly strong variations during our observing campaign. However, the line profiles are also highly variable, and is is therefore not clear that the changes in line strength are due to temperature changes only. Especially the short-term changes could be significantly affected by stellar wind effects (cf. e.g. Fig 7), while the long-term changes, where these effects are averaged out, are more likely to be dominated by temperature variations. Fig. 4 shows an example of the spectral changes in Crii lines. The Crii are also sharper in 2002, which is also visible when comparing with Fig. 7 of Sterken et al. (1991): they show a spectrum obtained in 1990, where the Crii lines have a strong blue-shifted component. Compare with Fig. 4 of Wolf et al. (1974), where the Crii lines are always strong. Also by comparing our equivalent width measurements with Table 6 of Wolf et al. (1974), we find much lower values for Crii and Tiii in most years. Only in 1999 and 2002, our values approach the values of Wolf et al. (1974). The weakening of Tiii during minimum has also been reported by Sterken et al. (1991).

Fig. 4 also shows strong variations of the H$`\beta `$ line, both in strength and line profile, between 1992 and 2002. The line strength variations are probably caused by a density decrease in the stellar wind, caused by a radius increase and/or a decrease of the stellar radius. The line profile shows that, while the terminal velocity is unchanged, the minimum of the absorption moved to smaller expansion velocities. This could be due to decreased optical depth in the H$`\beta `$ line or a decrease of the stellar wind velocity in the inner part of the wind.

The Hei$`\lambda `$5876 line is relatively symmetric and thus appears well-suited to study the photospheric variations of HD 160529. The radial velocity of this line is highly variable, but also the FWHM and equivalent width of the line change. We fitted Gaussian profiles to all spectra and so obtained not only the radial velocity but also the residual intensity and line width as a function of time. The results are shown in Figs. 2 and 5. It is obvious from this figure that the line is not only variable in radial velocity but also in line width and intensity. Radial velocity and equivalent width are not correlated. The equivalent width and radial velocity are highly variable within each observing season, i.e. on time scales of a few months. In addition, the average equivalent width seems to decrease from 1995 to 2002. This probably indicates a lower temperature in these years which is supported by e.g. Tiii line measurements. The H$`\alpha `$ strength also changed strongly in 1999–2002. The equivalent width was varying between about 14 and 17 Å from 1992 to 1997. In 1999 it was down to 10 Å and further decreased to 7.5 Å in 2001 and to 6.5 Å in 2002.

The relation of spectroscopic behaviour and light curve on short time scales can only be studied from the spectra obtained in 1992, when we have sufficient overlap in the spectroscopic and photometric data set (cf. Fig. 5). The sudden decrease of radial velocity seems to coincide with a peak of the visual brightness. The equivalent width of Hei$`\lambda `$5876 and also the FWHM have a peak at the same time. Note that this correlation of line strength with visual brightness on short time scales contrasts with the pronounced anti-correlation on long time scales (cf. Fig. 2). These short-term variations of the Hei$`\lambda `$5876 line could be due to pulsational variations or to stellar-wind effects. If they are due to pulsation, the variations of the equivalent width imply an increase of the temperature with brightness. If we interpret the variations as changes of the stellar wind, the observations would suggest a strengthening of the wind with brightness, with a phase-delay of the radial velocity.

Other lines which mainly originate in the photosphere are e.g. Mgii$`\lambda `$4481 and the Siii$`\lambda \lambda `$6347, 6371 lines.

If variations in the strength of the Tiii, Crii and Hei lines are due to photospheric temperature change of the star, we would expect the changes in Hei and the metal lines to be anti-correlated. Indeed, there is a weak evidence for this behaviour in some years, but the scatter is very large within individual years. In the long-term trend, the anti-correlation is clearly visible.

This is probably due to the fact that the metal lines form mostly in the stellar wind (see Figs. 6 and 7). Only at the end of the observing run, the lines appear close to the systemic velocity and with a line width that is compatible with $`v\mathrm{sin}i`$. In other years, the radial velocity, line width and line profile variations indicate a stellar wind origin for these lines. Line splitting is quite common and has already been observed by Wolf et al. (1974). In particular in the first years of the observations, the radial velocities are more negative than observed by Wolf et al. (1974) in maximum visual brightness.

### 3.3 Time series analysis

Lamers et al. (1998) have previously analyzed the photometric data of HD 160529. For two seasons, they derived a period of 45 and 55 days, respectively. They explained the variations by non-radial $`g`$-mode pulsations of low order $`l`$ and the period change by mode-switching.

We examined the photometric data and radial velocity measurements for periodicities with the ESO-Midas TSA package using the method of Scargle (1982). We limited the measurements to seasons with a significant number of measurements. For the photometry, these were the years 1991–1993 and for the spectroscopy the years 1992–1995. We restricted the search to periods below 120 days, since this is the approximate length of the observing runs. For the $`y`$-photometry we obtain the strongest peaks in the power spectrum at about 108, 83 and 68 days and for the radial velocity of Hei at 93 and 67 days. The equivalent width of this line shows two peaks at similar periods at 92 and 67 days and in addition two peaks at 126 and 82 days. We checked that the results do not depend critically on the choice of the data set. If data sets from individual years are removed or added, the detected periods change, but we always find peaks at similar periods, i.e. within a few days. This could indicate that the variations are not strictly periodic or not coherent over longer time scales. However, the small changes of the position of the peaks of the power spectrum are more likely due to the irregular sampling of our data.

It should also be mentioned that the radial velocity changes are definitely not sinusoidal. Although the time series analysis gives time scales of around 100 days, scales for the changes are of the order of 10 days. The time series analysis measures the time scale on which such sudden changes of the radial velocity repeat. Examples of the observed variations are shown in Figs. 2 and 5. Similar events have been observed in other years as well.

The similarity between the time scales of the photometric and spectroscopic variations suggests a physical connection. In 1992, when we obtained photometry and spectroscopy almost simultaneously, a peak in the light curve was observed to coincide with the a rapid change in radial velocity (Fig. 5). The equivalent width of Hei$`\lambda `$5876 also has a peak at the same time. In fact, the equivalent width seems to follow quite closely the visual light curve. Unfortunately, our observations do only allow to observe these relation in the 1992 observing seasons. Therefore, we can only conclude that the rapid changes in radial velocity may be related to the photometric variability.

### 3.4 Line profile variations

Many lines in the spectrum of HD 160529 show peculiar line profile variations. While this was already known from the work of Wolf et al. (1974), the time sampling of the spectra of these authors was not adequate to follow the changes in line profiles in detail.

We searched in particular for propagating features in the line profiles, in order to study the stellar wind of HD 160529. Most of the lines are dominated by absorption and the observed variability seems mainly due to additional absorption. In order to increase the contrast, we therefore examined the residuals of the profiles with respect to a “maximum flux” profile, which was constructed from the time series by computing the highest profile for each wavelength bin in each observing season. This procedure is very sensitive to noise. Therefore, a running median was applied along the $`t`$-axis before computing the maximum flux. The subtraction of a “maximum flux” profile is more appropriate than subtracting a mean profile if the variations are due to additional absorption on top of a stable profile. The subtraction of a mean profile always introduces emission-like features in the residuals.

A number of strong absorption lines do not show any emission component or a strong wind absorption and, therefore, superficially look like photospheric lines. An example of such metal lines are the Siii lines at 6347, 6371 Å. Profile variations as observed in 1993 are shown in Fig. 8.

Most of the stronger metallic lines, in particular Feii lines, show pronounced P Cyg profiles. These lines show peculiar variations of the line profiles. An example of the variations is shown in Fig. 10. Also in these profiles, most of the variability can be ascribed to additional absorption of variable strength, but with little change in radial velocity. The variability is strongest at very low velocities around 0 km s<sup>-1</sup>. This component appears directly on top of the emission component. Therefore the P Cyg profile changes to a splitted absorption profile (see Fig. 10, top) and later back to a P Cyg profile. Additional absorption at blue-shifted velocities of about $`50`$ km s<sup>-1</sup> is also observed, but the features at low and high expansion velocities do not seem to be connected. This second, blue-shifted component seems to move to higher expansion velocities, but the changes are quite small. In some cases, the absorption features even seem to move towards longer wavelengths. Features which are clearly moving to higher expansion velocities, i.e. towards shorter wavelengths, are not often observed. Therefore, while we have some evidence for structures which originate in an accelerating stellar wind, the strongest variations occur close to the systemic velocity. A cool photospheric absorption layer may be responsible for these variations.

The low-velocity feature is particularly strong in 1997, where it is clearly seen as a separate absorption feature. As an example we show in Fig. 9 the profile of the Feii$`\lambda `$4731 line in the various observing seasons. In this figure it can be seen that the line was clearly blue-shifted from 1991 to 1994 and later moved to longer wavelengths. The low-velocity feature could be related to the slightly red-shifted absorption features seen sometimes in LBVs (Wolf & Stahl 1990; Stahl 2001).

A comparison of the line profile variations of lines which have a photospheric appearance (e.g. Siii$`\lambda `$6347 in Fig. 8) and lines which clearly show a P Cyg profile (e.g. Feii$`\lambda `$6248 in Fig. 10) shows that the variations are very similar, although the mean profiles of Siii$`\lambda `$6347 and Feii$`\lambda `$6248 are very different. It should be mentioned, however, that also the lines with pure absorption profiles show clear variations from year to year in radial velocity and line width.

### 3.5 The resonance lines Nai, Caii and Ki

The NaiD lines show, apart from a strong interstellar component, clearly also a stellar contribution: a weak emission component and multiple variable absorption components. The variable absorption components are evidence for accelerating expanding shells. The maximum expansion velocity is typically around $`150`$ km s<sup>-1</sup> and reaches up to $`200`$ km s<sup>-1</sup>. Discrete components are frequently observed in the velocity range from $`60`$ to $`150`$ km s<sup>-1</sup>.

The changes in radial velocity of the components during one observing season are relatively small (up to about about 15 km s<sup>-1</sup>, but typically below 10 km s<sup>-1</sup>) The changes from year to year are so large (cf. Fig. 11), that it is not clear how the components are connected from year to year. Therefore, the development of the velocities on longer time-scales cannot reliably be determined from our observations. Similar components may be present at lower velocities as well, but cannot be seen because of the strong interstellar contribution, which saturates the absorption at velocities smaller than about $`50`$ km s<sup>-1</sup>. Acceleration of the components in general is very small. In 1992, no significant acceleration was seen. In 1993 and 1995, an accelerating component was seen, most clearly in 1993 (cf. Fig. 12). Even in this year, we observed an acceleration of about 0.12 km s<sup>-1</sup> d<sup>-1</sup> at velocities from $`75`$ to $`90`$ km s<sup>-1</sup>. Another component at higher velocities ($`115`$ km s<sup>-1</sup>) was stable in radial velocity. Both components first weakened and later strengthened again. In 1995, three components were observed ($`125`$, $`77`$ and $`103`$ to $`114`$) and one was accelerating. In 1994, two components were seen. Both components were very stable in radial velocity. The component at large expansion velocity ($`120`$ km s<sup>-1</sup>) was weakening, while the component at smaller velocity ($`75`$ km s<sup>-1</sup>) was getting stronger with time.

For CaiiH, K much less data is available. In particular, these lines where not covered in the years 1992–1994, where we have obtained the most extended time series. The yearly mean spectra of the CaiiK line are shown in Fig. 13. In cases where both the Caii and the Nai lines are available, we find different velocities in both groups of lines. The blue-shifted components are much stronger in Caii and therefore in most cases not as clearly separated. For the un-shifted component, which is mainly of interstellar origin, the Nai lines are stronger than Caii. In Caii, the interstellar contribution is split in two components at $`26`$ and $`4`$ km s<sup>-1</sup>. These components are only resolved in the Ucles and Feros spectra. This splitting is also marginally visible in the Nai lines.

Overall, the picture is not clear. There is no discernible long-term trend in the velocity or acceleration of the components. From an expanding stellar wind ($`R=200`$ R, $`v_{\mathrm{}}=200`$ km s<sup>-1</sup>, $`\beta =4`$), the expected acceleration is about 1 km s<sup>-1</sup> d<sup>-1</sup>, i.e. significantly larger than what is observed. This implies that the line-splitting is probably not due to radially expanding shells moving with the stellar wind. This behaviour is similar to what is observed for discrete absorption components (DACs) in OB supergiants (cf. e.g. Fullerton et al. 1997). The origin of these features and their slow acceleration is still not understood.

The Ki line at 7699 Å is another interstellar absorption line. This line has only been observed from 1995 to 2002. Due to the low abundance of K, the line is much weaker than Nai or Caii. No blue-shifted components are seen. The line has two absorption components: a strong main component at $`2`$ and a weaker component at $`29`$ km s<sup>-1</sup>. The components are clearly resolved only in the Feros spectra. From the two Feros spectra it appears that the fainter, blue component is slightly variable and therefore at least partly of stellar origin.

## 4 Stellar parameters

### 4.1 Distance

The distance of HD 160529 is important for the derivation of stellar parameters. Unfortunately, the distance is very uncertain: Sterken et al. (1991) estimate a distance of 2.5 kpc based on assigning to HD 160529 the same absolute visual magnitude as the LBV HD 269662 (=R 110) in the LMC ($`M_\mathrm{V}=8.9`$). R 110 is spectroscopically very similar to HD 160529 and also has a comparable photometric amplitude. In addition, an interstellar extinction of $`A_\mathrm{V}`$ = 3.4 mag ($`E_{\mathrm{B}\mathrm{V}}=1.1`$) is assumed. The empirical relation between luminosity and amplitude of LBVs (Wolf 1989) would give a lower bolometric luminosity of about $`M_{\mathrm{Bol}}=8.3`$, if we assume an amplitude of 0.5 mag. This would also yield a slightly smaller distance estimate of about 1.9 kpc. If we assume that HD 160529 has an absolute visual brightness of $`M_{\mathrm{Bol}}=9.5`$, similar to HD 33579, the brightest A hypergiant of the Large Magellanic Cloud, the distance would be 3.5 kpc.

A kinematic distance estimate is difficult for several reasons: the systemic velocity is difficult to determine because of variable radial velocities and line profile variations. In addition, the line-of-sight to HD 160529 is close to the galactic center and kinematic distances therefore have large errors. In the following, we will assume a distance of 2.5 kpc but note that this distance has a rather large uncertainty of about 30%.

### 4.2 Temperature changes

The temperature of HD 160529 is another uncertain parameter. Line strengths, derived e.g. from Siii, Hei and Feii or Tiii give contradictory results, when they are compared to normal supergiants or analyzed with static model atmospheres. This is certainly due to temperature stratification in the expanding atmosphere and stellar wind contributions in most absorption lines.

The only metal which appears in different ionization stages is Mg. The Mgi/Mgii line ratio can be used to estimate the spectral type. The Mgi lines $`\lambda \lambda `$3832, 3838, 5173, 5183 and the Mgii lines $`\lambda \lambda `$4481, 7877 were used. The results are summarized in Table 2. For comparison, the A0Ia star HD 92207, the A2Ia star HD 100262 and the F3Ia star HD 74180 were also measured. It should be noted that the line profiles of the measured lines also change with time. Therefore, the measured line ratios do not necessarily reflect changes in the effective temperature only. It is likely, however, that most of the variations of the line ratios reflect temperature changes.

The near-IR also contains a large number of Ni lines, mainly in the wavelength range around 7 440 and 8 700 Å. These lines are very strong, significantly stronger than in any of the late B/early-A supergiants observed by Przybilla & Butler (2001). For the spectra taken between 1995 and 1999, where the equivalent widths of these lines are relatively stable, we obtain for the lines $`\lambda \lambda `$7424, 7442, 7468 mean equivalent widths of 325, 417 and 533 mÅ, respectively. In 2002 the lines are still stronger, with equivalent widths of 423, 529 and 614 mÅ, respectively. For HD 100262 (A2Ia), we measure 87, 141 and 196 mÅ, respectively. The latter numbers are in very good agreement with the A2Ia star HD 111613 investigated by Przybilla & Butler (2001). The much stronger lines for HD 160529 may indicate an enhanced abundance of N, but could also be due to non-LTE effects.

Lines from neutral oxygen are strong in HD 160529. We measured an equivalent width of between 2.15 and 2.40 Å for the near-IR triplet $`\lambda \lambda `$7771-5. For A stars, the strength of this feature increases with decreasing $`\mathrm{log}g`$ (Przybilla et al. 2000).

In principle, the temperature can also be derived from the continuum energy distribution. The optical energy distribution is not very temperature sensitive, so the UV energy distribution is needed. In the case of HD 160529 this is difficult since the interstellar reddening is high and quite uncertain. However, low dispersion IUE spectra obtained in 1992 and 1979, i.e. during the visual minimum and maximum phase, respectively, show that the flux at about 2 500 Å is lower by about a factor 2 in visual maximum. This convincingly demonstrates that the temperature was significantly lower at visual maximum.

We tried to use synthetic spectra computed with the hydrostatic model atmosphere code Tlusty (Hubeny & Lanz 1995; Lanz & Hubeny 2001) to determine the temperature of HD 160529 and its changes. However, it turned out that the line spectrum of HD 160529 cannot be explained by hydrostatic models. Models with solar abundance which fit the Hei$`\lambda `$5876 line in minimum require about 11 000 K, while fits of metal lines such as, e.g., Mgii$`\lambda `$4481 would require about 9 000 K. Increasing the He abundance reduces the discrepancy only slightly. We interpret this result as a strong indication of stratification effects in the expanding atmosphere. Non-LTE models which include the atmospheric expansion and line blanketing are required. While such models are now available for WR stars (Hillier & Miller 1998), they are so far not yet available for cooler stars like HD 160529.

The changes of temperature and radius over the variability cycle can be constrained by assuming that the bolometric magnitude of HD 160529 is constant. The change of 0.5 magnitudes in the visual brightness then reflects a change in the bolometric correction. Assuming a minimum temperature of 8 000 K and a corresponding radius of 330 R, we derive a maximum temperature of 12 000 K and a corresponding radius of 150 R. The change of the temperature is necessary for a corresponding change in bolometric correction. This conclusion follows assuming a black-body distribution or if we use empirical bolometric corrections (Schmidt-Kaler 1982). The large change in radius is required by the constancy of the bolometric luminosity. The temperatures correspond to spectral type of about B8 in visual minimum and A9 in visual maximum (Schmidt-Kaler 1982), in agreement with the spectroscopic evidence.

### 4.3 Systemic velocity and rotational velocity

The systemic velocity of HD 160529 is difficult to determine since all absorption lines are influenced to some extent by the stellar wind, their radial velocities are variable and different for different groups of lines. We therefore selected emission lines which appear little disturbed by absorption to estimate the systemic velocity. The Feii lines of multiplets 40 and 46 appear well suited. From these lines we derive a radial velocity of $`29\pm 3`$ km s<sup>-1</sup>. We also used the forbidden line Nii$`\lambda `$6583 to estimate the systemic velocity. This line is flat-topped and probably forms at large distance from the star. The line center, as estimated from the bi-sector of the left and right edge, is at $`8`$ km s<sup>-1</sup>, if we use a rest wavelength of 6583.454 Å (Spyromilio 1995). Unfortunately, the line is faint and broad, so that the uncertainty of this measurement is large.

The average velocity measured from Hei$`\lambda `$5876 is about $`23`$ km s<sup>-1</sup> (see below), which is relatively close to this value. Other absorption lines such as the Crii multiplet 30, which appear as symmetric and narrow absorption lines in maximum, have velocities of around $`20`$ km s<sup>-1</sup>. The red Oi and Ni lines give similar results: they vary between $`27`$ and $`19`$ km s<sup>-1</sup> from 1995 to 2002. Because the radial velocity may be variable also for symmetric lines, it is not clear whether this velocity represents the true systemic velocity of HD 160529.

Because of the variable line profiles, the rotational velocity $`v\mathrm{sin}i`$ is also difficult to determine for HD 160529. Feros spectra taken in 1999 and 2002 show particularly narrow absorption lines. We used these spectra to estimate an upper limit for $`v\mathrm{sin}i`$. For this purpose, we fitted rotational profiles to some metal lines, taking into account the instrumental resolution of the spectrograph, but no other broadening mechanisms of the lines. In this way, we derive a $`v\mathrm{sin}i`$ value of about 45 km s<sup>-1</sup>. This has to be considered as an upper limit, since stellar wind effects and other broadening mechanisms are most likely important.

## 5 Modeling of H$`\alpha `$

Mass-loss rates derived from radio observations are only little model-dependent. HD 160529 is one of the few stellar objects which are bright enough for current radio telescopes. Leitherer et al. (1995) observed the star at 8.64 and 4.80 GHz and found that the measured flux is consistent with thermal emission from an optically thick expanding wind at constant velocity. The authors adopt the spectroscopic distance of $`d=2.5`$ kpc from Sterken et al. (1991), an electron temperature of $`T_\mathrm{e}`$ = 5 000 K for the radio emitting region and a terminal velocity of $`v_{\mathrm{}}`$ = 180 km s<sup>-1</sup> to derive a mass-loss rate of $`\dot{M}`$ = $`10^{4.87\pm 0.24}`$$`\mathrm{M}_{}`$yr<sup>-1</sup>.

From our observations (cf. Fig. 14), we find that the emission-line flux of HD 160529 was slowly decreasing during our observing period, i.e. with increasing visual brightness. The variations during each season are much smaller than this long-term trend.

In order to quantitatively investigate the mass-loss rate of HD 160529 versus time, we modeled the mean H$`\alpha `$ line of the years 1992 and 2002, which show the strongest and weakest H$`\alpha `$ emission, respectively. We used the multi-level line-transfer code described by Bastian (1982) and Stahl et al. (1983). A terminal velocity $`v_{\mathrm{}}`$ of 200 km s<sup>-1</sup> was used for both minimum and maximum. Note that stellar-wind theory predicts a decrease of $`v_{\mathrm{}}`$ with decreasing surface gravity, i.e. with increasing visual brightness of LBVs. This is not easily observed in HD 160529, since there are no lines in the optical with a well-defined blue absorption edge. The absorption edge in H$`\alpha `$ is variable, but this can also be explained by changes in the optical depth of the wind. A constant wind velocity cannot be excluded from our observations. We used a $`\beta `$-type velocity law with $`\beta `$=4. A velocity law with such a slow acceleration was found by Stahl et al. (2001) to fit the line profiles of AG Car in outburst. The stellar radius and effective temperature were estimated from the photometry and the line spectrum of the given year. Then we only varied $`\dot{M}`$ and $`v_{\mathrm{sys}}`$ until a satisfactory fit was achieved. Solar abundances were assumed. Possible clumping was not taken into account. The results are summarized in Tab. 3 and the fits are shown in Fig. 15.

We find that the mass-loss rate of HD 160529 is almost independent of temperature, i.e. visual brightness. Note that the derived mass-loss rate increases slightly from minimum to maximum, although the equivalent width decreases in this period. The reduced equivalent width is mainly due to the radius increase and – to a smaller extent – to the decrease in temperature. The radius increase causes (at constant mass-loss rate) a density decrease which will decrease the equivalent width of the Balmer lines. Therefore, if the radius increase is smaller than derived from the assumption of constant bolometric luminosity, the derived mass-loss rate could even slightly decrease from minimum to maximum. The derived mass-loss rate is in good agreement with the results of Leitherer et al. (1995), which are based on radio observations. We estimate that our mass-loss rates are accurate within about a factor of two. We therefore conclude that within the errors the mass-loss rate did not change significantly between minimum and maximum.

A more advanced modeling of the stellar wind of HD 160529 with a code as e.g. described by Hillier & Miller (1998) would certainly be worthwhile. However, results of such models for cool stars such as HD 160529 are not yet published.

## 6 Discussion and conclusions

For the interpretation of the short-term variability of HD 160529, the relevant time scales are the rotational time scale $`P_{\mathrm{rot}}`$, the stellar-wind time scale $`P_{\mathrm{wind}}`$ and the pulsation time scale $`P_{\mathrm{puls}}`$.

The time scales are defined by the following relations:

$`2\pi R_{}/v_{\mathrm{break}}<P_{\mathrm{rot}}<2\pi R_{}/v_{\mathrm{rot}}`$ (1)
$`P_{\mathrm{wind}}=R_{}/v_{\mathrm{}}`$ (2)
$`\mathrm{log}P_{\mathrm{puls}}=0.275M_{\mathrm{Bol}}3.918\mathrm{log}T_{\mathrm{eff}}+14.543`$ (3)
In Equation 1, the upper limit to $`P_{\mathrm{rot}}`$ is estimated from $`v\mathrm{sin}i`$ assuming $`\mathrm{sin}i=1`$. This value can also be slightly larger because the derived $`v\mathrm{sin}i`$ value is an upper limit. Note that it is expected that the rotational velocity is a function of phase. Because of the complex line profiles it could not be determined in light minimum and is therefore set constant here. Because of the acceleration of the wind to its terminal velocity in the extended stellar wind, the typical expansion time scale can be significantly longer than given by Equation 2. The pulsation period for the fundamental radial mode given by Equation 3 is the empirical fitting formula given by Lovy et al. (1984).

We used the following parameters from Sterken et al. (1991): $`M_{\mathrm{bol}}=8.9`$, $`M_{}=13\mathrm{M}_{}`$. In addition, we used $`v_{\mathrm{}}`$ = 200 km s<sup>-1</sup> and $`v_{\mathrm{rot}}`$ = 45 km s<sup>-1</sup> for both maximum and minimum phase. For the stellar radius and temperature we adopted 150 R/12 000 K and 330 R/8 000 K for minimum and maximum phase, respectively. With these parameters, we obtain the time scale estimates given in Table 4.

Within the uncertainties, the time scales $`P_{\mathrm{wind}}`$ and $`P_{\mathrm{puls}}`$ are compatible with the typical variation time scales of 50–100 days, while the rotational time scale is significantly longer. However, since rotational modulation can also lead to variations with an integer fraction of the rotational period, the time scale of the variations cannot be used to distinguish the possible mechanisms of the short-term variations in HD 160529. The simultaneous variation of the optical brightness and the Hei$`\lambda `$5876 line, however, is most easily explained by pulsations. Pulsations have also been suggested as a cause for the photometric micro-variations of LBVs (cf. e.g. Lamers et al. 1998).

The correlation of the equivalent width variations with visual brightness on short time scales is in marked contrast with the behaviour on long time scales. This strongly suggests that the physical mechanism for the variations on short time scales is different from the LBV-type variations.

It is rather surprising, that so far no strong dependence of the observed time scale on the visual brightness – and therefore the radius – has been detected. All plausible variability mechanisms would predict such a dependence. It should be mentioned, however, that Sterken et al. (1991) derived a period of 57 days in minimum while the derived value in maximum was 101 days – about a factor of two longer. It is not clear, however, if the difference is due to a real period change.

Most LBVs seem to increase their mass-loss rate in maximum. However, it is not clear that this is a general property of LBVs (Leitherer 1997). Vink & de Koter (2002) have shown that the mass-loss rate of LBVs versus temperature carries important diagnostic information. The mass-loss rate of AG Car versus phase has been studied by Stahl et al. (2001). Vink & de Koter (2002) compared these results with their predictions for line-driven stellar winds and found a good agreement. They also predict that the mass-loss rate of LBVs is a complicated function of temperature, surface gravity and abundances. In particular, the mass-loss rate can even decrease with decreasing temperature at temperatures below about 15 000 K, in particular at low masses. This predicted behaviour nicely fits the observed behaviour of HD 160529, which is an LBV with very low temperature and mass.

Although the mass-loss rates of LBVs can be predicted with some precision, the details of the mass-loss process are still not understood. In particular the complicated line-profile variations and the line splitting cannot be explained by radially expanding features. Rotational modulation may be the cause for these variations.

###### Acknowledgements.

We thank the European Southern Observatory (ESO) for the generous allocation of observing time and the staff at La Silla for the support in the installation of our Heros instrument at the telescope. We acknowledge the use of the SIMBAD database (CNRS data center, Strasbourg). This research has made use of NASA’s Astrophysics Data System (ADS) Abstract Service. This work was supported by the Deutsche Forschungsgemeinschaft (Wo 296, 9-1). TG would like to thank C. Leitherer and acknowledges support from STScI-DDRF grants. This project was supported by the Belgian Federal Office for Scientific, Technical and Cultural Affairs (IUAP P5/36), the Belgian Fund for Scientific Research (FWO), and by the Danish Natural Science Research Council, partly through the center for Ground-Based Observational Astronomy. We express our gratitude to T. Arentoft, A. Bruch, H.W. Duerbeck, H. Melief, M. Nolte and A. Visser for contributing some photometric measurements in the framework of LTPV. We thank S. Bagnuolo for taking some of the Feros spectra.


text:
# Thermal emissivity for finite three-dimensional photonic band gap crystals

## I INTRODUCTION

It has been shown that there exist periodic dielectric structures such that for a range of frequencies, no electromagnetic waves can propagate in any direction for all polarizations old (1993, 1994); Joannopoulos et al. (1995); Yablonovitch and Leung (1991); Sözüer and Dowling (1994). The properties of such an infinite crystal have been worked out in detail John and Wang (1990, 1991). Such work is typically limited to non-absorptive materials, since otherwise all waves become evanescent.

Much theoretical work has been done on the properties of finite one dimensional photonic band gap (PBG) crystals Bendickson and Dowling (1996), including recent calculations of the thermal emissivity of such one-dimensional structures Cornelius and Dowling (1999); Stimpson and Dowling (2001). The strong angular dependence of the gap effect with a one-dimensional structure has motivated successful experimental work with three-dimensional structures Lin et al. (2000). The three-dimensional model we will use here is an inverted opal structure, and with the use of a computer program published by Stefanou et al. Stefanou et al. (1998); and V. Yannapapas and Modinos (2000), the transmittance, reflectance, and absorbance of a finite number of layers of the crystal is calculated, as well as the complex band structure of the corresponding infinite crystal. As we shall show, our structure exhibits a broad band gap in the emissivity that has virtually no angular dependence, and as such, also suppresses the thermal emittance of the crystal in an omnidirectional fashion. We also predict emissivity enhancement at the photonic band edge.

For a given crystal structure, we can calculate $`𝒯`$ and $``$, the transmission and reflection coefficients of the crystal as a function of the frequency $`\omega `$ and the wave vector $`𝐤`$ of the incident wave. If we consider a crystal with real indices of refraction, then $`+𝒯=1`$ from conservation of energy. However, if there is an imaginary part to the any of the indices of refraction, there will be a non-zero absorbtion coefficient $`𝒜=1𝒯`$. Calculation of this quantity is useful because of a Kirchoff’s second law which states that a material’s thermal emmittance $``$ is proportional to its absorbance (and that they are equal for a blackbody). One can obtain the absorbance by simply multiplying Planck’s blackbody power spectrum by the absorbtion coefficient. So by modelling the absorbtion of several layers of a three-dimensional crystal coated on the surface of a thick, highly absorptive substrate, we are also modelling the thermal emissivity $`E`$ of the entire structure. The three-dimensional superlayer acts as a passive filter, which alters the emissivity of the underlying graybody substrate.

## II THEORY

### II.1 One-dimensional case

We employ software developed by Stefanou et al. Stefanou et al. (1998); and V. Yannapapas and Modinos (2000) that uses a generalization of the matrix-transfer techniques developed for one-dimensional crystals, such as thin-film dielectric stacks. In the one-dimensional case, the $`𝒯`$ and $``$ are $`\left|t\right|^2`$ and $`\left|r\right|^2`$ respectively, where t and r are the complex transmittance and reflectance coefficients, respectively. These coefficients are linearly related to one another by a $`2\times 2`$ matrix $`\widehat{𝐌}`$, where

$$\left[\begin{array}{c}1\\ r\end{array}\right]=\widehat{𝐌}\left[\begin{array}{c}t\\ 0\end{array}\right]$$
(1)
and

$$\widehat{𝐌}=\left[\begin{array}{cc}M_{11}& M_{12}\\ M_{21}& M_{22}\end{array}\right].$$
(2)
The column vectors on each side of Eq. (1) represent the fields outside each side of the crystal. Then $`t=1/M_{11}`$ and $`r=M_{21}/M_{11}`$. This matrix $`\widehat{𝐌}`$ can be constructed for a one-dimensional structure simply by multiplying the transfer matrices for the individual layers. For a homogenous plate, these can be readily computed by taking into account the the reflection and transmission at the interface at the edges of the plate with the Fresnel equations, along with the accumulated phase (and possible change in amplitude if the index of refraction has an imaginary component) from propagating though the material.

Because the phase accumulated from the wave propagating though a given layer depends on the portion of the wave-vector that is perpendicular to the interfaces of the plates, the frequency $`\omega `$ at which the gap is centered is highly dependent on the angle of incidence of the initial radiation. As a consequence, while the on-axis thermal emissivity at a given $`\omega `$ might be suppressed, off-axis emissivity might have a band-edge enhancement, as shown in earlier work Cornelius and Dowling (1999); Stimpson and Dowling (2001).

### II.2 Three-dimensional model

The model used by the Stefanou and co-workers Stefanou et al. (1998); and V. Yannapapas and Modinos (2000) also allows for a finite number of parallel planes of spheres of differing index of refraction from the host material in which they are embedded. For a given frequency $`\omega `$, the incident plane wave is expanded into spherical waves about an origin centered at one sphere, and then the scattered wave from a single sphere in one plane is calculated by satisfying the boundary conditions for Maxwell’s equations at the surface of a single sphere. The scattered wave for the entire infinite plane of spheres is just the sum of the scattered waves for a single sphere expanded about a translated origin, with a phase factor corresponding to the portion of the incident wave vector that is parallel to the plane of spheres. By summing this first-order scattering of all the spheres, *except* the sphere at the origin, and then expanding the result into spherical waves about the origin, we can get a better approximation of the true wave that is incident on the sphere at the origin.

Here, the restriction that each plane of spheres has the same two-dimensional periodicity allows for the difficult part of this expansion to be calculated only once per frequency, which is what makes this method much more numerically feasible than other approaches.

By solving the boundary conditions again, with the spherical wave expansion of the incident plane wave plus this first order scattered wave, we solve for a second-order scattering off the sphere at the origin. By symmetry, all the spheres of the infinite plane have this same scattered wave. Hence, the total scattered wave for the plane is again the sum of the scattered waves for a single sphere, expanded about a translated origin, with a phase factor. The end result is then expanded into plane waves and added to the incident wave to to solve for the total transmitted and reflected field. By repeating this process for each plane of spheres, the total field for the entire crystal is obtained.

## III RESULTS

The three-dimensional crystal chosen for our calculations is a diamond fcc inverted opal (Fig. 1) with the complex relative dielectric constant $`ϵ_1=12+i/10`$ for the space around the spheres and $`ϵ_2=1`$ inside the spheres, which are of radius $`0.30618621a`$. Four layers of the unit cell are stacked on a highly absorptive substrate ($`ϵ=12+7i`$) which is eight orders of magnitude thicker than the height of the total crystal, so that it can be considered essentially infintely thick. The emissivity dependence on normalized frequency and polar angle for both polarizations is given in Fig. 2.

Several important features are noted. Firstly, there is a large band-gap in the emissivity that is largely angle independent. Secondly, the emissivity of the graybody substrate goes exactly to unity at the band-edge resonance frequencies. This means that the emission is enhanced up to the pure blackbody limit rate at these frequencies. This observation is consistent with our previous one-dimensional calculations Cornelius and Dowling (1999); Stimpson and Dowling (2001), and recent experimental observations Lin et al. (2000). As noted in our previous work, around the band-edge resonances the three-dimensional photonic crystal superstrate acts as an antireflective coating, allowing all incident radiations at the frequency to be absorbed—which is the definition of a blackbody.

You will notice that a complete gap does not exist for this structure, as is the case for most inverted opals, as you can see emission is not suppressed for a particular angle for the center gap frequency of $`2.27c/a`$. While other structures might show more promise of having complete stop bands, opals are unique in that large crystals can be grown by self-assembly of individual spherical particles, which can then be doped out to provide the inverted opal structure we have presented.

A complementary set of plots for a one-dimensional crystal is given in Fig. 3. For the one-dimensional crystal, indices of refraction and thicknesses chosen such that the absorption gap is centered at the same frequency and of the same width as the three-dimensional crystal. This is achieved by setting $`ϵ_1=2.6`$, $`d_1=0.6a`$, $`ϵ_2=1.44`$, and $`d_2=0.81a`$ where $`d_1`$ and $`d_2`$ are the thicknesses of the the two layers and $`ϵ_1`$ and $`ϵ_2`$ are the corresponding relative dielectric constants. Here, $`a`$ is once again the length of the edge of the unit cell of the three-dimensional crystal, so all frequencies are comparable.

While the frequencies corresponding to suppressed emission in the one-dimensional crystal are highly dependent on angle, the emissive gap effect in the three-dimensional crystal is nearly independent of angle. Though, as mentioned before, a band structure cannot be calculated for a finite, absorptive structure, one can calculate the band structure for the infinitely extended crystal with the absorption removed by setting $`\text{Im}(ϵ_1)=0`$ and removing the backplane. However, the existence of the gap effect in our configuration is highly dependent on the boundary condition of the thick absorptive backplane and the presence of absorbtion in the crystal itself, and no gap effect is present if they are removed. This geometry does exhibit gap behavior at a different frequency, though, if $`ϵ_1`$ is changed appropriately (Fig. 4).

## IV CONCLUSION

We have calculated the thermal emissivity of a three-dimensional PBG crystal using extensions of tested numerical techniques developed for 1D crystals. Our matrix-transfer approach is conceptually much simpler than the traditional QED method of solving the non-trivial boundary conditions involved, and also not nearly as computationally intensive as the Finite Difference Time Domain (FDTD) simulations used by many groups for solving Maxwell’s equations in PBG crystals. Our technique also takes into account the essentially finite nature of the structure, crucial for experimental work where often it is possible to fabricate only a few layers of crystal. Our results agree qualitatively with recent experiments by S. Lin Lin et al. (2000).

## ACKNOWLEDGEMENTS

This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a grant from the National Aeronautics and Space Administration. Additional support was provided by the National Reconnaissance Office, the National Security Agency, the Advanced Research and Development Activity, the Office of Naval Research, and the Defense Advanced Research Projects Agency.


text:
# Scaling exponent of the maximum growth probability in diffusion-limited aggregation

## Abstract

An early (and influential) scaling relation in the multifractal theory of Diffusion Limited Aggregation (DLA) is the Turkevich-Scher conjecture that relates the exponent $`\alpha _{min}`$ that characterizes the “hottest” region of the harmonic measure and the fractal dimension D of the cluster, i.e. $`D=1+\alpha _{min}`$. Due to lack of accurate direct measurements of both $`D`$ and $`\alpha _{min}`$ this conjecture could never be put to serious test. Using the method of iterated conformal maps $`D`$ was recently determined as $`D=1.713\pm 0.003`$. In this Letter we determine $`\alpha _{min}`$ accurately, with the result $`\alpha _{min}=0.665\pm 0.004`$. We thus conclude that the Turkevich-Scher conjecture is incorrect for DLA.

Multifractal measures are normalized distributions lying upon fractal sets. As such, they present rich scaling properties that have attracted considerable attention in the last two decades. In this Letter we address the harmonic measure of Diffusion Limited Aggregates 81WS , which is the probability measure for a random walker coming from infinity to hit the boundary of the fractal cluster. This was one of the earliest multifractal measures to be studied in the physics literature 86HMP , but the elucidation of its properties was made difficult by the extreme variation of the probability to hit the tips of a DLA versus hitting the deep fjords. Thus the understanding of its scaling properties has been a long standing issue. These scaling properties are conveniently studied using the notion of generalized dimensions $`D_q`$, and the associated $`f(\alpha )`$ function 83HP ; 86HJKSP . The simplest definition of the generalized dimensions is in terms of a uniform covering of the boundary of a DLA cluster with boxes of size $`\mathrm{}`$, and measuring the probability for a random walker coming from infinity to hit a piece of boundary which belongs to the $`i`$’th box. Denoting this probability by $`P_i(\mathrm{})`$, one considers 83HP

$$D_q\underset{\mathrm{}0}{lim}\frac{1}{q1}\frac{\mathrm{log}_iP_i^q(\mathrm{})}{\mathrm{log}\mathrm{}},$$
(1)
where the index $`i`$ runs over all the boxes that contain a piece of the boundary. The limit $`D_0lim_{q0^+}D_q`$ is the fractal, or box dimension of the cluster. $`D_1lim_{q1^+}D_q`$ and $`D_2`$ are the well known information and correlation dimensions respectively 56BR ; 82Far ; 83GP . It is well established by now 86HJKSP  that the existence of an interesting spectrum of values $`D_q`$ is related to the probabilities $`P_i(\mathrm{})`$ having a spectrum of “singularities” in the sense that $`P_i(\mathrm{})\mathrm{}^\alpha `$ with $`\alpha `$ taking on values from a range $`\alpha _{\mathrm{min}}\alpha \alpha _{\mathrm{max}}`$. The frequency of observation of a particular value of $`\alpha `$ is determined by the function $`f(\alpha )`$ where (with $`\tau (q)(q1)D_q`$)

$$f(\alpha )=\alpha q(\alpha )\tau \left(q\left(\alpha \right)\right),\frac{\tau (q)}{q}=\alpha (q).$$
(2)
Of particular interest are the values of the minimal and maximal values, $`\alpha _{min}`$ and $`\alpha _{max}`$, relating to the largest and smallest growth probabilities, respectively. The maximal value $`\alpha _{max}`$ was a subject of a long controversy that was settled only recently (cf. 01DJLMP ; 02JLMP  and references therein). The issue of $`\alpha _{min}`$ appears to be one of the last of the multifractal properties of DLA that has resisted settling. This is the subject of this Letter.

Consider DLA clusters containing $`n`$ particles of radius $`\sqrt{\lambda _0}`$, and denote the radius of the minimal circle that contains the cluster as $`R_n`$. An incoming random walker from infinity has some probability to hit any of the existing particles of the cluster. Denote the maximal of these probabilities as $`p_{max}`$. The average of this probabilities over many clusters of $`n`$ particles appears to scale as

$$p_{max}\left(\frac{\sqrt{\lambda }_0}{R_n}\right)^{\alpha _{min}}n^{\alpha _{min}/D},$$
(3)
where for the last step we have used the obvious scaling law $`n(R_n/\sqrt{\lambda _0})^D`$. Turkevich and Scher have made the plausible assumption that the position of the cluster particle associated with $`p_{max}`$ is at the outermost tip of the cluster. Thus, a scaling relation can be derived by stating that upon adding a new particle to the cluster, $`R_n`$ will grow by one unit $`\sqrt{\lambda _0}`$ with probability $`p_{max}`$, or will not grow at all with probability $`1p_{max}`$. Then

$$\frac{dR}{dn}\sqrt{\lambda _0}p_{max}n^{1/D1},$$
(4)
where again the last step stems from the definition of the fractal dimension. Equating the RHS of (3) and (4) we get the Turkevich-Scher conjecture

$$D=1+\alpha _{min}.$$
(5)
We will show here that this conjecture is incorrect simply because the position of maximal probability is not at the outermost tip of the DLA cluster. In fact, one can introduce in analogy to Eq. (3) a scaling law for the probability to hit the actual tip of the cluster (the particle which is furthest away from the origin), i.e.

$$p_{tip}\left(\frac{\sqrt{\lambda }_0}{R_n}\right)^{\alpha _{tip}}n^{\alpha _{tip}/D}.$$
(6)
A scaling law

$$D=1+\alpha _{tip},$$
(7)
is then a tautology. We will show that for DLA $`\alpha _{tip}>\alpha _{min}`$.

To achieve accurate estimates of $`\alpha _{min}`$ (and in passing of $`\alpha _{tip}`$) we resort to the method of iterated conformal maps that was shown to be extremely useful for dealing with DLA and related growth processes. The method was amply described before, so we just remind the reader that it is based on compositions of fundamental conformal maps $`\varphi _{\lambda ,\theta }`$ which map the exterior of the unit circle to its exterior, except for a little bump at $`e^{i\theta }`$ of linear size proportional to $`\sqrt{\lambda }`$. The composition of these mappings is analogous to the aggregation of random walkers in the off-lattice DLA model. We shall here use the mapping introduced in 98HL  which produces two square root singularities which we refer to as the branch cuts, and the tip of the bump which we refer to as the micro tip. The dynamics is given by

$$\mathrm{\Phi }^{(n)}(w)=\mathrm{\Phi }^{(n1)}(\varphi _{\lambda _n,\theta _n}(w)).$$
(8)
where $`\mathrm{\Phi }^{(n)}`$ maps the exterior of the unit circle to the exterior of the cluster of $`n`$ bumps. The size of the $`n`$’th bump is controlled by the parameter $`\lambda _n`$ and in order to achieve particles of fixed size we have that, to leading order,

$$\lambda _n=\frac{\lambda _0}{|\mathrm{\Phi }_{}^{(n1)}{}_{}{}^{}(e^{i\theta _n})|^2}.$$
(9)
Using the iterated conformal maps it is very easy to keep track of where the maximum growth probability is located, and where the outermost tip is, as more particles are added. Let us assume that at the ($`n`$-1)’th growth step the site with the largest probability is located at the angle $`\theta _{max}`$ on the unit circle, i.e. for all $`\theta `$

$$\frac{1}{|\mathrm{\Phi }_{}^{(n1)}{}_{}{}^{}(e^{i\theta _{max}})|}\frac{1}{|\mathrm{\Phi }_{}^{(n1)}{}_{}{}^{}(e^{i\theta })|}$$
(10)
When we add a new bump in the $`n`$’th growth step the position of maximal probability may not change (up to reparameterization of the angle $`\theta _{max}`$), or move to the new bump. We can easily find the reparameterized angle and determine the new position from

$$p_{max,n}=\mathrm{max}\{\frac{1}{|\mathrm{\Phi }_{}^{(n)}{}_{}{}^{}(\varphi _{\lambda _n,\theta _n}^1(e^{i\theta _{max}}))|},\frac{1}{|\mathrm{\Phi }_{}^{(n)}{}_{}{}^{}(e^{i\theta _n})|}\}.$$
(11)
If $`p_{max,n}`$ is located at $`\theta _n`$ we put $`\theta _{max}=\theta _n`$ in the $`(n+1)`$’th growth step. Similarly we track the position $`|z|_{max}`$ on the cluster by finding the value $`\theta _{tip}`$ which assign the maximal value of $`|\mathrm{\Phi }^{(n)}(e^{i\theta })|`$. We compute $`p_{tip}`$ there as $`1/|\mathrm{\Phi }_{}^{(n)}{}_{}{}^{}|`$.

A direct measurement of $`\alpha _{min}`$ and $`\alpha _{tip}`$ is displayed in Fig. 1. From the direct measurement we find $`\alpha _{min}0.681`$ while $`\alpha _{tip}0.713`$. Clearly the latter is in agreement with (7) while the former is in disagreement with (5) (taking as a datum the result of 00DLP , $`D=1.713\pm 0.005`$).

The direct measurement, while correct, cannot guarantee that very slow convergence of the power laws as a function of $`n`$ may somehow hide an asymptotic identity of $`\alpha _{min}`$ and $`\alpha _{tip}`$. To remove this worry we adopt now the scaling function technique of 00DLP  to achieve an accurate determination of $`\alpha _{min}`$. In this approach one acknowledges that Eq. (3) may be realized only asymptotically for high values of $`n`$. For low and medium values of $`n`$, $`p_{max}`$, which is a function of the discrete $`n`$ and of $`\lambda _0`$, is in fact a scaling function of a single scaling variable,

$`p_{max}=f_{\delta ,\beta }(x),`$
$`x={\displaystyle \frac{1}{\sqrt{\lambda }_0}}(n+\delta )^\beta ,`$ (12)
where we have denoted $`\beta =\alpha _{min}/D`$. The difference with Eq. (3) is that $`f_{\delta ,\beta }`$ is in general not a linear function of its argument, except at exceedingly small values of $`x`$, when $`n`$ is very large. In Fig. 2 we demonstrate the existence of the scaling function and the excellent data collapse achieved using it.

In the upper panel we plot $`f_{\delta ,\beta }(\lambda _0,n)`$ for five values of $`\lambda _0`$ and $`1n250`$. In the lower panel the same data are collapsed using the single scaling variable. We draw the reader’s attention to the following two observations: (i) the data collapse is available immediately, even for the smallest values of $`n`$ 00DLP , and (ii) the scaling function is not linear throughout the range explored here. Thus the scaling law (3) is not obeyed yet for values of $`n`$ of the order of a few hundreds. The set of parameters $`\beta `$ and $`\delta `$ which give the best data collapse in the lower panel are $`\beta =.389`$ and $`\delta =.505`$. These parameters are used in the lower panel and give the estimate $`\alpha _{min}=0.666`$ when assuming that the fractal dimension is $`D=1.713`$.

An even more accurate determination of $`\alpha _{min}`$ is achieved next. Taking the data collapse as an evidence for the existence of a scaling function, we conclude that for any two pairs of numbers $`(n,\lambda _0)`$ and $`(\overline{n},\overline{\lambda }_0)`$ that satisfy the equation

$$\frac{1}{\sqrt{\lambda _0}}(n+\delta )^{\alpha _{min}/D}=\frac{1}{\sqrt{\stackrel{~}{\lambda }_0}}(\stackrel{~}{n}+\delta )^{\alpha _{min}/D},$$
(13)
it follows that

$$f_{\delta ,\beta }(\lambda _0,n)=f_{\delta ,\beta }(\stackrel{~}{\lambda }_0,\stackrel{~}{n}).$$
(14)
These equations offer a calculational procedure. We find $`p_{max}`$ for a given $`n`$ and $`\lambda _0`$, and then for another value $`\overline{n}`$ seek the value $`\overline{\lambda }_0`$ for which $`p_{max}`$ is the same. From Eq. (13) we then deduce that

$$\alpha _{min}=\frac{1}{2}D\frac{\mathrm{log}\lambda _0\mathrm{log}\stackrel{~}{\lambda }_0}{\mathrm{log}(n+\delta )\mathrm{log}(\stackrel{~}{n}+\delta )}$$
(15)
In Fig. 3 we present the results of such a calculation with $`\overline{n}=n+1`$, and $`1n250`$. Since $`\delta `$ is not known a priori, we used the value $`\delta =0.505`$ which was extracted from the data collapse in Fig. 2. We checked the sensitivity to $`\delta `$ by bracketing the results with $`\delta =0`$ and $`\delta =1`$ respectively. The data in Fig. 3 correspond to $`\delta =0.505`$. Fitting the data with a cubic polynomial and extrapolating to $`x0`$ we get the value $`\alpha _{min}0.665`$. On the other hand if we repeat the procedure using the values of $`0\delta 1`$ we are able to bracket the estimate in the interval

$$0.662<\alpha _{min}<0.669.$$
(16)
We thus conclude the analysis with the estimate $`\alpha _{min}=0.665\pm 0.004`$.

Finally, we explain why the Turkevich-Scher conjecture (5) fails. The reason is that the points corresponding to $`p_{max}`$ and $`p_{tip}`$ are not at all the same in typical DLA. In Fig. 4 we present the calculated value of $`R_n`$, computed from the position of largest $`|z|`$ on the cluster, compared with the position corresponding to the maximal harmonic measure. We see that the position of maximal probability fluctuates wildly, and the fluctuations do not appear to go down with the increase in the cluster size. The loss of the conjecture $`(\text{5})`$ means that there is no clear connection between the spectrum of singularities $`f(\alpha )`$ and the fractal dimension of DLA. As said above, the relation (7) is a tautology once the existence of the scaling law (6) has been established 02BS2 . Since the value of $`\alpha _{tip}`$ has nothing to do with the edge of the $`\alpha `$ spectrum, it appears as hard to determine it from first principles as to determine the dimension $`D`$ itself.