text:
# X-Ray Spectral Variability of Extreme BL Lac AGN H1426+428

## 1 Introduction

Blazars, including flat spectrum radio quasars (FSRQs) and BL Lacs, represent the most extreme examples of active galactic nuclei (AGN) due to their rapid variability and the extreme energies to which their spectra extend. One of the most exciting developments in extragalactic astrophysics during the last decade has been the detection of more than 65 AGN in the energy range from 30 MeV - 30 GeV (Hartman et al., 1999) by the EGRET experiment on the Compton Gamma Ray Observatory (CGRO). The detection of some AGN (e.g. Mrk 421, Mrk 501, 1ES 2344+514, PKS 2155-304 and most recently H1426+428 and 1ES1959+650) above 300 GeV by ground-based air cherenkov telescopes (ACTs) has reinvigorated the study of these objects (Catanese and Weekes, 1999; Horan et al., 2002; Holder et al., 2003). The dominant radiation from blazars is widely believed to arise from relativistic jets that are viewed at small angles to their axes (Blandford and Koenigl, 1979). EGRET and Whipple observations of high fluxes above 1 GeV and 300 GeV (respectively), the short-term variability of the gamma-ray emission (Kniffen et al., 1993; Mattox et al., 1993; Quinn et al., 1996; Buckley et al., 1996), and constraints on gamma-ray opacity due to pair-production with low energy photons are all pieces of evidence leading to the conclusion that the gamma-ray emission is also beamed and that jets are involved.

The broadband radiation spectrum of blazars consists of two parts: a synchrotron spectrum that spans radio to optical-ultraviolet wavelengths (and to X-rays for high-peaked objects) and a high-energy part that can extend from X-rays to gamma rays. Observationally, the spectra appear to have two distinct bumps when plotted in a $`\nu `$$`f_\nu `$ representation (Fossati et al., 1998; Ghisellini et al., 1998). The emission in the first peak of the spectral energy distribution (SED) is generally considered to be synchrotron emission from relativistic electrons. The most widely invoked models which attempt to explain the higher energy emission of the second SED peak fall into the category of leptonic models. These leptonic models posit that the X-ray to gamma ray emission is produced by inverse Compton scattering of lower energy photons by relativistic electrons within a narrow jet. The low energy photon fields could arise from synchrotron continuum photons within the jet (e.g. Koenigl (1981)), or they could arise from ambient photons from the accretion disk which enter the jet directly (e.g. Dermer, Schlickeiser, & Mastichiadis (1992)) or after scattering or reprocessing (e.g. Sikora, Begelman, & Rees (1994)). In addition to these leptonic models, hadronic models also attempt to explain the second bump of the SED. These models involve proton-initiated cascades (e.g. Mannheim (1993)) and/or proton synchrotron radiation (Muecke and Protheroe, 2001; Aharonian, 2000), as well as synchrotron emission from secondary muons and pions(Muecke et al., 2003).

The gamma-ray and simultaneous X-ray observations of blazars have placed strong constraints on these models (von Montigny et al., 1995; Macomb et al., 1995; Buckley et al., 1996) by constraining the properties of the emission regions (size, Doppler boosting factor, magnetic field, electron cooling time, total injection energy, …) and by comparing flaring timescales at different wavelengths; but so far the observations cannot definitively reject either hadronic or leptonic models of the high energy emission. X-ray observations are particularly important for observations of high-peaked BL Lacs since the peak of the synchrotron component lies in the band covered by X-ray detectors (e.g. Beppo-SAX observations of Mrk 501 (Pian et al., 1998) and H1426+428 (Costamante et al., 2001)). The determination of the X-ray peak is highly constraining to the suite of models for the emission from these objects, as are the timescales of the X-ray variability.

Further constraints to emission mechanisms may be obtained by exploring the characteristics of spectral hardness-intensity diagrams (HIDs) in the X-ray energy band, with the purpose of searching for hysteresis patterns. The presence of hysteresis loops and the handedness of such loops can be used to constrain the cooling and acceleration/injection timescales of the relativistic electron population (Böttcher and Chiang, 2002; Li and Kusunose, 2000; Kirk, Rieger, and Mastichiadis, 1998). Spectral hysteresis has been observed for several high-peaked BL Lac objects, including Mrk 421 and PKS 2155-304 (Takahashi et al., 1996; Kataoka et al., 2000; Zhang, 2002). Another set of observations of Mrk 421 using XMM reported no significant hysteresis (Sembay et al., 2002), however this result appears to be inconclusive due to the lack of coverage throughout the entire cycle of flux rise and subsequent decay coupled with the low amplitude variation of the one complete rise and decay.

H1426+428 is an X-ray selected high-peaked BL Lac with a relatively strong X-ray flux. It was first characterized as a BL Lac, due to its lack of prominent optical emission lines, by Remillard et al. (1989). Based upon observations by Beppo-SAX up to 100 keV, H1426+428 is an example of an extreme high energy peaked BL Lac (HBL) with the peak of the synchrotron emission lying above 100 keV (Costamante et al., 2001) during that observing campaign. This fact led Costamante et al. (2001) to speculate that H1426+428 is a prime candidate for TeV gamma-ray emission; a speculation recently confirmed by the Whipple collaboration (Horan et al., 2002; Petry et al., 2002). Characterizing the X-ray properties of this extreme object is important, especially due to the presence of this TeV emission in spite of a large redshift (z=0.129). Previous X-ray observations presented by Sambruna et al. (1997) and Madejski et al. (1992) using ASCA, ROSAT, and BBXRT (at different epochs separated by as much as 15 years) have shown evidence for spectral variability from one observation to the next, as well as evidence for spectral curvature and a spectral absorbtion feature at $``$0.6 keV. The results presented here represent a fairly continuous and deep exposure obtained for the purpose of exploring the details of the X-ray variability on multiple timescales.

## 2 Observations

X-ray observations were carried out between 7 March and 15 June of 2002. These observations were part of a multiwavelength campaign (Falcone et al., 2003), including radio, optical, X-ray, and gamma-ray energy bands. This report will describe the X-ray observations (multiwavelength observations will be reported on in a separate publication). X-ray observations were performed by PCA and HEXTE on the RXTE satellite. The results reported here are based on PCA data. The HEXTE data, which can only be analyzed on significantly larger timescales due to sensitivity constraints, will be described in the forthcoming multiwavelength paper. Typical RXTE exposures were between 1-4 hours per night for 3 weeks bracketing the new moon of each month during the 4 month long campaign. This resulted in $``$450 ksec of RXTE observations.

### 2.1 Instrument

The PCA instrument is composed of five proportional counter units (PCUs) and has a total collection area of 6500 cm2 when all counters are operating. However, during the observations discussed in this paper, only 2-3 of the proportional counters were on during normal operation. For nearly all of this campaign, the proportional counters known as PCU0 and PCU2 were collecting data, while the other PCUs were dormant. For consistency during data reduction, a cut was placed on the data to select only the PCU0 and PCU2 data. Data from these PCUs were extracted independently of one another and they were combined by jointly modelling the spectra at the time of analysis. The best statistics were obtained by using the 1st Xenon layer of each PCU and excluding layers two and three. Cuts were also made to eliminate data that included effects due to passage through the South Atlantic Anomaly, as well as data for which the electron rate was anomalously high. An elevation cut ($`>`$10$`\mathrm{deg}`$) was also imposed.

The extraction and processing of the data was done using ftools 5.2. A response matrix was calculated, using the latest available software (PCARSP V8.0), individually for each PCU with the same configuration used for the data extraction process. The weak source background model provided by RXTE/GOF was used for background subtraction. This model appeared stable and reproduced the observed spectrum at high energies, where the spectrum is expected to be background dominated. It is worth noting that the background model worked well for PCU0, which no longer has an operating propane veto layer. The model background for PCU0 converged to the observations at high energies and the rate variations in PCU0 mimicked those of PCU2 which has an operating propane veto layer.

### 2.2 Analysis

After extraction of the data, XSPEC 11.2.0 was used to analyze and apply models to the extracted spectra. The nominal energy range of the PCA is 2 to $``$60 keV, but since the response of the instrument is poorly characterized in the lowest energy bins, only data above 2.9 keV is used. Due to the combination of falling source spectra for H1426+42.8 and sensitivity as energy increases, there is no significant measurement above 24 keV. Based upon these two bounds, all analysis is confined to the region between 2.9 and 24 keV, and all other energy bins are excluded. The 3.8-5.4 keV energy region is also excluded from the analysis since there is a known Xenon absorbtion feature within this energy region that is not modelled well by the instrument response. The spectra are then fit to data that are binned into 1$``$day bins, as described below, and the fits are evaluated using a $`\chi ^2`$ analysis and an inspection of the residuals. Count rates in particular energy regions are obtained prior to application of spectral models in order to evaluate variability in the lightcurve. Flux is derived based upon the parameters of the model that fits the spectral data.

## 3 Results

### 3.1 Flux and Flux Variability

A lightcurve showing the X-ray rate as a function of time in different X-ray energy bands is shown in figure 1. This data has been binned into one day time intervals beginning with the first day of the campaign, and all error bars are 1$`\sigma `$ statistical errors. The 2.9-10 keV and the 10-24 keV energy bands have been evaluated independently of one another and plotted, along with the entire 2.9-24 keV band. The 2.9-24 keV flux observed by PCA varied between about 1.3x10-11 to 6.6x10-11 erg cm-2 s-1, and the mean 2.9-24 keV flux was 3.6x10-11 erg cm-2 s-1.
Variability is evident in all of the three bands. Based on this analysis, which uses bins of 1-day duration, variability studies are limited to the timescale range bracketed by 1-day and 4-months. There appears to be long term variability as the flux changes throughout the entire campaign with timescales in excess of two weeks, and there also appears to be variability on timescales as short as 1-2 days. While the variability was significant, especially on long timescales, there were no particularly large flares relative to those observed for other TeV emitting BL Lacs, such as Mrk 421, Mrk 501, and 1ES 1959+650.

### 3.2 Spectra and Spectral Variability

For each set of data in a one day time interval, a model of the spectra was derived. It was found that a power law with galactic absorbtion, column density $`N_H=1.4\times 10^{20}`$ (Stark et al., 1992), led to an excellent fit based upon a $`\chi ^2/dof`$ analysis. (In fact, rather than hovering at $``$1, the $`\chi ^2/dof`$ is actually consistently less than one, due to the fact that the background rate error bars, which are added in quadrature to the source rate error bars, are conservatively large. This occurs since the error bars on the background are assigned based on the counts in the observation, rather than the counts used to construct the background model, which is a sum of many observations with far more numerous counts.) The power law spectra are of the form:

$$f=[C(\frac{E}{1keV})^\mathrm{\Gamma }][e^{N_H\sigma (E)}]$$
(1)
where the units of the normalization constant C are photons $`keV^1cm^2s^1`$, $`\sigma (E)`$ is the photo-electric absorbtion cross section given by Morrison and McCammon (1983), and $`\mathrm{\Gamma }`$ is the power law spectral index. Allowing the absorbtion column density to vary as a free parameter did not significantly effect the fit so it was fixed to the galactic value quoted above. An example of a fit using this data is shown in figure 2.

The power law spectral index derived using the above procedure is shown as a function of time in figure 3, along with $`\chi ^2/dof`$ for the particular fit. All error bars are 1$`\sigma `$. The error bars for the 2.9-24 keV flux were calculated by applying simple error propagation to the 1$`\sigma `$ error bounds produced by XSpec for the power law spectral index and normalization constant. It can be seen from this plot that there was significant spectral variability during the campaign. The spectral index varies over the range from 1.46$`\pm `$0.05 to 2.03$`\pm `$0.09. The X-ray flux in the 2.9-24 keV band has also been derived based upon the spectra that were fit to the data. This flux is shown in the top panel of figure 3.

### 3.3 Variability Correlations and Hysteresis

Over the time frame of the entire campaign, there is no tight correlation between the X-ray flux and the power law spectral index. The correlation coefficient between the 2.9-24 keV flux and the spectral index has been calculated using the 1-day timescale data throughout the entire four-month campaign, and its value is 0.51. By simple inspection of the plots in figure 3, one can see that the spectral variability does not have a consistent relationship with the 2.9-24 keV flux variability over long timescales. A plot of the spectral index vs. flux is shown in figure 4. While, on average throughout the campaign, the flux increases as the spectrum becomes harder, there is a great deal of scatter on the HID for the entire campaign. This plot includes flux variations averaged over long timescales so it is certainly not reflective of the source characteristics during a single flare event. This lack of strong correlation does not preclude a relationship during individual flare timescales.

Hardness-intensity diagrams have also been produced from the data during individual periods of flux variation in order to search for evidence of hysteresis effects, which could be used to evaluate the cooling and acceleration timescales. Based upon the behavior of the flux, several time periods have been evaluated independently, as defined in table 1 and shown in figures 5 and 6. The time periods evaluated in figure 5 are those that have a rise and a subsequent decay in flux observed by RXTE-PCA. For three of these time periods, a well-defined spectral hysteresis is observed with the trace following a clockwise pattern in the representation plotted. For the fourth (region B), there is also a consistent hardness-flux relationship observed, which tends to harden the spectrum with increasing flux, but the handedness of the HID can not be determined within the resolution of the measurements.

The time periods evaluated in figure 6 are those which do not have a well-defined rise and decay, but do exhibit variable flux. Region E exhibits a fast decay, but there were no observations during the rise of the flare. The HID exhibits a softening spectra with the decrease in flux. Region F is characterized by a long rise in flux over a time period of $`>`$2 weeks, but there are no observations during the decay of flux. This HID diagram exhibits an overall hardening of the spectrum as the flux increases, but there is some scatter in the data points. Hysteresis properties cannot be determined due to the lack of data during the presumed decay of the flux. Region G is similar to region F in that it is a long rise of flux over a period of $``$2 weeks, but it does turn over and begin to decay in the last two days. This HID exhibits spectral hardening with increased flux. An overall clockwise hysteresis appears to be evident for region G, however this is not well determined due to two factors; the complete long-timescale flux variation is not observed throughout the entire decay, and the short-timescale flux variations make it difficult to evaluate the characteristics of the longer timescale variability. The two smaller HID loops associated with regions C and D flaring can also be seen in this HID. The effect of relatively short term variability combined with long-term variability is evident in the HID diagram for region G.

## 4 Discussion and Conclusions

This long campaign was designed to continuously monitor H1426+428 at X-ray energies (and at other wavelengths that will be reported on in a forthcoming paper) in an effort to explore variability at multiple timescales. In spite of the fact that there were no large flares, relative to those frequently observed for other TeV emitting AGN, the campaign successfully observed variability with doubling timescales ranging from $``$1 day to $`>`$2 weeks. In several cases, the short timescale flaring was observed on top of the long timescale gradual variations. This type of variability complicates studies of emission mechanisms and should be considered in general, especially for studies that include non-contemporaneous multiwavelength observations.

Throughout the entire campaign, a simple power law with galactic absorbtion was able to fit the data well between 2.9 and 24 keV. Prior observations performed by Sambruna et al. (1997) between 1985 and 1994 with ASCA and BBXRT led to similar spectra at some times, but at other times there were some notable differences. At some times, the spectra from Sambruna et al. (1997) could be fit by a single power law, but at other times a broken power law that became softer at higher energies was inferred. For the part of the spectrum between the break energy and $``$10 keV, the spectra had a typical power law index of $``$2.3 above a break energy of $``$2 keV for the observations described by Sambruna et al. (1997). This is somewhat softer than any of the spectra derived during the present campaign. The observations reported by Sambruna et al. (1997) all show 2-10 keV fluxes that fall within the range of fluxes observed during this campaign, whereas the spectra of these historical observations are softer than the more recent observations reported in this paper and in Costamante et al. (2001). This is most easily interpreted as time variability, rather than interpreting it as a contradiction between the two results, when consideration is given to the BeppoSAX observations of Costamante et al. (2001) which also leads to a different spectral index and when consideration is given to the observed variability during this four month campaign. It is also important to note that this has involved a comparison of different instruments with different band passes. Another interesting feature of the BBXRT observations reported by Sambruna et al. (1997) was the evidence for a spectral line at $``$0.6 keV, which would imply the presence of absorbing material. No spectral lines were observed in the 2.9 to 24 keV band covered by the observations presented in this paper (PCA can not observe down to 0.6 keV where the spectral feature was previously observed).

The power law spectral index was observed to vary throughout the range bounded by 1.46$`\pm `$0.05 and 2.03$`\pm `$0.09. If one interprets this as being due to a shift in the location of the first peak of the spectral energy distribution, then it implies that the peak is sometimes in excess of 100 keV, while at other times it falls into the 2.9-24 keV region observed during this campaign. While past observations of H1426+428 have implied a high first peak in the spectral energy distribution (SED) (Costamante et al., 2001), there have also been observations from which a spectrum with $`\mathrm{\Gamma }`$$`>`$2 was derived (Sambruna et al., 1997). Based on this campaign it can be seen that significant spectral variability is present on multiple timescales. The average spectral index was generally harder than a flat $`\nu `$fν spectrum, which implies that the first peak of the average SED was probably in excess of $``$100 keV (PCA alone cannot derive the actual peak location for these cases, but past observations (Costamante et al., 2001) with BeppoSAX have shown the peak to be in excess of $``$100 keV when the spectral index is $``$1.9 at PCA energies). At other times the spectral index was flat ($`\mathrm{\Gamma }`$2), which implies that the first peak of the SED was in the observed 2.9-24 keV region. It is interesting to note that this movement in the first peak of the SED during this campaign is not accompanied by any measurable fluctuations in the TeV emission from this known TeV source (Falcone et al., 2003), but the probability of detecting such associated TeV fluctuations from such a weak source has yet to be evaluated. The overall spectral energy distribution and its implications will be discussed in more detail in a forthcoming article that will incorporate radio, optical, X-ray, and VHE gamma-ray observations carried out contemporaneously during this campaign.

The longest flux variability timescale observed was the long rise and short decay in June (region G). In spite of the fact that only a few data points were available for the decay phase, some hysteresis was evident. The clockwise loop approached closure, as shown in the third panel of figure 6. Although it is only marginally significant, it is interesting to note the small loop in the middle of the HID diagram, which coincides with the shorter timescale flux variability that occurred in the midst of the long rise in June. This behavior illustrates the need for sufficient sampling over complete flaring timescales. It is probably this effect of hysteresis within hysteresis due to multiple flares, as well as flares that are non-simultaneous, that causes the scatter on the overall HID plotted in figure 4. In spite of this scatter, it has been shown that there is a tighter relationship between flux and spectral hardness when each flux variability episode is studied independently on its proper timescale.

Spectral hysteresis has been observed during several periods of flux variability during this campaign. For an unambiguous observation of spectral hysteresis, it is necessary to observe during a large fraction of both the decay and the rise of a flare. For three of the four cases in which this condition was met during this campaign, spectral hysteresis loops tracing a clockwise pattern were observed (given the axis orientation shown in figure 5). In the fourth case (region B), no hysteresis and associated orientation could be observed, but the spectrum did harden as the flux increased during the rise and decay of flux. In the past, two other high-peaked BL Lac AGN that are known TeV emitters (Mrk 421 and PKS 2155-304) have exhibited similar clockwise loops in their hysteresis plots (Takahashi et al., 1996; Kataoka et al., 2000). As mentioned previously, Sembay et al. (2002) did not observe hysteresis with XMM observations of Mrk 421, however these observations did not sample the rise and fall of the flux variability so it is reasonable to expect a null result in that case. The one flare that was completely observed by Sembay et al. (2002) did display the characteristic spectral softening with decreasing flux, but no loop could be discerned from this low amplitude variation. Another observation of Mrk 421 using BeppoSAX (Zhang, 2002) showed evidence of spectral hysteresis with the orientation in the counter-clockwise direction. So, for at least one object, the hysteresis characteristics can change from flare to flare. Since the hysteresis probes the timescale of the cooling processes as a function of energy, it appears as though some flares have a cooling time that acts faster at higher energies, consistent with synchrotron cooling, while other flares exhibit cooling times that are longer. For the H1426+428 observations presented here, the generally clockwise hysteresis curves are consistent with the scenario in which the X-ray flux is dominated by synchrotron emission since the synchrotron cooling timescale is shorter at higher energies. In contrast, if the hysteresis loops were oriented in the opposite direction, one would expect the source to be in a regime in which the cooling and acceleration times were approximately equal, according to models of Kirk, Rieger, and Mastichiadis (1998).

By performing a study of the X-ray hysteresis, one can also explore the potential of various acceleration mechanisms to produce the TeV emission. Typically, clockwise hysteresis patterns (using the axis orientation shown in preceding figures) are interpreted as a soft time lag in X-ray emission, while counter-clockwise orientation is interpreted as a hard lag (Kirk, Rieger, and Mastichiadis, 1998; Kataoka, 2000). A soft lag at X-ray energies, which is interpreted as the energy dependence of the synchrotron cooling time, is characteristic of synchrotron self Compton models (Takahashi et al., 1996). A hard lag is not consistent with some standard one-zone synchrotron self Compton models, while models such as that of Kirk, Rieger, and Mastichiadis (1998), which incorporates a time-dependent propagating shock that accelerates electrons in the region of the shock front and allows for the case where $`\tau _{cooling}\tau _{accel}`$, can be consistent with a hard lag. The observations presented here are generally in agreement with standard models that predict a soft lag, but the lack of hysteresis from the flare in region B does not favor simple one-component models during that particular time period. However, the model of Li and Kusunose (2000) can reproduce hysteresis loops that are not clockwise (and presumably intermediate situations with no obvious loop orientation) by increasing the injection energy to the point where synchrotron self Compton losses begin to dominate the electron cooling. The TeV emission during this entire campaign was not in a high flaring state so firm conclusions regarding TeV acceleration during any short time period cannot be made.

Detailed modelling can be used to interpret the observed hysteresis patterns. Böttcher and Chiang (2002) have worked on such modelling of low-peaked BL Lacs, and they have found that external Compton and synchrotron self Compton mechanisms will produce different signatures at X-ray energies. However, this study was aimed at low-peaked BL Lacs, rather than high-peaked BL Lacs. We have found that a comparison to these models is inconclusive, probably due to the significant difference in the parameter space sampled by X-ray observations since low peaked objects are viewed in an energy region much higher with respect to the synchrotron peak, relative to high-peaked objects. Li and Kusunose (2000) did model the hysteresis characteristics of high peaked BL Lacs, but they confined themselves to only synchrotron self Compton mechanisms. They did find that the hysteresis characteristics are sensitive to total injection energy leading them to conclude that it would be difficult to draw any conclusions based upon hysteresis observations. However, for high-peaked BL Lac observations at energies at or below the synchrotron peak (such as the observations presented here), Li and Kusunose (2000) found that clockwise orientation could be expected when the injection energy was low and synchrotron losses were dominant, and the opposite orientation could be expected when the injection energy was high and inverse Compton emission was dominant. This supports the conclusion that for at least some of the flares observed during this campaign (those with clockwise orientation) the synchrotron losses were the dominant cooling mechanism. More detailed modelling, which is beyond the scope of the current paper, of X-ray spectral characteristics that can be used to differentiate between various high energy emission models for high-peaked BL Lacs is needed.

We acknowledge the helpful service of HEASARC and the RXTE guest observer facility, particularly the helpful interaction with A. Smale and K. Jahoda. This research is supported by grants from the U.S. Department of Energy.


format: text

text:
# Measurement of synchronization properties in systems with different excitation levels

## Abstract

We devised a measure based on the distributions of relative event timings of two coupled units. The measure dynamically evaluates temporal interdependencies between the two coupled units. Using this we show that even in the event of symmetrical coupling, non-identical units (having different control parameters) have small but persistent shifts in their event timings that change dramatically with the relative properties of their internal trajectories. We link that finding with the activity dependent synaptic modification in the brain showing that internally or externally driven excitation levels of the neuron may determine direction of information processing in a network.

It is hypothesized that the temporal structure of activity in the brain can play a crucial role in information processing von der Malsburg (1985). It was also established that relative timings between the spiking of coupled neurons may lead to long- and/or short- term synaptic modifications Bi and Poo (1998); McMahon and Barrionuevo (2002); Song and Abbott (2001). To be able to monitor successfully such spatio-temporal patterning in biological systems in general it is crucial to develop adequate methods to measure relative timings of the observed events. Any method designed for this purpose must be statistical, based on information obtained from many events but at the same time able to detect dynamical changes in the system. We have developed such a measure that is based on the dynamical monitoring of entropy changes in relative firing patterns between two coupled units. We will refer to it as a conditional entropy (CE) since it is calculated based on the relative timings of the events of both coupled units.

The relative timings measured for the two units are shown in Fig. 1. When the conditional entropy of the event timing pattern of unit 2 was calculated with respect to unit 1, we measured the timings of events that occurred at unit 2 with respect to the last event of unit 1. Conversely when the conditional entropy of unit 1 with respect to unit 2 was calculated, we measured the event timing of unit 1 with respect to the last event of unit 2. The probability distributions are dynamically updated (when the new event takes place) by increasing the bin of the distribution in which the latest timing falls by a fixed $`\mathrm{\Delta }P`$. Thus, if the timing falls in the bin $`I`$, $`P_I(t)=P_I(t1)+\mathrm{\Delta }P`$. The distributions were then renormalized. Thus, the $`\mathrm{\Delta }P`$ effectively determines to what degree the distribution will be skewed towards the timings of the newest events, so that the probability of the event timing after $`n`$ events declines to:

$$P_q(t)=\frac{1}{(1+\mathrm{\Delta }P)^n}P_q(tn).$$
(1)
The entropies are calculated from the normalized distributions, $`S=_IP_I\mathrm{ln}P_I`$, every time new events occur. We applied the measure to monitor relative timing patterns of two coupled non-identical Rössler oscillators Rössler (1976) and then two Hindmarsh-Rose models of the neurons Hindmarsh and Rose (1994). The coupled units differed through the value of their control parameters that influence dynamical properties of their trajectories. In light of the second example we will relate to those properties as excitation levels.

The event timings for the Rössler oscillator were defined as the times at which the oscillator’s trajectory (unit 1 or 2) crosses specified Poincare section ($`z=1`$). The equations for the coupled Rössler system are:

$$\begin{array}{c}\dot{x}_{1,2}=(z_{1,2}+y_{1,2})\hfill \\ \dot{y}_{1,2}=x_{1,2}+a_{1,2}y_{1,2}+\alpha (y_{2,1}y_{1,2})\hfill \\ \dot{z}_{1,2}=b+(x_{1,2}c)z_{1,2}\hfill \end{array}$$
(2)
where the subscript denotes the oscillator number; $`a_{1,2}`$, $`b=0.2`$, $`c=10.0`$ are the control parameters of the oscillators and $`\alpha =0.4`$ is the coupling. The values of parameters $`a_1`$ and $`a_2`$ were different for both units.

We found that the measure highlighted an important result. The conditional entropy of the oscillator with higher excitation level (a larger control parameter, $`a_i`$) is close to zero, whereas the entropy of the oscillator with lower excitation is significantly higher (Fig 2A). Moreover, a small change in the difference between the control parameters that changes their relative excitation levels causes the conditional entropies to change rapidly. We applied the measure to coupled oscillators having parameter $`a_{1,2}`$ span an extended range of values (Fig 2B); the observed phenomenon persisted over whole the range. Fig 2B plots grayscale map to represent the difference between the conditional entropies ($`S_1S_2`$) for the full range of the parameters of both units; white denotes minimal value whereas black denotes the maximal one. The difference values do not change significantly except in the neighborhood of the diagonal ($`a_1=a_2`$), where $`S_1=S_2=0`$. As the diagonal is crossed the values of the entropies switch and thus their difference takes opposite values. This is due to the fact that even though the coupling between the two units is symmetrical, there is asymmetry in the firing pattern that depends differentially on the relative excitation level of the coupled oscillators. That is, even if the firing seems to be coincident due to the coupling, the event timing of the unit that has higher excitation precedes narrowly that of the other unit. Fig 2C, D depicts this situation. We integrated the trajectories of the two oscillator for $`5000`$s; at $`t=2500`$ we changed the values of the control parameters so that at $`t<2500`$s, $`a_1<a_2`$ and at $`t>2500`$s $`a_1>a_2`$. We plot an example of the relative timing of the events of the two units for both regimes. Only at $`a_1=a_2`$, the two units synchronize with zero time lag, making their conditional entropies equal and zero. Fig. 3 shows the changes in the conditional entropies of the two units due to changes of relative values of the control parameter $`a`$. As the relative values of control parameter changes, the entropies change dramatically following the changes of the timing interdependencies of the two units: $`\alpha =0`$ in I and V sections - the two units are independent and both have high conditional entropies; $`\alpha =0.4`$ and the $`a_1=a_2`$ in section II - the conditional entropies of both units decrease to zero; when the units have different relative excitation levels (sections III and IV), their conditional entropies start changing rapidly at the onset of the sections, with the conditional entropy of the unit with higher excitation level tending to zero, and that of the other one being high.

We investigated the stability of the described phenomena for different values of coupling (Fig 4). As before, the control parameter of one unit was constant and that of the second unit was varied. When the coupling is equal to zero both oscillators are independent and their conditional entropies are high (Fig. 4A). However, the described timing interdependencies materialize even for small couplings. For a coupling as low as $`\alpha =0.2`$ the conditional entropies for both oscillators are significantly different (Fig. 4B). The conditional entropies $`S_1=S_2`$ and are equal to zero only when both oscillators have the same values of their control parameters (Fig 4C). Thus, this is the only time when the two oscillators mutually synchronize. At larger couplings the events of the oscillator with higher excitation precede within the narrow window the ones of the oscillator with lower excitation (Fig. 4D).

We also investigated the effects of noise on the relative values of conditional entropies. A term of the form, $`A\psi _{1,2}(t)`$ was added to the $`y`$-coordinate of each of the oscillators. The amplitude, $`A`$, varied from $`0`$ to $`10`$ and $`\psi (t)[1,1]`$ was a randomly generated, uniformly distributed, variable (different for both oscillators). The observed phenomenon was very robust with respect to noise (Fig. 5). The noise effectively destroyed the synchronization at the point where the control parameters of the oscillators were equal (Fig. 5B-D). This is due to the fact that at this point the noise induced jitter in the timings of the events causing the same oscillator to sometimes precede and sometimes follow the other one. However, the salient features of the behavior remain the same: the unit with the higher excitation level always fires in a narrow window before the one with the lower excitation level and thus the timing interdependencies remain unchanged when the control parameters are different. Thus, when the excitation levels of the oscillators are different, the conditional entropy of the unit with higher excitation is close to zero, whereas that of the other oscillator is high.

We show that the same phenomenon can be observed for a system of two interacting Hindmarsh-Rose models of the neuron. The equations for the coupled Hindmarsh-Rose neurons are:

$$\begin{array}{c}\dot{x}_{1,2}=y_{1,2}ax_{1,2}^3+bx_{1,2}^2z_{1,2}+I_{01,2}+\alpha (x_{2,1}x_{1,2})\hfill \\ \dot{y}_{1,2}=cdx_{1,2}^2y_{1,2}\hfill \\ \dot{z}_{1,2}=r\left[s(x_{1,2}x_0)z_{1,2}\right]\hfill \end{array}$$
(3)
where the subscript denotes the neuron number; $`a=1.0`$, $`b=3.0`$, $`c=1.0`$, $`d=5.0`$, $`r=0.006`$, $`s=4.0`$, and $`x_0=1.6`$ are the parameters of the model and $`\alpha =1.1`$ is the coupling strength. The parameter $`I_0`$ represents the amplitude of internal current to the neuron and is the control parameter of the system. We performed measurement of conditional entropy for two coupled neurons. Here, the event timings were the timings of spikes emitted by the neurons. The excitation levels of the two neurons varied (Fig 6). The conditional entropies changed in the same manner as in the case of the Rössler oscillators. That is the conditional entropy of the neuron with the higher excitation level was zero and that of the other neuron is significantly higher (Fig. 6A, B). This suggests that the behavior is quite general and not limited to purely theoretical models. As in the case of the Rössler oscillators, these dramatic changes in the CEs were due to the fact that the neuron with the higher excitation (i.e., higher $`I_0`$) fired a spike in a narrow window before the neuron with the lower excitation (Fig. 6C, D).

These results imply that a link can be made between levels of excitation and the dramatic changes in the spike timing interdependencies of the coupled neurons. The small shifts observed and measured may have big implications in connection to experimentally observed activity dependent synaptic modification. Experimentally measured synaptic Long Term Potentiation (LTP) and Long Term Depression (LTD) as well as short term changes are directly linked to the relative timings of the emitted spikes Bi and Poo (1998); McMahon and Barrionuevo (2002); Song and Abbott (2001). If the presynaptic neuron emits a spike in a narrow window before the postsynaptic neuron the connecting synapse (and in some cases other synapses) will be potentiated. Conversely, if the postsynaptic neuron fires systematically before the presynaptic neuron the related synapses will be depressed. We, on the other hand show that relative spike timings depend upon the relative level of excitation between the neurons in question. Thus we established a possible link between the relative levels of excitation of the neurons themselves and the pattern of synaptic modifications. Furthermore, we can postulate that since the synaptic modifications will depend differentially on the relative excitation levels of the coupled neurons, the direction of information processing will be directly linked to them: neurons having higher excitation levels will effectively drive the neurons with lower excitation levels, since this monosynaptic pathway will be potentiated. Since it is known that those levels of excitation can be controlled intrinsically by neurons (through expression of appropriate receptors) Desai and L.C. Rutherford (2002), or externally through the signal arriving at the given neuron, this implies that the directionality of information processing can be dynamically controlled by the information itself. We implement this hypothesis on a model of two coupled neurons, with non-symmetrical modifiable couplings. Changes in the coupling strengths, $`\alpha _{1,2}`$, of both neurons are linked to their respective conditional entropies. The general form of those changes is given by: $`\alpha _{1,2}(S)=\frac{\lambda }{1+\beta S_{1,2}^2}`$, where $`\lambda =1.1`$ and $`\beta =5.0`$ are constants. Initially the neurons have the same coupling, but different control parameters, $`I_{01}`$ and $`I_{02}`$. The coupling $`\alpha _1`$ (driven by the value of conditional entropy) from the neuron with higher excitation level got progressively stronger whereas that of the other one became significantly smaller until a steady state was reached, where $`\alpha _1`$ achieved its maximal value and $`\alpha _2`$ decreased to its minimal one. At t=3000ms the control parameter of the neurons was switched, so that the neuron with the previously lower excitation level now had a higher one. The change evoked rapid synaptic modifications (changes in the coupling strength $`\alpha _{1,2}`$) driven by the changing conditional entropies. When steady state was reached again the synapses effectively switched their efficacy (the values of $`\alpha _{1,2}`$ had reversed). Thus the neuron that was initially driven became a driver (Fig. 7A, B) - the information flow is reversed.

In conclusion, we have constructed a measure that monitors dynamically changing distributions of event timings of two coupled units. The use of this measure led to discovery of a relation between the event timings of two coupled units with different values of control parameters (i.e. excitation levels). The event timings of the unit with higher excitation level happen within a narrow window before those of the other unit. We have shown this phenomenon on two different systems: Rössler oscillators and Hindmarsh-Rose models of the neurons. We believe that the observed phenomenon together with the known experimental results of activity dependent synaptic modification may provide an important link between dynamic modifications of the network during the information processing and the content of the processed information itself. The obtained results have to be now investigated experimentally.