Alessandro D A M Spallicci
email: spallicci@cnrs-orleans.fr

Fabio Ragosta
email: fabio.ragosta@inaf.it

José A Helayël-Neto

Martín López-Corredoira

Supernovae and photon frequency shift induced by the Standard-Model Extension

Keywords: Standard-Model Extension, Photons, Light Propagation, Cosmology, Supernovae

We revisit for 714 SNeIa the discrepancy between the red-shift associated to the distance modulus µ and the spectroscopic red-shift. Previous work has shown that the total red-shift z might be a combination of the expansion red-shift zC and of a static, blue or red shift zLSV(r), r being the comoving distance. The latter is due to the energy non-conservation of the photon propagating through Electro-Magnetic (EM) background fields (host galaxy, intergalactic and Milky Way), under Lorentz(-Poincaré) Symmetry Violation (LSV), associated to the Standard-Model Extension (SME). The non-conservation stems from the vacuum expectation value of the vector and tensor LSV fields. For zero radiation Ω rad and curvature Ω k densities, and matter density Ωm = 0.28, the SN1a positions in the (µ, z) plan are recovered according to the different strengths, orientations, alignments and space-time dependencies of the EM fields and LSV components. The LSV vacuum energy may be thus tantamount to ΩΛ ≃ 0.7, but unrelated to an accelerated expansion. We present models with red or blue-shifts zLSV, below 10% of z. The ν frequency variation is below 10 -19 ∆ν/ν per m.

Introduction. At the infinitesimal scale, the Standard-Model (SM) unifies three of the four fundamental forces: Electro-Magnetic (EM), weak, and strong interactions. The SM guarantees Lorentz(-Poincaré) Symmetry (LoSy), stating that experimental results are independent of the time, position or velocity of the observer. Despite its remarkable successes, the SM does not incorporate the particles corresponding to a, yet to be experimentally found, dark universe. Indeed, it does not account for the accelerating expansion of the universe (dark energy) nor has a viable dark matter particle. Furthermore, the photon is the only free massless particle in the SM, after the attribution of some mass to the neutrino. Thus, a theory going beyond the SM has been put forward: the SM Extension (SME) [1,2]. The SM is assumed LoSy invariant up to certain energy scales beyond which a LoSy Violation (LSV) might occur. The SME allows testing the low-energy manifestations of the LSV.

Light propagation in the SME. In [3,4] the foundations of the analysis of light propagation in the SME were laid. Going beyond the SM, we encounter a massive photon, compatible with the actual upper limit [5] of 10 -54 kg, but as pointed out in [6], this limit arises chiefly from modelling rather than measurements. The effective, and frame dependent, mass is proportional to the value of the LSV parameters and, conversely to the de Broglie [7,[START_REF] De Broglie | La Mécanique Ondulatoire du Photon. Une Nouvelle Théorie de la Lumière[END_REF]-Proca [START_REF] Proca | [END_REF] formalism, the SME massive photon is gauge-invariant [3,4]. Moreover, in specific conditions, sub-and super-luminal velocities, imaginary and com-plex frequencies, birefringence appear, and, evidently in any condition, LSV anisotropy and inhomogeneity are present. The effects of the latter two are not dominant in the regimes we explored. The massive photon group velocity differs from c of a quantity proportional to the inverse of the frequency squared [START_REF] De Broglie | La Mécanique Ondulatoire du Photon. Une Nouvelle Théorie de la Lumière[END_REF]. Incidentally, such dependence has been analysed recently with the signals from Fast Radio Bursts (FRBs) [10][11][12][13][14][15]. Most remarkably, the energy-momentum tensor of a light-wave is not conserved in vacuum when light crosses an EM background field [4].

Non-conservation occurs also when the EM background field and the LSV perturbations are independent from space-time coordinates, but interestingly, for a space-time dependent LSV, the presence of an EM field is not necessary [16]. Indeed, the photon exchanges energy with the LSV and EM fields. An estimate of the energy change that light would undergo was given [16]. The wave-particle correspondence, even for a single photon [17], leads to consider that the light-wave energy nonconservation is translated into photon energy variation and thereby into a red or a blue shift. The energy variations, if losses, would translate into frequency damping.

The LoSy breaking factors are represented by a k AF α four-vector when the handedness of the Charge conjugation-Parity-Time reversal (CPT) symmetry is odd and by a k ανρσ F tensor when even. The physical reasons inducing non-conservation notwithstanding constant EM background and LSV breaking vector (the breaking tensor appears either under a derivative or coupled to a derivative of the EM background) lie in the Carroll-Field-Jackiw action [18] for the CPT odd handedness. The action contains a contribution such that even if the EM background (and the LSV perturbation) is constant, the corresponding four-potential is not. Thereby, there is an explicit x α coordinate dependence at the level of the Lagrangian. This determines a source of energymomentum non-conservation, according to the Noether theorem. Otherwise put, there is an exchange of energymomentum between the photon and the EM background.

The LSV frequency shift. For f αν representing the photon field, a ν the four-potential, F αν the EM background field, j ν an external current and the symbol * the dual field, the energy-momentum density flux θ α 0 variation is given by [4] (in SI units µ 0 = 4π × 10 -7 NA -2 )

∂ α θ α 0 = j i f i0 - 1 µ 0 ∂ α F αi f i0 + k AF α * F αi f i0 + 1 2 ∂ α k AF 0 * f αν a ν -1 4 ∂ 0 k αiκλ F f αν f κλ + ∂ α k αiκλ F F κλ f i0 + k αiκλ F (∂ α F κλ ) f i0 . (1)
The dimensions are Jm -4 . In absence of an external current and of a large scale electric field, but considering only the constant components of the galactic and intergalactic magnetic fields, and finally neglecting at optical frequencies the contribution of k αiκλ F , an estimate was given through a simplified expression along the line of TABLE I: LSV shift types. In the first, the frequency variation is proportional to the instantaneous frequency and to the distance; in the second to the emitted frequency and the distance; in the third only to the distance (we don't consider a fourth option proportional to the observed frequency and the distance). k 1,2 have the dimensions of Mpc -1 , k 3 of Mpc -1 s -1 . The positiveness of the distance r constraints z LSV/1 > -1 for k 1 < 0, and

-1 < z LSV/1 < 0 for k 1 > 0. Type 1 2 3 ∆ν k1νdr k2νedr k3dr νo νe exp k 1 r νe(1 + k2r) νe + k3r zLSV exp -k 1 r -1 - k2r 1 + k2r - k3r νe + k3r sight observer-source [16] ∂ α θ α 0 ≈ - 1 µ 0 k AF 0 * F 0i f i0 . (2)
The * F 0i represents the magnetic components of the dual EM background tensor field, f i0 the electric components of the EM photon tensor field, and k AF 0 is the time component of the breaking vector. Assuming that the energy variation of a light-wave corresponds to a frequency shift for a photon, we converted the energy variation of Eq.

(2) into a ∆ν, being ν the photon frequency.

After some steps, an exact expression for |∆ν| due to LSV was given [16]. Herein, we reproduce it qualitatively

|∆ν| LSV ∝ B f i0 t LB k AF 0 ̺ , (3)
where ̺, is a positive arbitrary safe margin, taking into account that the magnetic fields (of the galaxy of the source, inter-galactic, and the Milky Way) estimated at B = 5 × 10 -10 -5 × 10 -9 T each, and crossed by light from the source to us in a look-back time t LB , have likely different orientations and partly compensate their effects on the wave electric field f i0 . We have not considered any magnetic field at the source.

In the following, we refer to a general z LSV given by all the terms in Eq. ( 1).

The LSV as vacuum energy. According to [19], the LoSy breaking four-vector, k AF , and the rank-four tensor, k F , correspond, respectively, to the vacuum condensation of a vector and a tensor field in the context of string models. They describe part of the vacuum structure, which appears in the form of space-time anisotropies. Therefore, their presence at the right-hand side of Eq. ( 1) reveals that vacuum effects are responsible for the energy variation of light waves, which, in turn, correspond to a photon frequency shift. Otherwise put, vacuum anisotropies are the true responsible for changing the frequency of light emitted by astrophysical structures.

Superposing the shifts. We recall that z = ∆ν/ν o where ∆ν = ν e -ν o is the difference between the observed ν o and emitted ν e frequencies, or else z = ∆λ/λ e for the wavelengths.

We pose the following conjecture. Expansion causes λ e to stretch to λ c , that is λ c = (1 + z C )λ e . The wavelength λ c could be further stretched or shrunk for the LSV shift to

λ o = (1 + z LSV )λ c = (1 + z LSV )(1 + z C )λ e . But since λ o = (1 + z)λ e , we have 1 + z = (1 + z C )(1 + z LSV ); thus z = z C + z LSV + z C z LSV . (4)
The second order is not negligible for larger z C .

Behaviour of the LSV shift with distance. We model the behaviour of z LSV with distances in three different ways, Tab. I. Since z C stems from the universe expansion, we estimate z LSV on the basis of the comoving distance. The k i parameters, proportional to the Lemaître-Hubble-Humason constant H, herein 70 km/s per Mpc, take either positive (frequency increase) or negative (frequency decrease) values.

The LSV frequency shift from the SNIa data. A greater than expected SNIa distance modulus for a given redshift led to consider a cosmological constant or dark energy [20,21]. Later, suppositions of alternatives theories of gravity, see [22,[START_REF] López-Corredoira | Proc. Int. Conf. Cosmology on small scale[END_REF] for reviews, and of photons oscillating into axions producing a dimming effect for SNeIa [START_REF] Csáki | [END_REF]25] were made. It must be recalled that SNIa data reliability is still debated [26][27][28][29][30][31][32]. We challenge the meaning of the red-shift uniquely explained by expansion and conceive a complementary point of view.

For appreciating d L , as function of z C , we have computed for for Ω rad = Ω k = Ω Λ = 0

d L = 2c HΩ 2 m 2 -Ω m (1 -z C ) -(2 -Ω m ) 1 + Ω 2 m z C , (5
) using the Chebyshev polynomial expansion [33], which coefficients of are linear in Ω m .

We took a set of 714 SNeIa from the Union2 catalogue [34], devoted to the dark energy search, using data from the Hubble Space Telescope (HST) and combining multiple datasets. The catalogue provides the luminosity distance, d L , here in parsecs, related to the maximum of the light curve in the B-band

µ = m -M = 5 log d L -5 ( 6
)
where µ is the distance modulus, m and M are the apparent and absolute magnitudes, respectively. The Union2 catalogue provides also the error on µ, the red-shift z of the host galaxy and its error. The comoving distance r, computed from the Lemaître-Friedman-Robertson-Walker metric (FRWL), is given by [START_REF] Peebles | Principles of physical cosmology[END_REF][START_REF] Zhang | [END_REF]

r = 2c H ( √ z C + 1 -1) . ( 7
)
Our iterative computation starts from the cosmological parameters and compares the final simulated s and the observed value of the distance modulus

Ω k , Ω m (Ω Λ = 0) → d L → z c → r → z LSV → µ s ? = µ o ,
= Ω k = Ω Λ = 0, Ωm = 0.
28. The the error on k i is due to the Hierarchical Bayesian Model (HBM) approach, considering different source of uncertainties to render homogeneous the photometric data of the distance modulus from different datasets in Union 2 catalogue. There is no actual error propagation due to the fact that HBM takes into account the error considering it statistically implemented into the model analysis The intrinsic standard deviation σ for k 1 is 3 × 10 -6 . The simulation results provide a mean z LSV blue-shift for types 1 and 2 and a z LSV red-shift for type 3. where the z LSV shift may be positive or negative. Herein, we have held to the measured µ and z while decomposing the latter in two contributions z c and z LSV . Else, we could have considered the measured µ as the sum of two contributions µ(z C ) and µ(z LSV ) and looking for an agreement with spectroscopic data.

Modelling. For achieving the best fit, we employed a Hierarchical Bayesian Model (HBM) and determined the entire posterior probability density function of the these parameters. We engaged an HBM to retrieve from the distance distribution of the SNeIa all the possible values of k, for each of the three types in Tab. I. We assumed flat priors in the range [-100, 100] to avoid heavy constraints in exploring the domains of the parameters. The Bayesian analysis is very depending upon the choice of the priors, conversely to other kinds, as the Gaussian approaches [27][28][29][30][31]. The flat priors allow an even probability distribution and thus unbiased estimates.

Unfortunately, the HBM is heavily influenced by the error on the measured values, namely on the distance modulus leading to a statistical non-informative error on the model parameters. From the Gaussian shape of the posterior distribution of k i , it is evinced that a non-zero probability exists for having both positive and negative values of k i , due to the measurement errors. We decided to use HBM due to the multiple origin of the distance modulus estimates and differences in data reduction schemes, thus avoiding the influence of some systematic error in the estimates of the parameters. The HBM improves the precision of the estimates of the model parameters, even though the accuracy is rather coarse.

Results and discussion. Table II provides, for Ω rad = Ω k = Ω Λ = 0 and Ω m = 0.28, the k i values and the associated rms for the three types of k i associated to z LSV . The latter consists of a mean blue-shift for types 1 and 2 and of a red-shift for type 3. Figure 1 shows the interpolating curve well fitting µ and z data, for type 1. Figure 2 shows the distribution of possible values of z LSV for a given z, again for type 1. We have well interpolating µ(z) curves also for types 2 and 3 of z LSV shifts. The rms values of types 1 and 2 are similar, but the frequency variation proportionality to the instantaneous (type 1) rather than emitted (type 2) frequency determines a 'classic' exponential decay or growth.

The absolute value of z LSV increases with the distance. If blue, the drawback is to explain the photon gaining energy travelling from the source to us, which nevertheless may well occur as a net result of the different contributions of the LSV (vector and tensor) and EM fields along the path. The advantage is a clear understanding of why we measure a lower spectroscopic z, as z LSV is subtracted from the real z. If red, the drawback is that z LSV must be comparable or smaller than the observation error [37][38][39]. The advantage is the correspondence to an intuitive dissipation effect taking place along the photon path. Remarkably, for Ω m ≥ 0.75, k 1 turns negative and the z LSV becomes a red-shift as in model 3.

The peculiarity of our approach is that a single mech-anism explains all the positions of the SNeIa in the (µ, z) plan, including the outliers. This shift may be red for some SNeIa and blue for others, depending on cosmology parameters, as well as the orientation of the components of the LSV vector or tensor and of the EM fields (host galaxy, intergalactic medium, Milky Way). Evidence of dark energy anisotropy has been shown recently [32]. The LSV vacuum energy appears to act on light propagation and it is equivalent to Ω Λ ≈ 0.7, rather than igniting an accelerated expansion. Numerical applications were outlined [16]. Working them out when considering a SNIa at z = 0.5, t LB = 1.57 × 10 17 s, f i0 = 3.79 × 10 -9 V s m -2 , an average B = 2.75 × 10 -9 T, the proportionality factor in Eq. (3) becomes 5.87 × 10 47 . The parameter k AF 0 has a laboratory upper limit of 10 -10 m -1 but a more stringent, and less favourable for our study, astrophysical upper limit of 5.1 × 10 -28 m -1 . In this worst case, it is sufficient that ̺ ≥ 7.6 × 10 -8 , to get z LSV in the order of 10% of z.

Where Ω m was let free, the constraint Ω = 1 leads obviously to huge values of the matter density, as, incidentally, dark energy emerges at high red-shift.

Perspectives. We have applied to SNeIa a new frequency shift in vacuum occurring to a photon in a LSV scenario, accompanied by an EM field. The experimental and observational limits on LSV and magnetic fields are fully compatible with our findings.

Obviously, the LSV shift is generally applicable and it is not limited to the SNIa case. In future work, leaving aside the SME, we will show that a frequency shift is also produced by a generalised non-linear electromagnetism, encompassing the formulations of Born-Infeld and Euler-Kockel-Heisenberg. The found limits on z LSV will be applicable to z NL for the non-linear electro-magnetism.

The XIX century Maxwellian linear electro-magnetism and Einsteinian non-linear gravitation have been well tested. This has not impeded the proposition of alternative formulations of gravity. The lack of experimental proofs on the dark universe and the successes of general relativity prompt to revisit astrophysical observations, largely based on light signals, with non-Maxwellian electro-magnetism, opening the door to radically new interpretations. In future work, we will show that massive photons may fake time dilation for not-moving sources. Such effect, though, is very marginal for SNeIa.

The z LSV shift is a fraction below 10% of the measured red-shift and thus below 2 × 10 -19 ∆ν/ν per m. It would be desirable to test frequency invariance in vacuum [40][41][42][43] with a ground or space interferometer.

3 . 5 × 10 - 3 2. 6 × 10 - 5 ± 2 . 8 × 10 - 3 < 1 FIG. 1 :FIG. 2 :

 35103610528103112 FIG. 1: For type 1 of LSV shift, Tab. I, Ω rad = Ω k = Ω Λ = 0 and Ωm = 0.28, µ, in mag units, versus z, Eq. (4).

TABLE II :

 II Simulation results for Ω rad

Acknowledgments. FR thanks the Large Synoptic Survey Telescope (LSST) Data Science Fellowship Program, funded by LSST Corporation, the Brinson, Moore and National Science (NSF) Foundations, the Adler Planetarium, and the Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) at Northwestern University; he acknowledges the NSF CyberTraining Grant #1829740. Discussions with L. Barbosa de Lima, M. A. Lopes Capri, V. E. Rodino Lemes, R. de Oliveira Santos (Rio de Janeiro) and S. Capozziello (Napoli) are acknowledged.


Ann-Sofie Bak Nielsen

Gemma H Janssen

Golam Shaifullah

Joris P W Verbiest

David J Champion

Grégory Desvignes

Lucas Guillemot

Ramesh Karuppusamy

Michael Kramer

Andrew G Lyne

Andrea Possenti

Ben W Stappers

Cees Bassa

Ismaël Cognard

Kuo Liu

Gilles Theureau

Timing stability of three black widow pulsars

Keywords: binaries: close -pulsars: general -pulsars: individual: PSR J0023+0923, J2214+3000, J2234+0944

We study the timing stability of three black widow pulsars, both in terms of their long-term spin evolution and their shorter-term orbital stability. The erratic timing behaviour and radio eclipses of the first two black widow pulsar systems discovered (PSRs B1957+20 and J2051-0827) was assumed to be representative for this class of pulsars. With several new black widow systems added to this population in the last decade, there are now several systems known that do not show these typical orbital variations or radio eclipses. We present timing solutions using 7-8 yrs of observations from four of the European Pulsar Timing Array telescopes for PSRs J0023+0923, J2214+3000 and J2234+0944, and confirm that two of these systems do not show any significant orbital variability over our observing time span, both in terms of secular or orbital parameters. The third pulsar PSR J0023+0923 shows orbital variability and we discuss the implications for the timing solution. Our results from the longterm timing of these pulsars provide several new or improved parameters compared to earlier works. We discuss our results regarding the stability of these pulsars, and the stability of the class of black widow pulsars in general, in the context of the binary parameters, and discuss the potential of the Roche-lobe filling factor of the companion star being an indicator for stability of these systems.

INTRODUCTION

The first radio pulsar was discovered in 1967 by Jocelyn Bell Burnell (Hewish et al. 1968) and since then more than 2700 radio pulsars have been found (Manchester et al. 2005). This sample includes both isolated pulsars and pulsars in binary systems, and slow and millisecond pulsars (MSPs). MSPs are old pulsars that have gone through the recy-E-mail: nielsen@strw.leidenuniv.nl cling process (Alpar et al. 1982;Bhattacharya & van den Heuvel 1991). MSPs are spun up in X-ray binary systems through accretion, which transfers angular momentum onto the pulsar. During the X-ray pulsar phase, the magnetic field strength of the pulsar decreases (Chen et al. 2013). A special group of MSPs, the so-called 'spiders', with very low companion masses and tight orbits (P b <24 hr), are the 'redback' pulsars (RBPs), which have companion star masses of M2 0.1 -0.4 M , and the 'black widow' pulsars (BWPs) with M2 0.1 M (Chen et al. 2013;[START_REF] Roberts | Neutron Stars and Pulsars: Challenges and Opportunities after 80 years[END_REF]. The low mass companions of BWPs are often thought to be semidegenerate stars, similar to brown dwarfs (see section 4.4 in Lazaridis et al. 2011). The companions of RBPs are probably non-degenerate stars, as suggested by a few cases where the optical counterpart was discovered (Ferraro et al. 2001;[START_REF] Roberts | Neutron Stars and Pulsars: Challenges and Opportunities after 80 years[END_REF]Al Noori et al. 2018). The fate of the companion star in the BWP systems is unclear. However, one postulated scenario is that the companion could be fully ablated, creating an isolated MSP. Alternatively, tidal effects could play a significant role in eventually destroying the companion (Stappers et al. 1998;Chen et al. 2013).

Since the discoveries of PSRs B1957+20 and J2051-0827 (Fruchter et al. 1988;Stappers et al. 1996b), the populations of BWPs and RBPs have increased greatly due to recent surveys, using either radio observations which were optimised to find MSPs, or targeted radio surveys which are guided by high energy point sources such as those seen with the gamma-ray observatory Fermi. More than 40 MSPs were found using such surveys [START_REF] Roberts | Neutron Stars and Pulsars: Challenges and Opportunities after 80 years[END_REF]. The two original BWP systems show eclipses of the radio emission from the pulsar. The eclipses are due to material being ablated off of the companion star by a strong pulsar wind (Fruchter et al. 1988;Stappers et al. 1996b;Chen et al. 2013). The presence of this matter in the binary, often makes the BWPs complicated systems to use for high precision pulsar timing experiments. The companion stars are often bloated, and their unstable and in-homogeneous mass distributions exerts small torques on the pulsar, which in turn causes small changes in the orbit. Such orbital variations increase the number of parameters that are required in timing models, and mean that it is not possible to predict the timing behaviour and thus not possible to model the orbital variations for longer timescales. The orbital variability could also reduce the sensitivity of the timing solution to other signals of interest, such as for example gravitational-wave (GW) signals.

By observing multiple stable pulsars in different parts of the sky, an array of pulsars, and cross-correlating their timing residuals, it should be possible to obtain the sensitivity required to measure the expected quadrupolar signature of a nHz GW. This is the general idea behind the pulsar timing arrays (PTAs) (Foster & Backer 1990;Tiburzi 2018). With new BWPs being discovered, Bochenek et al. (2015) considered the feasibility of using BWPs in PTAs, by testing if fitting for multiple orbital frequency derivatives significantly reduced the sensitivity to the GW signal. They used simulated data sets of five different pulsars to constrain how sensitive a timing solution would be to a GW signal, if several orbital-period derivatives were fitted. They found that the sensitivity to the GW signals was not reduced significantly, partly due to the comparatively long periods of the GW signals and the typical BWP orbital period of 2 to 20 hr. Recently a few BWPs have been included in PTAs, e.g. PSR J0610-2100 is already used by the EPTA (Desvignes et al. 2016). There are three other BWPs that show promising timing precision and which are candidates for use in a PTA: PSRs J0023+0923, J2214+3000, and J2234+0944. All three are currently observed by the North American Nanohertz Observatory for Gravitational Waves (NANOGrav, Arzoumanian et al. 2018) and the European Pulsar Timing Array (EPTA, this work). PSR J0023+0923 was discovered with the Green Bank Telescope (GBT) in a survey of unassociated Fermi γ-ray sources (Hessels et al. 2011). It is an MSP with a spin period of ∼3 ms, an orbital period of about 3.3 hr and a minimum companion mass of 0.017 M . The optical counterpart to PSR J0023+0923 was discovered by Breton et al. (2013) with the Gemini North telescope. By combining measurements of the filling factor with the estimated distance from the dispersion measure value they inferred that the companion star never fills its Roche lobe. Their results show that, when taking the small filling factor into account, the companion star could be as small as 0.05 R . PSR J0023+0923 has also been observed in X-rays, but no pulsations were detected (Ransom et al. 2011).

Both PSRs J2214+3000 and J2234+0944 were found in radio searches in the directions of unassociated Fermi-LAT sources, with the GBT and the Parkes radio telescope respectively (Ransom et al. 2011;Ray et al. 2012). PSR J2214+3000 was discovered to be a BWP system, with a low-mass companion (∼0.014 M ) and an orbital period of about 9.9 hr and a spin period of 3.12 ms. The companion star of PSR J2214+3000 was detected in the optical by Schroeder & Halpern (2014). PSR J2234+0944 has an orbital period of about 10.1 hr, with a companion star of ∼0.015 M , and a spin period of ∼3.36 ms (Ray et al. 2012). Phase resolved γ-ray pulsations have been detected for all three pulsars1 (Abdo et al. 2013).

In this paper we present updated timing solutions for the three BWPs J0023+0923, J2214+3000, and J2234+0944. In Section 2 we present the observations of the three BWPs, in Section 3 we present our timing solutions and in Section 4 we discuss the stability of these pulsars and compare their properties with other BWPs. The data used in this paper is available online2 .

OBSERVATIONS

The data in this paper consists of times of arrival (TOAs), collected and observed in the context of the European Pulsar Timing Array (EPTA). The EPTA is a collaboration of scientists, the observational part of which is supported by the data collected by five different telescopes around Europe: the Westerbork Synthesis Radio Telescope (WSRT) in the Netherlands, the Effelsberg 100-meter Radio Telescope (EFF) in Germany, the Lovell Radio Telescope (JBO) in the UK, the Nançay Radio Telescope (NRT) in France, and the Sardinia Radio Telescope (SRT) in Italy (Desvignes et al. 2016). As the SRT has only just started observing on a regular basis, data from this telescope are not used in this paper. The three data sets for PSR J0023+0923, J2214+3000 and J2234+0944 span about 7.5, 8 and 7 years of data, respectively, see Fig. 1. As described in Table 1 the data used are from several pulsar recording instruments (Karuppusamy et al. 2008;Desvignes et al. 2011;20092010201120122013201420152016 et al. 2016;Guillemot et al. 2016;Bassa et al. 2016).

Lazarus

The data were collected using dual-polarization receivers, and backends providing the capabilty of coherent dedispersion (except for data from the Lovell digital filterbank (DFB) backend, which employs incoherent dedispersion), thus recording four coherently dedispersed timeseries (one for each Stokes parameter), integrated modulo the spin period (i.e., 'online folding') into subintegrations of between 10-60 seconds, depending on the telescope and date of observation (see Table 1).

Specific parts of the data that were affected by radio frequency interference (RFI) and not detected or cleaned in the preprocessing stages were cleared using standard PSRCHIVE tools (van Straten et al. 2012), where it should be noted, that the Lovell ROACH backend uses a real-time spectral kurtosis (Nita & Gary 2010) RFI removal method as implemented in dspsr in the pre-processing of the data. The cleaned data were then averaged in frequency and time to maximise the signal-to-noise ratio (S/N). For PSR J0023+0923, some observations were taken with a total observing length of several hours, to cover the full orbital period. We split observations that were more than 1 hr long into smaller subsets of 30 min each, to avoid smearing as a result of a non-optimal timing solution. Due to its short orbital period, data for PSR J0023+0923 were averaged in time several times in an iterative process where in each iteration the timing model was improved. This was needed to optimally adjust the orbital delays when time-averaging the data.

Standard templates were constructed for each pulsar by fitting von Mises functions to the total-intensity profiles formed by adding a large number of the averaged observations from the previous step. The TOA for each averaged observation was measured by cross correlation with the respective standard template in the Fourier domain (Taylor 1992) and uncertainties of the TOAs were estimated using the Fourier domain Markov chain (FDM) Monte Carlo algorithm of the pat tool (van Straten et al. 2012).

Since the interstellar medium (ISM) is not only magnetised and ionised but also turbulent and inhomogeneous, this modulates the observed flux density as a function of time and observing frequency (Rickett 1969). This effect, known as interstellar scintillation, leads to the pulsar emission being observed across the wide bandwidths of the receivers employed at varying intensities. This is particularly significant for the NUPPI backend at NRT where data is recorded with 512 MHz of bandwidth around a centre frequency of 1484 MHz, leading to the 'patchy' dynamic spectra shown in Fig. 2. Integration along the frequency axis leads to TOAs which are referred to the centre of the band rather than the frequency at which the maximum intensity occurs. Standard profiles which do not account for this frequency-dependent feature therefore lead to biased estimates of the TOAs. To account for scintillation, the NUPPI data for PSRs J0023+0923 and J2214+3000 were split in two bands and into four bands for PSR J2234+0944.

The TOAs were measured at each observatory with reference to a local clock, which was synchronised using GPS to the Terrestrial Time standard, that is derived from the "Temps Atomique International" time standard, TT(TAI). The TOAs were corrected to the Solar-System barycentre using the NASA-JPL DE421 planetary ephemerides (Folkner et al. 2009). The Tempo2 software package was used (Hobbs et al. 2006) to fit the timing solutions presented in Table 2. Data sets generated by different telescopes were phase-aligned by using constant offsets (JUMPS; following the procedure described in Verbiest et al. 2016) and a few clock offsets were applied where necessary (see also Desvignes et al. 2016 andLazarus et al. 2016) 3 .

Error-scaling factors (called T2EFACs) were calculated for each receiver-frequency-band combination using the timing model derived and taking the square root of the reduced χ 2 from fits to the TOAs for each such combination alone. The T2EFAC value corrects the statistical weights of the data by re-scaling the error bars to account for probable under estimation at the template matching stage.

RESULTS

We obtained stable timing solutions for the three pulsars J0023+0923, J2214+3000 and J2234+0944 over time spans of 7-8 years, using data from four of the EPTA telescopes. The best-fit timing solutions for these three pulsars are shown in Table 2, and in Fig. 3 the timing residuals for the entire time-span are shown. Table 1. Telescope (backends), centre frequency, bandwidth, number of TOAs for the different pulsars, range of the error-scaling factor (T2EFAC) and the subintegration length. The four telescopes used are part of the EPTA and they are the Westerbork telescope (WSRT), Effelsberg, Nançay, and Lovell at Jodrell Bank. * The data from Nançay, NUPPI, are affected by scintillation, and the frequency bands have thus been split into smaller bands, of 256 MHz for PSR J0023+0923 and J2214+3000. For PSR J2234+0944 we have split the band in four, so we have a bandwidth of 128 MHz. The given T2EFAC range is for the range used for all three pulsars. Furthermore, the sub-integration length is given in the last coloumn. Timing stability of black widow pulsars 5

New parameter measurements

For low eccentricities the location of periastron is not well measured in the TOAs and a strong covariance between the longitude of periastron (ω) and the epoch of periastron passage (T0) occurs. This leads to very large uncertainties in the estimate of χ 2 of these parameters (see Lange et al. 2001). This is avoided by using the ELL1 model in tempo2.

When using the ELL1 model, the Laplace parameters, EPS1=e sin(ω) and EPS2=e cos(ω), instead of the eccentricity, e, and ω, are included (Hobbs et al. 2006). Furthermore TASC is included in the ELL1 model rather than T0. We measure EPS1 and EPS2, which are orthogonal components of the eccentricity. The eccentricity is given by e= EPS1 2 + EPS2 2 and the eccentricities of the pulsars are thus nominally e=2.1(25)×10 -5 , e=1.5(9)×10 -5 , e=9.6(33)×10 -6 respectively for PSRs J0023+0923, J2214+3000 and J2234+0944. Given the low significance, the eccentricities are interpreted as upper limits.

In this paper we also introduce and measure a firstorder derivative of the dispersion measure ( Ḋ M = dDM/dt) for PSR J2214+3000, of 3.6(4)×10 -4 cm -3 pc yr -1 . This accounts for the DM variations we show in Fig. 4 (for more details see Section 3.2). We further find a very marginal detection of the first-order derivative of the semi-major axis, of -8.4(58)×10 -15 for PSR J2214+3000. The semi-major axis marginal upper limit measurement of PSR J2214+3000 is not used in the final timing solution presented in Table 2. In addition we derived a significant value of Ṗb for PSR J0023+0923, of 4.7(23)×10 -12 , from the measured orbital frecuency (FB0) and first order orbital frequency derivative (FB1). Ṗb is slightly larger than previous measurement (Arzoumanian et al. 2018), however, this measurement is still within the error on the parameter found in this paper. In general, our results are consistent to the results found in Arzoumanian et al. (2018).

Dispersion Measure variations

The turbulent nature of the ionised ISM, which causes the scintillation described in Section 2, along with significant transverse velocities characteristic to many pulsars implies that the DM along the line of sight can vary significantly over a range of timescales. The effects of changes in the DM on the arrival times of the pulses is inversely proportional to the observing frequency squared, and are therefore more easily detectable at lower frequencies. DM variations are only detectable when the low-frequency data are included, as the DM uncertainty in the high-frequency data alone is too large to provide any sensitivity to DM variations. If we look at the combined low-and high-frequency WSRT data of PSR J2214+3000 (see Fig. 4), we detect DM variations. The DM measurements in Fig. 4 were obtained by fitting the DM to non-overlapping segments of data that were around 100-200 days long, to gain sufficient frequency coverage. In a separate analysis of the low-frequency WSRT data, we detected a first-order DM derivative of 3.6(4)×10 -4 cm -3 pc yr -1 . For PSRs J0023+0923 and J2234+0944 we did not detect any DM variations. We furthermore tested whether the orbital phase was dependent on DM, see Section 4.2.

DISCUSSION

In the following section we will discuss some of the features that the three sources show and relate these to the long-term stability of the three BWPs. Furthermore we will compare the properties of stable and unstable BWPs.

Orbital Stability

In this paper, we classify a BW pulsar as stable if it has a single phase-connected timing solution, with only one orbital period derivative of fixed sign (see also Freire et al. (2017)). However 'instability' counter to this definition would not disqualify a BW pulsar as a PTA source. A more applicable definition for that case is provided in Bochenek et al. (2015) and we do not elaborate on that any further in this article. We clarify that the instability as defined here is only linked to variability in the orbital period and not a red-noise process driven instability.

By this definition, the first BWP found, PSR B1957+20, is unstable over long time spans, since it was not possible to create a timing solution spanning more than a few months to about a year (Arzoumanian et al. 1994). The pulsar also shows radio eclipses due to the bloated companion and material surrounding the companion. The second BWP found, PSR J2051-0827, shows similar behaviour, including timing solutions that were only stable for short periods of time, of about 3 yrs (Lazaridis et al. 2011;Shaifullah et al. 2016). Contrary to the first two BWPs discovered, there is one BWP system PSR J0610-2100, which shows a stable timing solution and is included in the EPTA and IPTA (Desvignes et al. 2016;Verbiest et al. 2016).

The three BWPs analysed in this paper are very different from those early discoveries. None of these three BWPs show radio eclipses and it is possible to find a coherent timing solution for each of the three systems that is valid at least over a 7-8 year period. However, PSR J0023+0923 does show significant orbital period derivatives (see table 2), and all orbital frequency (FB) parameters are non-zero, which seems to indicate that there are orbital variations present similar to what is found for e.g. PSR J2051-0827 (Shaifullah et al. 2016). However, all three pulsars have reasonably low timing residual rms (∼ 2.1 -4.3µs) and all are good PTA sources.

Eclipses

We examined our data sets for the three pulsars for eclipses and did not find any evidence for radio eclipses for any of the three pulsars at any of our observed frequencies. See the right hand panels on Fig. 3, which presents averages over many orbits. These plots indicate that there is not enough material around the companion to create the eclipses, as is also suggested for PSR J0023+0923 from X-ray observations by Gentile et al. (2014). The lack of radio eclipses can be an indication that the companion stars are not filling their Roche lobes (see Section 4.3). To further test if we see any indications of eclipses, we studied the DM as a function of orbital period. If there was any indication of additional ionized material along the line of sight at around orbital phase 0.25 (i.e. when the companion is at the minimum distance from us and the MSPs is on the opposite part of the orbit), then it would be evident when looking at the DM as a function of the orbital phase (Stappers et al. 1996a). In Fig 5, we do not see any significant DM variation around orbital phase 0.25; however, we note that the orbital period is poorly sampled for all three pulsars in this paper. The poorly sampled orbital periods are the cause of the outliers in Fig. 5. Other evidence for eclipses not being present is seen in the mass function of the three black widows in this paper (see Table 2). As mentioned by both Freire (2005) and Guillemot et al. (2019) there seems to be a correlation between eclipses, or the lack thereof, and the mass function of the pulsars. The higher values of mass functions in general seem to belong to pulsars showing eclipses and the low mass functions, as we also see for the three pulsars in this paper, are often coinciding with BWs without eclipses (see table 2 in Guillemot et al. ( 2019)).

Roche-lobe filling factors

The companion stars for the first-known two BWPs, PSRs B1957+20 and J2051-0827, are suggested to be Roche-lobe filling, with a Roche-lobe filling factor (RLFF) of about 0.9 (the filling factors are determined from fitting the optical light curve of the companion stars) (Stappers et al. 2001;Reynolds et al. 2007;van Kerkwijk et al. 2011). This means that the outer layers of the companion star are less

Residuals [ms]

WSRT (345MHz) WSRT (1380MHz) Effelsberg ( 1347MHz) The timing residuals of PSRs J0023+0923, J2234+0944 and J2214+3000 after fitting the timing models presented in Table 2. The residuals show stability over the entire time span. Right panels: The same residuals as a function of orbital phase measured from the ascending node. There is no evident systematic variation of the residuals as a function of orbital phase due to delays in the companion atmospheres. This indicates that there are no eclipses, whose presence would be indicated by a clear excess delay around an orbital phase of 0.25.

Effelsberg (2627MHz) NRT, NUPPI (1484MHz) NRT, BON (1408MHz) Lovell, dfb (1400MHz) Lovell, roach (1532MHz)
gravitationally bound to the star and it is thus easier to expel material from, or ablate, the companion star. This is one of the essential details of the geometry of the BWPs, as it is the ablated material that causes the radio eclipses (van Kerkwijk et al. 2011). Furthermore the BWP J1810+1744 similarly has a RLFF of about 0.8. This system also shows radio eclipses and instabilities (Breton et al. 2013;Polzin et al. 2018).

The BWP J2256-1024 was found by Breton et al. (2013) to have a relatively low RLFF of ∼0.4, but was reported by Gentile et al. (2014) to show orbital variability. Recent studies of another BWP, PSR J0636+5128 (Kaplan et al. 2018), show that this pulsar has a large RLFF, of about 0.75, but it does not show eclipses. The missing eclipses are probably due to the relatively low inclination angle of ∼24 • . This system seems to be a good candidate for orbital variations (Bak Nielsen in prep.).

Of the three pulsars discussed in this paper, the RLFF has only been measured for the companion star of PSR J0023+0923, and it is ∼0.3 (Breton et al. 2013). If the Roche lobe of a BWP system is not filled, eclipses would not be expected, which coincides with what is observed for PSRs J0023+0924, J2214+3000 and J2234+0944. However, we do find that orbital variations are present for PSR J0023+0923. Thus it is evident that while RLFFs seem to be closely related to the inherent stability of a binary system, RLFF alone is not a good discriminator of timing solution stability.

Additional intrinsic or observational parameters at play

BWPs are still a rather heterogeneous group of objects and various factors impact on their observational properties as well as on the intrinsic process at play in these systems.

Alongside the RLFF, the irradiation of the companion star by the pulsar and the distance between the pulsar and companion star are factors that have the potential to affect the orbital stability. For example the irradiation affects how heated and ablated the companion star is and thus how much it will potentially fill its Roche lobe (Stappers et al. 2001;van Kerkwijk et al. 2011;Breton et al. 2013). The distance to the BWPs affects the measured value of the RLFF, because it is determined by lightcurve modelling and is thus dependent on the apparent luminosity, and could thus affect our discussion in Section 4.3 (Breton et al. 2013). The orbital inclination of the system has no effect on the orbital stability, however, it does influence whether or not we see eclipses in these systems. The inclinations of PSRs J0023+0923, J2051-0827 and B1957+20 are similar, respectively i ∼58 • , i ∼40 • and i ∼65 • , inferred from light curve modelling. The first one of these systems does not show eclipses and the other two do (Breton et al. 2013;Stappers et al. 2001;van Kerkwijk et al. 2011). However, one system where it is likely that the inclination angle does matter is J0636+5128, which has an inclination of ∼24 • , and shows no eclipses, even though it has a similar RLFF as PSRs J2051-0827 and B1957+20 (Kaplan et al. 2018).

PSR J2234+0944 AS A POSSIBLE TRANSITIONAL SYSTEM

Transitional millisecond pulsars are a relatively new class of binary pulsars, that show changes between emitting in radio and X-rays, and are believed to switch states between having an active accretion disc and no accretion disc (Archibald et al. 2009). Since the first transitional millisecond pulsar, PSR J1023+0038, was discovered to switch states, two more have been found, and all three systems are RBPs (Papitto et al. 2013;Bassa et al. 2014;Roy et al. 2015). The companion stars of BWPs and RBPs are very different, as RBPs are thought to have non-degenerate companions with a higher mass than the semi-degenerate companions of BWPs (Chen et al. 2013). However, variations in their orbit and the presence of radio eclipses in some of these systems are characteristics shared between the BWPs and the transitional systems (Fruchter et al. 1988;Stappers et al. 1996b;Archibald et al. 2009;Papitto et al. 2013;Chen et al. 2013;Roy et al. 2015). Torres et al. (2017) suggested that PSR J2234+0944 could be a potential transitional system, based on variations in the γ-ray spectrum, which is similar to what is seen in the transitional millisecond pulsar J1227-4853. However, PSR J2234+0944 has a timing solution that is stable over 8 yrs, with no obvious irregularities or eclipses. Given the data, this system appears to be different from the variable transitional systems identified so far.

We studied three BWPs: PSRs J0023+0923, J2234+0944, and J2214+3000, all of which show different behaviour than the first two BWPs discovered, PSRs B1957+20 and J2051-0827, that are counted as unstable systems. We present timing solutions for PSRs J0023+0923, J2214+3000 and J2234+0944 over a data span of 7-8 yrs. PSR J0023+0923 shows orbital variations, which can be modelled with up to a third order orbital frequency derivative (see Table 2). That seems to counter the argument that a low value of the Roche Lobe Filling Factor is an indicator for stable BWPs. From the orbital frequency derivatives we derived an orbital period derivative for PSR J0023+0923, of 3.9(24)×10 -12 The other two BWPs in this study, PSRs J2214+3000 and J2234+0944, show no orbital variations. It is possible that this is due to their companion star not filling its Roche lobe. This would make the companion stars more gravitationally tightly bound to themself, resulting in a cleaner system with less material to interfere with, or to affect the radio signal and the orbital parameters. Since the three pulsars show stability, or only minimal instability, over long timescales it is possible to use the systems discussed here as part of a pulsar timing array, and all three pulsars have been added to the EPTA source list. It remains to be seen if the lack of orbital stability in some of these gravitationally tightly bound systems is tied to specific high-energy signatures. Comparing with more BWPs and other pulsar systems will provide a clearer picture of what the most important factors are for determining why a BWP system has a stable or unstable orbital solution.

Figure 2 .

 2 Figure 2. Dynamic spectra obtained with Nançay and the NUPPI backend; the scintillation in time and frequency for PSR J2234+0944 is clearly seen, and it affects the timing because the profile shape changes with frequency. The MJDs represented in this plot are 55879, 55888, 55902, 55914, 56046 and 56060. The dark area corresponds to RFI masked data and the bright areas correspond to high intensity scintels.

Figure 3 .

 3 Figure3. Left panels: The timing residuals of PSRs J0023+0923, J2234+0944 and J2214+3000 after fitting the timing models presented in Table2. The residuals show stability over the entire time span. Right panels: The same residuals as a function of orbital phase measured from the ascending node. There is no evident systematic variation of the residuals as a function of orbital phase due to delays in the companion atmospheres. This indicates that there are no eclipses, whose presence would be indicated by a clear excess delay around an orbital phase of 0.25.

Figure 5 .

 5 Figure 5. Dispersion measure as a function of orbital phase. There appears to be no variation around the orbital phase of 0.25 where an eclipse would be expected. The plot shows measurements from the low frequency WSRT data. The blue points show the orbit sampled in steps of 0.1 P b while the red points are sampled in steps of 0.2 P b . The blue points sample about 3-11 TOAs and the red points sample around 6-20 TOAs. The grey line denotes the average DM value.

  Measured DM variations as a function of time, derived from the combined WSRT low-and high-frequency data of PSR J2214+3000. A clear gradient in DM can be detected during approximately the first three years of our timing campaign. The red line indicates DM(t)=DM+(t-t 0 )× Ḋ M. Ḋ M was measured only by a fit to the low frequency WSRT data. ∆ DM was measured relative to the DM quoted in table 2 and was measured from fits of the DM to 100-200 day long segments.

		2010	2011		2012	2013		2014	2015
		0.0012
		0.0008
		0.0004
	DM [cm 3 pc]	0.0004 0.0000
		0.0008
		0.0012
		0.0016
		55250	55500	55750	56000 Time [MJD] 56250	56500	56750	57000
	Figure 4.

https://confluence.slac.stanford.edu/display/GLAMCOG/ Public+List+of+LAT-Detected+Gamma-Ray+Pulsars

http://www.epta.eu.org/ MNRAS 000,1-10 (2020)

Throughout this work Tempo2 version 2014.2.1 was used.MNRAS 000,1-10 (2020)

MNRAS 000,1-10 (2020)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000,1-10 (2020)

ACKNOWLEDGEMENTS The European Pulsar Timing Array (EPTA: www.epta.eu.org) is a collaboration of European institutes to work towards the direct detection of low-frequency gravitational waves and the running of the Large European Array for Pulsars (LEAP). The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from The Netherlands Foundation for Scientific Research NWO. Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope is supported by a consolidated grant from the STFC in the UK. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS) and with the Université d'Orléans. We acknowledge financial support from the 'Gravitation, Références, Astronomie, Métrologie' (GRAM) national programme of CNRS/INSU, France. Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg. GJ and GS acknowledge support from the Netherlands Organisation for Scientific Research NWO (TOP2.614.001.602). GD gratefully acknowledges support from European Research Council (ERC) Synergy Grant "BlackHoleCam" Grant Agreement Number 610058. The authors also thank Stefan Os lowski for contributions to observations used in this paper.


Mohammed Abutaha
email: mohammad.abu-taha@univ-nantes.fr

Safwan El Assad
email: safwan.elassad@univ-nantes.fr

Mousa Farajallah
email: mousa.farajallah@univ-nantes.fr

Audrey Queudet
email: audrey.queudet@univ-nantes.fr

Olivier Deforge
email: olivier.deforges@insa-rennes.f

Chaos-based Cryptosystems using Dependent Diffusion : An Overview

Keywords: chaos-based cryptosystems, dependent confusion-diffusion, speed performance, security analysis

Based on some important properties of chaos, such as ergodicity, quasi-randomness, and high sensitivity to the secret key, chaos is a hot research field in secured communication, and recently a variety of chaos-based cryptosystems have been proposed for achieving the confidentiality of transmitted images over public channels. Most of chaos-based encryption algorithms are based on the Fridrich structure, which use separate confusiondiffusion layers. In this paper, we give an overview of some chaos-based block cryptosystems, including our cryptosystem, using dependent diffusion, in which the confusion process and diffusion process are performed sequentially on each pixel of the plain plain image. This kind of cryptosystems are more efficient than the traditional confusion-diffusion architecture of Shannon. A comparison study of efficiency in terms of speed performance and of robustness against cryptanalysis is done.

I. INTRODUCTION

Since two decades, many chaos-based image encryption schemes have been proposed. The typical structure of these encryption schemes is based on separate confusion-diffusion layers as that is done in Fridrich architecture [START_REF] Fridrich | Symmetric ciphers based on two-dimensional chaotic maps[END_REF], [START_REF] Li | Improving security of a chaotic encryption approach[END_REF], [START_REF] Chen | A symmetric image encryption scheme based on 3d chaotic cat maps[END_REF]. El Assad et al., [START_REF] Assad | Chaos-based block ciphers: An overview[END_REF] have given an overview of main chaos-based cryptosystems, using Fridrich's structure. Recently, many chaosbased cryptosystems, more efficient in terms of time consuming and of resistance against cryptanalysis, than the previous ones have been investigated [START_REF] Yang | A fast image encryption and authentication scheme based on chaotic maps[END_REF], [START_REF] Wang | A new chaos-based fast image encryption algorithm[END_REF], [START_REF] Zhang | An image encryption scheme using reverse 2-dimensional chaotic map and dependent diffusion[END_REF], [START_REF] Farajallah | Fast and secure chaosbased cryptosystem for images[END_REF], and [START_REF] Farajallah | Chaos based crypto and joint crypto-compression systems for images and videos[END_REF]. These cryptosystems are based on dependent confusion-diffusion layers, generally achieved by substitution-diffusion processes that use dynamic keys supplied by a chaotic sequence. The chaotic sequence is produced by a chaotic generator which is the heart of any chaos-based cryptosystem and so a big part of the efficiency of the system depends on it. Compared to the conventional cryptographic algorithms (3DES, AES), chaos-based cryptosystems have several advantages such as: a very high security level, more flexibility, more modularity, a low power consuming, and easily implemented, which make them more suitable for large scale-data encryption, such as images and videos. In this paper we present an overview of chaos based cryptosystems, including our cryptosystem, using dependent confusion-diffusion. The rest of the paper is organized as follows. In section II, we describe four known Fig. 1: General structure of chaos-based cryptosystems chaos based cryptosystems of the literature. Section III, gives a comparative results of the speed performance and the security analysis of these cryptosystems, before concluding.

II. PRINCIPAL CHAOS-BASED CRYPTOSYSTEMS USING DEPENDENT DIFFUSION OF THE LITERATURE

All chaos-based and non chaos cryptosystems must achieve the confusion and diffusion effects. The confusion effect is measured by how much a change in the secret key affect the ciphered message. The diffusion effect is measured by how much a change in the plain message affect the ciphered message. In the literature there are mainly two types of chaosbased cryptosystems. The structure of the first type, as we can see in figure 1, is composed of two layers: a confusion layer followed by a diffusion layer that work separately. The confusion process is applied rc times on the block (or on the whole image), then the diffusion process is applied rd times on the output of the confusion process, and finally, the two processes are repeated r times. Both layers required image-scanning (for rc = rd = r =1). Most of chaos-based cryptosystems of first type are considered insecure upon chosen/known plain text attacks. El Assad et al., [START_REF] Assad | Chaos-based block ciphers: An overview[END_REF] gave in their paper an overview of main chaos-based cryptosystems of first type.

The structure of the second type of cryptosystems is similar to the structure of the first type of cryptosystems, but the confusion and diffusion processes are performed sequentially on each pixel of the plain block or plain image as shown in figure2. This type of cryptosystems are more efficient, in Fig. 2: Dependent diffusion structure of chaos-based cryptosystems Fig. 3: Yang et al scheme terms of security and speed performance, than the first type of cryptosystems. Indeed, first, the diffusion process at the pixel level is governed by the confusion process, second, a single scan of plain image pixels is needed to perform the confusion and diffusion effects. In the following we will describe the main chaos-based cryptosystems of the second type including our cryptosystem. Yang et al., [START_REF] Yang | A fast image encryption and authentication scheme based on chaotic maps[END_REF], used a permutation operation, as a confusion layer, achieved by a modified standard map to avoid the problem of permutation of the corner pixel (s = 0, t = 0), while using a logistic map as a diffusion layer. In addition a keyed hash function is used to generate a 128bit hash value from both the plain image and the secret hash keys. The hash value plays the role of the key for encryption and decryption while the secret hash keys are used to authenticate the decrypted image. As we can see from the encryption scheme of figure 3, a one bit change in the plain image changing the dynamic keys of the processes of confusion and diffusion and then the encrypted image become completely different from the previous encrypted one. So, the immunity against knownchosen plain text attacks is easily obtained from one round. However for the decryption scheme it is necessary to transmit the secret encryption-decryption key for each new encrypted image. In their paper, Wang et al [START_REF] Wang | A new chaos-based fast image encryption algorithm[END_REF] introduced the idea of mixing the two layers of permutation and diffusion into a single layer of dependent permutation-diffusion. As a result, one image scanning is required instead of two scanning stages, to accelerate the encryption algorithm. The main steps of Wang et.al., cryptosystem [START_REF] Wang | A new chaos-based fast image encryption algorithm[END_REF] are summarized in the following:

1) Division step: The image is divided into a number of blocks (num), each one is 64 pixels.

num = L × P 64 (1)
Where L and P are the height and the width of the image, respectively. 2) First generation step: In this step, the used keys K 0 , K 1 , ...K 15 are generated, and each one is an 8-bit number. 3) Second generation step: 64 pseudorandom values are generated from the spatio-temporal chaotic map [START_REF] Kaneko | Spatiotemporal chaos in one-and two-dimensional coupled map lattices[END_REF]. These pseudorandom values are defined below as Φ(i, j). 4) Diffusion step: The pixel values inside each block are modified based on the following equation:

G t (i, j) = cycl{X, Y } X = [P t (i, j) ⊕ Φ(i, j) + C t (i -1, j)]Mod G Y = LSB 3 (C t (i -1, (j -1)mod 8) ⊕ Φ(i, j)) (2)
Where: G t (i, j) is the number of levels in the gray image (here is 256 levels) and is performed for all pixels in the t th block. P t (i, j) and C t (i, j) are the plain and the ciphered pixels at i and j positions, respectively. 5) Moving block positions: A block at the position t is moved to a new position according to the following equation: 0) is a generated value from the used chaotic map. If the t new is a non-visited place, then the block is moved to this place, otherwise, the t new is incremented until the new value point to a non-visited place. As not all pixels inside a given block are permuted, this step is necessary to increase the security of this model to the statistical analysis attacks.

t new = X(0) × num (3) X(

6) Exchanging lattice values:

The pixels C(7, i) and C(7, (i+d)mod 8) are exchanged, where d is calculated as: [START_REF] Zhang | An image encryption scheme using reverse 2-dimensional chaotic map and dependent diffusion[END_REF] where:

d = LSB 3 (C(7, 0)) (4) 7) First Block pixel exchange: C 0 (0, s) is exchanged with C k l (7,
s = LSB 3 (k l ) k l = [K0⊕K1⊕K2⊕...K15]×num 256
Steps, 3 → 6 are performed for all blocks, and repeated r rounds until the required security level is reached. In Zhang's paper [START_REF] Zhang | An image encryption scheme using reverse 2-dimensional chaotic map and dependent diffusion[END_REF], two cryptosystems were designed based on Fridrich's architecture. The first one consists of a dependent diffusion layer based on the reverse 2-D cat map. The second algorithm presents new mapping from a pseudo-random position to another pseudo-random one for the confusion effect. The diffusion layer in the cryptosystems is based on the logistic map. In these versions, Zhang tried to achieve the confusion and the diffusion effects sequentially. Then, the effect of one ciphered pixel is transferred to the next one and so on. From this idea, only two rounds (in the first version) and one round (in the second version) of the diffusion-confusion process are/is needed instead of many rounds of separated confusion and diffusion processes used in the traditional structures such as Fridrich cryptosystem and other cryptosystems. In the following, our work is directed to the first Zhang cryptosystem. The mathematical model of the first Zhang algorithm is:(Enc=the encryption process).

Enc = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ x y = 1 p i qi p i q i + 1 x y ModN ciph(x, y) = arr(x , y ) ⊕ f (t) t = ciph(x, y) (5)
The ) from the whole image. The values of r x j and r y j are calculated using equation ( 6) and equation [START_REF] Zhang | An image encryption scheme using reverse 2-dimensional chaotic map and dependent diffusion[END_REF], where j is a counter ranging from 0 to n -1 encryption rounds. 2) Array exchanges: the second step is to exchange the first byte arr(0, 0) with the random byte from the previous step arr(r x j , r y j ). 3) Dependent diffusion: then the cryptosystem goes to the dependent diffusion layer for m rounds (m = 2 in the Zhang algorithm case), which also includes three stages. a) New position estimation: in the dependent diffusion layer the first step is to calculate the new byte position (x , y ) from the old byte position (x, y) using equation [START_REF] Farajallah | Fast and secure chaosbased cryptosystem for images[END_REF].

x y = 1 p i qi p i q i + 1 x y ModN ( 8
)
where N is the size of the square test image. p i and q i are calculated using the following equations:

p i = (SQ 1 (2000 + i) × 10 9
)mod 512 (9)

q i = (SQ 2 (2000 + i) × 10 9 )mod 512 (10)
The variable i in the last equations is a counter ranging from 0 to m -1. The two sequences SQ 1

and SQ 2 as can be seen in equations (6, 7, 9, and 10) are calculated using the following equation:

f (x n ) = α × x n-1 (1 -x n-1 ) (11)
With the initial values x -1 =0.12345678912345 for SQ 1 , and x -1 =0.67856746347633 for SQ 2 . The value of α is set to 3.99999. b) Calculation of the local ciphered pixel: the next step of the dependent diffusion layer is to calculate the ciph(x, y) value using the following equation:

ciph(x, y) = arr(x , y ) ⊕ f (t) (12)
where

f (t) = [α( t 1000 ) × [1 - t 1000
] × 1000]mod 256 (13) c) Update of t: the last step of the dependent diffusion layer is to change the value of the t variable using equation [START_REF] Mar | New analysis methods on strict avalanche criterion of s-boxes[END_REF].

t = ciph(x, y) ( 14
)
The initial value t 0 is defined by the following equation:

t 0 = [4×key d ×(1-key d )×1000]mod 256. ( 15
)
Where the initial value of the key d =0.33456434300001. Farajallah et al., 2015 [START_REF] Farajallah | Fast and secure chaosbased cryptosystem for images[END_REF], and Farajallah [START_REF] Farajallah | Chaos based crypto and joint crypto-compression systems for images and videos[END_REF] first they partially cryptanalysis the first algorithm of Zhang (described above) based on the following partial cryptanalysis equation.

arr k k = ciph k ⊕ f (ciph k-1 ) (16)
where k = x × N + y and k = x × N + y (see equation ( 8)), arr k k is the input pixel of the last dependent diffusion round (m) in the last encryption round (n) and ciph k is the ciphered pixel. As the function f (t) is known, then equation ( 16) can be used to remove the diffusion effect of the last (m and n) rounds from the ciphered pixels. This allows recovery of a permuted version of the previous ciphered image. This removal gives the attacker the possibility to: 1) Decrease the dynamic key space of the whole cryptosystem.

2) Perform partial cryptanalysis of the Zhang cryptosystem

for (n = 1, m = 1) and (n = 1, m = 2). 3) Decrease the UACI and NPCR values significantly. Then, based on the previous analysis they proposed an efficient cryptosystem that overcome the weaknesses of Zhang cryptosystem while keeping a very high speed compared to the main chaos-based cryptosystem of the literature. The encryption side of the proposed cryptosystem is given in Figure 5), for the first block. Each pixel from the plain block (p 0 (k)) is XOR-ed with the initial byte (iv(k)) from the initial vector (IV ), then the output is XOR-ed with the discrete logistic map output to carry out the diffusion process. Then, the 8 least significant bits resulting from the diffusion process (LSB 8 (y 0 (k)) are relocated using the modified 2-D cat map to obtain the ciphered pixel at the new position (c 0 (k n )). It is important to note that the input of the discrete logistic map is based on the previous ciphered pixel (since c 0 (k n ) = LSB 8 (y 0 (k)) and the input of the discrete logistic map is 32 bits and the ciphered pixel is 8 bits. That is why the cryptosystem takes (y 0 (k -1)) before the LSB 8 function and not after. For the first encrypted byte, the input of the discrete logistic map is Kd m , and this value is re-initialized every new encryption round. Because the c 0 (k n ) is only a part of the logistic map input, it is impossible to recover y 0 (k -1) from c 0 (k n ) only. The encryption of the next blocks is almost the same. Each pixel from the plain block (p l (k)) is XORed with ciphered byte from the previous block at the same position(i.e., c l-1 (k) to achieve the CBC mode). Then the rest of the operations are the same as in the first encryption block. The modified cat-map is given by the following equation Fig. 5: Encryption structure of the proposed cryptosystem

i n j n = Mod 1 u v 1 + uv i j + ri + rj rj , M M ( 17
)
(17) which is a one-to-one function, which means that each point of the square matrix can be transferred to exactly one unique point. So, instead of exchanging the values at the new position (i n , j n ) with the old one (i, j), we use a transfer operation because of its speed compared to the swap operation that is usually used. The block size b s is M 2 (M is the square root of the block size in our proposed cryptosystem).

The system parameters u, v, ri and rj are in the range of [0, M -1]. The structure of the dynamic keys which are produced by the chaotic generator [START_REF] Assad | Generator of chaotic sequences and corresponding generating system[END_REF] during the permutation process is:

Kp = Kp 0 Kp 1 Kp 2 • • • Kp r-1 Kp m = u m v m ri m rj m (18)
The modulo operation of equation ( 17) makes it a noninvertible equation. But it is still a reversible one. Thus, in the decryption part of the proposed cryptosystem, the reverse layer is also achieved by equation ( 17).

The Logistic map is a non-linear chaotic discrete function that produces random sequences. In the proposed cryptosystem, the logistic map is used as a diffusion function to achieve the diffusion effect, by transferring the effect from one byte in the block to other bytes in the same block. This structure makes the proposed cryptosystem highly sensitive to the plaintext. The mathematical model of the discrete logistic map is:

X k+1 = ⎧ ⎨ ⎩ X k ×(2 N -X k ) 2 N -2 ifX k = 3 × 2 N -2 , 2 N 2 N -1 ifX k = 3 × 2 N -2 , 2 N (19)
where X k+1 is the new value calculated from the previous one X k . N is the number of bits representing the integer output of the discrete logistic map, which is equal to 32 bits. From the figure 5, we can write the encryption mathematical model as:

c l (k n ) = LSB 8 [y l (k)] (20)
y l (k) = p l (k) ⊕ s l-1 (k) ⊕ f (y l (k -1)) (21)
where y l (k) is a 32-bit variable, p l (k), s l-1 (k) are 8-bit variables and f is the logistic map. The following remarks should be considered: 1) During the encryption, equation (21) should be evaluated before equation (20), for each byte of a block and for all blocks.

2) The input of the logistic map for k = 0 is kd m when l = 0 and it is y l-1 (b s -1) for l > 0. 3) For k > 0 and for all l, the input of the logistic map is the result of equation ( 21) and not the previous output (see equation 19). Note that: k = i × M + j k n = i n × M + j n i n and j n are calculated using equation ( 17). The sequence s l-1 (k) is given by the following equation:

s l-1 (k) = iv(k) if l = 0 c l-1 (k) if l > 0 (22)
where

l = 0, 1, 2, ...b n -1, k = 0, 1, 2, ...b s -1, IV = {iv(0), iv(1), iv(2), ..., iv(b s -1)}, b s is block size in bytes b n = image size block size = L×C×P bs
, is the number of blocks. with, L, C, and P are the number of lines, the number of columns, and the number of planes of the image respectively.

III. TIME PERFORMANCE AND SECURITY ANALYSIS

A. Time performance

The time performance is determined by evaluating the running speed that can be measured by: the average encryption/decryption times, the average encryption throughput, and the average needed number of cycles to encrypt one byte. The encryption throughput (ET) and the number of cycles, which are required to encrypt (or decrypt) one byte, are defined as:

ET = Image Size (Byte) Encryption T ime (second) (23)
Number of cycles per Byte = CPU Speed (Hertz)

ET (Byte) (24)
The last equation permits to compare the running speed of different cryptosystems working on different platforms. In table I, we give the average encryption/decryption times for Lena image at different sizes with different cryptosystems(for our cryptosystem the used block size is 1024 byte). In table II we give the average throughput and the average number of cycles. The proposed cryptosystem is at least 3 times faster that zhang and wang cryptosytems.

B. security analysis

usually A cryptanalyst tries to break the cipher without knowing the secret key, and this with several levels of difficulties based on the available resources. Chosen plain-text attack is the easiest one for the attacker: The attacker has access to the system without knowing the secret keys. Then, he has the possibility to choose a set of plain-text messages and to encrypt them. If a cryptosystem can resist to chosen plain-text attack, then it can resist to all other attacks such as:cipher text only, known plain-text and chosen cipher attacks. The chosen plain-text attack can be realised by the plaintext sensitivity attack or differential attacks introduced by Eli Biham and Adi Shamir [START_REF] Biham | Differential cryptanalysis of des-like cryptosystems[END_REF]. Most researchers use two security parameters to measure the resistance of any chaosbased cryptosystem for plain-text sensitivity attacks. These parameters are: the Number of Pixel Change Rate (NPCR) and the Unified Average Changing Intensity (UACI); they are given by the following equations, respectively:

NP CR = 1 L × C × P × P p=1 L i=1 C j=1 D(i, j, p) × 100% (25)
where

D(i, j, p) = 0, if C 1 (i, j, p) = C 2 (i, j, p) 1, if C 1 (i, j, p) = C 2 (i, j, p) (26) UACI = 1 L × C × P × 255 × P p=1 L i=1 C j=1 |C 1 (i, j, p) -C 2 (i, j, p)| × 100% (27)
In the previous equations, i, j and p are the row, column, and plane indexes of the image, respectively. L, C and P are, respectively, the length, width, and plane sizes of the image. The optimal NPCR value is 99.61%, and the optimal UACI value is 33.46% [START_REF] Wu | Npcr and uaci randomness tests for image encryption[END_REF]. The previous metrics are necessary but not sufficient to ensure that the proposed cryptosystem is resistant against plain-text sensitivity attacks. Then, a new metric measurement, the Hamming Distance (HD) is used to quantify the avalanche effect. The Avalanche effect is achieved for any block cipher, when a small change (for example, flipping a single bit) in either the plain-text or the secret key, causes a drastic change in the cipher-text (e.g., half of the output bits are flipped), [START_REF] Mar | New analysis methods on strict avalanche criterion of s-boxes[END_REF]. Therefore, this evaluation test is used to measure the resistance of any cryptosystem to the plain-text and the key sensitivity attacks. The HD is defined by:

HD(C 1 , C 2 ) = 1 |Ib| |Ib| K=1 (C 1 (K) ⊕ C 2 (K)) (28)
where |Ib| = L × C × P × 8, is the size of the image in bits.

The optimum HD value is 50%. A good block cipher should produce an HD close to 50% (probability of bit changes, which means that a one bit difference in plain-image will make every bit of the corresponding cipher-image change with a probability of a half. In table III, we show the obtained comparative results of NPCR, UACI and HD for Boat bmb image of size 256 × 256 × 3.

Statistical attack

To resist the statistical attacks, the histogram of the ciphered image must be uniform and the correlation coefficient should be close to zero . In figure 6a, 6b we show the plain image for boat and its ciphered one, and in the figure 6c, 6d we give the histograms for the plain/cipher images. As we can see the histogram of the ciphered image seems to be uniform, to be sure we made the Chi square test with the following parameters: alpha=0.05, and number of classes equal to 256. The obtained experimental value is 255.12 which is less than the theoretical one that equal 293. This means that the histogram is uniform. Figure 7 shows the correlation curves of the adjacent pixels in the horizontal direction for the plain image and its ciphered one. The values of their corresponding correlation coefficient are 0.96606 and 0.0085. All these results imply a high security level of the proposed cryptosystem.

IV. CONCLUSION

In this paper, we presented an overview of chaos-based cryptosystems, including our cryptosystem, that are based on dependent confusion-diffusion processes, and we compared their speed performance and their robustness against cryptanalysis. The obtained results show that, our cryptosystem is faster than the others while having a very high security level.

Fig. 4 :

 4 Fig. 4: Zhang image encryption cryptosystem architecture

  r x j = (SQ 1 (2000 + 100 + j) × 10 9 ) mod (512) (6) r y j = (SQ 2 (2000 + 100 + j) × 10 9) mod (512)[START_REF] Zhang | An image encryption scheme using reverse 2-dimensional chaotic map and dependent diffusion[END_REF]

Fig. 6 :

 6 Fig. 6: Histogram of the boat plain image and its ciphered image

Fig. 7 :

 7 Fig. 7: correlation of the boat plain image and its ciphered image

TABLE I :

 I Average encryption/decryption times for different cryptosystems for Lena image with different sizes.

	Cryptosystem	256 × 256 × 3	512 × 512 × 3	1024 × 1024 × 3
	Proposed Cryptosystem 2.04/2.68	8.08/10.57	31.85/41.83
	Zhang	7.5/7.5	3 0 /30	120/120
	Wang	7.79/8.39	31.16/33.54	124.64/134.16

TABLE II :

 II Average Encryption throughput and average Number of cycles for to encrypt/decrypt one byte.

	Cryptosystem	ET in MBps	Number of cycles per byte
	Proposed Cryptosystem 92.975/70.879	29.87/40.57
	Zhang	25/25	122.07/122.07
	Wang	24.06/22.35	122.85/132.24

TABLE III :

 III The NPCR, UACI and HD

	Cryptosystem	NPCR	UACI	HD
	Proposed Cryptosystem 99.615	33.463	0.500030
	Zhang	99.633	33.34	-
	wang	99.96	33.33	-

ACKNOWLEDGEMENT

This work is supported by the European Celtic-Plus project 4KREPROSYS -4K ultraHD TV wireless REmote PROduction SYStems, 2015 -