Photon Counting for Studying Faint Astronomical Variable Signals in Optical Band

Abstract

Although physics of neutron stars as pulsars together with their emission mechanisms leave discussions open, such objects represent the best targets to be deeply investigated by photon counting through the high-speed photometry technique. In this scenario, the capability of devices based on the silicon photomultiplier technology allows detecting single photons with remarkable time resolutions (few tens of nanoseconds). Whenever performing (optical) ground-based observations of variable sources, time of arrivals of incoming photons must be corrected because Earth’s reference frame system is not inertial. Time corrections provide time of arrivals to be moved to the Solar System Barycentre inertial reference frame. If the pulsar belongs to a binary system, further corrective terms, due to the orbital motion of the companion star, have to be taken into account. In this chapter, we report experimental results obtained from observations performed on two different variable sources, the isolated Crab pulsar and Hz Her/Her X-1 binary system, with a very fast custom astronomical photometer.

1. Introduction

Pulsars are highly magnetised, rotating neutron stars (NSs) (i.e. stars at the final state of their evolution) with an associated typical mass of about 1.4 M_⊙ (solar masses) as stated by the Chandrasekhar limit¹. These objects emit electromagnetic radiation, detected as periodic signals (i.e. pulses), preferentially in radio band (especially if they are isolated). In some cases, energies associated to these pulsations are spread all over the electromagnetic spectrum. Being possible to find them either isolated or bounded in binary systems, different emission mechanisms take place. In this chapter, we only focus on optical ground-based observations performed on such sources.

Since the Earth’s reference frame system is not inertial, the clock does not tick at a constant rate. It is a differential quantity due to both the Earth’s rotation and revolution around the Sun. Hence, timing of periodic signals constitutes a fundamental problem that must be taken into account when performing data analysis in order to obtain corrected measurements (i.e. spin periods, orbital parameters, rotational light curves). For this reason, the Solar System Barycentre (SSB) is chosen as the origin of a new reference frame system considered as inertial in a very good approximation. Several corrective timing factors have to be included to perform a complete analysis of isolated pulsars. The equation converting photon ToAs from the terrestrial system to the SSB can be written as follows:

$$t_{\mathit{SSB}}=t_{\mathit{obs}}+t_{\mathit{clk}}+∆R_{\odot }+∆E_{\odot }-∆S_{\odot }-∆\mathit{DM}$$E1

where t_obs is the observed ToAs of photons, t_clk is the set of clock corrections through which universal time coordinated (UTC) is converted to Barycentric Dynamical Time (BDT) passing through terrestrial time (TT), ∆R_⊙ is the Rømer delay, ∆E_⊙ is the Einstein delay, ∆S_⊙ is the Shapiro delay, and ∆DM is the delay due to the dispersion measure [1]. All these corrective terms are presented and explained in Section 2.

Pulsars can be also found to be bounded in binary systems, either low-mass X-ray binary (LMXB) or high mass X-ray binary (HMXB) depending on the mass of the companion star. In the first case, the mass of the latter does not exceed 1 M_⊙, while in the second case masses involved are larger than 5 M_⊙ [2]. One more class has been recently introduced for such systems called intermediate mass X-ray binary (IMXB), in which range of masses varies within the range (2 − 3) M_⊙.

The equation that permits to compute the true emission time must be completed with the following term:

$$t_{\mathit{em}}=t_{\mathit{obs}}-\frac{d}{c}-\frac{Z\left(t_{\mathit{em}}\right)}{c}$$E2

where t_obs is ToAs of observed photons, d is the distance to the system, c is the speed of light, and Z(t_em) is the distance between the so-called line of nodes and the pulsar itself (see Section 2 for details). The term d/c can be neglected being a constant.

In this scenario, a high-speed photometer working in optical band (ranged from 320 nm up to 900 nm), called silicon fast astronomical photometer (SiFAP) [3, 4, 5] was conceived and realised at ‘Sapienza University of Rome’. It is comprised of two channels: the first one is dedicated to study the science target, while the second one is devoted to monitor a reference star in the field of view (FoV). For both the two channels, a custom electronic chain capable to integrate the signal coming from the source in time windows down to 20 μs has been developed. Concerning the channel observing the target, an independent electronic system was added in order to tag the ToA of each single incoming photon with a time resolution of 25 ns.

A global positioning system (GPS) unit provides a reference time marker through the 1 pulse per second (1 PPS) signal, linked to UTC, with 25 ns resolution at 50% of the rising edge of the pulse itself. This signal drives two light-emitting diodes (LEDs) to have two optical markers superimposed on the data; the first one is at the beginning and the second one is at the end of the acquisition.

SiFAP was successfully mounted at 3.58 m Telescopio Nazionale Galileo (TNG) and 1.52 m Cassini telescope. Exciting results were obtained for the isolated Crab pulsar at TNG and Hz Her/Her X-1 binary system using Cassini telescope, respectively. By using two different data analysis approaches, Fourier and epoch folding search (EFS), a very good agreement between the observed spin period of the Crab pulsar and the one expected from ephemeris provided by Jodrell Bank Observatory (JBO) was obtained. The rotational light curve, showing the expected shape, was reconstructed for the Crab pulsar. Concerning Hz Her/Her X-1 system, its orbital period was found to be compatible with that extrapolated from the available literature, and the associated (simplified) orbital light curve was also successfully reconstructed.

2. The timing problem for isolated pulsars

2.1. Clock corrections

As stated earlier, the clock does not tick at a constant rate in the Earth’s reference frame system suffering from the effects due to both the Earth’s rotation and its revolution around the Sun. Time reference system is thus required to be constant. This can be achieved by considering some time conversions.

The first one is used to convert UTC into International Atomic Time (TAI). In this time domain, 1 s is defined as the time required for a ¹³³Cs atom at the ground state to complete exactly 9 192 631 770 oscillations [6]. It is given by:

$$\mathit{TAI}=\mathit{UTC}+N_{\mathit{ls}}+10s$$E3

where N_ls is the number of leap seconds. If the difference between UTC and UT1 (Universal Time, also known as astronomical time or solar time, referred to the Earth’s rotation) reaches 0.6 s, leap seconds are added in order to maintain this difference below 0.9 s. Until now (August, 2017), the number of leap seconds is 27, although 10 s more have to be considered because of historical reasons [7].

The second one in turn converts TAI into geocentric TT, using the following equation:

$$\mathit{TT}=\mathit{TAI}+32.184s$$E4

where 32.184 s is an offset arising from historical issues.

The last time conversion term consists and allows TT to be converted into BDT (i.e. the time one would have when photons are collected from the SSB). An approximate formula, composed of 791 coefficients, developed by Fairhead and Breatagnon in 1990 [8], permits to convert TT into BDT taking into account the effects due to the time dilation. Such a formula can be expressed as:

$$\begin{array}{l}BDT\approx TT+{\displaystyle \sum }_{i=1}^{478}{A}_i\sin \left({\omega }_{A_i}T+{\varphi }_{A_i}\right)+T{\displaystyle \sum }_{i=1}^{205}{B}_i\sin \left({\omega }_{B_i}T+{\varphi }_{B_i}\right)+\\ +T^2{\displaystyle \sum }_{i=1}^{85}{C}_i\sin \left({\omega }_{C_i}T+{\varphi }_{C_i}\right)+T^3{\displaystyle \sum }_{i=1}^{20}{D}_i\sin \left({\omega }_{D_i}T+{\varphi }_{D_i}\right)+T^4{\displaystyle \sum }_{i=1}^3{E}_i\sin \left({\omega }_{E_i}T+{\varphi }_{E_i}\right)\end{array}$$E5

where A_i, B_i, C_i, D_i, and E_i are coefficients expressed in microsecond, ω_Ai, ω_Bi, ω_Ci, ω_Di, and ω_Ei are angular velocities expressed in rad/10³ years. ϕ_Ai, ϕ_Bi, ϕ_Ci, ϕ_Di, and ϕ_Ei are angular phases expressed in radians, and T is the number of Julian centuries since January 1, 2000 (i.e. the beginning of the J2000 epoch). T can be written as:

$$T=\frac{\mathit{JD}-2451545}{36525}$$E6

where JD (in TT scale) is the time expressed in Julian Date² [8]. An accuracy of ~1 ns is achieved computing BDT through Eq. (5). More details describing the computation of the 791 coefficients can be found in Ref. [8].

2.2. The Rømer delay: A geometrical correction

The Danish astronomer Rømer was the first to study the geometrical correction caused by the variation of the path between the telescope collecting data and the SSB because of the Earth’s rotation and movement around the Sun. Taking into account these two effects, the classical light-travel time is affected by a delay, called Rømer delay.

According to Figure 1, it is possible to compute the vector r_B − O pointing from the SSB towards the observatory site as the difference r_B − O = r_E − O − r_B − E [9] if the vector r_B − E pointing from the SSB towards the Geocentre and the vector r_E − O pointing from Geocentre towards observatory site are known. In addition to the geometrical correction, parallax effects should be included whenever observing targets close to the solar system. In our case, paths of incoming photons can be considered as parallel (i.e. not affected by the parallax) in the whole solar system in a good approximation because we are studying point sources located at infinity. Hence, considering such an approximation, the equation for the geometric time correction (i.e. the Rømer delay, ∆R_⊙) can be expressed as:

$$∆R_{\odot }=\frac{{\mathit{r}}_{B-O}·\widehat{\mathit{n}}}{c}$$E7

where $\widehat{\mathit{n}}$ is the unitary vector for the incoming photons and c is the speed of light.

Figure 1.

Schematic representation of the Rømer delay.

2.3. The Einstein and Shapiro delays: Relativistic corrections

The time-varying gravitational potential and the Doppler shifts experienced by the observatory clock cause the so-called Einstein delay (∆E_⊙). This delay is the combined effect of two different terms: the gravitational redshift³ and the time dilation due to motions of the Earth and other bodies.

The first effect can be corrected through Eq. (5), considering the clock to tick at the SSB. The second effect can be computed and corrected by using ephemeris (position and velocity) of the observatory site. The analytic equation describing the Einstein delay is given by:

$$∆E_{\odot }=\frac{{\mathit{r}}_{E-O}·{\mathit{v}}_{\mathit{\text{Earth}}}}{c^2}$$E8

where v_Earth is the Earth’s velocity with respect to the SSB.

When the light coming from a distant source travels close to massive objects, its path is not straight anymore but curved because of their gravitational fields. Because a curved line is longer than a straight path, the light will need more time to complete its travelling. The Shapiro delay (∆S_⊙) caused by the space–time curvature around massive objects can be thus written as:

$$∆S_{\odot }=-\frac{2G{M}_{\mathit{\text{Body}}}}{c^3}\mathit{ln}\left(1+\mathit{cos}\theta \right)$$E9

where G is the gravitational constant, M_Body is the mass of the current object of which the gravitational field is considered, and θ is the angle between the pulsar and the Earth as seen from the Sun [11].

2.4. The dispersion measure delay

The time delay due to the propagation of a radiation with a given frequency with respect to one of infinite frequencies along a path of length d from the pulsar to the Earth is called Dispersion Measure (ΔDM) delay. The Dispersion Measure is expressed as follows:

$$\mathrm{\Delta }\mathit{DM}=\frac{1}{c}{\int }_0^d\left(1+\frac{f_p^2}{2f^2}\right)\mathit{dl}-\frac{d}{c}=\mathfrak{D}\frac{\mathit{DM}}{f^2}$$E10

where f_p is the frequency of a partially ionised plasma crossed by a signal with a given frequency f, and $\mathfrak{D}$ is the Dispersion Constant. The Dispersion Constant is defined as:

$$\mathfrak{D}\equiv \frac{e^2}{2{\mathit{\pi m}}_{e}c}$$E11

where e is the fundamental charge, and m_e is mass of the electron.

DM is the so-called Dispersion Measure, computed through the following integral:

$$\mathit{DM}={\int }_0^d{n}_e\mathit{dl}$$E12

where n_e is the Galactic electron density distribution.

According to Eq. (10), DM is equal to zero in X-ray band being proportional to the inverse of the square of the signal’s frequency. Moreover, it can be neglected for observations performed in optical band, while it must be taken into account for those in radio band [1].

3. Pulsars in binary systems

As a pulsar is not isolated but bounded to a companion star in a binary system, further corrective terms must be considered together with those explained in Section 2. The time related to the emission of radiation towards the observer varies cyclically depending on the orbital phase in which the pulsar is. It is immediately clear that for binary systems Doppler effect (concerning only the motion of the pulsar along the line of sight) must be considered. Starting from the Kepler’s laws, a binary orbit can be fully described by seven parameters, five of which provide to define both the shape and orientation of the orbit, and the remaining two determine how the considered body moves along its orbit.

According to Figure 2, for data analysis of binary systems, it is sufficient to consider the following:

• the semi-major axis (a), defined as a = [q/(1 + q)]a_sep, where q = M₂/M₁ which is the masses of the companion star and compact object, respectively, and a_sep is their orbital separation;

• the eccentricity (e) of the orbit;

• the inclination angle (i) as the angle between the orbit plane and the plane of the sky;

• the argument of periastron (ω) orienting the ellipse in the orbital plane. It is defined as the angle measured from the ascending node to the periastron;

• the epoch (T_asc) at which the object is passing on the ascending node of the orbit. Sometimes this parameter is replaced by the epoch $\left(T_{\frac{\pi }{2}}\right)$ of the passage at the superior conjunction (shifted by an angle of 90^° with respect to the ascending node);

• the mean anomaly (m) specifies the current position of the body at a given time t. It is expressed as m = 2π(t − T_asc)/P_orb, where P_orb is the orbital period of the binary system.

Figure 2.

Orbital elements for a binary system. The line of nodes is the intersection between the orbital and sky planes.

The distance Z(t) between the so-called line of nodes and the body (see Figure 2), when projecting the orbit onto the plane of the sky, can be expressed with a first-order approximation as:

$$Z\left(t\right)=a\mathit{sin}i\left[\mathit{sin}\left(m+\omega \right)+\frac{1}{2}e\mathit{sin}\left(2m+\omega \right)-\frac{3}{2}e\mathit{sin}\omega \right]$$E13

As stated earlier, Z(t) is measured starting from the line of nodes, and it is considered as positive when the pulsar is farther than the line of nodes considering observer’s reference direction. In this scenario, ToAs of observed photons are delayed with respect to emitted ones. Time of true emission can be computed through Eq. (2).

When a binary system is observed, the Doppler shift affecting the observed spin frequency (ν_obs) of the pulsar must be considered. Being the projected velocity of the pulsar $\dot{Z}\left(t\right)$, such an effect can be quantified as:

$${\nu }_{\mathit{obs}}-\nu \left(t\right)=\nu \left(t\right)\frac{\dot{Z}\left(t\right)}{c}=u\frac{2\mathit{\pi \nu }\left(t\right)}{P_{\mathit{orb}}}\left[-\mathit{cos}\left(m+\omega \right)-e\mathit{cos}\left(2m+\omega \right)\right]$$E14

where u = a sin i/c. Hence, if the orbital Doppler shift is not accounted for, it greatly limits the total amount of time over which pulsars can be observed and affects the data obtaining then a distorted folded pulse profile. Moreover, a maximum exposure time such that even a frequency residual δν does not distort the folded pulse profile can be defined. For a given time bin ∆t, the Doppler effect causes a variation on the frequency that can be evaluated as:

$${\mathit{\delta \nu }}_{\mathit{\text{Dopp}}}\simeq {\dot{\nu }}_{\mathit{\text{Dopp}}}∆t\le \mathit{u\nu }{\left(\frac{2\pi }{P_{\mathit{orb}}}\right)}^2$$E15

In order to satisfy this condition, we must have that:

$${\mathit{\delta \nu }}_{\mathit{\text{Dopp}}}\le {\mathit{\delta \nu }}_{\mathit{min}}$$E16

where δν_min = 1/(nT_exp), with n the number of phase bins which the folded spin period has been divided into, and T_exp is the total exposure time. In this way, Eq. (15) can be rewritten as:

$$∆t_{\mathit{max}}\le {\left(\frac{1}{\mathit{nu\nu }}\right)}^{\frac{1}{2}}\frac{P_{\mathit{orb}}}{2\pi }$$E17

The knowledge of the relevant orbital parameters is surely needed to reconstruct corrected light curves. Starting from Eq. (2), and considering circular orbits (i.e. e = 0), the delay due to the Doppler effect on the true emission time of photons can be thus computed as:

$$t_{\mathit{em}}=t_{\mathit{obs}}-u\mathit{sin}\left\{\frac{2\pi }{P_{\mathit{orb}}}\left[\left(t_{\mathit{obs}}-T_{\mathit{asc}}\right)-\frac{1}{2}{\dot{P}}_{\mathit{orb}}\frac{{\left(t_{\mathit{obs}}-T_{\mathit{asc}}\right)}^2}{P_{\mathit{orb}}}\right]\right\}$$E18

$$t_{\mathit{em}}=t_{\mathit{obs}}-u\mathit{cos}\left\{\frac{2\pi }{P_{\mathit{orb}}}\left[\left(t_{\mathit{obs}}-T_{\frac{\pi }{2}}\right)-\frac{1}{2}{\dot{P}}_{\mathit{orb}}\frac{{\left(t_{\mathit{obs}}-T_{\frac{\pi }{2}}\right)}^2}{P_{\mathit{orb}}}\right]\right\}$$E19

where both t_obs and T_asc $\left(T_{\frac{\pi }{2}}\right)$ must be in the same clock reference system (e.g. either UTC or BDT).

Orbital parameters that are mandatory to correct ToAs of photons coming from binary systems can be obtained in two ways. The first one consists in retrieving already existing ephemeris available in the literature. Unfortunately, problems could occur if tabulated parameters are too old and thus are needed to be updated.

In particular, in order to have more reasonable values of both epoch of ascending node and orbital period with its time derivative, if these two parameters are known at a specific reference time, they can be propagated [12] at the epoch of observation as:

$$T_{\mathit{asc}}=T_{\mathit{asc}}^0+N{P}_{\mathit{orb}}^0+\frac{1}{2}{N}^2{P}_{\mathit{orb}}^0{\dot{P}}_{\mathit{orb}}^0$$E20

$$T_{\frac{\pi }{2}}=T_{\frac{\pi }{2}}^0+N{P}_{\mathit{orb}}^0+\frac{1}{2}{N}^2{P}_{\mathit{orb}}^0{\dot{P}}_{\mathit{orb}}^0$$E21

$$P_{\mathit{orb}}=P_{\mathit{orb}}^0+{\dot{P}}_{\mathit{orb}}^0∆t$$E22

where $T_{\mathit{asc}}^0\left(T_{\frac{\pi }{2}}^0\right)$, $P_{\mathit{orb}}^0$, and ${\dot{P}}_{\mathit{orb}}^0$ are the epoch of ascending node (superior conjunction), the orbital period, and its first time derivative tabulated at a given reference time. ∆t is the time interval spanning from $T_{\mathit{asc}}^0\left(T_{\frac{\pi }{2}}^0\right)$ up to the observation date, and N is the nearest integer of the ratio between ∆t and $P_{\mathit{orb}}^0$ [13]. The second way to get ephemeris is performing several observations in order to have many spin frequencies as a function of time.

Relativistic corrections (i.e. Einstein and Shapiro delays) due to the gravitational field of the pulsar’s companion can be neglected for optical observations, not affecting the results. In particular, the Einstein delay is proportional to the eccentricity of the system and thus close to zero for small values of the eccentricity.

Although the Shapiro delay can contribute for tens of microseconds in the worst cases, it is not always possible to estimate its direct contribution. First, its contribution can be separated from that of the orbital motion only for nearly edge-on systems (70^° < i ≤ 90^°) being dependent on the inclination angle.

Second, since the knowledge of both the mass of the companion star and that of the pulsar are sometimes unknown or known with a large relative error, time distortions could be even introduced when applying such a correction. Moreover, astrometric effects like parallax and proper motion (taken into account when considering coordinates of the object) can be also neglected because their contribution is much smaller than the effect caused by geometric issues (i.e. the Doppler shift).

4. SiFAP: high-speed photometry of variable sources

As already mentioned in Section 1, SiFAP was developed at the Department of Physics of ‘Sapienza University of Rome’ to study variable sources since 2009. During subsequent years, the instrument was improved to achieve both better photometric sensitivity and time resolution. The present version of SiFAP is comprised of two channels, Channel 0 observing the science target and Channel 2 monitoring a reference star in the FoV by using their own dedicated sensor. The sky background signal is acquired by the same sensors moving the telescope away from the two objects.

The top-level diagram of the architecture of SiFAP is shown in Figure 3. Detectors used, especially selected from the production batch for this application, are based on the Silicon Photo Multiplier (SiPM) technology. They are called multi pixel photon counters (MPPCs) provided by Hamamatsu Photonics⁴.

Figure 3.

Top-level block diagram of the architecture of SiFAP.

Main characteristics of MPPCs (C13366 series) used by SiFAP are as follows:

• a double-stage Peltier cell cooling the sensor down to −20°C;

• an active area of 1.3 · 1.3 mm2;

• a squared pixel size of 50 μm;

• a typical dark count rate of ~2500 cps (count per second);

• a photon detection efficiency (PDE) (the percentage of detected photons with respect to incident ones) covering the spectral range (320 – 900) nm. It is peaked at a wavelength of 450 nm (blue) with a value of ~40%.

The whole system lies in a single unit, which contains two main custom blocks named photon counting gate (PCG) and photon temporal tagging (PTT), respectively. PCG is capable to collect data coming from the sensors in fixed integration time windows down to 20 μs, thanks to the handshake between a field-programmable gate array (FPGA) and a micro-processor (μP). PTT is an embedded system dedicated to observe only the science target (Channel 0) and developed to allow single photon tagging with a time resolution of 25 ns. The electronic chain of MPPC sensors provides counting photons with longer integration time windows within the range (1 – 100) ms.

A further custom unit, named MS (Master Sync), was realised to both remotely control and synchronise PCG and PTT units. MS is composed of a dedicated micro-controller (μC) connected to a PC through a RS232 serial interface.

The GPS unit allows to link UTC to the 1PPS signal used to drive an electrical-optical transducer (i.e. LED) in order to have an optical temporal marker superimposed on the data at the beginning and the end of the acquisition. More detailed technical descriptions of SiFAP instrument are available in [3, 4, 5].

5. Data analysis techniques for variable signals

5.1. Discrete Fourier analysis

When observing pulsars searching for any regular pulsations (i.e. variations of flux) hidden into the signal, many numerical analysis techniques can be used. One of the most used methods is the Fourier analysis, which provides a transition from the time domain to the frequency domain. The continuous function that permits such a transition is the so-called Fourier transform, which is expressed as:

$$H\left({\nu }_k\right)={\int }_{-\infty }^{\infty }h\left(t\right)e^{2\mathit{\pi i}{\nu }_{k}t}\mathit{dt}$$E23

where H(ν_k) is the amplitude of the transformed signal h(t), and ν_k is the sampling frequency.

For M data samples, the Fourier transform can be computed on M independent outputs only at a frequency value of:

$${\nu }_k=\frac{k}{\mathit{M\Delta t}}$$E24

with k =  − M/2, …, M/2.

According to the Nyquist’s theorem, the maximum achievable frequency (ν_max) is given by:

$${\nu }_{\mathit{max}}=\frac{1}{2\mathit{\Delta t}}$$E25

so that frequencies higher than ν_max cannot be resolved. Because observational data are not continuous but sampled in discrete time intervals (Δt), the Fourier transform for M data samples can be rewritten in its discrete version as:

$$H\left({\nu }_k\right)=\mathrm{\Delta }t\sum _{l=0}^{M-1}{h}_l{e}^{\frac{2\mathit{\text{πikl}}}{M}}$$E26

which represents the so-called discrete Fourier Transform (DFT). The accuracy (Δν) on frequency peaks achievable by such a method strongly depends on the total exposure time (T_exp) as:

$$\mathrm{\Delta }\nu =\frac{1}{T_{\mathit{exp}}}$$E27

systematically affecting the estimates of frequencies. Unfortunately, since observations require the number of data samples to be very large in order to have both more statistics on faint signals and time accuracy, it would be appropriate using a fast algorithm (e.g. Cooley-Tukey) to perform the Fourier transform (FFT, Fast Fourier Transform).

The Fourier technique is particularly suitable in astronomy to detect periodicities. Being these periodic signals usually unknown and faint, their detection is realised by computing the power spectrum representing the frequency distribution of the squares of the Fourier transform coefficients.

In fact, the Parseval’s theorem states that:

$${\sum }_l{h}_l^2=\frac{1}{M}{\sum }_k{\left[H\left({\nu }_k\right)\right]}^2$$E28

This means that the squared modulus of a function in the time domain is equal to the sum of the squares of its projection in the frequency domain. Hence, using the Parseval’s theorem, the power spectrum can be written as:

$$P\left(\nu \right)=\frac{2}{N_{\gamma }}{\left|H\left(\nu \right)\right|}^2$$E29

where N_γ is the total number of photons collected.

5.2. Epoch folding search

Unlike the Fourier analysis, EFS is a more refined method and is directly performed in the time domain. The standard approach to EFS consists in taking a dataset of a given total exposure time T_exp and defining a reasonable target period (P^⋆), which can be determined by a preliminary FFT analysis although not always there are clear features in it, and anyway its accuracy on frequency peaks is limited by the signal time bin.

Starting from P^⋆, a set of equispaced ‘trial’ periods (P_i) is created with a given time resolution. It is clear that P_i ∈ [P_min, P_max] with P_min < P^⋆ < P_max. In turn, each trial period is divided into n time bins (often a phase is used instead of a time, making n phase bins be within the range from 0 up to 1). Hence, it is possible to map the whole dataset with such time (or phase) bins, producing a corresponding number of folded curves.

In this context, it can be easily verified that in absence of pulsations (or any secular trend), counts in each phase bin of folded curves produced by the set of trial periods are Poisson distributed with mean and variance best estimated by the mean number of counts per bin. Since the number of events in each phase bin is usually large, the number of counts (x_i) in the i-th bin is normally distributed with the mean equal to the variance ${\mu }_{\mathit{exp}}={\sigma }_{\mathit{exp}}^2$. The statistic can be thus expressed as:

$$S=\sum _{i=1}^n\frac{{\left(x_i-{\mu }_{\mathit{exp}}\right)}^2}{{\mu }_{\mathit{exp}}}$$E30

where S represents a χ² distribution with n − 1 DoF (Degrees of Freedom). Under the hypothesis of absence of pulsation, one expects that S ≈ n − 1. Hence, if S is much greater than its expected value, the statistics belongs no longer to a χ² distribution, showing that a nonuniform behaviour (periodic) is present in the acquired data.

The original procedure used for calculating S was provided by Leahy in 1983 [14] adopting the following parameters:

• R = x_tot/T_exp, where R is the total counting rate, x_tot is the total number of valid acquired events, and T_exp is the total exposure time;

• R_i = x_i/T_i, where R_i is the counting rate at the i-th bin, x_i is the number of counts in the i-th bin, and T_i is the time duration of i-th bin (it can differ bin-to-bin because of possible gaps in the acquisition);

• μ_exp = μ_exp, i = RT_i, where μ_exp is the expected count rate which cannot be the same for each bin as T_i varies;

• ${\sigma }_i^2=R/T_i$, where ${\sigma }_i^2$ is the variance of i-th bin.

Starting from Eq. (30), and after having defined these parameters, S can be expressed as:

$$\sum _{i=1}^n\frac{{\left(R_i{T}_i-R{T}_i\right)}^2}{R{T}_i}=\sum _{i=1}^n\frac{{\left(R_i-R\right)}^2}{{\sigma }_i^2}$$E31

A first threshold on the accuracy of the spin period of a pulsar can be estimated by guessing how much two periods must differ in order to provide two different profiles within statistical fluctuations. Another relation is based on the assumption that a shift of one phase bin when folding the light curve to a pulse profile with n phase bins has a significant influence of the pulse shape.

This assumption leads to write:

$$\delta P_{\mathit{min}}=\frac{P^2}{2n{T}_{\mathit{exp}}}$$E32

According to Monte Carlo simulations done by Leahy himself in 1987 [15], the accuracy on the best folding period found by EFS can be expressed as:

$${\sigma }_P=\frac{P^2}{2T_{\mathit{exp}}}{\left[\frac{S_{\mathit{max}}}{n-1}-1\right]}^{-0.63}$$E33

where S_max is the maximum value of S (corresponding to the best folding period).

6. Observational results

In this section, we report results of data analysis of two different pulsar signals acquired with SiFAP at 3.58 m TNG⁵ (Telescopio Nazionale Galileo, Observatorio Astronomico Roque de Los Muchachos, La Palma, Canary Islands, Spain) and 1.52 m Cassini telescope⁶ of Bologna Astronomical Observatory (Loiano, Bologna, Italy).

6.1. SiFAP at TNG: the Crab pulsar

The first science case we report here is the well-known isolated Crab pulsar (PSR B0531 + 21). The Crab pulsar is a relatively young NS, which is in the central region of Crab Nebula⁷ (see Figure 4), the remnant of supernova SN 1054, and it was the first pulsar to be connected with a supernova remnant [16]. Isolated pulsars are in general associated to radio emissions, but some of them (like Crab pulsar) show also optical counterparts.

Figure 4.

Top: X-ray and optical combined images of Crab Nebula. The emitted radiation is well visible through magnetic poles. Bottom: optical image of Crab Nebula taken using Cassini telescope. The pulsar is in the central region of the Nebula and is indicated by two orthogonal segments.

This source was chosen as a benchmark because of being one of the most widely studied variable sources in every band of the electromagnetic spectrum. Moreover, this object provides a very high signal-to-noise (S/N) ratio, allowing safe detection of the pulsed signal, even without performing post-processing data analysis if the primary mirror of the telescope is enough large (3 m class or larger). In addition, it is very easy to obtain a large amount of data, which allow making comparisons among different data observations for such an object. Main properties of Crab pulsar are summarised in Table 1:

Crab pulsar
Constellation		Taurus
Right ascension [Epoch J2000]	α	05^h34^m31^s.97232
Declination [Epoch J2000]	δ	+22^°00^′52^′′.069
Apparent magnitude	Vmag	16.5 mag
Spin period @ (02/15/2014)	Ps	0.033692938448829(12) s
Spin period first derivative	${\dot{P}}_s$	4.2(1) · 10^−13 s/s
Distance	d	2.2(5) kpc

Table 1.

Main properties of Crab pulsar.

1σ uncertainties affecting the last digit(s) are presented within parentheses

It is worth to emphasise that the spin period of the pulsar slows down by about 38 ns per day due to large amounts of energy carried away in the pulsar wind, although the second-time derivative of the spin period must be kept in mind to be not constant. Such an effect is thought to be due to rotational instabilities (glitches) of NS, which lead to a strong evidence for the existence of a fluid component inside it. During a glitch, in fact, a small sudden increase in rotation rate is observed in the sudden early arrival of pulses [1]. The light curve of the Crab pulsar is reported in Figure 5 shows its typical double-peaked shape.

Figure 5.

The Crab pulse profile showing the optical light curve (o) in the range (600 − 750) nm, the average radio light curve at 1380 MHz (r), and a single giant pulse at 1357.5 MHz (gr). Two periods are shown for clarity. Phase 0 corresponds to the arrival at the SSB of the peak radio pulse. The optical light curve was divided into 5000 phase bins. Image taken from http://www.ing.iac.es/PR/SH/SH2003/triffid.html.

The first taller and thinner peak is believed to be due to the intrinsic emission of the pulsar, while the second one, lower and larger, should be caused by the re-emission of the first peak by the surrounding Nebula. Despite the Crab pulsar has a faint Optical counterpart, it is very bright in X-ray band, and the flux density and spectrum are known to be constant, with the exception of the pulsar itself. The pulsar provides a strong periodic signal that is usually used to check the timing of X-ray detectors.

In X-ray astronomy, this source is sometimes used as a flux density calibrator. In fact, crab and millicrab units were introduced. In particular, a millicrab corresponds to a flux density of about 2.4 · 10¹¹ erg/s/cm² (2.4 · 10¹⁴ W/m²) within the range from 2 to 10 keV, for a ‘crab-like’ X-ray spectrum. Very few X-ray sources ever exceed one crab in brightness.

SiFAP observed the Crab pulsar at TNG on 26-02-2014 [17] for about 2400 s. After having performed both Fourier and EFS techniques, we reconstructed its rotational (spinning) light curve, shown in Figure 6.

Figure 6.

Top: Crab pulsar light curve obtained from data collected with SiFAP by using the Fourier approach. Bottom: Crab pulsar light curve obtained from data collected with SiFAP by using the EFS approach. The phase shift is due to the slight difference between the two computed spin periods.

We obtained two slightly different spin periods for the Crab pulsar. FFT analysis provided a spin period equal to (0.0336927957 ± 0.0000000014) s, while the one obtained with EFS was (0.0336929420 ± 0.0000000050) s, causing a phase shift as well. These two results had to be compared to the one extrapolated from JBO⁸ database at the same observation date. The difference between the Crab pulsar spin period estimated from JBO and ours computed using FFT is ~140 ns. Such a difference was reduced to ~3 ns in the case of EFS, demonstrating that this method is more refined and robust.

It is worth to emphasise that several previous observations were performed using the Cassini telescope in order to optimise both the observational and analysis strategies.

6.2. SiFAP at Cassini telescope: Hz her/her X-1 binary system

Hz Her/Her X-1, also known as 4 U 1656 + 35, is a moderately strong X-ray binary source first studied by the Uhuru satellite. Such a system is classified as Intermediate Mass X-ray Binary (IMXB) being composed of a pulsar (Hercules X-1, spinning with a period of about 1.24 s) accreting mass from an A7 star (Hz Herculis) [18].

Hz Herculis was discovered in 1936, and classified as a variable star, while its collapsed companion was identified in 1972, thanks to the X-ray emission of the latter. Orbital parameters obtained from X-ray observations (i.e. X-ray ephemeris) and main features of Hz Her/Her X-1 [12, 19] are summarised in Table 2.

Binary system
Right Ascension [Epoch J2000]	α	16^h57^m50^s.5
Declination [Epoch J2000]	δ	+35^°20^′52^′′
Distance	d	6.6(4) kpc
Orbital period	Porb	1.700167202(1) d
Orbital period first derivative	${\dot{P}}_{\mathit{orb}}$	−1.33(7) · 10^−8 d/yr
Projected orbital radius	a sin i/c	13.1831(4) lt-s
Superior conjunction	$T_{\frac{\pi }{2}}$	54345.558195(80) MJD
Inclination angle	i	83(4)^°
Eccentricity	e	<1.3 · 10^−4
Her X-1
Mass	MX − 1	1.5(3) M⊙
Spin period	Ps	1.237744750(60) s
Spin period first Derivative	${\dot{P}}_s$	−1.778(56) · 10^−13 s/s
Magnetic field	B	4.1 · 10^12G
X-ray luminosity	Lx	3 · 10^37 erg/s
Hz Her
Mass	MHz	2.3(3) M⊙
Radius	RHz	4.2(2) R⊙
Apparent magnitude	Vmag	13.2 − 14.7 mag

Table 2.

Orbital parameters and main properties of Hz Her/Her X-1 binary system.

1σ uncertainties affecting the last digit(s) are presented within parentheses

The possibility to observe periodic occultations between the two stars was favoured by the inclination angle, which is close to 90^° (edge-on systems). Hz Herculis is included in the eclipse variable stars because its optical variability is due to the darkening by the collapsed companion (Her X-1). The magnitude variation of Hz Herculis in B-band is reported in Figure 7.

Figure 7.

B-band orbital light curve of Hz Herculis ([19] and reference therein).

A strong modulation of about 35 d in X-ray intensity, which is believed to be due to the occultation of NS by the accretion disk, is also shown by such a system. This variation has a characteristic shape, with two maxima called main-on and short-on as shown in Figure 8. The main-on, during which the maximum X-ray intensity can be observed, is followed by a 2 – 3 d lasting rapid drop of intensity called turn-off. The second maximum (called short-on or short-state) is about 78 d long, starts around ϕ₃₅ = 0.6, where ϕ₃₅ represents the 35 d phase.

Figure 8.

Light curve of the variability over 35 d, due to an accretion disk surrounding Her X-1 ([19] and reference therein).

The flux during the short-state has a maximum level of (10 – 20)% of the main-on intensity. In addition, it is worth to be noted that the optical signal is very much fainter than the X-ray one, with a pulsed fraction down to 0.1%.

Her X-1 generates a well-distinguishable X-ray pulse due to its rotation. Such a high-energetic radiation heats Hz Herculis atmosphere affecting its luminosity with a periodic cadence and making the spin signal to be detected in the optical band, thanks to thermal re-emission. It is also important to stress that the spin period is not monotonic but shows an evidence for spinning-up, in contrast with the monotonic spinning-down of the Crab pulsar. Such behaviour could be explained because of the inflow of matter towards NS, causing acceleration on its rotation.

SiFAP observed Hz Her/Her X-1 from August 25–28, 2016 at Cassini telescope. Four observations lasting about 2.5 h each were performed. Having only four data acquisitions (i.e. four data points and thus poor statistics) available, it was not possible to perform neither FFT nor EFS analysis.

In a first approximation, the expected orbital light curve of Hz Herculis, illustrated in Figure 7, can be considered to be well described by a sinusoid, although its shape is more complex. The orbital light curve reconstructed by analysing all data belonging to the four runs is reported in Figure 9. The plot shows four data points representing the instrumental magnitudes computed by averaging the total count of each observation as a function of time expressed in MJD unit.

Figure 9.

Light curve of the orbital period of Hz Her. Four data points represent the instrumental magnitude computed by averaging the total count of each observation as a function of time.

The value of the orbital period obtained by sinusoidally fitting data was found to be equal to (1.70017 ± 0.00058) d, in an optimal agreement with that expected one, differing by only 3.6 · 10⁻⁶ d. It is worth considering that this result could not be directly compared with that represented in Figure 7 because observations done with SiFAP were performed without using any filter. The resulting plot thus shows a profile integrated all over the optical band, the region of the electromagnetic spectrum where PDE of MPPC spans.

We also tried to reconstruct the rotational light curve without any success. This was due to two main factors. The first one concerned the poor knowledge about optical ephemeris of the system. It is known that the optical counterpart of Hz Her/Her X-1 system arises from the reprocessing of X-ray radiation. Therefore, X-ray ephemeris reported in Table 2 was not suitable to correct ToAs for the orbital motion of the system. In this way, we were not able to merge efficiently acquired data in order to increase the statistics and thus the S/N ratio. In fact, putting together observations which have not been previously corrected for the Doppler effect would smear pulsations out. This happens because spin periods (frequencies) are varying during the orbital period.

The second constraint was due to the physical condition of the binary system itself. Unfortunately, we estimated that our observations were performed when the system was in a not optimal orbital phase for being observed. Moreover, the periodic precession of the disk (35 d modulation) played an unfavourable role because the system was in the state with a very low intensity.

7. Conclusions

In this chapter, we showed how astronomical variable sources could be studied by photon counting through the (optical) high-speed photometry technique with ground-based telescopes and instrumentation. The most difficult effort to obtain correct results is analysing data taking into account several perturbation factors such as the Earth’s not inertial reference frame system and orbital motion of binary systems. Although all these effects could deteriorate measurements, anyway very good and promising results can be obtained by applying above-described time corrections.

In this scenario, the SiFAP instrument gave its contribution with interesting results, despite more specific and deeper knowledge on this kind of variable sources can be achieved, thanks to multiwavelength (even simultaneous) measurements. In fact, the opportunity to observe such targets in more than one electromagnetic band is very useful to understand both the physics and emission mechanisms and their possible correlation.

Acknowledgments

The authors acknowledge the research project Protocol C26A15YCJ4 funded by the Department of Physics of the ‘Sapienza University of Rome’ for the financial support.

Moreover, the authors would remind that results obtained by SiFAP on the Crab pulsar are based on observations made with the Italian 3.58 m Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.

In addition, the authors would also like to thank the 1.52 m Cassini telescope, which is run by INAF-Osservatorio Astronomico di Bologna at Loiano (Bologna), supported measurements on the binary system Hz Her/Her X-1 performed with SiFAP.