11institutetext: INAF-Osservatorio Astronomico di Palermo, Piazza del
Parlamento 1, 90134 Palermo, Italy
11email<EMAIL_ADDRESS>22institutetext: Signal Processing Group
Dpt. of Signal Theory and Communications, Universidad Carlos III de Madrid,
Avda. de la Universidad 30 28911 Leganés, Madrid, Spain

# X-ray variability of HD 189733 across eight years of XMM-Newton observations

I. Pillitteri 11 G. Micela 11 A. Maggio 11 S. Sciortino 11 J. Lopez-Santiago
22

(Received ; )

The characterization of exoplanets, their formation, evolution, and chemical
changes is tightly linked to our knowledge of their host stars. In particular,
stellar X-rays and UV emission have a strong impact on the dynamical and
chemical evolution of planetary atmospheres. We analyzed 25 XMM-Newton
observations encompassing about eight years and totaling about 958 ks in order
to study the X-ray emission of HD 189733 A. We find that the corona of HD
189733 A has an average temperature of 0.4 keV and it is only during flares
that the mean temperature increases to 0.9 keV. Apart from the flares, there
is no significant change in the flux and hardness of the coronal emission on a
timescale of several months to years. Thus, we conclude that there is no
detectable activity cycle on such timescales. We identified the flares and
built their energy distribution. The number of flares observed around the
phases of the planetary eclipses is not statistically different from the
number of flares during transit phases. However, we do find a hint of a
difference in the flare-energy distributions, as the flares observed around
the planetary eclipses tend to be more energetic than the flares observed
around the primary transits of the planet. We modeled the distribution of the
number of flares per day with a power law, showing that it is steeper than the
one observed in the Sun and in other Main Sequence stars. The steepness hints
at a significant fraction of undetected micro-flares. Altogether, the plasma
temperatures below 1 keV observed during the flares, along with the slightly
larger fraction of energetic flares seen at the secondary transits highlight
the peculiarity of the corona of HD 189733 A and points to star-planet
interaction as the plausible origin of part of its X-ray emission. However,
more observational and modeling efforts are required to confirm or disprove
this scenario.

###### Key Words.:

stars:HD189733 – stars:flare – stars:activity – stars: planetary systems –
X-rays:stars

##  1 Introduction

More than 4500 extra-solar planets are already known thanks to spaceborne
missions such as Kepler (Borucki et al., 2006) and TESS (Ricker et al., 2015)
or ground-based photometric and spectroscopic surveys. However, the fraction
of exoplanets that can be characterized with respect to their masses, density,
and chemical composition is biased towards planets orbiting very close to
their hosts. Transmission spectroscopy to infer the chemical composition of
the planetary atmospheres is more feasible in close-in planets than in
exoplanets in wider orbits. In the future, the next JWST observations (Gardner
et al., 2006) and Ariel (Tinetti et al., 2018) will enable detailed chemical
studies of the atmospheres of a significant fraction of known exoplanets. A
full characterization of these atmospheres requires a multi-wavelength
approach to reveal the high energy flux and space weather conditions where
exoplanets are in orbit.

To properly investigate the conditions for the formation and the evolution of
planetary atmospheres, the effects that originate from close star-planet
separations must be taken into account. These effects are related to the
intense flux originating from the star especially at high energies, namely UV
and X-rays bands, the stellar wind and the magnetized medium where such
planets orbit. Models of photo-chemical evolution (Cecchi-Pestellini et al.,
2009, Locci et al., 2022) and evaporation (Kubyshkina et al., 2018) need to be
constrained by the dose of XUV fluxes to which the planets are exposed. In
this context, a detailed characterization of the stellar coronae, their X-ray
spectra, and their time variability in form of flares and coronal mass
ejections (CMEs) is of crucial importance. The present paper addresses the
X-ray characterization of the corona of HD 189733, which is one of the best-
studied exoplanet host. Our goal is to provide the most updated description of
the X-ray emission of this nearby system.

HD 189733 is a wide binary system composed of a K1 type star and a M star
(Bakos et al., 2006). The stellar companion, HD 189733 B, is an M4 type star
in an 3200 yr period orbit at about 216 AU of projected separation from the
primary. HD 189733A hosts a transiting hot Jupiter in a 2.2 days orbit (Bouchy
et al., 2005). Extensive observations in the UV and optical/IR bands have
enabled the inference that the planet of HD 189733 A is inflated and actively
evaporating (Lecavelier des Etangs et al., 2012; Bourrier et al., 2013, 2020).

X-rays from HD 189733A have been investigated in detail with Chandra and XMM-
Newton (Pillitteri et al., 2010, 2011, 2014; Poppenhaeger & Wolk, 2014;
Pillitteri et al., 2015). HD 189733 A has an X-ray luminosity of order of 1028
erg s-1and a corona with mean temperature around 5 MK, thus brighter and
hotter than the solar corona. HD 189733 B has been detected in Chandra
observations (Poppenhaeger et al., 2013) with a quiescent X-ray luminosity of
$L_{\mathrm{X}}\sim 4.7\times 10^{26}$ erg s-1and a coronal temperature of
about 3.5 MK.

The ages of HD 189733 A and HD 189733 B derived from their levels of X-ray and
chromospheric activity are discrepant with HD 189733 A looking younger ($\leq
1$ Gyr)) than HD 189733 B (few Gyr). Pillitteri et al. (2010) and Poppenhaeger
& Wolk (2014) suggested that the star has been spun up from the exchange of
angular momentum between the planet and its star.

HD 189733 A is the first star for which a planetary transit in X-rays has been
reported based on the stacked light curve of six Chandra observations
(Poppenhaeger et al., 2013). The authors suggested that the radius of HD
189733 Ab in X-rays is larger than in the optical band due to the opaque
envelope of the inflated atmosphere surrounding the planet.

In this paper, we present a comprehensive analysis of the observations of HD
189733 made with XMM-Newton for a total of $\sim$958 ks. In particular, we
study the coronal variability and the occurrence of flares, their duration,
shape, and energy. We compare these characteristics to those of other variable
coronal sources to seek any difference or analogy with flares of single stars
without hot Jupiters. The structure of the paper is as follows. In Section 2,
we present the archival observations and their analysis. In section 3, we
present our results and in Section 4, we discuss the results and present our
conclusions.

##  2 Observations and data analysis

HD 189733 has been observed 25 times with XMM-Newton during primary (phase
$\phi=0$) and secondary ($\phi=0.5$) planetary transits. Table 1 lists the
available observations and their duration. The observations span about eight
years between 2007 and 2015 for a total exposure time of 1012 ks (with science
exposures totaling 958 ks). Figure 1 shows the range of observed phases during
each exposure. About 20 out of 25 observations were taken at the primary
transits of the planet, the remaining five were taken at or just after
secondary transits. About 400 ks accumulated in 16 observations are part of a
large program (P.I. Wheatley) that observed planetary transits of HD 189733A
b. The typical length of the exposures is of order of 35-40 ks, however the
fourth observation was more than 60 ks long. Twenty of the observations were
acquired with the Thin1 filter, while five were obtained with the Medium
filter. The Thin1 filter has a slightly better sensitivity than the Medium
filter at very low energies ($<0.5$ keV), however, it could be more prone to
high background episodes and UV leaks.

The five observations carried with the Medium filter were acquired in
FullWindow mode (field of view $\sim 35\arcmin$), while the remaining 20
observations that used the Thin1 filter were obtained in SmallWindow mode,
where only a fraction of the chips were used to imaging the region around the
target. The faster timing allowed by the SmallWindow mode reduces pile-up
effects, however this is irrelevant in the case of HD 189733 A, since the
source count rate is not exceptionally high (as compared to other very bright
astrophysical sources).

Table 1:  Log of the XMM-Newton observations.

| ObsId | RA_NOM | Dec_NOM | Start_UTC | Duration (ks) | P.I. | Filter | Planetary Phases
---|---|---|---|---|---|---|---|---
1 | 0506070201 | 20:00:43.70 | +22:42:39.0 | 2007-04-17 14:06:31.000 | 54.9 | Wheatley | Thin1 | 0.84–1.13
2 | 0600970201 | 20:00:43.70 | +22:42:39.0 | 2009-05-18 21:15:54.000 | 37.3 | Wolk | Medium | 0.45–0.63
3 | 0672390201 | 20:00:43.70 | +22:42:39.0 | 2011-04-30 23:14:20.000 | 39.1 | Pillitteri | Medium | 0.41–0.6
4 | 0690890201 | 20:00:43.70 | +22:42:39.1 | 2012-05-07 18:24:32.000 | 61.5 | Pillitteri | Medium | 0.45–0.76
5 | 0692290201 | 20:00:43.70 | +22 42 35.8 | 2013-05-09 20:16:00.000 | 39.2 | Wheatley | Thin1 | 0.90–1.10
6 | 0692290301 | 20:00:43.71 | +22 42 35.8 | 2013-11-03 07:54:13.000 | 37.1 | Wheatley | Thin1 | 0.90–1.09
7 | 0692290401 | 20:00:43.71 | +22 42 35.8 | 2013-11-21 00:58:40.000 | 42.0 | Wheatley | Thin1 | 0.89–1.09
8 | 0744980201 | 20:00:43.71 | +22 42 35.3 | 2014-04-05 05:05:20.000 | 48.0 | Wheatley | Thin1 | 0.81–1.05
9 | 0744980301 | 20:00:43.71 | +22 42 35.3 | 2014-05-02 01:22:25.000 | 33.7 | Wheatley | Thin1 | 0.91–1.08
10 | 0744980401 | 20:00:43.71 | +22 42 35.3 | 2014-05-13 01:55:22.000 | 40.0 | Wheatley | Thin1 | 0.88–1.08
11 | 0744980501 | 20:00:43.71 | +22 42 35.3 | 2014-05-15 09:57:00.000 | 32.0 | Wheatley | Thin1 | 0.93–1.09
12 | 0744980601 | 20:00:43.71 | +22 42 35.3 | 2014-05-17 14:21:12.000 | 32.0 | Wheatley | Thin1 | 0.92–1.07
13 | 0744980801 | 20:00:43.71 | +22 42 35.3 | 2014-10-17 16:08:26.000 | 37.0 | Wheatley | Thin1 | 0.87–1.06
14 | 0744980901 | 20:00:43.71 | +22 42 35.3 | 2014-10-19 20:38:36.000 | 34.4 | Wheatley | Thin1 | 0.91–1.1
15 | 0744981001 | 20:00:43.71 | +22 42 35.3 | 2014-10-22 01:39:14.000 | 39.9 | Wheatley | Thin1 | 0.9–1.07
16 | 0744981101 | 20:00:43.71 | +22 42 35.3 | 2014-10-24 06:15:47.000 | 39.0 | Wheatley | Thin1 | 0.9–1.09
17 | 0744981301 | 20:00:43.71 | +22 42 35.3 | 2014-11-08 20:16:34.000 | 34.6 | Wheatley | Thin1 | 0.88–1.08
18 | 0744981201 | 20:00:43.71 | +22 42 35.3 | 2014-11-11 00:37:26.000 | 44.0 | Wheatley | Thin1 | 0.89–1.11
19 | 0744981401 | 20:00:43.71 | +22 42 35.3 | 2014-11-13 06:46:05.000 | 32.0 | Wheatley | Thin1 | 0.91–1.08
20 | 0744980701 | 20:00:43.71 | +22 42 35.3 | 2014-11-15 09:48:00.000 | 39.0 | Wheatley | Thin1 | 0.91–1.06
21 | 0744981501 | 20:00:43.71 | +22 42 35.3 | 2015-04-13 02:37:23.000 | 45.4 | Wheatley | Thin1 | 0.89–1.12
22 | 0744981601 | 20:00:43.71 | +22 42 35.3 | 2015-04-17 12:34:55.000 | 41.0 | Wheatley | Thin1 | 0.88–1.09
23 | 0744981701 | 20:00:43.71 | +22 42 35.3 | 2015-04-19 19:06:26.000 | 38.0 | Wheatley | Thin1 | 0.91–1.09
24 | 0748391401 | 20:00:43.69 | +22:42:39.1 | 2015-04-03 03:30:49.000 | 47.4 | SCHARTEL (PS) | Medium | 0.41–0.64
25 | 0748391501 | 20:00:43.69 | +22:42:39.1 | 2015-04-23 06:17:15.000 | 44.0 | SCHARTEL (PS) | Medium | 0.51–0.69

Note: The observations are listed by increasing ObsID. The range of orbital
phases of the planet are also listed. The last two observations were obtained
with a request of Director Discretionary Time (DDT, P.I. Schartel).

Figure 1:  Planetary phases observed as a function of MJD. Phases of the
primary and secondary transits are marked by horizontal lines for sake of
clarity.

The datasets were downloaded from the XMM-Newton web
archive111ttp://nxsa.esac.esa.int/nxsa-web/#search and reduced with SAS v18.0
in order to obtain tables of events calibrated in sky and detector positions,
arrival times, pulse energies, and quality grades. The selection of events and
their filtering in appropriate energy bands and quality grade has been
performed through the SAS task evselect. The same task allows us to accumulate
light curves and spectra for the EPIC datasets.

###  2.1 Contamination from nearby sources

A circular region of a 30″radius used to extract the events associated with an
X-ray source imaged with XMM-Newton-pn contains about 80% of the total events
emitted by the source. HD 189733 has two nearby X-ray sources within about
12″: the M type star HD 189733 B, and an X-ray source that is likely of extra-
galactic nature and perhaps associated with an AGN (hereafter source C)
(Pillitteri et al., 2010; Poppenhaeger et al., 2013). Figure 2 shows the
positions in one of the MOS1 images. We display here the MOS because allows
for a better resolution over the pn. For the analysis of light curves and
spectra, we used the data from pn because they offer a higher count statistics
over MOS. The contribution of events from these sources to the light curves
and spectra of HD 189733 A must thus be taken into account. For HD 189733 B
and source C we extracted the events and light curves in two 5″regions around
their centroids. Source C has a hard spectrum (Pillitteri et al., 2010, 2011)
and, from an examination of its light curves, it does not seem to vary
significantly in time. These two conditions allow us to minimize the
contamination from source C by simply restricting our analysis to the band
$0.3-2.0$ keV and excising the core of the counts centroid up to 8″. Previous
analyses of the X-ray spectra of HD 189733 A have demonstrated that most of
the coronal emission of HD 189733 A is concentrated at energies below 1 keV
(Pillitteri et al., 2014).

The operative way to minimize the contamination of the light curves of HD
189733 A is thus summarized as follows:

  * •

i) select pn events in a 30″region centered on HD 189733 A and in the band
0.3-2.0 keV;

  * •

ii) exclude the events in two 8″circular regions centered on the positions of
HD 189733 B and source C;

  * •

iii) select the events in a semi-circular region opposite; to and not
encompassing the positions of HD 189733 B and source C. This is a control set
for the variability of HD 189733 B, (cf. Fig. 3);

  * •

iv) select events from a 40″circular region for background subtraction to be
applied to the light curves obtained in ii and iii

Light curves with binning of 600 s obtained from events selected at points i
and ii have been used for the analysis of variability of HD 189733 A and HD
189733 B, respectively. The choice of the binning interval is justified by the
subsequent use of an algorithm to detect variability and flares, which assumes
that in each bin, the counts follow an underlying Gaussian distribution (see
Sect. 3.2.2). A bin of 600 s gives on average about 100–150 counts per bin,
which fulfills the requirement of Gaussianity for all the observations, even
during quiescent phases. When HD 189733 B was quiescent, its flux was very low
and similar to the local background flux. However, in the case of flares from
HD 189733 B, the light curves of HD 189733 A obtained at i+ii also show an
excess due to the flares from HD 189733 B. The flux of HD 189733 B can
overwhelm the flux of HD 189733 A at the position of the centroid of HD 189733
A even when excising the events in a 8″region as in ii. Hence, the flares of
HD 189733 B can be mistakenly attributed to HD 189733 A. To check any flaring
of HD 189733 B we used the light curves obtained in point iii in order to
determine the intervals with flares from HD 189733 B in the light curves
obtained from points i+ii. The light curves of the events in a semicircular
region of 30″not comprising HD 189733 B and source C (point iii) served as a
control sample to be compared to the light curves obtained at points i and ii.
Figure 3 shows both the selection region in i through iii,  along with three
examples where the light curves of the primary are contaminated by flares of
the secondary.

While the background can be variable in our XMM-Newton observations, we did
not remove high background intervals. We preferred to keep uninterrupted time
series within each observation and subtract the scaled background count rates
to obtain the net source count rates. However, we did not consider the
intervals of strong background variability for the subsequent analysis of
flares (see Sect. 3).

Figure 2:  MOS1 image of HD 189733 during observation 0744981201. The
brightest source is HD 189733 A, followed by HD 189733 B and source C (likely
a background AGN). The circles indicate regions of size 5″. The large circle
has radius 30″. HD 189733 B was flaring during this observation while during
quiescence was barely detected.

Figure 3:  Discriminating among different sources of variability. Left panel:
pn image and selection of events of HD 189733 A in a semicircular region
during observation 0744981201 in the band 0.3-2.0 keV. The circles indicate
regions of size 5″. The large circle has radius 30″. Right panel: Three
examples of light curves where flaring of HD 189733 A or HD 189733 B is
detected. The gray light curves are from events in a 30″region minus two
8″regions around HD 189733 B and source C (selection i + ii ). The time bins
are 600 s each. The black light curves are accumulated from the events in the
semicircular region opposite to HD 189733 B and source C as the figure in the
left panel. Green light curves are made from the events of HD 189733 B in two
5″regions. The background light curves scaled to the area of selection i + ii
are also shown. In the top panel, a flare is visible both in the gray and
green curves, but not in the black curve; this allows us to conclude that the
flaring source was HD 189733 B and contaminated the selection in i + ii but
not the selection in the semicircular region on the opposite side. In the
middle panel, a flare is detected both in gray and black curves, this flare is
genuinely due to HD 189733 A. In the bottom panel three flares are detected,
the first one, seen in gray and black curves, attributed to HD 189733 A; while
the second one is visible in gray and green curves and attributed to HD 189733
B. Another small flare of HD 189733 A is visible after the flare of HD 189733
B.

##  3 Results

Here, we introduce the variability analysis of HD 189733 B and the variability
of HD 189733 A.

###  3.1 Flares of the M star HD 189733 B

As mentioned in the introduction, HD 189733 B is a M4 star with a relatively
low level of quiescent emission, We identified the flares by eye as the main
episodes of variability, characterized by a steep rise followed by a slower
decay. We thus identified four bright flares from HD 189733 B in the 25 XMM-
Newton observations. In the remaining 21 observations, HD 189733 B was
quiescent and its flux was barely in excess with respect to the background
level.

Figure 4 shows the four pn light curves during which HD 189733 B was flaring.
The flare in ObsId 0744980201 is heavily affected by background variability,
so we do not discuss it in detail here. Its duration was about 5 ks with a
rise time of about $900-1200$ s, the peak intensity was about a factor of four
with respect to the quiescent rate, and it is the faintest among the flares
detected in HD 189733 B.

During observation 0744981201, we observed the brightest flare from HD 189733
B, with a ratio of the peak-to- quiescent rate equal to $\sim 30$. The rise
phase of the flare was about $1200\pm 300$ s. The decay phase is characterized
by an initial steep decrease of duration of about 8000 s (e-folding time $\sim
3500$ s), followed by a much flatter decay that covers the remaining exposure
time (about 16 ks). The star did not return to its initial quiescent rate,
rather it kept emitting a flux about 3.5 times higher than the quiescent flux
before the flare. Small impulses are visible during this post-flare phase –
this flickering perhaps occurred in the same coronal region after the main
ignition resulting thus in an overall flux higher than before the flare.

Flares detected in ObsIds 0692290301 and 0744980701 are very similar, with
rise times $t_{r}\sim 1200-1500$ s and decay encompassing about 5000 s. The
flare in ObsId 0692290301 has been cited in Bourrier et al. (2020). These two
flares were less bright than the flare in 0744981201, reaching a peak rate of
28 and 13 times the quiescent rate, respectively. All the three flares seem to
show a secondary peak after the initial fast decay hinting a second ignition
or a triggered ignition in nearby coronal loops after the first event.

The spectra of the quiescent phases are well described with a thermal
component at $0.3\pm 0.1$ keV with a flux of $4-6\times 10^{-14}$ erg s-1 cm-2
in $0.3-2.0$ keV band corresponding to a luminosity of $2-3\times 10^{27}$ erg
s-1. These values are a factor of four to five higher than the X-ray
luminosity derived by Poppenhaeger et al. (2013) with Chandra. During the
flares a second thermal component is found at $1.3-1.5$ keV in addition to the
cool one. The flux rises by a factor of 10 or more and the luminosity reaches
$\sim 3\times 10^{28}$ erg s-1. Using some diagnostics based on peak
temperatures and decay time (Serio et al., 1991; Reale, 2007), we can infer
flaring loop sizes of order of $5\times 10^{10}$ cm, which is comparable to
the radius of the star. The shape and temperatures of these flares resemble
those observed in the Sun and other solar analogs that are due to mechanisms
of coronal heating and are related to magnetic reconnection. The energy
released by these flares is of order of $1.2-1.4\times 10^{32}$ erg and the
cumulative energy is more than $4\times 10^{32}$ erg released on a total
exposure time of 958 ks. Two of these flares were observed four days apart
(ObsId 0744980701 and 0744981201). We speculate that they could have been
originated in the same region – one that is more active than the rest of the
corona.

Figure 4:  pn light curves of HD 189733 B with flares. The time bins are 600
s. The light curves of the star (black curve) and the scaled background (gray
curve) are shown.

###  3.2 Emission and variability of HD 189733 A

####  3.2.1 Hardness ratio and mean coronal temperature of HD 189733 A

We constructed the light curves of HD 189733 A in 0.3-1.0 keV and 1.0-2.0 keV
bands with binnings of 600 s in order to calculate an hardness ratio (HR:
$(H-S)/(H+S)$) (Fig. 10). From the HR values we derived a mean coronal
temperature of the star by means of a library of synthetic spectra calculated
with XSPEC software. The spectra were modeled with a single thermal APEC
component with temperatures in a range between 0.1 and 2.0 keV. From the
synthetic spectra convolved with the spectral response and effective area of
XMM-Newton pn, we derived counts and HRs in the bands $0.3-1.0$ keV (soft) and
$1.0-2.0$ keV (hard) as for the real data. We used the effective area files
calculated for the Thin1 and the Medium filters and used the appropriate set
of simulated spectra for each observation to calculate the expected HR values
corresponding to each value of kT of the spectra. Then we interpolated the
simulated HRs onto the observed HRs to obtain the corresponding values of kT.
Figure 5 shows the values of HR during the 25 observations as well as the
relationship between HR and kT of the model and the values of temperatures as
a function of the XMM-Newton exposure time. The differences between the Thin1
and Medium filters are minimal (less than 1%). The intervals where the flares
of HD 189733 B occurred were removed from the analysis. The kT light curve
indicates that almost all the coronal emission of HD 189733 A is always below
1 keV. The relationship between HR and kT becomes insensitive for temperature
above 1.4 keV but it works well for HD 189733 A because of the softness of its
spectrum. Overall, we do not recognize cyclic variations of the hardness of
the coronal spectrum of HD 189733 A over a timescale of eight years or fewer,
which could arise from an activity cycle on such a timescale.

Figure 6 shows the histogram of the values of the mean temperatures (kT) of
the corona of HD 189733 A derived from HR. The distribution of kT has a median
of 0.4 keV corresponding to about 4.4 MK, the $25\%-75\%$ quantile range is in
$0.3-0.5$ keV and the maximum temperature is $\sim 0.9$ keV (10.8 MK).

From the count rates and the mean coronal kT derived from HR, we obtained
fluxes and luminosities by adopting a conversion factor (CF) between count
rates and fluxes in the 0.3–2.0 keV band. We used PIMMS to convert pn Thin1
(Medium) filter count rates to fluxes with a grid of optically thin plasma
thermal APEC models spanning kT in 0.1–1.2 keV. We set N${}_{H}=0$ cm-2 as
this is compatible with the gas absorption measured from the spectra of HD
189733 A and its nearby distance (Pillitteri et al., 2010), and the overall
metal abundance was chosen equal to $Z=Z\odot$.222 This choice is supported by
the results of the best fit of the full exposure spectrum of observation nr.
15 where there was negligible flaring activity. From a best fit to a model
with two thermal APEC components we obtained (90% confidence ranges quoted in
braces): $kT_{1}=0.25\ (0.22-0.28)$ keV, $kT_{2}=0.8\ (0.76-0.91)$ keV, and
$Z=0.15Z_{\odot}\ (0.11-0.19)$. We thus derived a relationship between kT and
CF for observations made with either Thin1 or Medium filters. The CF values
range between $1.08\times 10^{-12}$ and $1.31\times 10^{-12}$ erg cnt-1 cm-2.
We interpolated these relationships to obtain the CFs for each time bin of
each observation and relative to Thin1 or Medium filter. We also corrected the
count rates for systematic offsets like the fraction of PSF comprised in the
extraction region of the events and dead time due to the read out of the
frames with the SAS task epiclcorr. From the fluxes, we finally obtained
luminosities, assuming a distance of 19.775 pc (plx$=50.57\pm 0.03$ mas, Gaia
Collaboration, 2020). A synoptic view of the X-ray luminosity for all the XMM-
Newton observations is given in Fig. 11.

Figure 5:  Hardness ratio and coronal temperature. Top panel: Time series of
the hardness ratio (HR) calculated in the bands $0.3-1.0$ keV (soft) and
$1.0-2.0$ keV (hard) as a function of the exposure time (bins of 600 s).
Vertical lines indicate the durations of each XMM-Newton observation. The
average uncertainty of HR is about 0.09 (quantile range 25%–75%: 0.08–0.11).
Middle panel: Relationship between HR and temperature (kT) derived with a
library of synthetic spectra and XSPEC software for the case of Thin1 filter
(solid line) and Medium filter (dashed line). Bottom panel: Time series of kT
expressed in keV (average uncertainty is 0.09 keV, range 25%–75%: 0.07–0.12
keV).

Figure 6:  Histogram of the values of the mean temperature (kT) inferred from
the hardness ratio values.

####  3.2.2 Flares of HD 189733 A

Once we identified the main flares of HD 189733 B, we went on to identify the
flares of HD 189733 A. We first discarded the four time intervals where HD
189733 B was flaring and thus the predominant source of variability in the
light curves. For identifying the flares of HD 189733 A and their lengths, we
adopted a semi automatic method. We first determined a quiescent rate for each
light curve corresponding to the 10% quantile of the distribution of the count
rates in that curve. This approach excludes the bins where the count rate
drops to low values due to telemetry issues (as in observation 2 of Fig. 9 at
$\sim 23$ ks). The median quiescent luminosity of HD 189733 A in the $0.3-2.0$
keV band is about $4.5\times 10^{27}$ erg s-1($25\%-75\%$ range $4.3-5.3\times
10^{27}$ erg s-1).

We used the ChangePoint package (version 2.2.2) within R language (Killick &
Eckley, 2014) to detect changes of the count rates and to segment the light
curves in intervals of statistically constant count rates. In our case we
applied the ChangePoint package with the PELT algorithm to the background
subtracted count rates of HD 189733 A and to the light curves of the
background. The penalty parameter ($p$) is crucial for a coarser or finer
segmentation of the light curve and needs to be properly set. A higher $p$
results in fewer time segments, while a lower $p$ produces more segments. The
optimal value of $p=0.006$ was assessed with a set of simulations of a
constant source with count rates normally distributed around a mean value
($r_{q}$), equal to the quiescent count rates that we determined in each light
curve and with standard deviation ($\sigma_{q}$) equal to the median of the
errors of the net count rates (after propagating the errors due to background
subtraction). Thus, a value of $p$ was chosen such that ChangePoint detects at
most one break per light curve on average under the expectation that a
constant (noisy) count rate produces only one time segment equal to the full
exposure.

As explained in Sect. 2, the bin size of light curves was chosen in a way to
guarantee enough counts per bin. We tested the use of smaller bin sizes of 100
s and 300 s with Changepoint. The results were that such finer binning
produced a few more short time segments due to noise without additional
information on the variability.

Using simulated light curves, we quantified also the sensitivity to small
flares that could be confused with small statistical variability and were thus
not detected. To do so, we injected a number of flares in the light curves of
constant source described above. The flare profile was modeled with two
exponentials for the rise and the decay phases. The analytical form is:

$F(t)=\cases{P}\,e^{(t-t_{0})/\tau_{r}}&\quad\mathrm{for\ t<t_{0}}\\\
P\,e^{-(t-t_{0})/\tau_{d}}\quad\mathrm{for\ t>t_{0},}{}$

where t0 sets the time of the peak, $\tau_{r}$ and $\tau_{d}$ are the typical
rise and decay times, respectively, and P is the peak normalization. To reduce
the parameter space, we chose 1 ks for the rise time and 5 ks for the decay
time. These values are reasonable choices that match well the observed
characteristic rise and decay phases. The only free parameters of the flare
model were thus the time of the peak ($t_{0}$) and intensity of the peak ($P$)
in units of sigma ($\sigma_{q}$).

The number of flares to inject was determined from a power law relationship
$dN/dE/day\propto E^{-2}$ (see Sect. 4). We simulated flares with energies
between $2\times 10^{30}-1.2\times 10^{31}$ erg, which is the range where we
observe a flattening of the flare distribution (cf. Fig. 8) and that
corresponds to peaks in the range $0.5-1\sigma_{q}$ above the quiescent rates.
The expected number of flares per day in this range of energy is about 16. We
scaled this number for the duration of each observation and, thus, the number
of injected flares was 5–11, depending on the length of each light curve. The
ChangePoint algorithm applied to the simulated light curves resulted in about
three time segments per light curve due to small flares detected at the
penalty threshold chosen for the real data. These elevated segments are
generated from a single flare or, more frequently, from the blend of several
small flares.

Once we identified the time segments with ChangePoint, we grouped them into
flares defined as the union of contiguous segments starting with a rising or a
peaking segment and followed by one or more decay segments. In observation no.
2, the dip during the decay is due to a telemetry glitch; in this case, we
considered the decay lasting from peak to the end of the observation. A
reference level is given by the quiescent rate, $r_{q}$, defined above. In
some cases, the average rate is significantly above the quiescent level in a
block more than 5 ks long, but it cannot be clearly associated to a single
flare. In these cases, we checked whether the flickering could be due to
blended small flares, as has been observed in several simulations. In those
cases, and only when the average count rate in one block was higher than
1$\sigma_{q}$ above the quiescent level $r_{q}$, we tried to identify small
flares subdividing the elevated interval into small ones. As a guidance, we
checked the HR that tends to increase during the flares because of plasma
heating. We applied the ChangePoint algorithm to the background curves as well
and we discarded the time intervals when the background was very high from
searching for flares (see Fig. 9)

The number of flares identified in this way was 60 – with 47 occurring in
observations taken around primary transits of HD 189733 Ab and another 13
occurring in observations taken at secondary transits. The rates of flares per
ks were very similar: $0.063\pm 0.008$ ks-1 for the total sample, $0.065\pm
0.01$ ks-1 for the observations at primary transits and $0.054\pm 0.015$ ks-1
for the observations at secondary transits. The rates of flares per day were
$5.4\pm 0.7d^{-1}$, $5.7\pm 0.8d^{-1}$, $4.7\pm 1.3$, respectively, for the
full sample and the two sub-samples observed at the primary and secondary
transits.

We multiplied the flares luminosities for each bin time (600 s) after
subtracting the quiescent level and summed the energies per bin to obtain the
total flare energies. In the next section, we discuss the properties of the
flares, their energies, and the comparison between them and the flares
observed in the Sun and solar type stars.

During transits, the planet could obscure a surface spot onto the star with
projected size similar to that of the planet. Indeed, the available TESS light
curve shows that some of the profiles of the optical transits are distorted
due to the planet passing over a stellar spot (Colombo, private
communication). Moreover, it seems that the same spot lasted for more than one
orbital period ($\sim 2.2$ days), so that the planet could have passed in
front of the same region more than once. Photospheric spots are the anchors of
magnetic loops and coronal active regions. The time required for the planet to
transit across an active region of the size of Jupiter is about 1000 s. The
possibility that the obscured region would flare during that interval and that
the same flare would be unnoticed is, however, reputed to be low. Given the
count statistics of the X-ray light curves and the chosen bin time, the event
is hardly detectable with a high significance. Moreover, the patchy and
variable stellar surface in X-rays makes these kinds of obscuration even less
detectable.

##  4 Discussion and conclusions

In this paper, we analyze, in a homogeneous way, 25 observations of the HD
189733 system obtained with XMM-Newton between the years 2007-2015. This set
of observations is quite unique in the context of stars hosting planets
because they accumulated almost 1 Ms, devoted to observations of the same star
in X-ray at two opposite phases of the orbit of its planet. The present study
offers thus a comprehensive characterization of the variability across the
lifetime of eight years for HD 189733 A.

In Table 2, we report the list of identified flares with their duration,
energy, total luminosity, and flux. We also report the flux received at the
surface of HD 189733 Ab. When considering all the flares, the total flaring
time was 423.6 ks out of 958 ks for a flare percentage of 44% of the exposure
time. The total time of the flares that occurred around the secondary transits
amounts to about 123.6 ks or $\approx 51.5$% of the observing time (240 ks).
The corresponding time at the primary transits was 300 ks or 41.8% of the
observing time (718 ks). In the observed range of flare energies, the flares
detected at secondary transits are more frequent than the flares detected at
primary transits.

Figure 7:  Empirical cumulative distribution functions (ECDFs) of the energies
of flares detected around primary and secondary transits. The arrow marks the
largest separation between the two ECDFs at $\log E\approx 31.63$.

The range of flare energies spans about 1.2 orders of magnitude between $\log
L_{X}\approx 30.6$ erg and $\log L_{X}\approx 32.2$ erg. The medians of flare
energies observed around primary and secondary transits are slightly different
($\log E_{X}\approx 31.31$ and 31.58, respectively) but the ranges of energies
are similar in both sub-samples. A Kolmogorov-Smirnoff test applied to the two
samples indicates a probability of $\sim 13$% that the energies of flares seen
at the secondary transit are drawn from the same underlying distribution of
energies of flares detected at the primary transits. From this, we conclude
that there is a slight, but not statistically significant hint that flares at
the secondary transits were more energetic. In Fig. 7, we show the empirical
cumulative distribution functions (ECDFs) of the two samples of flare
energies. The greatest difference between the two ECDFs is seen at $\log
E\approx 31.6$. The peak luminosities and the duration of the flares of the
two samples do not show a trend, meaning that it is likely a combination of
flares duration and peak luminosities to result in the slight difference of
energies.

Figure 8 shows the rate of flares per day with energies above a given
threshold for the total sample, the flares detected near the planetary
transits and those detected near or after the secondary transits. The
frequency of X-ray flares with respect to their energies in the Sun and in
Main Sequence stars is well described by a power law relationship:

$\frac{dN}{dE_{X}}\propto E_{X}^{-\alpha}$,

where the index $\alpha$ is of order of 2 (Hudson, 1991; Kashyap et al., 2002;
Drake et al., 2013). The log-log distribution shown in Fig. 8 allows us to
estimate the slope $\gamma$ of the linear relationship between the number of
flares and their energies:

$\log N(E>E_{0})=\beta+\gamma\log E_{X}$,

with the slope $\gamma$ being related to $\alpha$ by the relationship
$\gamma=1-\alpha$.

The bias due to completeness of unidentified small flares makes the curves in
Fig. 8 flatten below a given threshold ($E_{th}$). The value of $E_{th}$
depends on the sensitivity of the instruments and the intensity of the flares.
Above $E_{th}$ the cumulative distribution is deemed complete and not affected
by unidentified smaller flares, thus, it is possible to fit the tail of the
distribution without any completeness bias. We calculated the slope at
different energy thresholds and checked that for values $\log E_{X}>31.2,$ the
slope does not change significantly. We thus considered only the 43/60 flares
with energies above $\log E_{th}=31.2$ that roughly correspond to flares with
peak luminosities between $5\times 10^{27}$ and $2.8\times 10^{28}$ erg s-1or
about one to six times the quiescent luminosities. In Sect. 3, we gave a
statistical estimate of the number of small flares we could have missed. Their
number amounts to about seven per light curve and 175 for the total
observation time (11.1 days), or about ten small flares per day with energies
$30.5<\log E<31.2$. However, in the case of a larger power law index
$\alpha\sim-2.5,$ the number of undetected small flares with $30.5<\log
E<31.2$ can be up to 40 per day. In principle, our simulation of injected
flares could be used to retrieve the completeness limit as a function of the
flares energy. However, the flare model depends on a few choices we made about
the shape and the rise and decay characteristic times, which we kept fixed to
reduce the space parameter. The flare energies were obtained from the
integration of the flare profiles, but any given value of total energy could
correspond to flares with different profiles and duration. For this reason, we
preferred to infer the completeness limit directly from the data.

Figure 8:  Cumulative distributions of flares per day with energies above a
given threshold ($\log N(E>E_{0})$). Dashed line indicates the full sample of
flares ($N_{\mathrm{flares}}=60$); thin solid line is relative to the flares
observed at planetary transits ($N_{\mathrm{flares}}=47$); thick solid line is
relative to the flares observed at secondary transits
($N_{\mathrm{flares}}=13$). The lines refer to a linear fit of each sample in
the form $\log N(E_{X}>E_{0})\propto\beta\log E_{X}$. The range with energies
$\log E_{X}<31.2$ not used for the fit to the slopes is shaded in gray color.

For the full sample of flares, we obtained: $\alpha\sim 2.45\pm 0.06$; for the
flares observed around the primary transits: $\alpha\sim 2.53\pm 0.05$; and
for the flares observed at the secondary transit phases: $\alpha\sim 2.04\pm
0.14$. Overall, the modeling of the distribution of flare energies with a
power law is a good approximation, however, the index of the power law that we
find for HD 189733 A is steeper than the analog power law indices found is the
Sun and in other solar-type stars.

The subsample of flares detected at the secondary transits shows a flatter
energy distribution than the corresponding distribution at the primary
transits. The respective power law indices differs by about $2-3\sigma$ ,
although the uncertainties in the slopes could be underestimated given the
correlation of the errors among the bins.

We go on to consider whether the flares observed in HD 189733 A could be
associated with coronal mass ejections (CMEs). For the Sun, Yashiro et al.
(2004) observed that flares of classes higher than X2 (corresponding to fluxes
of 0.2 erg s-1 cm-2) are always associated with CMEs. We calculated the
average luminosities ($0.3-2.0$ keV) during flares and converted these to the
GOES band ($1-8\AA$), with a scaling factor calculated assuming an APEC model
with $kT=0.86$ keV, which is appropriate for the peak temperatures we inferred
from HR, assuming $N_{H}=0$ cm-2 and $Z=Z_{\odot}$. The scaling factor was
thus 0.1122. From the luminosities, we derived the fluxes at 1 AU and compared
these to the GOES scale. Using Yashiro’s classification, a hypothetical Earth
around HD 189733 A would have experienced flares of classes between X1 and
X20, more than 70% of them associated with CMEs.

Drake et al. (2013) give an estimate of the mass rate of CMEs as a function of
the GOES luminosity. The average luminosities of the flares of HD 189733 A
converted in the GOES band give values above $6\times 10^{28}$ erg. From Fig.
3 in Drake et al. (2013), we infer that the mass rates ejected in the supposed
CMEs are in excess of $10^{-12}M_{\odot}yr^{-1}$, with most of them peaking at
$10^{-11}M_{\odot}$ yr-1 and reaching up to $10^{-10}M_{\odot}yr^{-1}$ in the
most powerful flares, or about $6-600\times 10^{13}$ g s-1. Drake et al.
(2013) give also the relationship between X-ray luminosity of flares and
kinetic energy loss rate of the CME. From the flare luminosities of HD 189733
A, we infer a kinetic energy in excess of $10^{-4}$ L⊙ and a velocity of
$v_{CME}\approx 1000$ km s-1.

Since the average X-ray activity of the stars declines with age, it is
expected that also the most energetic flares become less frequent (see
Johnstone et al., 2021). Audard et al. (2000) derived an empirical
relationship between the number of flares per day with energies above
$10^{32}$ erg and the average X-ray luminosities of the stars in the band
0.01–10 keV. After converting the luminosities in 0.0–10 keV, we detected five
flares with energies above $10^{32}$ erg in 11.1 days or
N(E¿$10^{32}\mathrm{erg})\approx 0.45\,\mathrm{d}^{-1}$ flares. The
relationship described by Audard et al. (2000) for an average luminosity of HD
189733 A of $6\times 10^{28}$ erg s-1predicts about 0.46 flares per day, which
is perfectly matched with the observed rate. Two out of the five flares with
E¿$10^{32}$ erg were detected at secondary transits, again hinting that
energetic flares occur preferentially at those phases. The finding that more
energetic flares occur more frequently after the planetary eclipse is in
agreement with prior results (Pillitteri et al., 2011, 2015). The reason for
this could be some kind of star-planet interaction due to the escaping
planetary gas trapped in the global magnetic field around the star. The
interplay between ionized gas and magnetic field could generate reconnection
events leading to energetic flares (cf. Cohen et al., 2011). In support of
this scenario, the analysis of the flare observed in ObsID 0690890201
(Pillitteri et al., 2014) led us to conclude that the flaring structure had a
size of a few stellar radii, which is compatible with a magnetic structure
elongated towards the planet – produced by the dynamics of the ionized gas
escaping the planet and launched towards the star (Colombo et al. 2021 in
prep., see also the MHD modeling in Matsakos et al., 2015). Using Eq. 14 of
Serio et al. (1991) and by assuming a typical decay time $\tau_{d}=5$ ks and a
peak temperature of 0.9 keV, we estimated an average loop length of $\approx
1.6$ stellar radii. For longer lasting flares with $\tau_{d}\sim 10$ ks, the
loop would be twice in length with the top of the loop distant about 1 stellar
radius above the stellar surface. HD 189733 was one of the first cases for
which star-planet interaction was suggested from the varying emission of Ca II
H&K lines (Shkolnik et al., 2008). Pillitteri et al. (2014) derived a loop
length of about four stellar radii from the analysis of an oscillatory pattern
detected in a flare light curve. The height of this loop was close to 25% of
the separation between the planet and the star and it is above the star’s co-
rotation radius. The latter indicates there is a mechanism that maintained the
loop stable. The authors suggested that the star and the planet are
magnetically connected. This hypothesis was later supported by MHD simulations
and FUV observations (Matsakos et al., 2015; Pillitteri et al., 2015). The
latter authors speculated that the star could accrete gas evaporating from the
planet. The impact of the stream onto the stellar surface would cause a hot
spot that is approximately 90 degrees ahead of the sub-planetary point.
Pillitteri et al. (2015) also showed that the star exhibited FUV variability
in phase with the planet motion. The hot spot would become visible during the
secondary transit and be hidden to the observer during the primary transit.
The stellar activity in the hot spot is higher in terms of flare occurrence
and intensity but does not cause a perceptible enhancement in the quiescent
emission of the star, possibly due to the reduced extension of the spot.
Further MHD simulations show that the stream from the planet to the star is
dynamically varying due to the interaction with the stellar wind (Colombo et
al., 2021). In a concurrent scenario (not necessarily an alternative to the
first one), the magnetized environment between the planet and the star traps
the out-flowing ionized gas from HD 189733Ab, but the planetary orbital motion
stretches the resulting magnetic field and triggers the flares (cf. Lanza,
2013). These speculations are, however, biased by the fact that we have
observed a relatively narrow portion of planetary orbital phases and that a
continuous and repeated monitoring covering at least half orbit between the
two opposite phases has not been obtained. Further monitoring at the secondary
transit phases is desirable to match the unbalanced exposure times at the two
opposite phases and to better determine any difference of flare energies when
the planet is facing the two different sides of the star.

The average temperature of the corona of HD 189733A is estimated at around 0.4
keV (10%-90% range in 0.27-0.58 keV). It is only during flares that the
temperature rises up to about 0.9 keV. The flares observed near primary
transits and those observed at secondary transits show similar temperatures.
The lack of temperatures in excess of 1 keV is not due to a lack of
sensitivity of our analysis (which can, in principle, detect temperatures up
to $\sim 1.5$ keV). The corona of HD 189733 A remained relatively cold even
during the most energetic flares, while, for example, HD 189733 B showed
temperatures on the order of $1.3-1.5$ keV during the four flares that we
observed and that released energies $E\approx 10^{32}$ erg.

From the monitoring of chromospheric lines as activity indicators and in
X-ray, a number of stars are known to possess activity cycles even among
planet host stars. The ubiquity of such cycles are not established and depends
on several factors such as the stellar age. Saar & Brandenburg (1999) and
Böhm-Vitense (2007) showed that when comparing rotational periods with the
lengths of activity cycles, stars can be divided into two sequences: a steeper
sequence formed by young and active stars and a flatter sequence of inactive
and older stars. The Sun fits in between the two sequences. HD 189733 A does
not seem to possess an X-ray activity cycle, at least when considering a
timescale of eight years. HD 189733 A has a rotational period of 11.7 days, as
an active star (mimicking an age of less than 1 Gyr), it belongs to the steep
sequence of active and young stars and, thus, it would have an activity cycle
of $\geq 12$ years.

One peculiar behavior characterizing HD 189733 A is the fact that its
variability in X-rays is totally different from the variability in the optical
band. TESS observed HD 189733 A (TIC 256364928) during sector 41 for about 27
days and about 11 transits were recorded. The TESS light curve does not
clearly show short timescale variability due to flares. As remarked in Sect.
3, the planet obscured some surface spots but yet the star was relatively
inactive, behaving instead like an old quiet star in spite of its level of
X-ray activity and luminosity. We speculate that the stark difference between
the activity and the relatively cold flares observed in X-rays and the optical
band could be related to the angular momentum transfer from the planet to the
star and its influence on the mechanism of magnetic dynamo that produces the
corona. Furthermore, the magnetic interaction between the star and the planet
could produce a significant part of the observed X-ray emission and
variability. Lanza (2013) quantifies the power emitted by two possible
mechanisms of loop stretching due to magnetic interaction of the star+planet
magnetospheres or of magnetic loops anchored between a star and its hot
Jupiter. In particular, the second mechanism could easily produce a power of
order of $10^{28}$ erg s-1(similar to the average X-ray luminosity of HD
189733 A) for a planetary magnetic field of $\sim 3$ G at the pole and
separation equal to that of HD 189733 Ab.

We analyzed the variability of HD 189733 B, taking into account the close
angular separation with the primary HD 189733 A in order to carefully
distinguish its time variability in X-ray. This mainly consisted of four
bright flares detected in 958 ks, which correspond to a rate of one flare
every 2.8 days. On average, the flares released energies in excess of
$10^{32}$ erg. The quiescent luminosity of HD 189733 B is about $\log
L_{X}\sim 27.27$ in $0.3-2.0$ keV, which is about a factor four higher than
found in Chandra observations (Poppenhaeger et al., 2013). The luminosity is
still consistent with that of a dM star with an age of a it is Gyr; however,
in contrast to the age that would otherwise be inferred from the X-ray
activity of HD 189733 A, which is less than 1 Gyr. The age discrepancy between
the two components of the HD 189733 system has been discussed in literature
and explained as the result of spin up on the part of the star with the hot
Jupiter, accompanied by the transfer of angular momentum from the orbital
motion of the planet to the rotation of the star (Pillitteri et al., 2010;
Poppenhaeger & Wolk, 2014). The flare temperatures of HD 189733 B were in
excess of 1.3 keV – in stark contrast to the flare temperatures of the flares
of HD 189733 A, which remained below 1 keV.

###### Acknowledgements.

The authors acknowledge the anonymous referee for their useful comments and
suggestions. IP, AM and GM acknowledge financial support from the ASI-INAF
agreement n.2018-16-HH.0, and from the ARIEL ASI-INAF agreement n.2021-5-HH.0.
Based on observations obtained with XMM-Newton, an ESA science mission with
instruments and contributions directly funded by ESA Member States and NASA

## References

  * Audard et al. (2000) Audard, M., Güdel, M., Drake, J. J., & Kashyap, V. L. 2000, ApJ, 541, 396
  * Bakos et al. (2006) Bakos, G. Á., Pál, A., Latham, D. W., Noyes, R. W., & Stefanik, R. P. 2006, ApJ, 641, L57
  * Böhm-Vitense (2007) Böhm-Vitense, E. 2007, ApJ, 657, 486
  * Borucki et al. (2006) Borucki, W. J., Koch, D., Basri, G., et al. 2006, ISSI Scientific Reports Series, 6, 207
  * Bouchy et al. (2005) Bouchy, F., Udry, S., Mayor, M., et al. 2005, A&A, 444, L15
  * Bourrier et al. (2013) Bourrier, V., Lecavelier des Etangs, A., Dupuy, H., et al. 2013, A&A, 551, A63
  * Bourrier et al. (2020) Bourrier, V., Wheatley, P. J., Lecavelier des Etangs, A., et al. 2020, MNRAS, 493, 559
  * Cecchi-Pestellini et al. (2009) Cecchi-Pestellini, C., Ciaravella, A., Micela, G., & Penz, T. 2009, A&A, 496, 863
  * Cohen et al. (2011) Cohen, O., Kashyap, V. L., Drake, J. J., et al. 2011, ApJ, 733, 67
  * Colombo et al. (2021) Colombo, S., Pillitteri, I., Orlando, S., & Micela, G. 2021, arXiv e-prints, arXiv:2111.11110
  * Drake et al. (2013) Drake, J. J., Cohen, O., Yashiro, S., & Gopalswamy, N. 2013, ApJ, 764, 170
  * Gaia Collaboration (2020) Gaia Collaboration. 2020, VizieR Online Data Catalog, I/350
  * Gardner et al. (2006) Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485
  * Hudson (1991) Hudson, H. S. 1991, Sol. Phys., 133, 357
  * Johnstone et al. (2021) Johnstone, C. P., Bartel, M., & Güdel, M. 2021, A&A, 649, A96
  * Kashyap et al. (2002) Kashyap, V. L., Drake, J. J., Güdel, M., & Audard, M. 2002, ApJ, 580, 1118
  * Killick & Eckley (2014) Killick, R. & Eckley, I. A. 2014, Journal of Statistical Software, 58, 1
  * Kubyshkina et al. (2018) Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018, ApJ, 866, L18
  * Lanza (2013) Lanza, A. F. 2013, A&A, 557, A31
  * Lecavelier des Etangs et al. (2012) Lecavelier des Etangs, A., Bourrier, V., Wheatley, P. J., et al. 2012, A&A, 543, L4
  * Locci et al. (2022) Locci, D., Petralia, A., Micela, G., et al. 2022, , 3, 1
  * Matsakos et al. (2015) Matsakos, T., Uribe, A., & Königl, A. 2015, A&A, 578, A6
  * Pillitteri et al. (2011) Pillitteri, I., Günther, H. M., Wolk, S. J., Kashyap, V. L., & Cohen, O. 2011, ApJ, 741, L18
  * Pillitteri et al. (2015) Pillitteri, I., Maggio, A., Micela, G., et al. 2015, ApJ, 805, 52
  * Pillitteri et al. (2010) Pillitteri, I., Wolk, S. J., Cohen, O., et al. 2010, ApJ, 722, 1216
  * Pillitteri et al. (2014) Pillitteri, I., Wolk, S. J., Lopez-Santiago, J., et al. 2014, ApJ, 785, 145
  * Poppenhaeger et al. (2013) Poppenhaeger, K., Schmitt, J. H. M. M., & Wolk, S. J. 2013, ApJ, 773, 62
  * Poppenhaeger & Wolk (2014) Poppenhaeger, K. & Wolk, S. J. 2014, A&A, 565, L1
  * Reale (2007) Reale, F. 2007, A&A, 471, 271
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Saar & Brandenburg (1999) Saar, S. H. & Brandenburg, A. 1999, ApJ, 524, 295
  * Serio et al. (1991) Serio, S., Reale, F., Jakimiec, J., Sylwester, B., & Sylwester, J. 1991, A&A, 241, 197
  * Shkolnik et al. (2008) Shkolnik, E., Bohlender, D. A., Walker, G. A. H., & Collier Cameron, A. 2008, ApJ, 676, 628
  * Tinetti et al. (2018) Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Experimental Astronomy, 46, 135
  * Yashiro et al. (2004) Yashiro, S., Gopalswamy, N., Michalek, G., et al. 2004, Journal of Geophysical Research (Space Physics), 109, A07105

##  Appendix A Table of flares

Table 2:  List of identified flares.

Nr. Obs. | Duration | $\log E_{X}$ | $\log L_{X}$ | $\log f_{X}$ | $\log f_{X,Pl}$ | Nr. Obs. | Duration | $\log E_{X}$ | $\log L_{X}$ | $\log f_{X}$ | $\log f_{X,Pl}$
---|---|---|---|---|---|---|---|---|---|---|---
| ks | erg | erg s-1 | erg s-1 cm-2 | erg s-1 cm-2 |  | ks | erg | erg s-1 | erg s-1 cm-2 | erg s-1 cm-2
1 | 2.4 | 30.89 | 28.11 | -12.56 | 3.68 | 13 | 5.4 | 31.24 | 28.46 | -12.21 | 4.03
1 | 3.6 | 31.00 | 28.23 | -12.44 | 3.79 | 13 | 4.2 | 31.26 | 28.48 | -12.19 | 4.05
1 | 6 | 31.21 | 28.44 | -12.23 | 4.00 | 14 | 4.8 | 31.15 | 28.38 | -12.3 | 3.94
2 | 18.6 | 32.04 | 29.26 | -11.41 | 4.83 | 15 | 8.4 | 31.51 | 28.73 | -11.94 | 4.3
3 | 7.8 | 31.75 | 28.97 | -11.7 | 4.54 | 16 | 6.0 | 31.58 | 28.80 | -11.87 | 4.37
3 | 6.6 | 31.41 | 28.63 | -12.04 | 4.20 | 16 | 13.2 | 31.74 | 28.96 | -11.71 | 4.53
4 | 6 | 31.21 | 28.44 | -12.23 | 4.00 | 16 | 5.4 | 31.43 | 28.65 | -12.02 | 4.22
4 | 8.4 | 31.29 | 28.51 | -12.16 | 4.08 | 17 | 10.8 | 31.75 | 28.97 | -11.7 | 4.54
4 | 22.8 | 31.92 | 29.14 | -11.53 | 4.71 | 19 | 5.4 | 31.36 | 28.58 | -12.09 | 4.15
5 | 9 | 31.57 | 28.79 | -11.88 | 4.36 | 19 | 1.8 | 30.62 | 27.84 | -12.83 | 3.41
5 | 5.4 | 31.40 | 28.62 | -12.05 | 4.19 | 20 | 3.0 | 30.84 | 28.06 | -12.61 | 3.63
5 | 7.8 | 31.16 | 28.38 | -12.29 | 3.95 | 20 | 3.0 | 30.87 | 28.09 | -12.58 | 3.66
5 | 7.8 | 31.41 | 28.63 | -12.04 | 4.20 | 20 | 14.4 | 31.81 | 29.03 | -11.64 | 4.6
6 | 2.4 | 30.96 | 28.18 | -12.49 | 3.75 | 21 | 10.2 | 31.83 | 29.05 | -11.62 | 4.62
7 | 4.8 | 31.33 | 28.55 | -12.12 | 4.12 | 21 | 8.4 | 31.56 | 28.78 | -11.89 | 4.35
7 | 2.4 | 30.95 | 28.17 | -12.5 | 3.74 | 21 | 5.4 | 31.22 | 28.44 | -12.23 | 4.01
7 | 1.8 | 30.66 | 27.88 | -12.79 | 3.45 | 22 | 6.0 | 31.70 | 28.92 | -11.75 | 4.49
7 | 5.4 | 31.41 | 28.64 | -12.03 | 4.20 | 22 | 18.0 | 31.93 | 29.15 | -11.52 | 4.72
8 | 6 | 31.31 | 28.53 | -12.14 | 4.10 | 23 | 5.4 | 31.19 | 28.41 | -12.26 | 3.98
8 | 8.4 | 31.55 | 28.77 | -11.9 | 4.34 | 23 | 4.2 | 31.25 | 28.48 | -12.19 | 4.04
9 | 2.4 | 30.99 | 28.22 | -12.45 | 3.78 | 23 | 2.4 | 31.05 | 28.27 | -12.4 | 3.84
9 | 6 | 31.31 | 28.53 | -12.14 | 4.10 | 23 | 6.0 | 31.57 | 28.79 | -11.88 | 4.36
9 | 2.4 | 30.84 | 28.06 | -12.61 | 3.63 | 23 | 2.4 | 30.93 | 28.16 | -12.51 | 3.72
10 | 6.6 | 31.45 | 28.67 | -12 | 4.24 | 24 | 6.0 | 31.78 | 29.00 | -11.67 | 4.57
10 | 4.2 | 31.22 | 28.45 | -12.22 | 4.01 | 24 | 8.4 | 31.68 | 28.90 | -11.77 | 4.47
10 | 6 | 31.37 | 28.59 | -12.08 | 4.16 | 24 | 3.0 | 31.08 | 28.30 | -12.37 | 3.87
11 | 4.8 | 31.34 | 28.57 | -12.1 | 4.13 | 24 | 10.8 | 31.58 | 28.80 | -11.87 | 4.37
11 | 4.8 | 31.27 | 28.49 | -12.18 | 4.06 | 25 | 4.2 | 31.23 | 28.45 | -12.22 | 4.02
11 | 13.2 | 31.95 | 29.17 | -11.5 | 4.74 | 25 | 15.0 | 31.81 | 29.04 | -11.63 | 4.6
12 | 22.2 | 32.08 | 29.30 | -11.37 | 4.87 | 25 | 6.0 | 31.08 | 28.31 | -12.36 | 3.87

Note: We report the observation number (cf. Table 1), the duration of the
flares, their energy in the band 0.3-2.0 keV, the total emitted luminosities
and fluxes in the same band, and the fluxes received at the mean separation of
HD 189733 Ab.

##  Appendix B Light curves

In this appendix, we display the pn light curves of HD 189733 A (Fig. 9)
obtained in 25 XMM-Newton observations, the HR light curves derived from the
pn count rates in the bands $0.3-1.0$ keV and $1.0-2.0$ keV (Fig. 10), and the
light curves of the X-ray luminosity of HD 189733 A in the band $0.3-2.0$ keV
(Fig. 11).

Figure 9:  Light curves of HD 189733 A obtained with pn in the band $0.3-2.0$
keV. The points with error bars are the net count rates of HD 189733 A for the
events in the selection operated at points i and ii of Sect. 2. The horizontal
red segments mark the intervals identified with the Changepoint algorithm.
Similarly, the light red segments represent the intervals of constant
background. In intervals of high background variability (gray shaded areas),
we did not search for flares in HD 189733 A. The flares are shaded with
different colors in each panel. The horizontal solid blue lines represent the
quiescent rates and the $\pm 1sigma$ range. Panels contoured with a green
border refers to observations at the secondary transits. Panels marked with a
red border indicate where variability of HD 189733 B was prominent.

Figure 10:  Solid black lines with error bars: light curves of hardness ratio
of HD 189733 A obtained from pn events for each observation. Red and blue
curves: Soft (0.3-1.0 keV) and hard (1.0-2.0) band light curves scaled between
zero and the maximum of soft rate. As in Fig. 9, the panels marked with a red
border indicate observations where flares of HD 189733 B were detected. In
these curves, we discarded the intervals relative to the flares of HD 189733
B.

Figure 11:  Synoptic plot of the logarithm of the X-ray luminosity of HD
189733 A in the band 0.3-2.0 keV observed in the XMM-Newton observations. Gaps
in the light curves correspond to intervals where flares from HD 189733 B
affected the estimated of the count rate of HD 189733 A.


# Generating Function For Network Delay

A.M. Sukhov, N.Yu. Kuznetsova, A.K. Pervitsky and A.A. Galtsev  Samara State
Aerospace University, Moskovskoe sh., 34, Samara, 443086, Russia
e-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper correspondence between experimental data for packet delay and
two theoretical types of distribution is investigated. Statistical tests have
shown that only exponential distribution can be used for the description of
packet delays in global network. Precision experimental data to within
microseconds are gathered by means of the RIPE Test Box. Statistical
verification of hypothesis has shown that distribution parameters remain
constants during 500 second intervals at least. In paper cumulative
distribution function and generating function for packet delay in network are
in an explicit form written down, the algorithm of search of parameters of
distribution is resulted.

##  I Introduction

The special area of the control theory, named networked control systems in
which transfers as environment of operating signals were used computer
networks, has arisen in the late nineties of the XX-th century [18].
Originally, as the network environment of control systems local networks [8]
which differ high-speed data transfer and in the minimum percent of packet
loss were used.

The networked control systems in which as the handle environment the global
network Internet is used, is extremely complicated because of random character
of distribution of packet delay and their big absolute values [10, 15, 17].
However till now, results of the advanced network researches are not used in
the control theory and algorithms on their basis are not created. The present
project assumes introduction of new network decisions in networked control
systems.

Except classical problems of the control theory, there is requirement of
management and monitoring of network processes. For example, for video
transfer over networks TCP/IP networks the important parameter is the
available bandwidth between two hosts [4]. For problems of routing it is
necessary to know throughput between routers. Special network emulators are
applied to modeling of the majority of network processes, for example,
INET/OMNET ++ [16].

Now numerous programs which imitate transmission of packets through TCP/IP
protocol network are created. In a basis of operation of all emulators
position that the type of delay distribution is unknown is occupied. The
purpose of our research consists in that not only to define delay type, but
also to find generating function for the traffic emulators.

For the decision of problems of the networked control systems on the basis of
stack TCP/IP it is more convenient to use the process approach to the control
theory, based on idea of existence of some universal functions of control. The
purpose of our research is the finding of this function for network
components. In the modern theory of computer networks there were many
utilities working with a network delay, there is a progress in studying and
modeling of transmitting of packets. Our problem consists in trying to
describe process of a network delay of management packets and to show ways of
practical calculation of all parameters entering into corresponding
distribution functions [9].

By transmission of control signals through TCP/IP network, the separate
packets of the controlling data flow transferring the information are supplied
non-uniformly, and the part of packets in general is lost by transmission on a
network and does not reach a target. For rise of efficiency of control
algorithms it is necessary to reduce to a minimum of packets delay and their
variation, and also percent of packet loss. Similar algorithms are used for
transmission voice and video streams, in grid systems, at control of robust
systems, in network computer games, etc.

At first it would be desirable to result the brightest research on a
distribution type for network delay. To understand, about what there is a
speech in described papers, will give definitions of notations used in them:

  * •

Round-trip time (RTT) time is the time required for a packet to travel from
the testing host to a remote computer that receives the packet and retransmits
it back to the source.

  * •

The One-Way Delay (OWD) value is calculated between two synchronized points A
and B of an IP network, and it is the time in seconds that a packet spends in
travelling across the IP network from A to B.

In particular, Elteto and Molnar [6] have spent measurements of round-trip
delay in the Ericsson Corporate Network, complex analysis of the received data
has allowed to build the supposition about distribution type for network
delay. The main finding of their research is that the round-trip delay can be
well approximated by a truncated normal distribution.

Konstantina Papagiannaki et al [12] in the research have measured and have
analyzed packet delay between two adjacent routers in the core network. On the
basis of the received measurements, they have made the supposition about the
factors influencing occurrence of delay, and very big delays which cannot be
explained in the way of batch processing in routers on algorithm FIFO have
been noticed.

Recently, fulfilling a series of operations on measurement of an available
bandwidth [13], we have installed that for a type definition of delay
distribution we should research only a variable part of delay while its most
part remains constant. This fact also has served as a starting point of our
operation.

##  II Premises for model

In 1999 Downey [5] for the first time has detected linear dependence of the
minimum possible round trip time on the size of transferred packets. In 2004
precise experiments by Choi et al [1], Hohn et al [11] proved that the minimum
fixed delay component $D^{fixed}(W)$ for a packet of size $W$ is a linear (or
precisely, an affine) function of its size,

$D^{fixed}(W)=W\sum_{i=1}^{h}1/C_{i}+\sum_{i=1}^{h}\delta_{i}$ (1)

where $C_{i}$ is each link of capacity of $h$ hops and $\delta_{i}$ is
propagation delay. To validate this assumption, they check the minimum delay
of packets of the same size for three path, and plot the minimum delay against
the packet size.

Let $D(W)$ represents the one way delay (point-to-point delay) of a packet.
Here we refer to it as the minimum path transit time for the given packet size
$W$, denoted by $D^{fixed}(W)=\min D(W)$. With the fixed delay component
$D^{fixed}(W)$ identified, we can now subtract it from the point-to-point
delay of each packet to study the variable delay component $d^{var}$. The
variable delay component of the packet, $d^{var}$, is given by

$D(W)=D^{fixed}(W)+d^{var}$ (2)

Computed minimal delay $D^{fixed}(W)$ is

$D^{fixed}(W)=D_{min}+W/C,$ (3)

where $C=(W_{2}-W_{1})/(D_{2}-D_{1})$ is end-to-end available bandwidth and it
is searched, comparing average time of packet delay of the different sizes
$W_{2}$ and $W_{1}$ [13]. Here

$D_{min}=\lim_{W\rightarrow 0}D^{fixed}(W)$ (4)

The value $D_{min}$ is related to the distance between the sites (i.e.
propagation delay) and per-packet router processing time at each hop along the
path between the sites [2, 3]. This value represents as the minimum delay
$D_{min}$ for which the very small package can be transmitted on a network
from one point in another.

The minimal delay [13] of datagram transmission $D_{min}$ may be calculated as

$D_{min}=\frac{W_{2}D^{fixed}(W_{1})-W_{1}D^{fixed}(W_{2})}{W_{2}-W_{1}}$ (5)

This value as well as the methods of its measurement has a important
significance in applied tasks of control theory [18]. The second significant
question of networking control theory is the distribution type for variable
delay component $d^{var}$ which should be studied in rest of paper. To know
the expression for this parameter we may easy calculate the duration of buffer
for streaming application on receiving side.

##  III Experimental search

In order to determine distribution type for a variable delay component
$d^{var}$ we should run considerable quantity of measurements between various
hosts in the Internet. The basic problem of experimental testing is the
precise of delay measurements that is necessary for accurate result. Similar
measurements demand, at least, microsecond precision for delay measurements;
we are reaching such accuracy with help of RIPE Test Box mechanism [7, 14]. In
order to prepare the experiments three Test Boxes have been installed in
Moscow, Samara and Rostov on Don during 2006-2008 years in framework of RFBR
grant 06-07-89074. Each RIPE Test Box represents a server under management of
an FreeBSD operating system with the GPS receiver connected to it.

Characteristic times of investigated processes (a packet delay, jitter) have
the order from 10 $ms$ to 1 $sec$, therefore is enough to use system hours of
a RIPE Test Box for their reliable measurement. At the first stage experiment
between tt01.ripe.net (RIPE NCC at AMS-IX, Amsterdam), tt143.ripe.net (Samara,
SSAU), tt17.ripe.net (Bologna) and tt74.ripe.net (Melbourne) have been made
which include precision measurement of packet delay with accuracy 2-12 $\mu
s$. Testing results are available via telnet to corresponding RIPE Test Box on
port 9142. It is important to come and write down simultaneously the data on
both ends of the investigated connection.

On the basis of the received data set it is easy to construct a cumulative
distribution function for network delay $D$:

$F(D)=P(x\leq D)$ (6)

For initial comparison truncated normal and exponential distributions have
been chosen, expressions for which are written down.

For truncated normal distribution it is possible to select following
approximation:

$F(D)=\left\\{\begin{IEEEeqnarraybox}[]{[}][c]{rl}0,&D<D_{min};\\\
\frac{\sqrt{2/\pi}}{\sigma}\int\limits_{D_{min}}^{D}\exp\left\\{-\frac{(x-D_{min})^{2}}{2\sigma^{2}}\right\\}dx,&D\geq
D_{min}\end{IEEEeqnarraybox}\right.$ (7)

where

$\sigma=D_{av}-D_{min}$ (8)

is the difference between average network delay $D_{av}(W)=\mathbb{E}[D(W)]$
and minimum delay $D_{min}(W)$.

It should be noted that all statistical data has been gathered by us for the
fixed size of a packets $W$. By default for RIPE Test Box it equals to 100
bytes. In Section V we update a cumulative distribution function $F(D,W)$
taking into account the packets size $W$.

The alternative type of allocation which will be checked on correspondence is
an exponential distribution, expression for which is written below.

$F(D)=\left\\{\begin{IEEEeqnarraybox}[]{[}][c]{rl}0,&D<D_{min}\\\
1-\exp\left\\{-\lambda(D-D_{min})\right\\},&D\geq
D_{min}\end{IEEEeqnarraybox}\right.$ (9)

where

$\lambda=1/(D_{av}-D_{min})$ (10)

is reciprocal to the difference between average network delay
$D_{av}(W)=\mathbb{E}[D(W)]$ and minimum delay $D_{min}(W)$.

For initial check of conformity to distribution type two methods will be used:
calculation of Pearson correlation coefficient and a graphic method. We will
designate as $K_{nor}$ correlation coefficient between experimental and normal
distributions then $K_{exp}$ is correlation coefficient between experimental
and exponential distributions.

The obtained data is shown in Table I, where the column host corresponds to a
direction between two RIPE Test Boxes, and the column $W$ specifies in the
size of a testing packet.

N | host |  $W$ (bytes) | $K_{nor}$ | $K_{exp}$
---|---|---|---|---
1 | bolonia |  |  |
| tt01$\Rightarrow$tt17 | 100 | 0.76 | 0.97
2 | samara |  |  |
| tt01$\Rightarrow$tt143 | 100 | 0.87 | 0.98
3 | samara |  |  |
| tt01$\Rightarrow$tt143 | 1024 | 0.99 | 0.99
4 | melburn |  |  |
| tt01$\Rightarrow$tt74 | 100 | 0.66 | 0.97
TABLE I: Precise measurements

Except correlation coefficients it is possible to compare graphics
representation of cumulative distribution functions (CDF), showing all three
functions on one plot. On the uniform graphics (see Figures 1, 2) dash line
selects an experimental curve, dot-dash curve corresponds to normal
allocation. In dot style the exponential distribution is painted.

The plots constructed in Figures 1 and 2, represent dependence of CDF on delay
of a packet on a site from Amsterdam to Samara (tt01$\Rightarrow$tt143). The
first plot describes testing of a network by packages in the 100 bytes size,
the second plot corresponds to packages in the size of 1024 bytes. Time on
axis is measured in milliseconds.

Figure 1: Experimental (dash), normal (dash-dot) and exponential (dot) CDFs,
precise testing. Direction: tt01$\Rightarrow$tt143, $W=$100 bytes Figure 2:
Experimental (dash), normal (dash-dot) and exponential (dot) CDFs, precise
testing. Direction: tt01$\Rightarrow$tt143, $W=$1024 bytes

All experiments resulted above testify that the preferable type of
distribution describing packet delay in a global network is an exponential
distribution. Thus, as have shown our researches, the random variable of
packet delay between two network points is arranged on an exponential low with
the parameter calculated from experimental values under the Equation (10).

However, not each investigator who is engaged in the control theory has the
equipment for the precision measurements, similarly RIPE Test Boxes. Therefore
in this part it would be desirable to present technique which uses the data of
well-known utilities and doesn’t demand the expensive equipment.

For testing it is possible to use the utility ping as it is the most
widespread resource for verification of connection quality in TCP/IP networks.
Let’s mark only that this utility measures round-trip time, instead of one way
delay.

The data received with help of ping has the millisecond precision that is
exact enough to judge delay distribution. The utility ping allows testing
connections between points AIST - New Zealand (tt47.ripe.net), Volgatelekom -
Australia (tt74.ripe.net) and SSAU-Melbourn (tt74.ripe.net). As remote hosts
were used servers of RIPE measurement system, AIST, Volgatelecom (VT),
Infolada and SSAU is local Internet Service Providers from Samara region,
Russia. Processing the obtained data on the above described algorithm, we have
received the results presented in the Table II.

N | host |  $W$ (bytes) | $K_{nor}$ | $K_{exp}$
---|---|---|---|---
1 | AIST$\Rightarrow$ |  |  |
| New Zeland | 32 | 0.94 | 0.95
2 | Volgatelecom$\Rightarrow$ |  |  |
| Australia | 32 | 0.96 | 0.98
3 | SSAU$\Rightarrow$ |  |  |
| Melburn | 64 | 0.66 | 0.97
4 | Infolada$\Rightarrow$ |  |  |
| Athens | 32 | 0.98 | 0.98
TABLE II: ping measurements

The evident illustration is resulted in definition of distribution type in
Figure 3-5.

Figure 3: Experimental (dash), normal (dash-dot) and exponential (dot) CDFs,
precise testing. Direction: Samara$\Rightarrow$Holland, $W=$32 bytes Figure 4:
Experimental (dash), normal (dash-dot) and exponential (dot) CDFs, precise
testing. Direction: Infolada$\Rightarrow$Athens, $W=$32 bytes Figure 5:
Experimental (dash), normal (dash-dot) and exponential (dot) CDFs, precise
testing. Direction: SSAU$\Rightarrow$Australia, $W=$1064 bytes

It should be noted that the utility ping allows finding automatically values
of variables $D_{av}$ and $D_{min}$ (see Eqns. (7) and (9)) which completely
define the distribution form, both normal and exponential types. It is enough
to give sequence from 20 packets to obtain the given values with split-hair
accuracy, sufficient for the description of processes of the control theory.

##  IV Statistical hypothesis testing

The verification executed in the previous section about conformity of
distribution type has preliminary character and isn’t strict. In this section
for check of distribution type Pearson’s chi-square test will be used.

For processing some data sets from 2000 to 2500 values the delays collected
with interval in 2 seconds have been used. This data divides into intervals in
50, 100, 250, 500, 1000 and 2000 values and was tested on Pearson. Testing
results could be found in Tables III, IV.

At construction of these Tables following designations were used:

  * •

Dimension of observations $N$ (number of measurements)

  * •

$n$ is the number of cells. All observations $N$ are divided among $n$ cells
according Sturges’ rule $n=(1+3.22\lg N)+1$

  * •

$t$ is the value of the test-statistic

  * •

$\chi^{2}_{0.95,n-1}$ is the theoretical value of threshold of hypothesis
acceptance

If $t<\chi^{2}_{0.95,n-1}$ then hypothesis about corresponding type of
distribution is accepted, differently hypothesis is rejected.

In the beginning we will check up on conformity to exponential distribution
the data, reception by means of RIPE Test Box, see Table III.

$N$ | 50 | 100 | 200 | 250 | 500 | 1000 | 2000
---|---|---|---|---|---|---|---
$n$ | 14 | 17 | 19 | 20 | 22 | 24 | 27
$\chi^{2}_{0.95,n-1}$ | 22.36 | 26.30 | 28.87 | 30.14 | 32.67 | 35.17 | 38.89
$t$ | 21.29 | 25.35 | 22.77 | 23.10 | 134.31 | 547.16 | 978.98
Hypothesis acceptance | Yes | Yes | Yes | Yes | No | No | No
TABLE III: Verification of exponential distribution, Samara-Amsterdam,
(tt143$\Rightarrow$tt01), packet size 100 bytes

From Table III follows that within 500 second intervals (250 measurements) the
packet delay is distributed on the exponential law.

The data received with RIPE Test Box, was checked also on conformity to the
truncated normal distribution, see Table IV.

$N$ | 50 | 100 | 200 | 250 | 500 | 1000 | 2000
---|---|---|---|---|---|---|---
$n$ | 14 | 17 | 19 | 20 | 22 | 24 | 27
$\chi^{2}_{0.95,n-1}$ | 22.36 | 26.30 | 28.87 | 30.14 | 32.67 | 35.17 | 38.89
$t$ | 43.32 | 217.46 | 2906.47 | 6077.07 | $\infty$ | $\infty$ | $\infty$
Hypothesis acceptance | No | No | No | No | No | No | No
TABLE IV: Verification of truncated normal distribution, Samara-Amsterdam,
(tt143$\Rightarrow$tt01), packet size 1024 bytes

Pearson’s chi-square test allows rejecting hypothesis about the truncated
normal type of distribution for the description of process of packet delay.

##  V Distribution Type for Delay and Generating Function for Traffic Emulator

In real the Internet processes the size of transferred packages can vary,
therefore the cumulative distribution function should be updated. For each
size of a packet $W$ there is the minimum time $D^{fixed}(W)$ defined by the
Eq. (3).

Then, the final cumulative distribution function $F(D,W)$ is

$F(D,W)=\left\\{\begin{IEEEeqnarraybox}[]{[}][c]{rl}0,D<D_{min}+W/C&\\\
1-\exp\left\\{-\lambda(D-D_{min}-W/C)\right\\},&\\\ D\geq
D_{min}+W/C&\end{IEEEeqnarraybox}\right.$ (11)

where

$\lambda=1/(D_{av}-D_{min})$ (12)

is reciprocal to the difference between average network delay
$D_{av}(W)=\mathbb{E}[D(W)]$ and minimum delay $D_{min}(W)$ and $C$ is end-to-
end capacity.

It is important to notice that the data on the delay, received with use of
packages of various length, comprises the additional variable component caused
network jitter $j$ (delay variation) [4]. Therefore, the best the controlling
algorithm will form packages of the identical size. If we use the utility ping
for delay definition in it there is a special key for resizing of a testing
package ($-l$ in Windows, $-s$ in Linux).

Our model has one more wide application except tasks of the control theory.
Results of our operation can be applied at writing of emulators of the traffic
in a global network [16]. Till now it was considered that the type of delay
distribution is unknown, and traffic emulators used own generating functions
for delay generation. On the basis of the type of delay found us (see Eqn.
(11)) it is possible to write generating function

$D=D_{min}+W/C-(1/\lambda)\ln(1-F(D,W))$ (13)

In this equation content distribution function (CDF) $F(D,W)$ can be set the
generator of random numbers. The received numbers will give values of delay
for network packet of the different size. We will notice once again that in
real networks these values can be calculated according to the specilised
utility.

##  VI Conclusion

In the present work for the description of process of the packet delay in a
global networks it has been chosen exponential distribution. In comparison
with truncated normal distribution it has shown the best correlation with
experimental results. Unlike the truncated normal distribution it passes check
by statistical methods.

Experimental data were gathered by means of the precision RIPE measuring
system to within microseconds, and also by means of the standard utility ping.
This utility measures round-trip time to within milliseconds. During small
periods about several minutes when it is possible to consider conditions of
transmission on a TCP/IP network invariable, such approach gives correlations
from above 0.99. At change of network conditions the elementary ping testing
by a series from 20 packets will allow to change exponential distribution
parameters instantly.

We write in an explicit form of cumulative distribution functions for normal
and exponential distribution of delay. Generating function for packet delay
has been found in real networks which can be used in emulators of the traffic
of a global network.

In summary we would like to thank Leonid Fridman, the Professor from
University of Mexico for fruitful dialogues in which course the idea of this
article has taken shape. Also it would be desirable to thank all collective of
technical service RIPE ncc and especially Ruben van Staveren and Roman
Kalyakin for constant assistance in comprehension of subtleties of a measuring
infrastructure. We also would like to express the gratitude to the Wolfram
Research corporation, which the first has marked our preprint and has given us
licenses for the right of use of Mathematica

## References

  * [1] Choi, B.-Y., Moon, S., Zhang, Z.-L., Papagiannaki, K. and Diot, C.: Analysis of Point-To-Point Packet Delay In an Operational Network. In: Infocom 2004, Hong Kong, pp. 1797-1807 (2004)
  * [2] Cottrell, L., Matthews, W. and Logg C.: Tutorial on Internet Monitoring $\&$ PingER at SLAC. http://www.slac.stanford.edu/comp/net/wan-mon/tutorial.html
  * [3] Crovella, M.E. and Carter, R.L.: Dynamic Server Selection in the Internet. In: Proc. of the Third IEEE Workshop on the Architecture and Implementation of High Performance Communication Subsystems (1995)
  * [4] Dovrolis C., Ramanathan P., and Moore D., Packet-Dispersion Techniques and a Capacity-Estimation Methodology, IEEE/ACM TRANSACTIONS ON NETWORKING, VOL. 12, NO. 6, DECEMBER 2004, p. 963-977
  * [5] Downey A.B., Using Pathchar to estimate internet link characteristics, in Proc. ACM SICCOMM, Sept. 1999, pp. 222–223.
  * [6] Elteto, T., Molnar, S., On the distribution of round-trip delays in TCP/IP networks, in The Proceedings of the Local Computer Networks (LCN 99) Conference, IEEE, 1999, pp. 172–181
  * [7] Georgatos, F., Gruber, F., Karrenberg, D., Santcroos, M., Susanj, A., Uijterwaal, H. and Wilhelm R., Providing active measurements as a regular service for ISP’s. In: PAM2001
  * [8] Georges J.-P., Divoux T., and Rondeau E., Confronting the performances of a switched ethernet network with industrial constraints by using the network calculus, _International Journal of Communication Systems(IJCS)_ , vol. 18, no. 9, pp. 877–903, 2005
  * [9] Fridman E., Seuret A., and Richard J.-P., “Robust sampled-data stabilization of linear systems: an input delay approach,” _Automatica_ , vol. 40, pp. 1441–1446, 2004
  * [10] Hespanha J.P., Naghshtabrizi P., and Xu Y., “A survey of recent results in networked control systems,” _Proceedings of the IEEE_ , vol. 95, pp. 138–162, January 2007
  * [11] N. Hohn, D. Veitch, K. Papagiannaki and C. Diot, Bridging Router Performance And Queuing Theory, Proc. ACM SIGMETRICS, New York, USA, Jun 2004
  * [12] D. Papagiannaki, S. Moon, C. Fraleigh, P. Thiran, F. Tobagi, and C. Diot, Analysis of measured single-hop delay from an operational backbone network, in Proc. IEEE INFOCOM 2002, New York, New York, June 2002
  * [13] A.M. Sukhov, T.G. Sultanov, M.V. Strizhov, A.P. Platonov, Throughput metrics and packet delay in TCP/IP networks, RIPE59 Meeting, Lisbon, 2009; arXiv:0907.3710
  * [14] Ripe Test Box, http://ripe.net/projects/ttm/
  * [15] Y. Tipsuwan and M.-Y. Chow, ”Control methodologies in networked control systems,” Control Engineering Practice, vol. 11, pp. 1099-1111, 2003
  * [16] A. Varga, The OMNeT++ distrete event simulation system, Software on-line: http://whale.hit.bme.hu/omnetpp/, 1999
  * [17] Zampieri S., “Trends in networked control systems,” _Proceedings of the 17th World Congress, The International federation of Automatic Control_ , July 2008
  * [18] W. Zhang and M.S. Branicky and S.M. Phillips, Stability of Networked Control Systems, IEEE Control System Magazine, 21(1):84-99, February 2001