text:
# Variational method for finding periodic orbits in a general flow

## I Introduction

The periodic orbit theory of classical and quantum chaos Ruelle (1978); Gutzwiller (1990) is one of the major advances in the study of long-time behavior of chaotic dynamical systems. The theory expresses all long time averages over chaotic dynamics in terms of cycle expansions Artuso et al. (1990), sums over periodic orbits (cycles) ordered hierarchically according to the orbit length, stability, or action. If the symbolic dynamics is known, and the flow is hyperbolic, longer cycles are shadowed by the shorter ones, and cycle expansions converge exponentially or even super-exponentially with the cycle length Rugh (1992).

A variety of methods for determining all periodic orbits up to a given length have been devised and successfully implemented for low-dimensional systems Cvitanović et al. (2003); Viswanath ; Diakonos et al. (1998); Biham and Wenzel (1989); Davidchack and Lai (1999). For more complex dynamics, such as turbulent flows Frisch (1996), nonlinear waves Aranson and Kramer (2002), or quantum fields Biró et al. (1994); Cvitanović (2000) with high (or infinite) dimensional phase spaces and complicated dynamical behavior, many of the existing methods become unfeasible in practice. In the most computationally demanding calculation carried out so far, Kawahara and Kida Kawahara and Kida (2001) have found two periodic solutions in a $`15,422`$-dimensional discretization of a turbulent plane Couette flow. The topology of high-dimensional flows is hard to visualize, and even with a decent starting guess for the shape and location of a periodic orbit, methods like the Newton-Raphson method are likely to fail. In ref. Cvitanović and Lan (2003) we have argued that variational, cost-function minimization methods offer a robust alternative. Here we derive, implement and discuss in detail one such new variational method for finding periodic orbits in general flows, and specifically high-dimensional flows.

In essence, any numerical algorithm for finding periodic orbits is based on devising a new dynamical system which possesses the desired orbit as an attracting fixed point with a sizable basin of attraction. Beyond that, there is much freedom in constructing such system.

For example, the multipoint shooting method eliminates the long-time exponential instability of unstable orbits by splitting an orbit into a number of short segments, each with a controllable expansion rate. The multiple shooting combined with the Newton-Raphson method is an efficient tool for locating periodic orbits of maps Christiansen . A search for periodic orbits of a continuous time flow can be reduced to a multiple shooting search for periodic orbits of a set of maps by constructing a set of phase space Poincaré sections such that an orbit leaving one section reaches the next one in a qualitatively predictable manner, without traversing other sections along the way. In turbulent, high-dimensional flows such sequences of sections are hard to come by. One solution might be a large set of Poincaré sections, with the intervening flight segments short and controllable.

Here we follow a different strategy, and discard Poincaré sections altogether; we replace maps between spatially fixed Poincaré sections, by maps induced by discretizing the time evolution into small time steps. For sufficiently small time steps such maps are small deformations of identity. We distribute many points along a smooth loop $`L`$, our initial guess of a cycle location and its topological layout. If both the time steps and the loop deformations are taken infinitesimal, a partial differential equation governs the “Newton descent”, a fictitious time flow of a trial loop $`L`$ into a genuine cycle $`p`$, with exponential convergence in the fictitious time variable. We then use methods developed for solving PDEs to get the solution. Stated succinctly, the idea of our method is to make an informed rough guess of what the desired cycle looks like globally, and then use a variational method to drive the initial guess towards the exact solution. For robustness, we replace the guess of a single orbit point by a guess of an entire orbit. For numerical safety we replace the Newton-Raphson iteration by the “Newton descent”, a differential flow that minimizes a cost function computed as deviation of the approximate flow from the true flow along a smooth loop approximation to a cycle.

In sect. II we derive the partial differential equation which governs the evolution of an initial guess loop toward a cycle and the corresponding cost function. An extension of the method to Hamiltonian systems and systems with higher time derivatives is presented in sect. III. Simplifications due to symmetries and details of our numerical implementation of the method are discussed in sect. IV. In sect. V we test the method on the Hénon-Heiles system, the restricted three body problem, and a weakly turbulent Kuramoto-Sivashinsky system. We summarize our results and discuss possible improvements of the method in sect. VI.

## II The Newton descent method in loop space

### II.1 A variational equation for the loop evolution

A periodic orbit is a solution $`(x,T)`$, $`x^d`$, $`T`$ of the periodic orbit condition

$$f^T(x)=x,T>0$$
(1)
for a given flow or discrete time mapping $`xf^t(x)`$. Our goal is to determine periodic orbits of flows defined by first order ODEs

$$\frac{dx}{dt}=v(x),x^d,(x,v)𝐓$$
(2)
in many (even infinitely many) dimensions $`d`$. Here $``$ is the phase space (or state space) in which evolution takes place, $`𝐓`$ is the tangent bundle Arnol’d (1992), and the vector field $`v(x)`$ is assumed to be smooth (sufficiently differentiable) almost everywhere.

We make our initial guess at the shape and the location of a cycle $`p`$ by drawing a loop $`L`$, a smooth, differentiable closed curve $`\stackrel{~}{x}(s)L`$, where $`s`$ is a loop parameter. As the loop is periodic, we find it convenient to restrict $`s`$ to $`[0,2\pi ]`$, with the periodic condition $`\stackrel{~}{x}(s)=\stackrel{~}{x}(s+2\pi )`$. Assume that $`L`$ is close to the true cycle $`p`$, pick $`N`$ pairs of nearby points along the loop and along the cycle

$`\stackrel{~}{x}_n`$ $`=`$ $`\stackrel{~}{x}(s_n),0s_1<\mathrm{}<s_N<2\pi ,`$
$`x_n`$ $`=`$ $`x(t_n),0t_1<\mathrm{}<t_N<T_p,`$ (3)
and denote by $`\delta \stackrel{~}{x}_n`$ the deviation of a point $`x_n`$ on the periodic orbit $`p`$ from the nearby point $`\stackrel{~}{x}_n`$,

$$x_n=\stackrel{~}{x}_n+\delta \stackrel{~}{x}_n.$$

The deviations $`\delta \stackrel{~}{x}`$ are assumed small, vanishing as $`L`$ approaches $`p`$.

The orientation of the $`s`$-velocity vector tangent to the loop $`L`$

$$\stackrel{~}{v}(\stackrel{~}{x})=\frac{d\stackrel{~}{x}}{ds}$$

is intrinsic to the loop, but its magnitude depends on the (still to be specified) parametrization $`s`$ of the loop.

At each loop point $`\stackrel{~}{x}_nL`$ we thus have two vectors, the loop tangent $`\stackrel{~}{v}_n=\stackrel{~}{v}(\stackrel{~}{x}_n)`$ and the flow velocity $`v_n=v(\stackrel{~}{x}_n)`$. Our goal is to deform $`L`$ until the directions of $`\stackrel{~}{v}_n`$ and $`v_n`$ coincide for all $`n=1,\mathrm{},N`$, $`N\mathrm{}`$, that is $`L=p`$. To match their magnitude, we introduce a local time scaling factor

$$\lambda (s_n)\mathrm{\Delta }t_n/\mathrm{\Delta }s_n,$$
(4)
where $`\mathrm{\Delta }s_n=s_{n+1}s_n,n=1,\mathrm{},N1,\mathrm{\Delta }s_N=2\pi (s_Ns_1)`$, and likewise for $`\mathrm{\Delta }t_n`$. The scaling factor $`\lambda (s_n)`$ ensures that the loop increment $`\mathrm{\Delta }s_n`$ is proportional to its counterpart $`\mathrm{\Delta }t_n+\delta t_n`$ on the cycle when the loop $`L`$ is close to the cycle $`p`$, with $`\delta t_n0`$ as $`Lp`$.

Let $`x(t)=f^t(x)`$ be the state of the system at time $`t`$ obtained by integrating (2), and $`J(x,t)=dx(t)/dx(0)`$ be the corresponding Jacobian matrix obtained by integrating

$$\frac{dJ}{dt}=AJ,A_{ij}=\frac{v_i}{x_j},\text{with }J(x,0)=\mathrm{𝟙}.$$
(5)
Since the point $`x_n=\stackrel{~}{x}_n+\delta \stackrel{~}{x}_n`$ is on the cycle,

$$f^{\mathrm{\Delta }t_n+\delta t_n}(\stackrel{~}{x}_n+\delta \stackrel{~}{x}_n)=\stackrel{~}{x}_{n+1}+\delta \stackrel{~}{x}_{n+1}.$$
(6)
Linearization

$$f^{\delta t}(x)x+v(x)\delta t,f^t(x+\delta x)x(t)+J(x,t)\delta x,$$

of (6) about the loop point $`\stackrel{~}{x}_n`$ and the time interval $`\mathrm{\Delta }t_n`$ to the next cycle point leads to the multiple shooting Newton-Raphson equation, for any step size $`\mathrm{\Delta }t_n`$:

$$\delta \stackrel{~}{x}_{n+1}J(\stackrel{~}{x}_n,\mathrm{\Delta }t_n)\delta \stackrel{~}{x}_nv_{n+1}\delta t_n=f^{\mathrm{\Delta }t_n}(\stackrel{~}{x}_n)\stackrel{~}{x}_{n+1}.$$
(7)
Provided that the initial guess is sufficiently good, the Newton-Raphson iteration of (7) generates a sequence of loops $`L`$ with a decreasing cost function Cvitanović and Lan (2003)

$$F^2(\stackrel{~}{x})\frac{N}{(2\pi )^2}\underset{i=1}{\overset{N}{}}(f^{\mathrm{\Delta }t_n}(\stackrel{~}{x}_n)\stackrel{~}{x}_{n+1})^2,\stackrel{~}{x}_{N+1}=\stackrel{~}{x}_1.$$
(8)
The prefactor $`N/(2\pi )^2`$ makes the definition of $`F^2`$ consistent with (13) in the $`N\mathrm{}`$ limit. If the flow is locally strongly unstable, the neighborhood in which the linearization is valid could be so small that the full Newton step would overshoot, rendering $`F^2`$ bigger rather than smaller. In this case the step-reduced, damped Newton method is needed. As proved in ref. Keller (1992), under conditions satisfied here, $`F^2`$ decreases monotonically if appropriate step size is taken. If infinitesimal steps are taken, decrease of $`F^2`$ is ensured. We parametrize such continuous deformations of the loop by a fictitious time $`\tau `$.

We fix $`\mathrm{\Delta }s_n`$ and proceed by $`\delta \tau `$ each step of the iteration, that is, multiply the right hand side of (7) by $`\delta \tau `$. According to (4), the change of $`\mathrm{\Delta }t_n`$ with respect to $`\tau `$ is equal to $`\delta t_n=\frac{\lambda }{\tau }(s_n,\tau )\delta \tau \mathrm{\Delta }s_n`$. As $`\delta \stackrel{~}{x}_n=\frac{}{\tau }\stackrel{~}{x}(s_n,\tau )\delta \tau `$, dividing both sides of (7) by $`\delta \tau `$ yields

$`{\displaystyle \frac{d\stackrel{~}{x}_{n+1}}{d\tau }}J(\stackrel{~}{x}_n,\mathrm{\Delta }t_n){\displaystyle \frac{d\stackrel{~}{x}_n}{d\tau }}v_{n+1}{\displaystyle \frac{\lambda }{\tau }}(s_n,\tau )\mathrm{\Delta }s_n`$
$`=f^{\mathrm{\Delta }t_n}(\stackrel{~}{x}_n)\stackrel{~}{x}_{n+1}.`$ (9)
In the $`N\mathrm{}`$ limit, the stepsizes $`\mathrm{\Delta }s_n,\mathrm{\Delta }t_n=O(1/N)0`$, and we have

$`v_{n+1}`$ $``$ $`v_n,\stackrel{~}{x}_{n+1}\stackrel{~}{x}_n+\stackrel{~}{v}_n\mathrm{\Delta }s_n,`$
$`J(\stackrel{~}{x}_n,\mathrm{\Delta }t_n)`$ $``$ $`1+A(\stackrel{~}{x}_n)\mathrm{\Delta }t_n,f^{\mathrm{\Delta }t_n}(\stackrel{~}{x}_n)\stackrel{~}{x}_n+v_n\mathrm{\Delta }t_n.`$
Substituting into (9) and using the scaling relation (4), we obtain

$$\frac{^2\stackrel{~}{x}}{s\tau }\lambda A\frac{\stackrel{~}{x}}{\tau }v\frac{\lambda }{\tau }=\lambda v\stackrel{~}{v}.$$
(10)
This PDE, which describes the evolution of a loop $`L(\tau )`$ toward a periodic orbit $`p`$, is the central result of this paper. The family of loops so generated is parametrized by $`\stackrel{~}{x}=\stackrel{~}{x}(s,\tau )L(\tau )`$, where $`s`$ denotes the position along the loop, and the fictitious time $`\tau `$ parametrizes the deformation of the loop, see Fig. 1(a). We refer to this infinitesimal step version of the damped Newton-Raphson method as the “Newton descent”.

The important feature of this equation is that a decreasing cost functional exists. Rewriting (10) as

$$\frac{}{\tau }(\stackrel{~}{v}\lambda v)=(\stackrel{~}{v}\lambda v),$$
(11)
we have

$$\stackrel{~}{v}\lambda v=e^\tau (\stackrel{~}{v}\lambda v)|_{\tau =0},$$
(12)
so the fictitious time $`\tau `$ flow decreases the cost functional

$$F^2[\stackrel{~}{x}]=\frac{1}{2\pi }_{L(\tau )}𝑑s\left(\stackrel{~}{v}(\stackrel{~}{x})\lambda v(\stackrel{~}{x})\right)^2$$
(13)
monotonically as the loop evolves toward the cycle.

At each iteration step the differences of the loop tangent velocities and the dynamical flow velocities are reduced by the Newton descent. As $`\tau \mathrm{}`$, the fictitious time flow alignes the loop tangent $`\stackrel{~}{v}`$ with the dynamical flow vector $`\stackrel{~}{v}=\lambda v`$, and the loop $`\stackrel{~}{x}(s,\tau )L(\tau )`$, see Fig. 1(b), converges to a genuine periodic orbit $`p=L(\mathrm{})`$ of the dynamical flow $`\dot{x}=v(x)`$. Once the cycle $`p`$ is reached, by (4), $`\lambda (s,\mathrm{})=\frac{dt}{ds}(\stackrel{~}{x}(s,\mathrm{}))`$, and the cycle period is given by

$$T_p=_0^{2\pi }\lambda (\stackrel{~}{x}(s,\mathrm{}))𝑑s.$$

Of course, as at this stage we have already identified the cycle, we may pick instead an initial point on $`p`$ and calculate the period by a direct integration of the dynamical equations (2).

### II.2 Marginal directions and accumulation of loop points

Numerically, two perils lurk in a direct implementation of the Newton descent (10).

First, when a cycle is reached, it remains a cycle under a cyclic permutation of the representative points, so on the cycle the operator

$`\overline{A}={\displaystyle \frac{}{s}}\lambda A`$
has a marginal eigenvector $`v(\stackrel{~}{x}(s))`$ with eigenvalue $`0`$. If $`\lambda `$ is fixed, as the loop approaches the cycle, (10) approaches its limit

$$\overline{A}\frac{x}{\tau }=0.$$

Therefore, on the cycle, the operator $`\overline{A}^1`$ becomes singular and the numerical woes arise.

The second potential peril hides in the freedom of choosing the loop (re-)parametrization. Since $`s`$ is related to the time $`t`$ by the yet unspecified factor $`\lambda (s,\tau )`$, uneven distributions of the sampling points over the loop $`L`$ could arise, with the numerical discretization points $`\stackrel{~}{x}_n`$ clumping densely along some segments of $`L`$ and leaving big gaps elsewhere, thus degrading the numerical smoothness of the loop.

We remedy these difficulties by imposing constraints on (10). In our calculation for Kuramoto-Sivashinsky system of sect. V, the first difficulty is dealt with by introducing one Poincaré section, for example, by fixing one coordinate of one of the sampling points, $`\stackrel{~}{x}_1(s_1,\tau )=const`$. This breaks the translational invariance along the cycle. Other types of constraints might be better suited to a specific problem at hand. For example, we can demand that the average displacement of the sampling points along the loop vanishes, thus avoiding a spiraling descent towards the desired cycle.

We deal with the second potential difficulty by choosing a particularly simple loop parametrization. So far, the parametrization $`s`$ is arbitrary and there is much freedom in choosing the best one for our purposes. We pick $`s`$ and $`\tau `$independent constant scaling $`\lambda (s,\tau )=\lambda `$. With uniform grid size $`\mathrm{\Delta }s_n=\mathrm{\Delta }s`$ and fixed $`\lambda `$, the loop parameter $`s=t/\lambda `$ is proportional to time $`t`$, and the discretization (10) distributes the sampling points along the loop evenly in time. As the loop approaches a cycle, $`\frac{\stackrel{~}{x}}{\tau }`$ is numerically obtainable from (10), and on the cycle the period is given by $`T_p=2\pi \lambda `$.

Even though this paper focuses on searches for periodic orbits, the Newton descent is a general method. With appropriate modifications of boundary conditions and scaling of time, (10) can be adapted to determination of homoclinic or heteroclinic orbits between equilibrium points or periodic orbits of a flow, or more general boundary value problems. Applied to 2-point boundary value problems, Newton descent is similar to the quasilinearization Stoer and Bulirsch (1983) but has the advantage that the free parameter $`\lambda (s,\tau )`$ is available for adjusting scales in the problem and that searches can be restricted to phase space submanifolds of interest. A simple example of a restriction to a submanifold are searches for cycles of a given energy, constrained to the $`H(q,p)=E`$ energy shell in the phase space of a Hamiltonian system. Furthermore, as we shall show now, the symplectic structure of Hamilton’s equations greatly reduces the dimensionality of the submanifold that we need to consider.

## III Extensions of Newton descent

In classical mechanics particle trajectories are also solutions of a variational principle, the Hamilton’s variational principle. For example, one can determine a periodic orbit of a billiard by wrapping around a rubber band of a roughly correct topology, and then moving the points along the billiard walls until the length (that is, the action) of the rubber band is extremal (maximal or minimal under infinitesimal changes of the boundary points). Note that the extremization of action requires only $`D`$ configuration coordinate variations, not the full $`2D`$-dimensional phase space variations.

Can we exploit this property of the Newtonian mechanics to reduce the dimenionality of our variational calculations? The answer is yes, and easiest to understand in terms of the Hamilton’s variational principle which states that classical trajectories are extrema of the Hamilton’s principal function (or, for fixed energy $`E`$, the action $`S=R+Et`$)

$$R(q_1,t_1;q_0,t_0)=_{t_0}^{t_1}𝑑t(q(t),\dot{q}(t),t),$$

where $`(q,\dot{q},t)`$ is the Lagrangian. Given a loop $`L(\tau )`$ we can compute not only the tangent “velocity” vector $`\stackrel{~}{v}`$, but also the local loop curvature or “acceleration” vector

$$\stackrel{~}{a}=\frac{^2\stackrel{~}{x}}{s^2},$$

and indeed, as many $`s`$ derivatives as needed. Matching the dynamical acceleration $`a(\stackrel{~}{x})`$ (assumed to be functions of $`\stackrel{~}{x}`$ and $`v(\stackrel{~}{x})`$) with the loop “acceleration” $`\stackrel{~}{a}(\stackrel{~}{x})`$ results in a new cost function and the corresponding PDE (11) for the evolution of the loop

$$\frac{}{\tau }(\stackrel{~}{a}\lambda ^2a)=(\stackrel{~}{a}\lambda ^2a).$$

We use $`\lambda ^2`$ instead of $`\lambda `$ in order to keep the notation consistent with (4), that is $`t=\lambda s`$. Expressed in terms of the loop variables $`\stackrel{~}{x}(s)`$, the above equation becomes

$`{\displaystyle \frac{^3\stackrel{~}{x}}{^2s\tau }}\lambda {\displaystyle \frac{a}{v}}{\displaystyle \frac{^2\stackrel{~}{x}}{s\tau }}\lambda ^2{\displaystyle \frac{a}{\stackrel{~}{x}}}{\displaystyle \frac{\stackrel{~}{x}}{\tau }}+({\displaystyle \frac{a}{v}}{\displaystyle \frac{\stackrel{~}{x}}{s}}2\lambda a){\displaystyle \frac{\lambda }{\tau }}`$
$`=\lambda ^2a\stackrel{~}{a},`$ (14)
where $`v=\frac{\stackrel{~}{x}}{\lambda s}`$. Although (14) looks more complicated than (10), in numerical fictitious time integrations, we are rewarded by having to keep only half of the phase space variables.

More generally, if a differential equation has the form:

$$x^{(m)}=f(x,x^{(1)},\mathrm{},x^{(m1)}),$$
(15)
where $`x^{(k)}=\frac{d^kx}{dt^k},k=1,\mathrm{},m`$ and $`x^d`$, the same technique can be used to match the highest derivatives $`\lambda ^mx^{(m)}`$ and $`\stackrel{~}{x}^{(m)}`$,

$$\frac{}{\tau }(\stackrel{~}{x}^{(m)}\lambda ^mx^{(m)})=(\stackrel{~}{x}^{(m)}\lambda ^mx^{(m)}),$$

with $`\stackrel{~}{x}^{(m)}=\frac{^m}{s^m}\stackrel{~}{x}(s)`$ calculated directly from $`\stackrel{~}{x}(s)`$ on the loop by differentiation. In loop variables $`\stackrel{~}{x}(s)`$ we have,

$`{\displaystyle \frac{^{m+1}\stackrel{~}{x}}{s^m\tau }}\lambda ^m{\displaystyle \underset{k=0}{\overset{m}{}}}{\displaystyle \frac{f}{x^{(k)}}}{\displaystyle \frac{}{\tau }}{\displaystyle \frac{^k\stackrel{~}{x}}{\lambda ^ks^k}}m\lambda ^{m1}\stackrel{~}{x}^{(m)}{\displaystyle \frac{\lambda }{\tau }}`$
$`=\lambda ^mx^{(m)}\stackrel{~}{x}^{(m)},`$ (16)
where $`x=x^{(0)}`$ and $`\stackrel{~}{x}^{(k)}=\frac{^k\stackrel{~}{x}}{\lambda ^k^ks},k=1,\mathrm{},m1`$ are assumed. Conventionally, (15) is converted to a system of $`md`$ first order differential equations, whose discretized derivative (see (17) below) are banded matrices with band width of $`5md`$. Using (16), we only need $`d`$ equations for the same accuracy and the corresponding band width is $`(m+4)d`$. The computing load has been greatly reduced, the more so the larger $`m`$ is. Nevertheless, choice of a good initial loop guess and visualization of the dynamics are always aided by a plot of the orbit in the full $`md`$-dimensional phase space, where loops cannot self-intersect and topological features of the flow is exhibited more clearly.

## IV Implementation of Newton descent

As the loop points satisfy a periodic boundary condition, it is natural to employ truncated discrete Fast Fourier Transforms (FFT) in numerical integrations of (10). Since we are interested only in the final, stationary cycle $`p`$, the accuracy of the fictitious time integration is not crucial; all we have to ensure is the smoothness of the loop throughout the integration. The Euler integration with fairly large time steps $`\delta \tau `$ suffices. The computationally most onerous step in implementation of the Newton descent is the inversion of large matrix $`\overline{A}`$ in (10). When the dimension of the dynamical phase space of (2) is high, the inversion of $`\overline{A}`$ needed to get $`\frac{\stackrel{~}{x}}{\tau }`$ takes most of the integration time, making the evolution extremely slow. This problem is partially solved if the finite difference methods are used. The large matrix $`\overline{A}`$ then becomes sparse and the inversion can be done far more quickly.

### IV.1 Numerical implementation

In a discretization of a loop, numerical stability requires accurate discretization of loop derivatives such as

$$\stackrel{~}{v}_n\frac{\stackrel{~}{x}}{s}|_{\stackrel{~}{x}=\stackrel{~}{x}(s_n)}(\widehat{D}\stackrel{~}{x})_n.$$

In our numerical work we use the four-point approximation Brandenburg (2003),

$$\widehat{D}=\frac{1}{12h}\left(\begin{array}{ccccccccccc}0& 8& 1& & & & & & & 1& 8\\ 8& 0& 8& 1& & & & & & & 1\\ 1& 8& 0& 8& 1& & & & & & \\ & & & & & \mathrm{}& & & & & \\ & & & & & & 1& 8& 0& 8& 1\\ 1& & & & & & & 1& 8& 0& 8\\ 8& 1& & & & & & & 1& 8& 0\end{array}\right)$$
(17)
where $`h=2\pi /N`$. Here, each entry represents a $`[d\times d]`$ matrix, $`88\mathrm{𝟙}`$, etc., with blank spaces are filled with zeros. The two $`[2d\times 2d]`$ matrices

$$M_1=\left(\begin{array}{cc}\mathrm{𝟙}& 8\mathrm{𝟙}\\ 0& \mathrm{𝟙}\end{array}\right),M_2=\left(\begin{array}{cc}\mathrm{𝟙}& 0\\ 8\mathrm{𝟙}& \mathrm{𝟙}\end{array}\right),$$

located at the top-right and bottom-left corners take care of the periodic boundary condition.

The discretized version of (10) with a fictitious time Euler step $`\delta \tau `$ is

$$\left(\begin{array}{cc}\widehat{A}& \widehat{v}\\ \widehat{a}& 0\end{array}\right)\left(\begin{array}{c}\delta \widehat{x}\\ \delta \lambda \end{array}\right)=\delta \tau \left(\begin{array}{c}\lambda \widehat{v}\widehat{\stackrel{~}{v}}\\ 0\end{array}\right),$$
(18)
where

$`\widehat{A}`$ $`=`$ $`\widehat{D}+\mathrm{diag}[A_1,A_2,\mathrm{},A_N],`$
with $`A_n=A(\stackrel{~}{x}(s_n))`$ defined in (5), and

$`\widehat{v}`$ $`=`$ $`(v_1,v_2,\mathrm{},v_N),\text{ with }v_n=v(\stackrel{~}{x}(s_n)),`$
$`\widehat{\stackrel{~}{v}}`$ $`=`$ $`(\stackrel{~}{v}_1,\stackrel{~}{v}_2,\mathrm{},\stackrel{~}{v}_N),\text{ with }\stackrel{~}{v}_n=\stackrel{~}{v}(\stackrel{~}{x}(s_n)),`$
are the two vector fields that we want to match everywhere along the loop. $`\widehat{a}`$ is an $`Nd`$ dimensional row vector which imposes the constraint on the coordinate variations $`\delta \widehat{x}=(\delta \stackrel{~}{x}_1,\delta \stackrel{~}{x}_2,\mathrm{},\delta \stackrel{~}{x}_N)`$. The discretized Newton descent (18) is an infinitesimal time step variant of the multipoint (Poincaré section) shooting equation for flows Christiansen . In order to solve for the deformation of the loop coordinates and period, $`\delta \widehat{x}`$ and $`\delta \lambda `$, we need to invert the $`[(Nd+1)\times (Nd+1)]`$ matrix on the left hand side of (18).

In our numerical work, this matrix is inverted using the banded LU decomposition on the embedded band-diagonal matrix, and the Woodbury formula Press et al. (1992) on the cyclic, border terms. The LU decomposition takes most of the computational time and considerably slows down the fictitious time integration. We speed up the integration by a new inversion scheme which relies on the smoothness of the flow in the loop space. It goes as follows. Once we have the LU decomposition at one step, we use it to approximately invert the matrix in the next step, with accurate inversion achieved by the iterative approximate inversions Press et al. (1992). In our applications we find that a single LU decomposition can be used for many $`\delta \tau `$ evolution steps. The further we go, the more iterations at each step are needed to implement the inversion. After the number of such iterations exceeds some given fixed maximum number, we perform another LU decomposition and proceed as before. The number of integration steps following one decomposition is an indication of the smoothness of the evolution, and we adjust accordingly the integration step size $`\delta \tau `$: the greater the number, the bigger the step size. As the loop approaches a cycle, the evolution becomes so smooth that the step size can be brought all the way up to $`\delta \tau =1`$, the full undamped Newton-Raphson iteration step. In practice, one can start with a small but reasonable number of points, in order to get a coarse solution of relatively low accuracy. After achieving that, the refined guess loop can be constructed by interpolating more points, and proceed with for a more accurate calculation in which $`\delta \tau `$ can be set as large as the full Newton step $`\delta \tau =1`$, recovering the rapid quadratic convergence of the Newton-Raphson method.

It is essential that the smoothness of the loop is maintained throughout the calculation. We monitor the smoothness by checking the Fourier spectrum of $`\stackrel{~}{x}(,\tau )`$. An unstable difference scheme for loop derivatives might lead to unbounded sawtooth oscillations Thomas (1995). A heuristic local linear stability analysis (described in Roždestvenskiǐ and Janenko (1983)) indicates that our scheme is stable, and that the high frequency components do not generate instabilities.

### IV.2 Initialization and convergence

As in any other method, a qualitative understanding of the dynamics is a prerequisite to successful cycle searches. We start by numerical integration with the dynamical system (2). Numerical experiments reveal regions where a trajectory spends most of its life, giving us the first hunch as to how to initialize a loop. We take the FFT of some nearly recurred orbit segment and keep only the lowest frequency components. The inverse Fourier transform back to the phase space yields a smooth loop that we use as our initial guess. Since any generic orbit segment is not closed and might exhibit large gaps, the Gibbs phenomenon can take the initial loop so constructed quite far away from the region of interest. We deal with this problem by manually deforming the orbit segment into a closed loop before performing the FFT. Searching for longer cycles with multiple circuits requires more delicate initial conditions. The hope is that a few short cycles can help us establish an approximate symbolic dynamics, and guess for longer cycles can be constructed by cutting and glueing the short, known ones. For low dimensional systems, such methods yield quite good systematic initial guesses for longer cycles Christiansen et al. (1997).

An alternative way to initialize the search is by utilizing adiabatic deformations of dynamics, or the homotopy evolution Conley (1978). If the dynamical system (2) depends on a parameter $`\mu `$, short cycles might survive as $`\mu `$ varies passing through a family of dynamical systems, giving in the process birth to new cycles through sequences of bifurcations. Most short unstable cycles vary little for small changes of $`\mu `$. So, a cycle existing for parameter value $`\mu _1`$ can be chosen as the initial trial loop for a nearby cycle surviving a small change $`\mu _1\mu _2`$. In practice, one or two iterations often suffice to find the new cycle.

A good choice of the initial loop significantly expedites the computation, but there are more reasons why good initial loops are crucial. First of all, if we break the translational invariance by imposing a constraint such as $`\stackrel{~}{x}_1(s_1,\tau )=c`$, we have to make sure that both the initial loop and the desired cycle intersects this Poincaré plane. Hence, the initial loop cannot be wildly different from the desired cycle. Second, in view of (12), the loop always evolves towards a local minimum of the cost functional (13), with discretization points moving along the $`\stackrel{~}{v}\lambda v`$ fixed direction, determined by the initial condition. If the local minimum corresponds to a zero of the cost functional, we obtain a true cycle of the dynamical flow (2). However, if the value of cost functional is not equal to zero at the minimum while the gradient is zero, (18) yields a singular matrix $`\widehat{A}`$. In such cases the search has to be abandoned and restarted with a new initial loop guess. In the periodic orbit searches of sect. V starting with blind initial guesses (guesses unaided by a symbolic dynamics partition), such local minima were encountered in about $`30\%`$ of cases.

### IV.3 Symmetry considerations

The system under consideration often possesses certain symmetries. If this is the case, the symmetry should be both be feared for possible marginal eigen-directions, and be embraced as a guide to possible simplifications of the numerical calculation.

If the dynamical system equations (2) are invariant under a discrete symmetry, the concept of fundamental domain Cvitanović and Eckhardt (1993); Cvitanović et al. (2003) can be utilized to reduce the length of the initial loop when searching for a cycle of a given symmetry. In this case, we need discretize only an irreducible segment of the loop, decreasing significantly the dimensionality of the loop representation. Other parts of the loop are replicated by symmetry operations, with the full loop tiled by copies of the fundamental domain segment. The boundary conditions are not periodic any longer, but all that we need to do is modify the cyclic terms. Instead of using $`M_1`$ and $`M_2`$ in (17), we use $`M_1Q`$ and $`M_2Q^1`$, where $`Q`$ is the relevant symmetry operation that maps the fundamental segment to the neighbor that precedes it. In this way, a fraction of the points represent the cycle with the same accuracy, speeding up the search considerably.

If a continuous symmetry is present, it may complicate the situation at first glance but becomes something that we can take advantage of after careful checking. For example, for a Hamiltonian system unstable cycles may form continuous families Hénon (2001); Sepulchre and Mackay (1997), with one or more members of a family belonging to a given constant energy surface. In order to cope with the marginal eigendirection associated with such continous family, we search for a cycle on a particular energy surface by replacing the last row of equation (18) by an energy shell constraint Christiansen . We put one point of the loop, say $`\stackrel{~}{x}_2`$, on the constant energy surface $`H(\stackrel{~}{x})=E`$, and impose the constraint $`H(\stackrel{~}{x}_2)\delta \stackrel{~}{x}_2=0`$, so as to keep $`\stackrel{~}{x}_2`$ on the surface for all $`\tau `$. The integration of (10) then automatically brings all other loop points to the same energy surface. Alternatively, we can look for a cycle of given fixed period $`T`$ by fixing $`\lambda `$ and dropping the constraint in the bottom line of (18). These two approaches are conjugate to each other, both needed in applications. In most cases, they are equivalent. One exception is the harmonic oscillator for which the oscillations have identical period but different energy. Note that in both cases the translational invariance is restored, as we have discarded the Poincaré section condition of sect. II.2. As explained in Viswanath , this causes no trouble in numerical calculations.

## V Applications

We have checked that the iteration of (18) yields quickly and robustly the short unstable cycles for standard models of low-dimensional dissipative flows such as the Rössler system Rössler (1976). A more daunting challenge are searches for cycles in Hamiltonian flows, and searches for spatio-temporally periodic solutions of PDEs. In all numerical examples that follow, the convergence condition is $`F^2<10^5`$.

### V.1 Hénon-Heiles system and restricted three-body problem

First, we test the Hamiltonian version of the Newton descent derived in Sect. III by applying the method to two Hamiltonian systems, both with two degrees of freedom. In both cases, our initial loop guesses are rather arbitrary combinations of trigonometric functions. Nevertheless, the observed convergence is fast.

The Hénon-Heiles system Hénon and Heiles (1964) is a standard model in celestial mechanics, described by the Hamiltonian

$$H=\frac{1}{2}(p_x^2+p_y^2+x^2+y^2)+x^2y\frac{y^3}{3}.$$
(19)
It has a time reversal symmetry and a three-fold discrete spatial symmetry. Figure 2 shows a typical application of (14), with the Newton descent search restricted to the configuration space. The initial loop, Fig. 2(a), is a rather coarse initial guess. We fix arbitrarily the scaling $`\lambda =2.1`$, that is, we search for a cycle $`p`$ of the fixed period $`T_p=13.1947`$, with no constraint on the energy. Figure 2(b) shows the cycle found by the Newton descent, with energy $`E=0.1794`$, and the full discrete symmetry of the Hamiltonian. This cycle persists adiabatically for a small range of values of $`\lambda `$; with $`\lambda `$ changed much, the Newton descent takes the same initial loop into other cycles. Figure 2(c) verifies that the cost functional $`F^2`$ decreases exponentially with slope -2 throughout the $`\tau =[0,10]`$ integration interval, as predicted by (12). The points get more and more sparse as $`\tau `$ increases, because our numerical implementation adaptively chooses bigger and bigger step sizes $`\delta \tau `$.

In the Hénon-Heiles case, the accelerations $`a_x,a_y`$ depend only on the configuration variables $`x,y`$. More generally, the accelerations could also depend on $`\dot{x},\dot{y}`$. Consider as an example the equations of motion for the restricted three-body problem Szebehely (1967),

$`\ddot{x}`$ $`=`$ $`2\dot{y}+x(1\mu ){\displaystyle \frac{x+\mu }{r_1^3}}\mu {\displaystyle \frac{x1+\mu }{r_2^3}},`$
$`\ddot{y}`$ $`=`$ $`2\dot{x}+y(1\mu ){\displaystyle \frac{y}{r_1^3}}\mu {\displaystyle \frac{y}{r_2^3}},`$ (20)
where $`r_1=\sqrt{(x+\mu )^2+y^2},r_2=\sqrt{(x1+\mu )^2+y^2}`$. These equations describe the motion of a test particle in a rotating frame under the influence of the gravitational force of two heavy bodies with masses $`1`$ and $`\mu 1`$ fixed at $`(\mu ,0)`$ and $`(1\mu ,0)`$ in the $`(x,y)`$ coordinate frame. The stationary solutions of (20) are called the Lagrange points, corresponding to a circular motion of the test particle in phase with the rotation of the heavy bodies. The periodic solutions in the rotating frame correspond to periodic or quasi-periodic motion of the test particle in the inertial frame. Figure 3 shows an initial loop and the cycle to which it converges, in the rotating frame. Although the cycle looks simple, the Newton descent requires advancing in small $`\delta \tau `$ steps in order for the initial loop to converge to it.

In order to successfully apply the Hamiltonian version of the Newton descent (14), we have to ensure that the test particle keeps a finite distance from the origin. If a cycle passes very close to one of the heavy bodies, the acceleration can become so large that our scheme of uniformly distributing the loop points in time might fail to represent the loop faithfully. Another distribution scheme is required in this case, for example, making the density of points proportional to the magnitude of acceleration.

### V.2 Periodic orbits of Kuramoto-Sivashinsky system

The Kuramoto-Sivashinsky equation arises as an amplitude equation for interfacial instability in a variety of contexts Kuramoto and Tsuzuki (1976); Sivashinsky (1977). In 1-dimensional space, it reads

$$u_t=(u^2)_xu_{xx}\nu u_{xxxx},$$
(21)
where $`\nu `$ is a “super-viscosity” parameter which controls the rate of dissipation and $`(u^2)_x`$ is the nonlinear convection term. As $`\nu `$ decreases, the system undergoes a series of bifurcations, leading to increasingly turbulent, spatio-temporally chaotic dynamics.

If we impose the periodic boundary condition $`u(t,x+2\pi )=u(t,x)`$ and choose to study only the odd solutions $`u(x,t)=u(x,t)`$, $`u(x,t)`$ can be expanded in a discrete spatial Fourier series Christiansen et al. (1997),

$$u(x,t)=i\underset{k=\mathrm{}}{\overset{\mathrm{}}{}}a_k(t)e^{ikx},$$
(22)
where $`a_k=a_k`$. In terms of the Fourier components, PDE (21) becomes an infinite ladder of ODEs,

$$\dot{a_k}=(k^2\nu k^4)a_kk\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}a_ma_{km}.$$
(23)
In numerical simulations we work with the Galerkin truncations of the Fourier series since in the neighborhood of the strange attractor the magnitude of $`a_k`$ decreases very fast with $`k`$, high frequency modes playing a negligible role in the asymptotic dynamics. In this way Galerkin truncations reduce the dynamics to a finite but large number of ODEs. We work with $`d=32`$ dimensions in our numerical calculations. In ref. Christiansen et al. (1997), multipoint shooting has been successfully applied to obtain periodic orbits close to the onset of spatiotemporal chaos ($`\nu =0.03`$). In this regime, our method is so stable that big time steps $`\delta \tau `$ can be employed even at the initial guesses, leading to extremely fast convergence. We attribute this robustness to the simplicity of the structure of the attractor at high viscosity values.

The challenge comes with decreasing $`\nu `$, with the dynamics turning more and more turbulent. Already at $`\nu =0.015`$, the system is moderately turbulent and the phase space portraits of the dynamics reveal a complex labyrinth of “eddies” of different scales and orientations. While the highly unstable nature of orbits and intricate structure of the invariant set hinder applications of conventional cycle-search routines, in this setting our variational method shines through. We design rather arbitrary initial loops from numerical trajectory segments, and the calculation proceeds as before, except that now a small $`\delta \tau `$ has to be used initially to ensure numerical stability. Topologically different loops are very likely to result in different cycles, while some initial loop guesses my lead to local nonzero minima of the cost functional $`F^2`$. As explained in Sect. IV, in such cases the method diverges, and the search is restarted with a new inital loop guess.

Two initial loop guesses are displayed in Fig. 4, along with the two periodic orbits detected by the Newton descent. In discretization of the initial loops, each point has to be specified in all $`d`$ dimensions; here the coordinates $`\{a_1,a_2\}`$ are picked so that topological similarity between initial and final loops is visually easy to identify. Other projections from $`d=32`$ dimensions to subsets of 2 coordinates appear to make the identification harder, if not impossible. In both calculations, we molded segments of typical trajectories into smooth closed loops by the Fourier filtering method of Sect. IV. As the desired orbit becomes longer and more complex, more sampling points are needed to represent the loop. We use $`N=512`$ points to represent the loop in the (a)-(b) case and $`N=1024`$ points in the (c)-(d) case. The space-time evolution of $`u(x,t)`$ for these two unstable spatio-temporally periodic solutions is displayed in Fig. 5. As $`u(x,t)`$ is antisymmetric on $`[\pi ,\pi ]`$, it suffices to display the solutions on the $`x[0,\pi ]`$ interval.

## VI Discussion

In order to cope with the difficulty of finding periodic orbits in high-dimensional chaotic flows, we have devised the Newton descent method, an infinitesimal step variant of the damped Newton-Raphson method. Our main result is the PDE (10) which solves the variational problem of minimizing the cost functional (13). This equation describes the fictitious time $`\tau `$ flow in the space of loops which decreases the cost functional at uniform exponential rate (see (12)). Variants of the method are presented for special classes of systems, such as Hamiltonian systems. An efficient integration scheme for the PDE is devised and tested on the Kuramoto-Sivashinsky system, the Hénon-Heiles system and the restricted three-body problem.

Our method uses information from a large number of points in phase space, with the global topology of the desired cycle protected by insistence on smoothness and a uniform discretization of the loop. The method is quite robust in practice.

The numerical results presented here are only a proof of principle. We do not know to what periodic orbit the flow (10) will evolve for a given dynamical system and a given initial loop. Empirically, the flow goes toward the “nearest” periodic orbit, with the largest topological resemblance. Each particular application will still require much work in order to elucidate and enumerate relevant topological structures. The hope is that the short spatio-temporally periodic solutions revealed by the Newton descent searches will serve as the basic building blocks for systematic investigations of chaotic and perhaps even “turbulent” dynamics.

###### Acknowledgements.

We would like to thank Cristel Chandre for careful reading of the manuscript and numerous suggestions.


format: text

text:
# A V-Band Survey for Variable Galactic Nuclei in the Hubble Deep Field

## 1 INTRODUCTION

An important goal of modern cosmology is to understand the evolution of active galaxies and their relationship to normal galaxy evolution. A key piece of information needed to further this understanding is an accurate knowledge of the luminosity function (LF) for AGN over a wide range of absolute magnitudes and redshifts. The AGN LF is populated by quasars at the bright end (e.g. Hartwick & Schade 1990) and Seyfert galaxy nuclei, considered to be their intrinsically fainter counterparts, at the low luminosity end (Cheng $`etal.`$ 1985; Huchra & Burg 1992; Maiolino & Rieke 1995; Koehler $`etal.`$ 1997; Ulvestad & Ho 2001; Londish $`etal.`$ 2002).

Understanding how the faint end of the AGN LF evolves is of particular importance for determining the frequency and total space density of AGNs at earlier epochs. This has obvious implications for determining their total contribution to the X-ray, IR and UV backgrounds. Large numbers of low-luminosity AGN have been proposed to explain the ionization of the intergalactic medium at high redshift (Steidel & Sargent 1989), although recent observational and theoretical results differ on the amount of their true contribution (Barger 2003; Schirber & Bullock 2003). Additionally, the faint end is an important constraint on evolutionary models for the AGN LF such as pure luminosity and luminosity-dependent density evolution (e.g. Boyle $`etal.`$ 2000).

One method to identify AGN in imaging surveys is via their variable nature. QSOs and Seyfert galaxies have long been known as variable objects, with significant optical flux changes occurring on timescales of months to several years. Variability has been a successful method for identifying AGNs, primarily QSOs (Hawkins 1986; Koo $`etal.`$ 1986; Hook $`etal.`$ 1994). A survey for variable sources in SA57 revealed many optically extended, Seyfert-like galaxies (Bershady $`etal.`$ 1998) which generally had higher variability amplitudes than the more luminous QSOs. This result suggests that variability is a good technique for selecting intrinsically faint QSOs and Seyfert nuclei.

While variability and other selection techniques (spectroscopy, color selection, etc.) have been successful in finding QSOs out to high redshifts, the intrinsically fainter AGN (i.e. Seyfert galaxies) have been difficult to detect beyond the local universe. Ground-based variability surveys, like spectroscopic surveys, are quickly limited to AGN which dominate the host galaxy light at higher redshift. Recent X-ray surveys with Chandra and XMM-Newton have revealed hundreds of X-ray sources at z$``$1. Many of these are confirmed AGN but still many more are ambiguous in nature (e.g Barger $`etal.`$ 2002). In summary, the population of moderate-to-high redshift, low-luminosity AGN has not been well surveyed and the evolution of the faint end of the AGN LF is still poorly known.

In this study, we have searched for optical nuclear variability in HDF111HDF refers to the Hubble Deep Field North galaxies to identify the population of active galaxies out to z$``$1. The primary data sets analyzed for galactic variability are the original HDF observations (Williams $`etal.`$ 1996) obtained in December 1995 and a second V-band (F606W) image obtained almost exactly 5 years later in December of 2000. The main advantage of HST over ground-based surveys is the ability to do accurate photometry within smaller apertures than can be done from the ground thus allowing us to probe much lower AGN/host galaxy luminosity ratios.

We discuss the details of the image processing and photometry in Section 2, followed by a discussion of the selected varying galactic nuclei in Sections 3 and 4. In Section 5, we cross-reference our sources with various multi-wavelength surveys. A discussion of the optical spectra for the varying sources is presented in Section 6. In Section 7, we present the luminosity function for variable nuclei. We close with a summary of our conclusions in Section 8.

## 2 PHOTOMETRY

### 2.1 Image Processing

The original epoch HDF data in F606W consisted of 108 individual exposures with a total integration time of 113,050 seconds. These data were obtained with a total of 9 unique dithers covering both sub-pixel offsets and multiple pixels on a scale of about $`\pm `$ 10 WF pixels. The primary data processing technique was detailed in Gilliland, Nugent, and Phillips (1999) for creation of over-sampled, by a factor of four, summed images with resulting noise and resolution properties competitive with, and sometimes slightly better than those reported in Williams $`etal.`$ (1996). The mean epoch of the original F606W data was MJD 50076.9. Exposure times used for the original epoch ranged from 350 to 1500 seconds with a mean sky rate of 0.0758 e-/s.

The second epoch HDF data in F606W was taken with GO-8389 (PI Rodrigo Ibata) with a science goal distinct from that here, but one that also required precise differential measurements across the two epochs. These data consisted of 70 individual exposures each of 1300 seconds for a total exposure time of 91,000 seconds. The goal was to have the new data obtained at the same pointing as the original HDF, and at the same orientation, and with dithering over the same spatial scale, but with an improved level of sampling using 35 unique dithers. This goal was quite well met, thus resulting in quite similar data in the two F606W epochs. The mean epoch of the second F606W data was MJD 51882.9 implying a delta of 4.94 years. The mean sky rate for the newer epoch was 0.0958 e-/s. Processing of the second epoch data was identical to that of the first, using the same software as described in Gilliland, Nugent and Phillips (1999). To facilitate simple, robust analysis of data from the two epochs the registrations of the second epoch were evaluated using a grid of galaxies from the original epoch (stars could not be used since these show unique proper motions over five years); this results in second epoch over-sampled combined images that are accurately aligned to the first epoch with maximum errors of $``$0.05 pixels. Therefore, comparative analysis for change across epochs can be made using precisely the same pixels, mapping back to common physical domains. For both epochs, each individual frame was sky subtracted before combining, cosmic rays were eliminated (effected pixels assigned zero weight in the combination) and hot pixels subtracted.

The extreme stability provided by HST is an essential factor in being able to probe to unique sensitivity levels with these extensive observations. Limiting factors, to be discussed further below, arise from (1) a time-varying Charge Transfer Efficiency (CTE) that must be empirically corrected for, and (2) limited sub-pixel dithering in the first HDF epoch that results in minor mismatches noticeable for the objects with tight cores at high signal to noise.

### 2.2 Source Selection and Aperture Size

We created a catalog of HDF objects using DAOFIND in IRAF222 IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. to ensure a systematic selection of all sources in the V-band images. All objects above the 8$`\sigma `$ detection threshold were identified in the 1995 images (the deeper of the two) using a kernel size of 1.8$`\times `$ the semi-major axis of the Gaussian convolution kernel (1.3 pixels) to allow more sensitivity to extended objects. The catalog was then edited to remove a minority of objects that did not appear to be centered on galaxies (i.e. multiple detections within single galaxies). For very amorphous galaxies with no clear center, all bright “knots” appearing within the source were retained in the catalog.

For our photometry, we chose a small aperture consistent with the FWHM of stellar sources in the field. Stars have FWHM$``$5.5 pixels in the over-sampled images corresponding to $``$0.14$`^{\prime \prime }`$. We chose an aperture size of 6 pixels in diameter to be most sensitive to photometric changes in an unresolved nucleus. This approach differs from that taken in Sarajedini, Gilliland & Phillips (2000; hereafter SGP), where the aperture size was scaled to the size of the galaxy. Our goal here is to provide better sensitivity to variable active nuclei with less dilution from the underlying, non-varying host galaxy light.

Aperture photometry for all sources was performed on the direct images using the PHOT algorithm in IRAF. Because the second epoch image was aligned with the first, exactly the same pixels were used in both epochs. During image processing, the background was subtracted and the images were scaled to an exposure time of 6000s. This allowed the sky value to be held constant at zero for both epochs and all WF CCDs. The resulting nuclear magnitudes for each galaxy in the HDF were then adjusted by a zero-point offset to bring them into agreement with the integrated, isophotal galaxy V-band AB magnitudes of Williams $`etal.`$ (1996).

### 2.3 Charge Transfer Efficiency Losses And Photometric Errors

After performing aperture photometry, nuclear variability can then be determined by simply obtaining the difference in magnitude for each source over the two epochs. Figure 1 shows the magnitude difference between the 1995 and 2000 images as a function of magnitude within an aperture of r$`=`$3 pixels, the equivalent of 0.15$`^{\prime \prime }`$ in diameter. The solid line represents no difference in magnitude between the two epochs.

The main culprit producing the obvious offset from the solid line is charge transfer efficiency (CTE) losses. As in SGP, we find a clear correlation of the photometric differences with position on the CCD and magnitude. This is the well known CTE effect (Whitmore $`etal.`$ 1999; Biretta $`etal.`$ 2001) which causes targets far from the CCD readout amplifier to appear fainter than similar targets near the amplifier. For this reason, the effect is usually more significant in the Y-direction, (the direction of the readout down the CCD) than the X-direction (the direction of the readout of the shift register at the bottom of the CCD). Although the 1995–1997 I-band study did not show significant correlations in the X-direction, the longer time interval sampled here does show a significant correlation in the X-direction for all 3 WF CCDs. Additionally, we note that the slope of the relation is generally steeper in the Y-direction than had been determined in SGP. This is expected due to the longer 5 year time interval between epochs in this study.

We first attempted to fit the relationship between object position, magnitude difference and magnitude globally for all 3 WF CCDs as done in SGP. However, we found that this approach was not sufficient to equally remove the CTE effect from all WF CCDs. Differences in the slope of the relation between the CCDs were significant and could easily be determined for each CCD individually. We adopted a two-step approach to determine this relationship. First, we fit a linear surface to the X and Y object positions and the magnitude differences using SURFIT in IRAF for each CCD. While the slope of the fit was found to be significant with both the X and Y positions, the value of the slope was steeper with Y position in all 3 WF CCDs, consistent with what is expected for CTE losses. The slope of the relation between Y position and magnitude difference was found to be -2.6$`\times `$10-5, -1.3$`\times `$10-5 and -1.46$`\times `$10-5 mag/pixel for WF2, WF3 and WF4 respectively. The slope of the relation with X position, while significant, was found to be about an order of magnitude lower in all 3 WF CCDs. These fits were used to correct the magnitude difference values for position-dependent offsets. Secondly we fit the dependence of the magnitude difference with nuclear magnitude. A 3rd order fit was required to properly model the dependence in all 3 WF CCDs. After applying this correction to the data in Figure 1, the resulting photometric differences are shown in Figure 2 as a function of nuclear magnitude.

Finally, we determine the expected noise level for non-variable sources. As in SGP, we divided the full set of individual exposures into two sets composed simply of the odd-numbered frames and the even-numbered frames. Each set represents an “epoch” without any real time difference. These images effectively carry through the effects of object and sky Poisson noise, readout noise, and possible errors in the adopted sky zero points.

The magnitude of each galaxy was then measured in the odd and even data sets and difference in magnitude vs. the average magnitude within r$`=`$3 pixel apertures is shown in Figure 3a. Figure 3b is the RMS of the intensity differences shown in a) within unity magnitude intervals. The solid line represents a quadratic fit to the points which is the adopted 1$`\sigma `$ galaxy photometric error as a function of magnitude.

The precision of the photometry must be independently limited at the bright end of the distribution, since systematic effects not captured by the odd-even data set comparisons become limiting factors (e.g. the sometimes poor distribution of sub-pixel phase-space dithers in the first epoch data). In SGP, we adopted a floor of 1.2%. Here we adopt a more conservative value of 1.5% consistent with the use of much smaller apertures on average for bright sources in the current study than in SGP. As should be apparent, this is a reasonable limit for photometric precision to adopt, but in detail remains somewhat arbitrary.

## 3 SELECTION OF VARIABLES AND SURVEY COMPLETENESS

Figure 4a is the absolute value, CTE-corrected magnitude differences for sources in the HDF. The solid line represents the 3$`\sigma `$ limit for variability significance (3$`\times `$ the RMS indicated by the solid line in Figure 3b). Objects above this limit are selected as significant variables and are indicated with open hexagons. Figure 4b is the result of normalizing the magnitude difference by this solid line. Here the Y-axis indicates the level of significance of each source in units of $`\sigma `$. The X-axis extends at the faint end to Vnuc$`=`$29.0, which is the estimated photometric completeness magnitude limit for galaxy nuclei in this survey. Beyond this limit, the number counts for galaxies in the HDF begin to decrease.

We find sixteen galaxies whose nuclei have undergone at least a 3$`\sigma `$ variation over the 5 year time interval. These sources are listed in Table 1 with columns as follows: (1) Williams $`etal.`$ (1996) ID, (2) & (3) RA & DEC (J2000), (4) Redshift from the literature, (5) Spectral type based on photometry from Fernández-Soto, Lanzetta and Yahil (1999, hereafter: FLY), (6) Vnuc internal to r=3 pixel aperture, (7) Magnitude difference between 1995 and 2000; positive implies brighter in 1995, (8) Significance of change obtained when normalized by the expected error as a function of magnitude, (8) Bulge-to-Total (B/T) 2-dimensional model fits for the I(F814W) images from Marleau & Simard (1998). The last several columns relate to the source detection at other wavelengths discussed in Section 5.

There are two distinct completeness issues that effect our survey. The first is related to the incomplete time-sampling of the variable sources we wish to detect. Most variability surveys employ a method where the survey field is imaged several times over many years (e.g. Trevese $`etal.`$ 1994; Hawkins 2002). Depending on the quality of the data, these surveys have shown that virtually all known AGN will be found to vary if observed periodically over several years. Our study is limited by the fact that we have only two epochs with which to determine variability and therefore sample just two points on the lightcurve of a varying source. Because of this, we will be incomplete in our census of AGN since some varying sources could lie at magnitudes close to their original magnitude measured 5 years earlier and would thus go undetected in our survey.

We estimate our incompleteness due to undersampling of the lightcurve by using variability data for AGN obtained over many years. We randomly select points along the lightcurve separated by 5 years, the time interval sampled by our HDF images. Because long term variability surveys for low-luminosity AGN have not yet been published, we have conducted this test with two sets of QSO lightcurve data; a sample of PG quasars from Giveon $`etal.`$ (1999) and a sample of SA57 quasars from Trevese $`etal.`$ (1994). Giveon $`etal.`$ monitored 42 quasars over 7 years with a typical sampling interval of 40 days while Trevese $`etal.`$ monitored 64 over 15 years about once a year. Sampling these lightcurves every 5 years and assuming photometric errors typical of our HDF images, we can estimate the probability that these sources would be detected in our survey. We find that $``$75–80% would be detected as varying by at least 3$`\sigma `$ with only 2 observations separated by 5 years. Therefore, the results of our analysis are likely to underrepresent the true number of variable nuclei by $``$20–25%.

The second incompleteness issue results from the use of a fixed aperture to detect nuclear variability. We have chosen a small, fixed aperture to minimize dilution of the AGN light from the underlying host galaxy to be more sensitive to AGN varying within bright hosts. We set this aperture to the size of an unresolved point source in the HST images (0.15$`^{\prime \prime }`$ in diameter). The fixed aperture will include the same amount of light from an unresolved source regardless of its redshift. However, the aperture will contain a larger fraction of the underlying, resolved galaxy for higher redshift objects than for low-z sources. Because of this, the dilution of the nuclear light by the galaxy increases as a function of redshift and, consequently, as a function of nuclear magnitude. If we assume a disk-dominated, re=1 kpc galaxy containing an AGN that is 10% of the host galaxy light and varying by $`\mathrm{\Delta }`$m=0.3 magnitudes, the measured $`\mathrm{\Delta }`$m within our fixed aperture would be 0.15 magnitudes at z=0.2 but only 0.06 magnitudes at z=1. The observed total nuclear magnitude would drop by $``$3 magnitudes over this interval. The increasing dilution with redshift is compounded by the increasing photometric errors with magnitude. For an AGN that has varied by $`\mathrm{\Delta }`$m=0.3–0.2 mags, consistent with typical structure function values for AGN (Trevese $`etal.`$ 1994), and reasonable estimates for the galaxy magnitude, morphology, physical size and true AGN fraction, we estimate that incompleteness due to significant AGN dilution begins to effect our survey at magnitudes between Vnuc$``$26.5 to Vnuc$``$27.5.

As can be seen in Figure 4, all of the variables lie at magnitudes greater than Vnuc=27.3. While this may partially be due to a real astrophysical effect (see discussion in Section 4), this limiting magnitude is within the range of our incompleteness estimate above. Therefore, although we survey galaxy nuclei to Vnuc$``$29, our limit for detecting variable AGN at a uniformly varying level within galaxy nuclei is probably closer to Vnuc$``$27.5, with decreasing completeness fainter than Vnuc$``$26.5.

The 16 variable nuclei detected in our survey represent 2.2% of the galaxies brighter than Vnuc$`=`$29.0 or 7.4% of the galaxies brighter than Vnuc$`=`$27.5. To emphasize the significance of the variable sources, Figure 5 is a histogram of the sigma distribution for galaxies in the HDF. The X-axis is the absolute value of the normalized magnitude difference ($`\sigma `$) and the Y-axis is the natural logarithm of the number of sources in 0.25$`\sigma `$ bins. Errorbars represent the poisson statistical errors. Filled circles show the histogram for all 719 nuclei brighter than Vnuc$`=`$29.0 and open squares represent the 217 nuclei brighter than Vnuc=27.5. The curved lines are gaussian fits to the data within 2.5$`\sigma `$ (solid is the fit to data brighter than Vnuc$`=`$29.0 and dashed is the fit to data brighter than 27.5). In both distributions, the data are well fit by gaussians out to $``$3$`\sigma `$ and show a “tail” of significant variables extending to higher $`\sigma `$ values. In a normal distribution, we would expect a total of $``$5 sources to be greater than 3$`\sigma `$ brighter than Vnuc$`=`$29.0 with $``$1.5 expected to lie brighter than Vnuc$`=`$27.5.

Based on Gaussian statistics, we estimate that $``$1–2 of our variables (all detected brighter than Vnuc$`=`$27.5) are spurious, implying that the true number of variables is $``$14. If we also correct by 20–25% for incompleteness due to lightcurve sampling, the incompleteness corrected number of variables in the HDF is $``$18 representing 8.3% of the nuclei brighter than Vnuc$`=`$27.5.

One of the 16 variables detected in this survey, 2-251.0, was also identified as variable in the I-band variability survey of the HDF (SGP). Eight variables were detected in that survey with significant ($`>`$3$`\sigma `$) magnitude changes over the 2 year time interval (1995–1997). Object 2-251.0 is by far the brightest of these sources, having a nuclear V magnitude $``$13$`\times `$ brighter than the next brightest source. Of the 7 I-band variables not detected in the present survey, 5 are fainter than Vnuc=28.0, placing them below the expected sensitivity limits for this study. The other two I-band variables, 3-266.0 and 3-404.0, have Vnuc$``$27.0. In the present study, object 3-266.0 changed by 0.025 mag (1.1$`\sigma `$ significance) and object 3-404.0 changed by 0.052 mag (2.0$`\sigma `$ significance). The latter source is a marginally significant detection and would lie above the 3$`\sigma `$ limit if it had been brighter than $``$26.5 in the V-band. Object 3-266.0, at just over 1$`\sigma `$ significance, may be variable but at lower significance due to the 2-epoch statistical incompleteness effect discussed above.

## 4 PROPERTIES OF THE SELECTED VARIABLE NUCLEI AND HOST GALAXIES

Figure 6 contains the 1995 V-band images for galaxies with variable nuclei. The postage stamp images are 3$`^{\prime \prime }`$ on a side and are scaled to the same maximum pixel value. It can be seen from this figure and Table 1 that the selected galaxies cover a range of morphologies, magnitudes and colors. Their redshifts range from 0.09 to 1.8 with a median $`<`$z$`>`$$``$0.66. Many of the galaxies have significant bulge components based to the 2-dimensional model fits. We note, however, that several sources with B/T$`=`$0 also have high chi-square values ($`>`$1.25) which in some cases appears to be due to a poorly fit bulge component. In addition to the morphological B/T parameter, we also list the spectral type for the source based upon the photometry of FLY to determine photometric redshifts for galaxies in the HDF. The selected galaxies appear to be evenly distributed over all spectral classes with no particular preference towards early or late-type galaxies.

Figure 7 is a color magnitude diagram of the galaxy nuclei for sources in the HDF. The variable nuclei are indicated with larger solid triangles. The variables cover the full range of V–I colors and do not appear to be preferentially bluer or redder than the non-varying nuclei of galaxies in the HDF.

All of the variables lie at the bright end of the magnitude distribution. As previously discussed, incompleteness due to increasing AGN light dilution as a function of redshift is likely to contribute to this observation. However, an astrophysical effect may be at least partly responsible. AGN are generally found in brighter host galaxies than those that do not harbor AGN (Huchra & Burg 1992; Ho $`etal.`$ 1997; Hamilton $`etal.`$ 2002). If these sources are indeed AGN, we might expect to find them in the brighter galaxies. An additional factor is the inclusion of the AGN light in the nuclear magnitude. We can minimize the above mentioned incompleteness effect by considering only those sources brighter than Vnuc=26.5, where incompleteness is not expected to be significant. We find that 23% (5/22) of nuclei are variable at 21.8$`<`$Vnuc$`<`$25.5 while only 15% (6/39) are variable at 25.5$`<`$Vnuc$`<`$26.5. While the small numbers do not provide significance, this trend is consistent with the previous studies concerning AGN host galaxies.

## 5 COMPARISON WITH MULTI-WAVELENGTH HDF SURVEYS

The HDF has been surveyed with a variety of other telescopes and instruments providing a wonderful resource of multi-wavelength information for this particular region of the sky. We have cross-referenced our list of variable nuclei with published catalogs of X-ray, radio and far-infrared sources as well as other photometric and spectroscopic surveys used to search for nuclear activity.

AGN have long been recognized as X-ray sources, being associated with a wide range of activity levels from the brightest quasars to low-luminosity Seyferts. Even highly obscured AGN can often be detected at X-ray wavelengths. Some of the light obscured by dust at other wavelengths will be reprocessed and emitted in the mid-IR. For these reasons, the X-ray and mid-IR regimes are excellent wavelengths in which to detect and study AGN.

Recent deep surveys with the Chandra X-ray Observatory have resolved many X-ray sources in the HDF and provide information about the nature of the X-ray radiation. The 2Ms Chandra X-ray survey (Alexander $`etal.`$ 2003) has detected 20 sources in the original HDF and an additional 2 sources of lower significance which are listed in their supplementary catalog. Of these 22 detections, 18 are included in our variability study. Two were not included because their optical fluxes were too faint to be included in our survey and an additional pair fell too close to the edge of the CCD to obtain accurate nuclear photometry in both epochs. Of these 18 sources, we find seven that match (within 1.2$`^{\prime \prime }`$) the positions of our variable nuclei (see columns 10 and 11 of Table 1). The filled symbols in Figure 4 represent the 18 X-ray sources included in our survey and show the level of variability significance for each source.

In addition, we have found several variables coincident with source positions from the ISOCAM 15$`\mu `$m survey of the HDF (Aussel $`etal.`$ 1999) and the 1.4GHz radio survey (Richards $`etal.`$ 1999). These are listed in columns 12 through 14 of Table 1. Below, we discuss the variable galaxy nuclei that overlap with one or more of these multi-wavelength surveys.

2-251.0: This z=0.960 early-type spiral galaxy is one of only two broad-line AGN in the HDF (see Section 6) and is a bright X-ray source detected with Chandra in the soft, hard and ultrahard bands (Hornscheiemer $`etal.`$ 2000; Brandt $`etal.`$ 2001a – hereafter H00, B01a) with a photon index of 0.67. We detect it as an optical variable at 6.5$`\sigma `$ significance. As previously mentioned, it is the only object that was also detected in our I-band survey of the HDF (SGP). Spectroscopic detection of broad MgII (Phillips $`etal.`$ 1997), as well as 1.4GHz radio (Richards $`etal.`$ 1999a) and ISOCAM (Aussel $`etal.`$ 1999) detections for this source all corroborate the AGN nature of this galaxy. Spectral fitting of the X-ray data has revealed a large intrinsic column density which appears to be related to the AGN (B01a).

3-355.0: This elliptical galaxy at z=0.474 is a 3.2$`\sigma `$ variable. The original 166ks Chandra exposure did not detect this source in the hard band, but B01a report hard X-ray counts from the deeper 479.7ks exposure. The photon index determined from the 2Ms exposure is 1.8. This galaxy is also a radio source with a steep radio spectrum ($`\alpha `$=1.0). Its position is coincident with an ISO source from the Aussel $`etal.`$ (1999) supplementary table of lower significance detections. H00 note that the X-ray luminosity for this source is consistent with that expected from hot gas in an elliptical but the nuclear variability detected here suggests this galaxy may also host an AGN.

3-659.1: The variable source is a spiral galaxy at z=0.401 (Barger $`etal.`$ 2002) with a variability significance of 4$`\sigma `$. This optical galaxy is $``$1.2$`^{\prime \prime }`$ from the hardest X-ray source ($`\mathrm{\Gamma }`$=0.56) in the HDF. The X-ray position is only 0.22$`^{\prime \prime }`$ from the second reddest source in the HDF NICMOS survey (Dickinson $`etal.`$ 2000). Even though the optical galaxy is positioned further from the center of the X-ray emission, the detected variability suggests that perhaps some of the X-ray emission originates from the galaxy nucleus. There is also significant 1.4GHz radio emission (Richards $`etal.`$ 1999a) and 15$`\mu `$m emission detected with ISO (Aussel $`etal.`$ 1999). The radio position is again more closely aligned with the NICMOS source than the optical galaxy (see Figure 1 in Richards 1999b) but both the NICMOS and optical sources are within the positional uncertainties.

3-965.111111: This is the most significant variable in our survey (9.4$`\sigma `$). It was not detected in either the radio or mid-IR bands and was first detected in the 1Ms Chandra survey (Brandt $`etal.`$ 2001b). It was detected in both the soft and hard X-ray bands in the 2Ms survey and has a photon index of 0.91. It is a bright elliptical galaxy near the edge of the HDF at z=0.663.

4-752.1: This optically variable (3.1$`\sigma `$) X-ray source is also an FR I radio galaxy (Richards 2000) with extensive radio structure. The host galaxy is a red elliptical at z=1.05. This source has not been detected in the hard X-ray band and thus appears to fall into the low X-ray luminosity/soft X-ray spectrum grouping of very red objects (Hogg $`etal.`$ 2000). Because of its faintness, the photon index cannot be determined. B01a suggests that the X-ray emission is associated with the central AGN or the hot interstellar medium and X-ray binaries of the elliptical galaxy. It is not detected in the mid-infrared.

4-976.1: This bright spiral galaxy at z=0.089 is variable at the 3.6$`\sigma `$ significance level. B01a notes that the position of the X-ray source is not coincident with the galaxy nucleus (which is 1.1$`^{\prime \prime }`$ away) but does lie near ($``$0.14$`^{\prime \prime }`$) an off-nuclear bright spot. They suggest that this bright spot may be a background AGN, starbursting region, or “super Eddington” X-ray binary (B01a, H00). Our variability measurements have been made on the central nuclear region of this galaxy which, if hosting an AGN, may be responsible for at least part of the X-ray emission. The photon index for this source is 1.7. This galaxy is within 3.3$`^{\prime \prime }`$ of a mid-IR source (Aussel $`etal.`$ 1999).

3-386.111: This is a spiral galaxy at z=0.474 which is listed in the 2Ms X-ray survey supplementary catalog containing lower significance sources associated with optically bright objects. The nucleus has a variability significance of 3.8$`\sigma `$. This source was not detected in the hard X-ray band. It is associated with an ISO source listed in the supplementary ISO source catalog.

One additional variable nucleus, 4-254.0, matched the position of a lower significance X-ray source indicated in the 1.38Ms Chandra survey (Brandt $`etal.`$ 2002). This X-ray source, however, was not listed among the 2Ms detections. The object is an elliptical galaxy at z=0.901 which is also associated with an ISO detection in the supplementary catalog.

Of the remaining 8 variable nuclei, one (2-860.0) is associated with a mid-IR detection and another (3-943.0) is within $``$2$`^{\prime \prime }`$ of a 1.4GHz radio source. The other 6 are not associated with any of the multi-wavelength surveys discussed here.

Overall, 44% (7/16) of our variable galaxy nuclei have X-ray counterparts from the 2Ms Chandra survey (or 31% excluding the two matches with greatest positional offsets). Optical variables make up 39% of the X-ray sources. Of the brightest X-ray sources (full-band flux $``$2.4$`\times `$10-16 ergs/s) five out of nine (56%) are associated with optical variables. If we assume that the observed optical variability indicates the presence of an AGN, this is an important check on the nature of the X-ray emitting source and could potentially help discriminate between low-luminosity X-ray emitting AGN and other X-ray emitting phenomena such as supernova remnants or X-ray binaries. There does not appear to be a relation between photon index and the detection of optical variability. Five of the eight X-ray sources in the HDF with enough signal to determine photon indices are optical variables. The X-ray/optical variables make up two of the softer X-ray sources and three of the harder sources in this small sample.

Seven of our 16 variables are associated with mid-IR detections (44%) and five (31%) have radio emission at 1.4GHz. Variable galaxies make up 22% of the detected sources at 15$`\mu `$m. This is comparable to the portion of the mid-IR integrated light attributed to AGN based on the correlation of mid-IR sources with Chandra sources (Elbaz $`etal.`$ 1999).

Finally, we have compared our results with two photometric studies to identify QSOs in the HDF on the basis of multi-band colors. Jarvis & MacAlpine (1998) identified 12 high redshift (z$`>`$3.5) QSO candidates. Eleven of these sources were included in our variability survey but none showed significant variability amplitudes, having a median $`\sigma `$ of $``$0.6. Based on structure functions for AGN/QSOs (Trevese $`etal.`$ 1994; Hook $`etal.`$ 1994; Hawkins $`etal.`$ 2002) and taking into account time-dilation effects ($`\mathrm{\Delta }`$trest=$`\mathrm{\Delta }`$t/(1+z)), the average magnitude change expected in our survey for high-z QSOs is $``$0.1–0.2 magnitudes. Even sources dominated by an AGN component, as expected for these color-selected candidates, would be difficult to detect above the 3$`\sigma `$ significance threshold due to their faint apparent magnitudes (Vnuc$`>`$27.7). Nonetheless, an average magnitude change of only $``$0.05 is significantly less than the 0.1–0.2 average magnitude change predicted from QSO structure functions.

Conti $`etal.`$ (1999) also identify 20 compact sources having QSO-like colors and morphologies with estimated redshifts between z$``$1 and 5.5. The nuclear magnitudes for these objects lie at Vnuc$``$27–28. None were identified as variable above the 3$`\sigma `$ threshold in our survey. The average nuclear magnitude of this sample is $``$1 magnitude brighter than the Jarvis & MacAlpine QSOs. If these are variable QSOs/AGN-dominated galaxies, we would again expect an average magnitude change of $``$0.1–0.2. The average change in magnitude for the Conti $`etal.`$ sample over the 5 year interval is $``$0.03 magnitudes, with all sources lying well below the 2$`\sigma `$ significance limit.

## 6 OPTICAL SPECTRA

The extensive redshift surveys of the HDF from Cohen $`etal.`$ (1996; 2000), Phillips $`etal.`$ (1997) and Barger $`etal.`$ (2002) have revealed only two broad-line AGNs (BLAGNs) in the HDF proper. These are 2-251.0 (z=0.96 variable galaxy and X-ray source discussed above) and 4-852.12 (z=0.943 X-ray source). The measured variability significance of 4-852.12 is 0.35$`\sigma `$. It is possible that this source is variable but has been observed at two points in its light curve that are close to the same magnitude. Further monitoring of the HDF would be necessary to rule out optical variability for this BLAGN.

We have studied the available spectra for all of the galaxies hosting variable nuclei to determine if any of the sources show specific emission lines or line flux ratios indicative of Type 2 AGN. Of the 16 variables, optical spectra exist for 13. One of these is the BLAGN 2-251.0 already discussed, which shows broad MgII ($`\lambda `$2800) emission and absorption in its spectrum. Almost all of the remaining 12 show emission lines in their spectra. Nine out of 11 show OII($`\lambda `$3727) when in range, with several sources also displaying H$`\beta `$ and OIII($`\lambda `$5007). The X-ray and strong radio FRII galaxy 4-752.1, displays only strong absorption lines in its spectrum. Only one source, 4-254.0, shows weak NeV($`\lambda `$3426) emission, a line indicative of the presence of an AGN (Hall $`etal.`$ 2000). NeIII($`\lambda `$3869), also stronger in AGN than in starforming galaxies (Rola $`etal.`$ 1997), is seen weakly in 2 sources, 3-143.0 and 3-386.111111. These galaxies also show OIII and H$`\beta `$, but with flux ratios consistent with star formation rather than AGN activity (Veilleux & Osterbrock 1997). Only the spiral galaxy 4-976.1 is low enough redshift to reveal strong H$`\alpha `$, SII($`\lambda `$6713+6731) and OI($`\lambda `$6300) emission in the optical spectrum. The line flux ratios of OIII/H$`\beta `$, SII/H$`\alpha `$ and OI/H$`\alpha `$ all indicate that this source is near the border that divides starforming galaxies and AGN (Veilleux & Osterbrock 1997). Therefore, based on the optical spectra alone, objects 2-251.0 (BLAGN), possibly 4-976.1 (LINER/Seyfert 2), and possibly 4-254.0 (through the weak presence of NeV) show evidence of AGN.

## 7 LUMINOSITY FUNCTION FOR VARIABLE NUCLEI

Knowledge of the space density of faint AGN at higher redshifts is critical to our understanding of AGN evolution and the AGN phenomenon in general. Our sample of 16 variable sources spans a redshift range of z=0.1 to 1.8 with most sources lying between z=0.4 and 1.1. If we make the assumption that these variables are indeed AGN, it is interesting to compare their number density with local samples of AGN.

We must first estimate the luminosity of the AGN candidates in our survey. Due to the faint magnitudes and small angular sizes of many of the variables, 2-dimensional modeling to separate the components (disk+bulge+AGN) of the galaxies is difficult and would not yield consistent results for all of the selected galaxies. A simpler approach is to estimate the AGN luminosities by assuming that the entire nuclear flux is attributed to the AGN. This value is an overestimate in one sense since the nuclear light will also include light from the underlying galaxy, though the small nuclear aperture minimizes the contribution from the host. Based on the expected level of variability over this time interval predicted from AGN structure functions, we estimate the variable component for our AGN candidates to be no less than 15% of the total nuclear light. In the most extreme case, the true AGN magnitude could be $``$2 magnitudes fainter than that measured within the nuclear aperture. In another sense however, we underestimate the AGN flux since we have not applied aperture corrections to include light from the wings of the PSF. Our aperture of r$`=`$0.075$`^{\prime \prime }`$ encircles about 50% of the flux from a point source (Holtzman $`etal.`$ 1995). If the point source makes up the entire flux within the aperture, an aperture correction would increase the magnitude by up to $``$0.8 mag. If the point source is some fraction of the total nuclear light, the aperture correction will be smaller. Determining the true aperture correction is not possible without knowing the relative fluxes of the point source and host within the fixed aperture.

Figure 8 is the absolute magnitude based on the nuclear flux vs the source redshift (Ho=75 km/s/Mpc; qo=0.5). The absolute B magnitudes are calculated assuming a power-law index of $`\alpha `$=–1.0 resulting in a K-correction of zero and VF606W–BJ=–0.32. We find that most sources lie between -15$`>`$MB$`>`$-17 in the redshift range 0.4$`<`$z$`<`$1.1.

We conduct the Luminosity-Volume test, or V/Vmax test, of Schmidt (1968) to determine if significant evolution is present in our sample of variable nuclei. For this calculation, we determine for each source in our survey the maximum redshift at which it could exist and still be included in our survey. The volume of space enclosed by a sphere with a radius equal to this zmax is the value Vmax. The volume enclosed at the true redshift of the source is $`V`$. Spectroscopic redshifts have been determined for 14 of our objects. Photometric redshifts are estimated for one other (FLY). There is no spectroscopic or photometric redshift published for 2-456.22. Based on the V–I color for this object, we estimate its photometric redshift to be z=1.8 which corresponds to those of other HDF sources with similar colors.

To define zmax, we determine the limiting redshift for each source to be detected in our survey. The brightest limiting magnitude in our survey is defined by the 3$`\sigma `$ line in Figure 4a, which indicates the faintest magnitude at which a source varying by a given $`\mathrm{\Delta }`$m would be considered a significant variable. Below this magnitude, a source varying by $`\mathrm{\Delta }`$m would not be detected as variable. Therefore, zmax is the redshift where the apparent magnitude of a nucleus varying by $`\mathrm{\Delta }`$m falls below the 3$`\sigma `$ variability significance limit shown in Figure 4a.

The mean V/Vmax for the 16 nuclei in our survey is 0.67$`\pm `$0.07. For a population of objects uniformly distributed in space, $`<`$V/Vmax$`>`$=0.5. Our value of 0.67$`\pm `$0.07 indicates that sources are found more often at larger distances than at nearby distances, indicating some evolution in the population of galaxies with variable nuclei. Our survey incompleteness, as discussed in Section 3, affects our ability to detect sources at higher redshifts and fainter magnitudes. Thus, the value of $`<`$V/Vmax$`>`$ would actually be higher if corrected for this incompleteness. An increase in the number density would be consistent with the picture of AGN/QSO evolution at higher luminosities, where the population increases rapidly out to z$``$3.

Even though evolution may be present among our sample, the small number of sources spread over a range in redshift makes it difficult to construct the luminosity function in more than one redshift bin. To accurately compare the density of our sources to local samples, we construct the LF in one redshift bin at z=0.4–1.1. We choose these limits to minimize the redshift range of our sources while maximizing the number of sources in the bin. Twelve of our 16 sources lie within this redshift range, which excludes the 2 lowest and highest redshift galaxies. The mean redshift for these sources is $`<`$z$`>`$=0.69.

We calculate the luminosity function using a technique similar in nature to the V/Vmax test described above. The LF is created by summing 1/Vmax for sources in discrete magnitude bins. In this case, we are actually determining Va, or the accessible volume in which each source could be detected. This volume is limited not only by the magnitude limit for detection, but also by the lower and upper redshift limits of the shell in which we are summing. Figure 9 is the LF for HDF variable nuclei in the redshift range z=0.4–1.1. Incompleteness due to the fact that our variability selection is based on only 2 epochs is expected to produce an underestimate in the number of AGN by a factor of $``$1.3, or a decrease in log$`\varphi `$ of $``$0.12. For comparison, we plot the local Seyfert LFs of Huchra & Burg (1992; HB92) from the CfA survey and Ulvestad & Ho (2001; UH01) from the Palomar survey. Both local surveys consist of spectroscopically selected Seyfert 1s and 2s.

Our LF for variable sources at 0.4$`<`$z$`<`$1.1 extends to fainter absolute magnitudes than the local LFs, largely due to the fact that the local LFs are calculated for the total galaxy light (AGN+host) rather than the magnitude of the AGN alone. Nonetheless, many of the AGN candidates detected in our survey may be truly fainter than those of HB92. Based on the spectroscopic selection criterion employed by HB92, most of our sources would not be detected as AGN (see Section 6) in the HB92 survey. Therefore, we assume that the AGN component in our galaxies do not comprise as much of the total galaxy light as the HB92 sources and are likely to be intrinsically fainter.

The UH01 LF covers the same integrated galaxy magnitude range as HB92 but shows an overall higher density for local Seyferts. Ulvestad & Ho explain that this is due to the fact that their LF includes intrinsically fainter seyfert nuclei than does the HB92 sample. For the Seyfert 1 galaxies, Ho & Peng (2001) decompose the nuclei from the host galaxies in both local samples and confirm this assertion, finding a median nuclear magnitude of MB=-14.6 for UH01 versus -17.4 for the HB92 Seyfert 1 galaxies. From their study, we also find that the nuclear magnitudes of the local seyferts in both samples can be extremely faint when compared to the total galaxy magnitude, in some cases up to 10 magnitudes fainter. Therefore, an LF consisting of nuclear magnitudes for local Seyferts would be a much better comparison to our HDF AGN candidates.

Since this is not currently available, an alternative is to produce the LF of our sources using the total integrated galaxy magnitude rather than the nuclear magnitude alone (dashed LF). The integrated galaxy magnitudes for our variable sources cover the same range as the faint end of the local Seyfert LFs but have a number density 6–10$`\times `$ greater than the local Seyfert density of UH01. If the AGN in these two samples cover a similar magnitude range, this is evidence for a significant increase in number density from z=0 to $``$0.7. The magnitudes for the UH01 AGN determined in Ho & Peng (2001) for the Seyfert 1 population cover the range -9$`>`$MB$`>`$-22 mag with most lying within -12$`>`$MB$`>`$-20. Our nuclear aperture LF (open circles in Figure 9) extends from -14.5$`>`$MB$`>`$-19. In the most extreme case, based on the over and underestimates of the true AGN magnitudes discussed previously, the variable AGN in our survey could cover the magnitude range -12.5$`>`$MB$`>`$-20. This magnitude range is consistent with that estimated for the UH01 sample. It is therefore unlikely that the increase in number density between the UH01 LF and our LF is due to the inclusion of instrinsically fainter AGN as appears to be the case for the difference between the UH01 and HB92 LFs.

Finally, we note that the LF for the variable nuclei continues to rise at faint magnitudes, indicating that number counts have not begun to turn over at magnitudes even as faint as MB$``$-15. The effects of incompleteness in our survey due to AGN light dilution would cause us to be less complete at the faint end and corrections in this sense would likely increase rather than decrease our faintest LF bin. Most local AGN LFs, such as those shown in Figure 9, show a flattening toward fainter magnitudes. The local LF of Londish $`etal.`$ (2002) for type 1 AGN reveals a much flatter slope at the faint end as compared to the bright end and can be well fit with a two power-law function. Our increasing numbers at the faint end may be an evolutionary effect or might be attributed to the inclusion of a greater fraction of obscured AGN. More data are needed as well as better classification into AGN types to improve the statistics and verify this trend.

## 8 SUMMARY AND CONCLUSIONS

We have investigated variability within the Hubble Deep Field North with the aim of identifying low-luminosity AGN in the nuclei of field galaxies. We compare photometry of the nuclei of galaxies in the initial V-band image obtained in 1995 and a second image obtained in 2000 and find evidence for significant variability in 16 galaxies. Correcting for possible spurious detections and incompleteness results in a total of $``$18 sources which make up 8.3% of galaxy nuclei to Vnuc=27.5. The galaxies hosting variable nuclei range in redshift from 0.09$`<`$z$`<`$1.8 with most lying between 0.4$`<`$z$`<`$1.1. Nuclear colors of the variable nuclei are consistent with those of non-variable nuclei.

Seven of our variable sources are coincident (within $``$1.2$`^{\prime \prime }`$) with X-ray sources detected in the 2Ms Chandra exposure (Alexander $`etal.`$ 2003). We find that 39% of the 2Ms Chandra detections in the original HDF are significant variables. The optical variability observed in these X-ray sources provides important evidence to confirm the presence of an AGN and could potentially help discriminate between low-luminosity X-ray emitting AGN and other X-ray emitting phenomena such as supernova remnants or X-ray binaries. In addition, 44% of our sources are associated with mid-IR detections at 15$`\mu `$m and 31% are detected at 1.4GHz.

Optical spectra are currently available for 13 of the 16 variable sources. Object 2-251.0 is a broad-line AGN. Two others show some evidence of harboring an AGN. Object 4-976.1, a low redshift spiral, has emission line ratios indicative of a LINER/Seyfert 2 and 4-254.0 reveals weak NeV emission.

Based on the V/Vmax test, there is evidence for evolution within our sample of variable nuclei in the sense that more variables are found at higher redshifts. We estimate the luminosity function of variable AGN in the redshift range 0.4$`<`$z$`<`$1.1. This LF is a lower-limit since our variability selection technique, employing only 2 epochs for detection, could miss $``$20–25% of variable nuclei. The nuclear magnitudes for our variables extend to MB$``$-15.0 with a median magnitude of MB=-16.3. Computing our LF with total integrated galaxy magnitudes, we compare to the LF of Ulvestad & Ho (2001) for local Seyfert galaxies. If our AGN cover a similar magnitude range as those included in their sample, we find evidence for an increase in the number density of faint AGN by a factor of 6–10 from z$``$0 to 0.7. The LF for our selected variable nuclei continues to rise at faint magnitudes. Unlike the LFs of most QSOs and local Seyferts, we do not see signs of flattening at the faint end.

The results of this paper can be tested and confirmed with a larger sample of variable nuclei detected at these redshifts. A similar variability survey is being conducted for the Groth-Westphal Survey Strip (Sarajedini $`etal.`$ 2003) over 7 years and will produce a larger sample of variables to improve the statistical significance of any evolutionary trends. Additional future work includes the analysis of the currently ongoing GOODS HST Treasury program (Dickinson $`etal.`$ 2003), an ACS survey of the sky around and including the HDF. The 5-epoch survey, at intervals of $``$45 days, will allow for the investigation of short-term variability as well as providing confirmation and improved completeness for the variables detected in this study.

Support for proposal GO-8389.02-97A for VS and GO-8389.01-97A and AR-7984.01-96A for RG was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. We gratefully thank Rodrigo Ibata who was the PI for obtaining the second epoch V-band image of the HDF. A special thanks to Amy Barger and Judy Cohen for supplying the Keck spectra in digital form to be analyzed for these sources. Thanks also to David Koo and the referee for helpful comments and suggestions for improving this paper.