Monitoring of rotational period variations in magnetic chemically peculiar stars

TL;DR

This paper introduces a new method for monitoring rotational period variations in magnetic chemically peculiar stars, revealing real changes in their angular velocities and providing insights into their magnetic and rotational behaviors.

Contribution

The paper presents a novel, generally applicable period analysis method that efficiently detects gradual period changes using all available observational data with phase information.

Findings

Confirmed period variations in all studied stars.

Demonstrated that period changes reflect real angular velocity variations.

Highlighted the inhomogeneous nature of the star group and potential different causes for period variations.

Abstract

A majority part of magnetic chemically peculiar (mCP) stars of the upper main sequence exhibits strictly periodic light, magnetic, radio, and spectral variations that can be fully explained by the model of a rigidly rotating main-sequence star with persistent surface structures and stable global magnetic field frozen into the body of the star. Nevertheless, there is an inhomogeneous group consisting of a few mCP stars whose rotation periods vary on timescales of decades, while the shapes of their phase curves remain nearly unchanged. Alternations in the rotational period variations, proven in the case of some of them, offer new insight on this theoretically unpredicted phenomenon. We present a novel and generally applicable method of period analysis based on the simultaneous exploitation of all available observational data containing phase information. This phenomenological method can monitor gradual changes in the observed instantaneous period very efficiently and reliably. We present up to date results of the period monitoring of V901~Ori, CU~Vir, $\sigma$\,Ori\,E, and BS~Cir, known to be mCP stars changing their observed periods and discuss the physics of this unusual behaviour. To compare the period behavior of those stars, we treated their data with an orthogonal polynomial model, which was specifically developed for this purpose. We confirmed period variations in all stars and showed that they reflect real changes in the angular velocity of outer layers of the stars, fastened by their global magnetic fields. However, the nature of the observed rotational instabilities has remained elusive up to now. The discussed group of mCP stars is inhomogeneous to such extent that each of the stars may experience a different cause for its period variations.

Tables

Table 1. Table 1: The comparison of some physical characteristics of CU Vir, V901 Ori, σ 𝜎 \sigma  Ori E, and BS Cir  as mentioned in Sect.  5 and the results of their orthogonal polynomial model of the phase functions (see Sect.  4.1 ).

Name	CU Vir	V901 Ori	$\sigma$ Ori E	BS Cir
$T_{\mathrm{eff}}$ [K]	13 400	23 000	23 000	8 800
sp. type	B8p, SiCr	B2p, He-strong	B2ep, He-strong	A4p, CrEu
age [Myr]	90	3	1-2	510
${\tilde{M}}_0$ [d]	2 451 217.2306(3)	2 453 348.710(3)	2454296.5980(2)	2453941.3871(12)
$\tilde{P}$ [d]	0\fd520 700 95(15)	1\fd538 728 6(5)	1\fd190 836 98(6)	2\fd0428486(7)
${\tilde{P}}^{\prime }$	$2.796\left(4\right)\times {10}^{-9}$	$1.12\left(3\right)\times {10}^{-8}$	$3.08\left(3\right)\times {10}^{-9}$	$5.4\left(4\right)\times {10}^{-9}$
${\tilde{P}}^{\prime \prime }$[d-1]	$-7.3\left(3\right)\times {10}^{-14}$	$-3.05\left(25\right)\times {10}^{-12}$	[$3\left(1\right)\times {10}^{-13}$]	[$1\left(3\right)\times {10}^{-13}$]
${\tilde{P}}^{\prime \prime \prime }$ [d-2]	$-2.620\left(16\right)\times {10}^{-16}$	$\left[1\left(3\right)\times {10}^{-16}\right]$	[$-2\left(1\right)\times {10}^{-16}$]	[$4\left(4\right)\times {10}^{-16}$]
${\tilde{P}}^{\prime \prime \prime \prime }$ [d-3]	$2.67\left(15\right)\times {10}^{-20}$
${\tilde{P}}^{\prime \prime \prime \prime \prime }$ [d-4]	$[8\left(3\right)\times {10}^{-24}$]
$\tilde{P}/{\tilde{P}}^{\prime }$ [yr]	$5.1\times {10}^5$	$3.8\times {10}^5$	$1.1\times {10}^6$	$1.1\times {10}^6$
$n/n_{\mathrm{phot}}$	18617/17922	4394/3656	27656/27373	13591/13591
$\sigma$ [mag]	0.0056	0.0050	0.0036	0.018
$g$	276	73	56	57
$n/{\sigma }^2/{10}^6$	600	180	2100	42
obs. int.	1949-2016	1976-2016	1974-2017	1975-2014
$n_{\mathrm{rev}}$	47 225	9556	12 915	6514
ph. ampl.	0.62	$>0.51$	$>0.26$	$>0.11$
${\beta }_{20}$	5.734$\times {10}^7$	$9.111\times {10}^6$	$1.120\times {10}^6$	$3.256\times {10}^6$
${\beta }_{21}$	-1.002$\times {10}^4$	$-4.486\times {10}^3$	$-9.418\times {10}^3$	$-3.001\times {10}^3$
${\beta }_{30}$	6.765$\times {10}^{11}$	$8.060\times {10}^9$	$-2.284\times {10}^9$	$6.139\times {10}^9$
${\beta }_{31}$	1.5570$\times {10}^8$	$1.005\times {10}^7$	$2.198\times {10}^7$	$-6.706\times {10}^8$
${\beta }_{32}$	-2.1819$\times {10}^4$	$-5.371\times {10}^3$	$-7.378\times {10}^3$	$-4.887\times {10}^3$
${\beta }_{40}$	3.025$\times {10}^{15}$	$-4.001\times {10}^{13}$	$-8.834\times {10}^{12}$	$3.961\times {10}^{11}$
${\beta }_{41}$	3.5026$\times {10}^{12}$	$4.055\times {10}^{10}$	$7.873\times {10}^{10}$	$1.3518\times {10}^{10}$
${\beta }_{42}$	-1.979$\times {10}^7$	$1.035\times {10}^7$	$6.940\times {10}^6$	$-5.4553\times {10}^6$
${\beta }_{43}$	$-3.4066\times {10}^4$	$-6.283\times {10}^3$	$-6.4403\times {10}^3$	$1.120\times {10}^6$
${\beta }_{50}$	7.999$\times {10}^{20}$
${\beta }_{51}$	-3.1565$\times {10}^{16}$
${\beta }_{52}$	-3.1346$\times {10}^{13}$
${\beta }_{53}$	-9.4321$\times {10}^7$
${\beta }_{54}$	-3.9449$\times {10}^4$

Equations

$$\vartheta =E+\phi ;\phi =\mathrm{FP}(\vartheta );E=\mathrm{IP}(\vartheta ),$$

$$\frac{\mathrm{d}\vartheta (t)}{\mathrm{d}t}=\frac{1}{P(t)}\vartheta (M_0)=0;\vartheta ={\int }_{M_0}^t\frac{\mathrm{d}\tau }{P(\tau )};$$

$$\frac{\mathrm{d}t(\vartheta )}{\mathrm{d}\vartheta }=P(\vartheta );t(\vartheta )=M_0+{\int }_0^{\vartheta }P(\zeta )\mathrm{d}\zeta ,$$

$${\vartheta }_0(t)={\int }_{M_0}^t\frac{\mathrm{d}\tau }{P_0}=\frac{t-M_0}{P_0};t({\vartheta }_0)=M_0+P_0{\vartheta }_0;\frac{\mathrm{d}t}{\mathrm{d}{\vartheta }_0}=P_0.$$

$$P({\vartheta }_0)=P_0{\left[\frac{\mathrm{d}\vartheta ({\vartheta }_0)}{\mathrm{d}{\vartheta }_0}\right]}^{-1};P(\vartheta )=P_0\frac{\mathrm{d}{\vartheta }_0(\vartheta )}{\mathrm{d}\vartheta }.$$

$$\dot{\nu }(t)=-K{\nu }^q\Rightarrow \dot{P}(t)=K{P}^{2-q},$$

$$K=-\frac{\dot{\nu }}{{\nu }^q}=-\frac{{\dot{\nu }}_0}{{\nu }_0^q}=\frac{\dot{P}}{P^{2-q}}=\frac{{\dot{P}}_0}{P_0^{2-q}}$$

$$\dot{\nu }={\dot{\nu }}_0{\left(\frac{\nu }{{\nu }_0}\right)}^q,\frac{\mathrm{d}\nu }{\mathrm{d}{\vartheta }_1}=\frac{{\dot{\nu }}_0}{{\nu }_0}{\left(\frac{\nu }{{\nu }_0}\right)}^q,\ddot{\nu }=q{\left(\frac{{\dot{\nu }}_0}{{\nu }_0^q}\right)}^2{\nu }^{2q-1},\Rightarrow \frac{\nu \ddot{\nu }}{{\dot{\nu }}^2}=q,$$

$$\displaystyle\displaystyle\dot{P}=\dot{P}_{0}\left(\frac{P}{P_{0}}\right)^{2-q},\ \frac{\mathrm{d}P}{\mathrm{d}\vartheta_{1}}=\dot{P}_{0}\,P_{0}\left(\frac{P}{P_{0}}\right)^{2-q},\$$

$$\ddot{P}=(2-q){\left(\frac{\dot{P_0}}{P_0^{2-q}}\right)}^2{P}^{3-2q},\frac{P\ddot{P}}{{\dot{P}}^2}=\frac{P_0{\ddot{P}}_0}{{\dot{P}}_0^2}=2-q,$$

$$\vartheta ({\vartheta }_0)\doteq {\vartheta }_0-\frac{{\dot{P}}_0}{2}{\vartheta }_0^2+\frac{q{\dot{P}}_0^2}{6}{\vartheta }_0^3,{\vartheta }_0(\vartheta )\doteq \vartheta +\frac{{\dot{P}}_0}{2}{\vartheta }^2+\frac{(3-q){\dot{P}}_0^2}{6}{\vartheta }^3.$$

$$\vartheta ({\vartheta }_0)={\vartheta }_0-\Delta ({\vartheta }_0);{\vartheta }_0=\vartheta +\Delta ({\vartheta }_0)\doteq \vartheta +\Delta (\vartheta )+\frac{1}{2}\frac{\mathrm{d}{\Delta }^2}{\mathrm{d}\vartheta }+\dots ;$$

$$P({\vartheta }_0)=P_0{\left[\frac{\mathrm{d}\vartheta ({\vartheta }_0)}{\mathrm{d}{\vartheta }_0}\right]}^{-1}=\frac{P_0}{1-\frac{\mathrm{d}\Delta }{\mathrm{d}{\vartheta }_0}}\doteq P_0\left[1+\frac{\mathrm{d}\Delta }{\mathrm{d}{\vartheta }_0}+{\left(\frac{\mathrm{d}\Delta }{\mathrm{d}{\vartheta }_0}\right)}^2+\dots \right];$$

$$P(\vartheta )\doteq P_0\left[1+\frac{\mathrm{d}\Delta }{\mathrm{d}\vartheta }+\frac{1}{2}\frac{{\mathrm{d}}^2{\Delta }^2}{\mathrm{d}{\vartheta }^2}+\dots \right].$$

$$\vartheta ({\vartheta }_0)={\vartheta }_0-\Delta ({\vartheta }_0);{\vartheta }_0(\vartheta )=\vartheta +\Delta (\vartheta );\Theta (E)=M_0+P_0\left[E+\Delta (E)\right];$$

$$P({\vartheta }_0)\doteq P_0\left(1+\frac{\mathrm{d}\Delta }{\mathrm{d}{\vartheta }_0}\right)\doteq P(\vartheta )=P_0\left(1+\frac{\mathrm{d}\Delta }{\mathrm{d}\vartheta }\right)=P(E)=P_0\left(1+\frac{\delta \Delta }{\delta E}\right),$$

$$\Delta ({\vartheta }_0)=\frac{1}{2!}{\dot{P}}_0{\vartheta }_0^2+\frac{1}{3!}{P}_0{\ddot{P}}_0{\vartheta }_0^3+\dots +\frac{1}{(k+1)!}{P}_0^{k-1}\frac{{\mathrm{d}}^k{P}_0}{\mathrm{d}{t}^k}{\vartheta }_0^{k+1}+\dots ;$$

$$\vartheta ({\vartheta }_0)={\vartheta }_0-\Delta ;{\vartheta }_0(\vartheta )\doteq \vartheta +\Delta ;\Theta (E)\doteq M_0+P_0[E+\Delta (E)];$$

$$P({\vartheta }_0)=P_0\left(1+{\dot{P}}_0{\vartheta }_0+\frac{1}{2}{P}_0{\ddot{P}}_0{\vartheta }_0^2+\dots +\frac{1}{k!}{P}_0^{k-1}\frac{{\mathrm{d}}^k{P}_0}{\mathrm{d}{t}^k}{\vartheta }_0^k+\dots \right);$$

$$\dot{P}({\vartheta }_0)=\frac{\mathrm{d}P}{\mathrm{d}{\vartheta }_0}\frac{\mathrm{d}{\vartheta }_0}{\mathrm{d}t}={\dot{P}}_0+P_0{\ddot{P}}_0{\vartheta }_0+\dots +\frac{1}{(k-1)!}{P}_0^{k-1}\frac{{\mathrm{d}}^k{P}_0}{\mathrm{d}{t}^k}{\vartheta }_0^{k-1}+\dots$$

$$\ddot{P}({\vartheta }_0)=\frac{\mathrm{d}\dot{P}}{\mathrm{d}{\vartheta }_0}\frac{\mathrm{d}{\vartheta }_0}{\mathrm{d}t}={\ddot{P}}_0+\dots +\frac{1}{(k-2)!}{P}_0^{k-2}\frac{{\mathrm{d}}^k{P}_0}{\mathrm{d}{t}^k}{\vartheta }_0^{k-2}+\dots$$

$$t=M_0+P_0{\vartheta }_0;\varphi (t,{\vartheta }_0)=\frac{t-T_0}{\Pi }=\frac{P_0{\vartheta }_0+M_0-T_0}{\Pi };\Delta (\varphi )=\Delta ({\vartheta }_0);$$

$$\vartheta ={\vartheta }_0-\Delta ({\vartheta }_0);{\Theta }_0(E)=M_0+P_{0}E;\Theta (E)={\Theta }_0(E)+P_0\Delta (E);$$

$$P({\vartheta }_0)=P_0\left(1+\frac{\mathrm{d}\Delta }{\mathrm{d}{\vartheta }_0}\right)=P_0\left(1+\frac{\mathrm{d}\Delta }{\mathrm{d}\varphi }\frac{\mathrm{d}\varphi }{\mathrm{d}{\vartheta }_0}\right)=P_0\left(1+\frac{P_0}{\Pi }\frac{\mathrm{d}\Delta }{\mathrm{d}\varphi }\right).$$

$${\theta }_0={\psi }_0;{\theta }_1={\psi }_1-{\alpha }_{10}{\psi }_0;{\alpha }_{10}=\frac{\overline{{\psi }_1{\theta }_0}}{\overline{{\theta }_0^2}};$$

$${\theta }_2={\psi }_2-{\alpha }_{20}{\theta }_0-{\alpha }_{21}{\theta }_1;{\alpha }_{20}=\frac{\overline{{\psi }_2{\theta }_0}}{\overline{{\theta }_0^2}};{\alpha }_{21}=\frac{\overline{{\psi }_2{\theta }_1}}{\overline{{\theta }_1^2}};$$

$${\theta }_3={\psi }_3-{\alpha }_{30}{\theta }_0-{\alpha }_{31}{\theta }_1-{\alpha }_{32}{\theta }_2;{\alpha }_{30}=\frac{\overline{{\psi }_3{\theta }_0}}{\overline{{\theta }_0^2}};{\alpha }_{31}=\frac{\overline{{\psi }_3{\theta }_1}}{\overline{{\theta }_1^2}};{\alpha }_{32}=\frac{\overline{{\psi }_3{\theta }_2}}{\overline{{\theta }_2^2}};$$

$${\theta }_j={\psi }_j-\sum _{k=0}^{j-1}{\alpha }_{jk}{\theta }_k;{\alpha }_{jk}=\frac{\overline{{\psi }_j{\theta }_k}}{\overline{{\theta }_k^2}}.$$

$${\theta }_0={\psi }_0;{\theta }_1={\psi }_1+{\beta }_{10}{\psi }_0,{\beta }_{10}={\alpha }_{10};$$

$${\theta }_2={\psi }_2-{\beta }_{21}{\psi }_1-{\beta }_{20}{\psi }_0,{\beta }_{21}={\alpha }_{21},{\beta }_{20}={\alpha }_{20}-{\alpha }_{21}{\alpha }_{10};$$

$${\theta }_3={\psi }_3-{\beta }_{32}{\psi }_2-{\beta }_{31}{\psi }_1-{\beta }_{30}{\psi }_0;{\beta }_{32}={\alpha }_{32},{\beta }_{31}={\alpha }_{31}-{\alpha }_{32}{\alpha }_{21},$$

$${\beta }_{30}={\alpha }_{30}-{\alpha }_{31}{\alpha }_{10}-{\alpha }_{32}{\alpha }_{20}+{\alpha }_{32}{\alpha }_{21}{\alpha }_{10};$$

$${\theta }_4={\psi }_4-{\beta }_{43}{\psi }_3-{\beta }_{42}{\psi }_2-{\beta }_{41}{\psi }_1-{\beta }_{40}{\psi }_0;{\beta }_{43}={\alpha }_{43},{\beta }_{42}={\alpha }_{42}-{\alpha }_{43}{\alpha }_{32},$$

$${\beta }_{41}={\alpha }_{41}-{\alpha }_{42}{\alpha }_{21}-{\alpha }_{43}{\alpha }_{31}+{\alpha }_{43}{\alpha }_{32}{\alpha }_{21},$$

$${\beta }_{40}={\alpha }_{40}-{\alpha }_{43}{\alpha }_{30}-{\alpha }_{42}{\alpha }_{20}-{\alpha }_{41}{\alpha }_{10}+{\alpha }_{43}{\alpha }_{32}{\alpha }_{20}+{\alpha }_{42}{\alpha }_{21}{\alpha }_{10}-{\alpha }_{43}{\alpha }_{32}{\alpha }_{21}{\alpha }_{10},$$

$${\vartheta }_0=\frac{t-{\tilde{M}}_0}{P_0};\vartheta ={\vartheta }_0-A\left[\cos (2\pi \varphi )-\frac{\overline{{\vartheta }_0\cos (2\pi \varphi )}}{\overline{{\vartheta }_0^2}}{\vartheta }_0-\overline{\cos (2\pi \varphi )}\right].$$

$$m(\vartheta )=m_0+A_1\cos (2\pi {\phi }_m)+A_2\cos (4\pi {\phi }_m)+A_3\cos (6\pi {\phi }_m)+$$

$$A_4\left[2\sin (2\pi {\phi }_m)-\sin (4\pi {\phi }_m)\right]+A_5\left[3\sin (2\pi {\phi }_m)+6\sin (4\pi {\phi }_m)-5\sin (6\pi {\phi }_m)\right],$$

$$\mathrm{where}{\phi }_m=(\vartheta -{\phi }_{0m})-\mathrm{round}(\vartheta -{\phi }_{0m}).$$

$$m(\lambda ,\phi )=m_0(\lambda )+\sum _{j=1}^{n_{\mathrm{s}}}{A}_j(\lambda )\left\{\exp \left[1-\cosh \left(\frac{\Delta {\phi }_j}{d_j}\right)\right]-2.289d_j\right\};$$

$$\mathrm{where}\Delta {\phi }_j=(\vartheta -{\phi }_{0j})-\mathrm{round}(\vartheta -{\phi }_{0j}).$$

$$W_j(\vartheta )={\overline{W}}_j\left[1+A_j{f}_{\mathrm{ion}}(\vartheta )\right];\mathrm{where}$$

$$f_{\mathrm{ion}}(\vartheta )=\frac{{\sum }_{k=1}^3\left[{\beta }_{2\cdot k-1}\cos (2\pi k\vartheta )+{\beta }_{2\cdot k}\sin (2\pi k\vartheta )\right]}{\sqrt{{\sum }_{l=1}^6{\beta }_l^2}},$$

Full text

††volume: 46

1 Introduction

It has been known since the times of the first systematic spectral classification that about 10% of the upper main sequence stars show peculiar spectra with enhanced or extinct absorption lines of some chemical elements that indicate their overabundance/underabundance in respect to solar chemical composition. Some of those chemically peculiar (CP) stars have also been identified as variable stars exhibiting periodic moderate light variability with amplitudes up to one-tenth of the magnitude accompanied by variations in the intensities and profiles of spectral lines of some elements (as a rule overabundant ones). Spectropolarimetry of the variable CP stars has subsequently shown that the majority of variable CP stars have a strong, stable, nearly dipolar magnetic field, whose axis does not coincide with the rotational one.

The magnetic chemically peculiar (mCP) stars are divided into several subtypes, namely classical Ap and Bp stars with an overabundance of the iron peak elements, silicon, strontium or rare earths, and hot mCPs with underabundance/overabundance of helium (He-weak and He-strong stars). The overabundant chemical elements in their atmospheres usually concentrate into large spots, stable for decades or centuries. The unequal distribution of chemical elements on the surface influences the horizontal inhomogeneity of their atmospheric structure resulting into the incidence of extended photometric spots with uneven energy distributions in their spectra (Krtička et al., 2007; Shulyak et al., 2010; Krtička et al., 2012). As the star rotates, periodic variations in the brightness, spectrum, and magnetic field are observed. The period of the observed variations is equal to rotational period. Combining both present and archival observations of mCP stars collected over the past several decades, we can reconstruct their rotational evolution with unprecedented accuracy.

2 Models for period monitoring

During the last two decades, the research team around the Department of Theoretical Physics and Astrophysics of the Masaryk University in Brno has developed several versatile instruments for the analysis of the periodicity of more or less periodically variable objects, using the methods of phenomenological modelling. In this section, we briefly outline the period analysis technique apt for an investigation of the period stability of magnetic chemically peculiar stars.

2.1 Phase function and its models

Most of the variations in periodical variable stars are cyclic with an instantaneous period $P(t)$, which is strictly constant or slightly variable with time $t$. The period itself and its progression over time cannot be observed directly, but both can be derived through analysing time series of light changes or extremum timings. For that purpose we introduced (see Mikulášek et al., 2008; Mikulášek, 2015)) a monotonically rising phase function $\vartheta(t)$ as a sum of the epoch $E(t)$ and the common phase $\varphi(t)$ and its inversion function $t(\vartheta)$,

[TABLE]

where $\mathrm{IP}(x)$ is an operator rounding $x$ to the nearest integer less than or equal to $x$, while $\mathrm{FP}(x)=x-\mathrm{IP}(x)$. The phase function equals zero at the time of the initial epoch, denoted usually by $M_{0}$, hence $\vartheta(t=M_{0})=0$. Using the inversion function $t(\vartheta)$ we can find for any value of the phase function the corresponding time. The discrete form of this function $t(E)$ predicts the time of the zero phase appertained to the particular epoch $E$.

Functions $\vartheta(t)$ and $t(\vartheta)$ are tied with the instantaneous observed periods $P(t),\,P(\vartheta)$ by the following differential equations (for details see Mikulášek, 2015)

[TABLE]

First of all, we remind that the optimal option of the model of the phase function $\vartheta(t)$ (and its inversion function $t(\vartheta)$) is crucial for the whole period analysis. The parameters of the model are iteratively determined by the process of modeling of the observed behavior of the variable star. The time dependence of the instantaneous period $P(t)$ or $P(\vartheta)$ is then a function derived from $\vartheta(t)$ according to general relations given by Eq. (2). So, the function $P(t)$ strongly depends on the chosen model of the phase function and its adequacy.

The basic and at the same time the simplest model of the variability supposes that the observed period of variations is punctually constant $P(t)=P(\vartheta)=P(\mathit{\Theta})=P_{0}$. Using Eq. (2) we get the corresponding, linear phase function $\vartheta_{0}(t)$ and its inversion $t(\vartheta_{0})$, described by two parameters $(M_{0},P_{0})$:

[TABLE]

Because $\vartheta_{0}$ is a linear function of time, we can use it instead of time. Combining Eq. (2) and Eq. (3) we obtain the useful relations

[TABLE]

2.1.1 Power law models

Several mechanisms, as a steady angular momentum loss via stellar winds or magnetic breaking, subsequently change the instantaneous periods of astrophysical objects $P(t)$, or their frequencies $\nu(t)=1/P(t)$, according to the so-called power law, as follows

[TABLE]

where $q$ is the so-called deceleration parameter, symptomatic for the dominating mechanism of the period change and $K$ is a proportionality constant characteristic of a particular object. If $K>0$, the frequency declines while the period is rising.

[TABLE]

where $\dot{P}_{0}\ \mathrm{and}\ \ddot{P}_{0}$ are the first at $t=M_{0}$.

Because the secular changes are typically slow, we can approximate the phase function $\vartheta(\vartheta_{0})$ and its inversion $\vartheta_{0}(\vartheta)$ using the Maclaurin expansion up to the cubic term. After some algebra, we obtain the following relations,

[TABLE]

The physically interesting parameter $q$ is presented only in cubic terms with the connection to the square of the period progression $\dot{P}_{0}^{2}$, which is almost always negligible111V901 Ori shows, among other period changing mCPs, a record value of rotational braking $\dot{P}=1.0\times 10^{-8}$. The cubic correction of the $\mathit{\mathit{\Theta}(E)}$ in the time of maxima even during 100 years does not exceed 1 minute, while the uncertainty of its determination from observations is always larger than 10 minutes.. Consequently, we have to reconcile the impossibility to determine the value of the deceleration parameter from observations and to reveal the background of the mechanisms causing a constant rise of the period. However, there are objects with rotational variability displaying cubic and higher terms. Nevertheless, these variations in their periods have to be caused by mechanisms not obeying a power law.

2.1.2 Standard Mclaurin polynomials

Let us admit that the observed phase curves display systematical phase (or time) shifts versus their linear model prediction amounting to the $-\Delta(\vartheta_{0})$ or $P_{0}\,\Delta(\vartheta)$ difference. We can interpret this as a result of the inconstancy of the instantaneous period $P(\vartheta_{0})$ or $P(\vartheta)$.

[TABLE]

Then, we can simply write:

[TABLE]

where $\mathit{\Theta}(E)$ is a model prediction of the time of the zeroth phase $(\varphi=0)$ corresponding to the epoch $E$.

In particular, if we express $P(t)$ by means of the Mclaurin expansion where $\dot{P}_{0}\ \mathrm{and}\ \ddot{P}_{0}$ are the first and second time derivatives at the time of the origin $t=M_{0}$, which may be fixed. Then

[TABLE]

2.1.3 Cyclic variation of the period

The modulation amounting to $\Delta(\vartheta_{0})$ of the ideal linear phase function $\vartheta_{0}(t)$ may represent a curving of it, expressing a subsequent change of the period, or may be cyclic222Note that the possible variations in the observed period of mCP star variations need not necessarily mean changes in the rotational period itself. They can also be caused by the inconstant radial velocity of the star as the result of orbital motion in a stellar system. Orbital motion results in some undulation of the basic phase function causing the light-time effect (for details see Liška et al., 2016)., with a proper period $\mathit{\Pi}$, different from the basic, in our case rotational period $P(\vartheta_{0})$.

Let us assume that the difference $\Delta(\phi)$ is a periodic function of the time dependent variable $\phi(t)=(t-T_{0})/\mathit{\Pi}$, where $T_{0}$ is the origin of counting of cycles with the period $\mathit{\Pi}$. Then, we can proceed according to the following algebra

[TABLE]

If we are allowed to opt for the time of the origin of the cyclic phase function $\phi$, $T_{0}$, and the definition of the cyclic function $\Delta(\phi)$, it is advantageous to choose them so that: $\Delta(\vartheta_{=}0)=0$ and $\dot{\Delta}(\vartheta_{=}0)=0$, then $\vartheta(t=M_{0})=0$ and $P(t=M_{0})=P_{0}$. Please also refer to Sect. 4.2, where we apply the set of relations (2.1.3) to the simple sinusoidal modulation of the basic phase function.

2.1.4 Orthogonalized models

The most widely used models of phase function $\vartheta$ are expressed as linear combinations of a basic set of functions of $\vartheta_{0}$: $\left\{\psi_{0},\,\psi_{1},\,\psi_{2},\,\psi_{3},\ldots\right\}$. For several good reasons (see e. g. Mikulášek, 2007b; Mikulášek et al., 2008), it is advantageous to switch from that to another set of functions $\left\{\theta_{0},\,\theta_{1},\,\theta_{2},\,\theta_{3},\ldots\right\}$ which are mutually orthogonal in the space of measurements. This means that the weighted mean value of the products of each two uneven functions is equal to zero: $\overline{\theta_{j}\,\theta_{k}}=0,\ \mathrm{for}\ j\neq k$. The functions of this set can be created iteratively by means of the well-known Gram-Schmidt procedure as the linear combination of the former set $\left\{\psi_{0},\,\psi_{1},\,\psi_{2},\,\psi_{3},\ldots\right\}$ as follows

[TABLE]

Although this way of expressing of newly created orthogonal functions is suitable for computation, it would be more illustrative to rewrite the relations for particular orthogonal terms $\theta_{j}$ using solely the functions $\psi_{k}$, where $k=(0,1,2,\ldots,j)$.

[TABLE]

and so on. The general repetitive orthogonalization routine, as described by Eq. (2.1.4), can be applied to the standard polynomial set of $\vartheta_{0}$, as it has been shown in Sect. 4.1.

Similarly, we can also orthogonalize models of the phase function $\vartheta$ containing non-polynomial terms, including the cyclical ones (see Sect. 2.1.3). The simplest possibility is a model with the constant basic period $P_{0}$ and some phase cosine-like modulation as a function of the phase $\phi$ (see Eq. (2.1.3)) with the set of elementary functions $\{1,\ \vartheta_{0},\ \cos(2\pi\phi)\}$. Assuming that $\tilde{M}_{0}$ has been chosen so that $\overline{\vartheta_{0}}=0$, we arrive at the following model

[TABLE]

2.2 Phenomenological models of mCP phase curves

The development of the period $P(t)$ of a mCP star can be derived from the mutual phase shifts among phase curves obtained at different times, where $\Delta\varphi=\varphi(t_{2})-\varphi(t_{1})=\vartheta(t_{2})-\vartheta(t_{1})-[\vartheta_{0}(t_{2})-\vartheta_{0}(t_{1})]$. For the determination of the phase shifts, we use the fact that the phase curves of mCP stars remain invariable for many decades. To obtain a maximum interval of observations we are forced to combine all available phase curve of varying characteristics of the star. It is advantageous to use proper phenomenological models of phase curves here, describing them accurately, using a minimum of free parameters.

The following review of the phase curve models is only rough and incomplete because the observed mCP star variations are in reality extremely rich and diverse. Therefore, it is common to tailor the models differently for each star.

2.2.1 Light curves

The period analyses of mCP stars are based mainly on photometry. Monochromatic light curves of mCP stars are the result of the presence of dull photometric spots on the surface of rotating stars. The resulting light curves are smooth and relatively simple. A majority of monochromatic light curves can be well approximated by the harmonic polynomial of the second or third order (Mikulášek et al., 2007b; Jagelka & Mikulášek, 2015). It is advantageous to centre the light curve into the phase of light maximum $\varphi_{m}$, which uses to be sharper than the minimum/minima. The model of a monochromatic light curve then has maximally five phase dependent components – three symmetric and two antisymmetric ones:

[TABLE]

However, the shapes of the spectral energy distribution of individual photometric spots are often dissimilar, which results in the fact that light curves in various colours are not the same; they differ both in their amplitudes and forms. Fortunately, the usage of the PCA or APCA technique (Mikulášek et al., 2004; Mikulášek, 2007c) considerably diminishes the number of parameters indispensable for good description of light curves obtained in all colours.

There is also another, more physical approach to the modelling of mCP light curves in various colors. It assumes that light curves in different colors can be expressed as a linear combination of a few symmetric basic profiles with that centers at phases $\vartheta_{0j}$ and half-widths $d_{j}$:

[TABLE]

The above mentioned phenomenological model applies to almost all mCP stars observed from ground-based observatories (see e. g. Mikulášek et al., 2015a; Krtička et al., 2015; Jagelka & Mikulášek, 2015). Nevertheless, $\sigma$ Ori E  forms an exception. Its light curves in the optical region can not be explained by the model of a rotating star with photometric spots (Oksala et al., 2015a, b, and references therein).

Most of the photometric data on the here studied mCP stars are freely accessible through the database On line catalogue of photometric observations of magnetic chemically peculiar stars, http://astro.physics.muni.cz/mcpod/ (Mikulášek et al., 2007c; Janík et al., 2011).

2.2.2 Magnetic fields

Magnetic fields of mCP stars are structurally much simpler, and mainly much stronger, than the fields of cool stars (Donati & Landstreet, 2009). The large-scale strength and geometry of the magnetic fields are stable, in the rotating stellar reference frame, on time scales of many decades (e. g. Wade et al., 2000; Silvester et al., 2014). Magnetic field observations are then very reliable subsidiary source of the phase information on the rotation of mCP stars. A major part of practically exploited data is time-series of the observations of the mean longitudinal component $B_{\rm{eff}}$ of the total magnetic field. They were derived from measurements of circular polarization induced in magnetically-split spectral line $\sigma$ components due to the longitudinal Zeeman effect (see, e. g. Mathys, 1988; Donati & Landstreet, 2009).

The magnetic fields of most of the mCP stars are more or less dipole-like with the dipole axis tilted to the rotational axis. The phase curves of the mean longitudinal component of the field $B_{\rm{eff}}$ are almost sinusoidal. However, the current spectropolarimetric measurements of high quality show that a quadruple component of the field of many mCP stars has to be considered (Wade et al., 2000). Nevertheless, a majority of magnetic phase curves are single waves, which can be described by a harmonic polynomial of maximally second order. An example of an mCP star harboring an extraordinarily complex magnetic field is one of our target stars – V901 Ori (Thompson & Landstreet, 1985; Kochukhov et al., 2011), the phase curve of which is only adequately described by the symmetric harmonic polynomial of the third order (see Fig. 1).

2.2.3 Spectral line variations

One of the specific features of mCP stars is a very uneven horizontal distribution of chemical elements on the stellar surface, which applies to those elements that are overabundant in respect to the solar abundance. Chemical elements concentrate to vast spots whose projection on the visible disc of the star changes as the star rotates. Owing to this, the spectral absorption lines strongly variable in their intensity and profiles with the rotational period of the star. Changes in equivalent widths of spectral lines of mCP stars are so conspicuous that they were noticed and measured at the very beginning of the study of mCP stars (see e. g. Belopolsky, 1913; Farnsworth, 1932; Deutsch, 1952).

The analysis of phase variation of the line profiles of ions of chemical elements, combined with changes in their polarization by the techniques of magnetic Doppler imaging, enables a successful mapping of the distribution of spectral spots in the atmospheres (see e. g. Rice et al., 1989; Khokhlova et al., 2000; Lüftinger et al., 2010). The pioneer systematic studies of the variability of He I and Si II lines in upper main sequence stars are those of Pedersen (1979), which helped us very much in anchoring our period research in the eighties of the last century. Of a rather limited use are data on radial velocity variations because they depend strongly on the measured spectral line selection and the method of measurements. Relatively reliable are the radial velocity of the hydrogen Balmer lines or other strictly defined spectral features. Nevertheless, the basic phase information hides present and archival equivalent width measurements of selected lines of overabundant chemical elements.

The handling of data of equivalent width time series of spectral lines of various elements, based on spectrograms of very different quality and derived by different techniques, is very demanding and time-consuming. However, well-homogenized data sometimes bear crucial information on the behavior of the phase function in the past and sometimes also in the future333It was also in the case of V901 Ori, where the spectroscopic data helped to reveal that the spin-down is slowing and would alternate into a spin-up (Mikulášek et al., 2008). The present data analysis confirms that this predicted break down moment came to pass in 2009 (see Mikulášek et al., 2011a, and Sect. 4.3)..

Phase curves of equivalent widths of a particular spectral line use to be a smooth, single or double wave, which we could fit by the third order harmonic polynomial, described by Eq. (18). We can also use the following time-tested model for light curves given by Eq. (19), which can pinpoint the phases when the spectral spots pass the stellar meridian. When using measurements of various spectral lines of an ion, the situation is more complex. As a rule, the phase curves are similar, but we must consider their different intensities and unavoidable blends with lines of other elements that, in general, somewhat suppress the relative variations of a particular spectral line. Nevertheless, the following simple model can take the circumstances mentioned above into account and give acceptable results (see Fig. 2)

[TABLE]

where $W_{j}(\vartheta)$ is the model predicted equivalent width of the $j$-th line of the studied ion for a phase function $\vartheta$. $\overline{W}_{j}$ is a mean value of all equivalent widths of the $j$-th line, $A_{j}$ is a relative amplitude of the $j$-th line, and $f_{\rm{ion}}$ is a normalized phase curve belonging to a particular ion of the chemical element (typically, He I), the coefficients of which are the same for all studied lines of the relevant ion.

2.2.4 Times of radio-pulse peaks

The only type of mCP stars’ observations where, instead of the standard phase curves, we used the times of some phase located events are the moments of radio-pulse peaks of CU Vir, known as the first main sequence radio pulsar (see e. g. Trigilio et al., 2000). The model describing this observation is very simple; we need just one parameter – the phase of the center of the peak $\varphi_{\rm{0r}},\ y_{\rm p}=\mathit{\Theta}(E_{i})+\tilde{P}\,\varphi_{\rm{0r}}$, while the measured quantities $y_{i}$ are the given times of pulses.

2.2.5 General remarks and recommendations for the modeling of phase curves

The suitable option of a phase curve model used for the period analysis is a crucial point of variable star modeling. Usage of an inappropriate model may influence the reliability of the complete solution and lead us to faulty conclusions. The neuralgic point of our effort is the fact that typically we are forced to use very inhomogeneous measurements of various nature and quality with phase curves connected to each other only freely. Typically, the phases of extremes of light curves, magnetic field curves, and spectral variations may coincide, but need not (e. g. Silvester et al., 2015). We are allowed to discuss a possible correspondence among the mentioned phase curves and use them always a posteriori, never a priori of the period analysis.

It is also desirable to plan new observations of every kind (i.e. photometric, spectroscopic and spectropolarimetric ones) so that they would be obtained more or less simultaneously (in a time interval of a few years). Then we could better phase the models of all observations with the phase information in use.

We have to pay attention to the correct weighting of entered data because the employed $\chi^{2}$ regression requires it (see Sect. 2.3). If we do not know the individual uncertainties of the original data in advance or if they appear suspect, we have to estimate them iteratively from the scatter of residuals $\{\Delta y_{i}\}$ for appropriately defined data subsets. We recommend eliminating outlier influence using their modified uncertainties by the well-established method given in Mikulášek & Zejda (2013).

The models of phase curves should be tailored to the studied object, available data, and the purpose of fitting the data. In particular, the number of used free parameters should be restricted to as few as possible, but without any serious influence on the accuracy and reliability of the results. Unfortunately, the effort of using the optimal phenomenological models considerably encumbers automation of the computational process. The diversity of real phase curves of particular mCPs requires that they have to be solved individually.

2.3 Solution of models

Our knowledge of the development of the periods of the studied variable stars in time $P(t)$ is derived from the analysis of the course of the phase function $\vartheta(t)$ whose models were dealt with in detail in Sect. 2.1. The parameters of the adopted phase function models are calculated through minimizing $\chi^{2}$ quantity

[TABLE]

where $\{t_{i},\,y_{i},\,\sigma_{i}\}$ is the set of all available $n$ observations of phase dependant quantities $y_{i}$ with uncertainties (or modified uncertainties) $\sigma_{i}$, obtained at $t_{i}$. $y_{\rm p}\left[\vartheta(t_{i})\right]$ is then a model prediction calculated for the particular phase function $\vartheta(t_{i})$ at $t_{i}$. The possible models of photometric, spectropolarimetric and spectroscopic phase curves are briefly described in Sect. 2.2. We standardly assume that the shapes of phase curves are invariable (functions of the phase $\varphi(t)$, only), nevertheless, in justified cases, we can assume their slow progress in time (see Sect. 4.5).

The minimization of $\chi^{2}$ brings $g$ non-linear equations for $g$ unknown free parameters of all used phase curves and the phase function444The number of free parameters $g$ depends mainly on the number of ’material constants’ describing all the phase curves, which could be rather large. For example, for CU Vir with $n=18\,641$ observations we need $g=274$ free parameters, whereas only four of them are necessary for the determination of the phase function. For details, see Sect. 4.2.. The equations are solved simultaneously by the standard iterative Newton-Raphson method, described e. g. in Press et al. (2002); Hartkopf et al. (1989); Mikulášek et al. (2011a), among others.

With a good initial estimate of the parameter vectors, the iterations converge fairly quickly. All estimates of uncertainties of model parameters were computed using the general law of uncertainty propagation also assuming correlations amongst individual coefficients (see e. g. Bevington & Robinson, 2003; Mikulášek & Zejda, 2013). For models of phase functions of all the studied mCP stars, except for CU Vir  which exhibits cyclic changes of its period, we preferentially used the orthogonal polynomials (see Sect. 2.1.4). This allows us to use a simple law of error propagations for evaluating functions of orthogonal polynomial parameters $\tilde{M}_{0},\,\tilde{P},\,\tilde{P}^{\prime},\,\tilde{P}^{\prime\prime}$ (see Eq. (22)).

3 Magnetic CP stars without detectable period changes

As an example of a well observed strictly periodic star, we mention CQ UMa = HR 5153 = HD 119213. This SrCrEu Ap star displays prominent variation in the Strömgren $v$ band with antiphase changes in R and I bands (see Fig. 3a).

Mikulášek et al. (2009), and Mikulášek et al. (2011b) used 1365 observations collected from eleven various sources of photometric data that cover a time interval of 42 years (6262 revolutions of the star). The mean period: $P\!=\!2\fd 449\,912\,0(27)$ can then be derived with the accuracy of 0.23 s. The time derivative of the period is $\dot{P}=1(2)\times 10^{-9}=(3\pm 7)$ s cen*-1*, which means that period is stable, as for most other CP stars.

Mikulášek et al. (2016) also tested the constancy of period of the photometrically revealed mCP star KIC 6 950 556, analyzing 64 793 detrended Kepler observations of an accuracy of 0.13 mmag. They found the rotational period of the star $P_{0}=1\fd 511\,785\,08(4)$ d and the time derivative of the period $\dot{P}=1(2)\times 10^{-10}=(0.3\pm 0.6)$ s cen*-1*. We shall note that the accuracy of the period rate determination of $2\times 10^{-10}$ is quite sufficient to reveal all presently known mCP stars with variable periods.

Wanting to obtain the relevant statistics of the incidence of period unstable stars among the mCP stars population, we will accordingly test all known and suspected mCP stars observed by the Kepler mission.

4 Magnetic CP stars with proven period variations

A few mCP stars might display minor secular changes in the shape of their light curves (see e. g. Žižňovský et al., 2000), which can be attributed to the precession of magnetically distorted stars (Shore & Adelman, 1976; Pyper & Adelman, 2004). However, there is also a small subgroup of mCP stars that have phase stable curves, but exhibit variable rotation periods (Mikulášek et al., 2011b, 2014, 2017, and references therein). The constancy of their light curves on the scale of decades disqualify precession as the cause of the observed period changes (Mikulášek et al., 2008).

The following text is based mainly on the period analyses of four of the best-monitored mCP stars – CU Vir, V901 Ori, $\sigma$ Ori E, and BS Cir, known for their period variability.

4.1 Orthogonal polynomial model of the phase function

As we want to compare and discuss the period development in our four mCP stars, we will use the same model for the period monitoring of all objects. It seems that the most illustrative and versatile model is the set of orthogonalized Mclaurin polynomial (Sect. 2.1.2) up to a fifth order if needed. Then, we can write

[TABLE]

where $\tilde{M}_{0},\,\tilde{P},\,\tilde{P}^{\prime},\,\tilde{P}^{\prime\prime},\tilde{P}^{\prime\prime\prime},\tilde{P}^{\prime\prime\prime\prime}$ are parameters of the orthogonal polynomial fit which have a similar meaning as the analogous parameters $M_{0},\,P_{0},\,\dot{P}_{0},\,\ddot{P}_{0}$, $P_{0}^{\prime\prime\prime}$, and $P_{0}^{\prime\prime\prime\prime}$, introduced in the standard Mclaurin polynomial model (see Sect. 2.1.2). Nevertheless, they are not generally equal, as they have another meaning. While the parameters of the standard model are the parameters of the polynomial expansion at $t=M_{0}$, which is also the initial epoch in our counting system, this means that the real phase function $\vartheta$ is going through this origin and coincides with its Mclaurin model in the close vicinity of the origin. The time-like function $\vartheta_{0}$ then represents the tangent of the phase function at $t=M_{0}$. The orthogonal model parameters $\tilde{M}_{0},\,\tilde{P},\,\tilde{P}^{\prime},\,\tilde{P}^{\prime\prime},\tilde{P}^{\prime\prime\prime},\tilde{P}^{\prime\prime\prime\prime}$ refer to the whole time interval covered by data regarding their particular weights. Thus, they are something like the mean values of the standard parameters. The auxiliary time-like function $\tilde{\vartheta_{0}}$ is in our orthogonal model the first (linear) approximation of the real phase function $\vartheta$. This means that it generally does not pass through the origin at the origin $t=M_{0}$, as well as its parameter $\tilde{P}$ is not always equal to the instantaneous period at the origin $P_{0}$. Therefore, the auxiliary functions $\tilde{\vartheta_{0}}$ and $\vartheta_{0}$ need not to be interchanged.

The six functions of $\tilde{\vartheta_{0}}$ from the set: $\{\theta_{0},\ \theta_{1},\ldots,\,\theta_{5}\}$ are mutually orthogonal on the set of the observational data regarding their weights. So, it is valid that $\overline{\theta_{j}\,\theta_{k}}=0$ for each non-equal pair of $\theta_{j}$ functions.

The functions were created by the Gram-Schmidt procedure (discussed in Sect. 2.1.4) from the base of common polynomials $\{1,\,\tilde{\vartheta_{0}},\,\tilde{\vartheta_{0}}^{2},\,\tilde{\vartheta_{0}}^{3}\ldots\}$, assuming that the origin of our epoch counting for the linear approximation is chosen so that $\tilde{M}_{0}\cong\overline{t}\ \Rightarrow\ \alpha_{10}=0;\ \theta_{1}=\tilde{\vartheta_{0}}$. Then, we can write:

[TABLE]

It is useful to express the set of basic orthogonal functions $\{\theta_{0},\ \theta_{1},\ldots,\,\theta_{5}\}$ as a function of $\tilde{\vartheta_{0}}$, directly

[TABLE]

The computation of derivatives and the derivatives of $\theta_{j}$ is then trivial. Coefficients $\beta_{jk}$ are given by relations

[TABLE]

The numerical results and values of pertinent dimensionless coefficients $\beta_{jk}$ of the orthogonalization for particular stars are given in Table 1.

It may be worthwhile to transform parameters of the orthogonal polynomials back to the standard Mclaurin parameters, here for the polynomials of the fifth order. After some algebra, we obtain the following transforming relations

[TABLE]

4.2 CU Virginis – a silicon mCP star

The famous very fast-rotating silicon mCP star CU Vir (HD 124224, HR 5313) displays an intriguing period variation. It is a common hot Si-type mCP star with a mass 3 M*⊙* and a radius of 2 R*⊙* (Stȩpień, 1998) and T_{\mathrm{eff}}$$=13\,000 K, $\log g=4.0,\ v\,\sin i=160$ km s${}^{-1},\ i=30^{o}$, (Kuschnig et al., 1999). Its nearly dipolar magnetic field with the moderate pole strength of $B_{\rm{p}}=3.0$ kG is tilted towards the rotational axis by $\beta=74^{\circ}$, the axis inclination being $i\simeq 43^{\circ}$ (Trigilio et al., 2000). CU Vir is the first known main sequence star that shows variable radio emission, resembling the radio lighthouse of pulsars (Trigilio et al., 2008, 2011; Ravi et al., 2010). CU Vir also exhibits variations in light and intensities of spectral lines of He I, Si II, H I, and other ions. The nature of its variability in UV and optical regions has been studied by Krtička et al. (2012). CU Vir now belongs to the most frequently and broadly studied mCP stars.

Occasional rapid increases in its rotation period were reported and discussed several times. Pyper et al. (1997, 1998) discovered an abrupt increase of the period from $0\fd 5206778$ to $0\fd 52070854$ that occurred approximately in 1984 and Pyper & Adelman (2004) then discussed two possible scenarios of explanation of the observed O-C diagram, namely a continually changing period or two invariable periods. After 1998, Trigilio et al. (2008, 2011) observed another increase in the period of radio pulses of $\Delta P=1.12$ s by comparing it with the period determined by Pyper et al. (1998). Pyper et al. (2013) presented a very deep study updating the period development of the star based on a significant amount of excellent observational data, namely their precise 2820 Strömgren $\mathit{uvby}$ values obtained from the Four College Automated Photometric Telescope (FCAPT) in the period 1998-2012. Their main result was that the O-C data since 1993 are consistent with a constant period of 0\fd5207137, which is the longest of all periods referred in Pyper et al. (1997, 1998); Pyper & Adelman (2004).

This statement agrees with findings of Mikulášek et al. (2011a) who collected and analyzed all available observations of CU Vir containing phase information between 1949 and 2011. They proved that the shapes of all phase curves were constant during several decades, while the period was continually changing. The rotation period was gradually shortening until the year 1968, when it reached its minimum. The period then started increasing, reaching its local maximum. The following spin-up was recently undoubtedly confirmed by Krtička et al. (2017).

The present period analysis is also based on published FCAPT photometry from 1998-2012 by Pyper et al. (2013) and our measurements from 2011-6. The complete observational material represents now 19 641 individual measurements of CU Vir including photometric measurements in photometric bands from 200 to 753 nm (see in Krtička et al., 2012), as well as spectroscopic, spectropolarimetric and radiometric observations.

Let us assume, in accordance with Krtička et al. (2017), that the difference $\Delta(\phi)$ is a sinusoidal function of $\phi(t)=(t-T_{0})/\mathit{\Pi}$, where $T_{0}$ is the origin of counting of cycles with the period $\mathit{\Pi}$. Applying relations (2.1.3), we obtain

[TABLE]

where $A$ is a semiamplitude of the change $\mathit{\Theta}(E)-\mathit{\Theta}_{0}(E)$ with the minimum at $T_{0}$, the semiamplitude of the mean period undulation being $A_{\mathrm{P}}=2\,\pi AP_{0}/\mathit{\Pi}$. $M_{0}$ was chosen so that $\Delta(\tilde{\vartheta_{0}}=0)=\mathrm{d}\Delta/\mathrm{d}\vartheta_{0}=0$. Analysing all the available observational data of BS Cir, we found $M_{0}=2\,446\,604.4390$ (fixed), $P_{0}=0.520\,694\,04(3)\,\mathrm{d},\ T_{0}=2\,446\,604(13),\ \mathit{\Pi}=24\,110(150)\,\mathrm{d}=66.0\pm 0.4$ yr, $A=0.1611(5)$ d, and $A_{\mathrm{P}}=1.888$ s (see Fig. 5). The employed data cover more than one cycle of the proposed sinusoidal variations.

The most natural explanation for the CU Vir O-C diagram offers the light-time effect. This, however was refused by Pyper et al. (1997, 2013), and Mikulášek et al. (2011a). Therefore we have to find a different explanation, admitting the fact that the observed period variations are due to changes in the angular velocity of the star or its surface layers at least (see also Stȩpień, 1998).

A very ambitious explanation of the observed behavior of CU Vir was proposed by Krtička et al. (2017), introducing a novel mechanism of rotation oscillation as a consequence of internal waves spreading within a rotating magnetic star.

We can also describe the observed phase function of CU Vir by the orthogonal polynomial model of the fifth order (see Sect. 4.1), determined by the set of coefficients $\beta_{jk}$, given by Eq. (23) and (4.1) and 6 parameters of the model (all these results necessary for ephemeris calculations are listed in Table 1). $\tilde{M}_{0}=2\,451\,217.2306$, the mean weighted period is $\tilde{P}=0\fd 520\,700\,95(15)$. The parameters $\tilde{P}^{\prime},\ \tilde{P}^{\prime\prime},\ \tilde{P}^{\prime\prime\prime}$ and $\tilde{P}^{\prime\prime\prime\prime}$ are determined with good accuracy, only the last parameter, $\tilde{P}^{\prime\prime\prime\prime}=8(3)\times 10^{-24}$ d*-4*, is so uncertain that we can neglect it. Moreover, it influences only the very beginning of CU Vir measurements, which at the time were rather unreliable. Fig. 4 displays the contributions of particular terms of the orthogonal polynomial model of a high degree. Nevertheless, it seems that the validity of the sinusoidal model is more or less confirmed, even if we used the model not assuming any cyclicity of period variations.

4.3 V901 Orionis – a He-strong mCP star

V901 Orionis = HD 37 776 is a very young He-strong B2p mCP star with an effective temperature of 23 000 K (Cidale et al., 2007), residing in the emission nebula IC 432, with a global, extraordinarily strong $(B_{s}\approx 20$ kG), and complex magnetic field (Thompson & Landstreet, 1985; Kochukhov et al., 2011). The observed moderate light variations are caused by the spots of overabundant silicon and helium (Krtička et al., 2007).

Thirty years of accurate photometric and spectroscopic monitoring enabled us to reveal a continuous rotational deceleration (Mikulášek et al., 2007a; Mikulášek et al., 2008), increasing the period of about $1\fd 5387$ by a remarkable 18 s! Ruling out (a) a light-time effect in a binary star, (b) the precession of the star’s rotational axis, and (c) evolutionary effects as possible causes of the period change, we interpreted the deceleration in terms of the rotational braking of the outer stellar layers as caused by the angular momentum loss in the stellar magnetosphere.

However, this cannot explain the discrepancy between the spin-down time, $\tau=P/\dot{P}=3\times 10^{5}$ yr, and the star’s age of one million years, or older (Mikulášek et al., 2008, 2011b, 2011a). The interpretation of the rotation period evolution by a simple angular momentum loss was also questioned by the negative value of the second derivative of the period: $\ddot{P}=-29(13)\times 10^{-13}\,\mathrm{d}^{-1}$, which indicated that braking could soon change into acceleration (Mikulášek et al., 2008). New precise measurements prove this conclusion without any doubts.

Presently, we have at our disposal 3656 photometric, 663 spectroscopic, and 75 magnetic measurements covering more or less evenly a time interval of 40 years. All the data were simultaneously modeled with the quartic (fourth order) orthogonal polynomial phase function model $\vartheta(\tilde{\vartheta_{0}})$, described in Sect. 4.1

[TABLE]

where the parameters of the model $\tilde{M}_{0},\ \tilde{P},\ \tilde{P}^{\prime},$ $\tilde{P}^{\prime\prime}$, and $\tilde{P}^{\prime\prime\prime}$ are listed in Table 1, including coefficients $\beta_{jk}$, necessary for ephemeris calculations.

The parameters of the orthogonal model are $\tilde{M}_{0}=2\,453\,348.710(3)$, $\tilde{P}=1\fd 538\,728\,6(5)$, $\tilde{P}^{\prime}=1.12(3)\times 10^{-8}$, $\tilde{P}^{\prime\prime}=-3.05(25)\times 10^{-12}\,\mathrm{d}^{-1}$, and $\tilde{P}^{\prime\prime\prime}=1(3)\times 10^{-16}\,\mathrm{d}^{-2}$. It seems that the last term of the expansion of the phase shift $\Delta(\tilde{\vartheta_{0}})$ (with $\theta_{4}$) can now be completely neglected555We should note that the uncertainty of the quartic term is quickly and proportionally diminishing to $n^{-1/2}\Omega^{-4}\sim\Omega^{-9/2}$ [see (31)], where $n$ is the total number of observations and $\Omega$ is the time of observations of the star. Consequently, in 15 years of constant monitoring, it will drop to $10^{-12}$ d*-2*, which should be enough for solving the question whether the variations of the period are cyclic or not..

According to equations (22) we can predict the course of the phase function $\vartheta(\tilde{\vartheta_{0}})$ the moments of the zeroth phase $\mathit{\Theta}(E)$ and that computed by the linear approximation $\mathit{\Theta}_{0}(E)$ for any epoch $E$ as follows

[TABLE]

as well as the courses of the instantaneous period $P(\tilde{\vartheta_{0}})$ and its derivatives $\dot{P}(\tilde{\vartheta_{0}}),\ \dot{P}(\tilde{\vartheta_{0}})$. Using them we can conclude that the deceleration of the stellar rotation switched in 2009.7(1.0) to acceleration instead (see Figs. 6 and 7).

It is likely that the monitored part of the O-C curve is only a segment of a cyclic curve, but we are not yet able to determine its parameters, like as amplitude and period $\mathit{\Pi}$. We can only estimate that the period must be longer than observed in hundred years. Then the reasons for the period change could be the same as that of CU Vir, but they necessarily need not be (Krtička et al., 2017).

V901 Ori has been followed up for the last forty years as a photometric, spectroscopic and magnetic variable, with each type of variations giving us information about the location of another type of surface structures. Photometric spots, accountable for light variations are determined namely by the distribution of silicon on the star’s surface (Krtička et al., 2007), while spectroscopic variability was studied solely using changes of the He I spectral line intensity (Mikulášek et al., 2008). This informs us about the location of spots with an overabundant helium, which do not coincide with silicon spots (Khokhlova et al., 2000). Observations of the effective magnetic field open up the possibility to study the magnetic geometry of the star (Thompson & Landstreet, 1985; Kochukhov et al., 2011).

If the mentioned structures have been moving mutually during the last decades, we should detect another trend in residuals of the three types of observations. Fig. 7 clearly proves that there were no such motion, so we can conclude that outer layers of the star rotate as a solid body. That fact confirmed the idea that stellar atmospheres of mCP stars with their strong global magnetic fields are horizontally stabilized due to this property.

4.4 $\sigma$ Ori E - a hybrid of mCP and Be star

$\sigma$ Ori E = HD 37479 = V1030 Ori is a hybrid of a classical He-strong mCP star and a Be star with strong stellar winds. It is an extremely young and massive star with a period of photometric and spectral changes $P=1\fd 19801(1)$ (Hesser et al., 1977; Hunger et al., 1989). The spin-down of the star has been noticed by Reiners et al. (2000) and confirmed by Oksala & Townsend (2007). Townsend et al. (2010) then revealed the period to be linearly increasing at the rate $\dot{P}=2.89(21)\times 10^{-9}=0.091(7)$ s yr*-1* using the method of the analysis of five times of primary minimum. The authors explained the observed lengthening of the period by magnetic braking through strong stellar winds. The established spin-down time $\tau=P/\dot{P}=1.4\times 10^{6}$ years was considered to be compatible with the estimated age of the star (Mikulášek et al., 2011b).

We considerably broadened the extent of the observations used by Townsend et al. (2010), added observations from the MOST satellite (courtesy of Townsend et al., 2013), some archival photometric observations, and our own $u$ observations from the winter 2016/17, combining them with our own and archival observations of equivalent width measurements of strong helium lines. In total, we used 27 656 individual measurements, 27 373 being photometric ones, covering a time interval of 43 years (1974–2017). Light curves (see Fig. 8) were modeled by special periodic functions originally developed for eclipsing binaries (Mikulášek, 2015), assuming asymmetry of their minima. Using the whole mentioned material and a more sophisticated method of the phase shifts we were able to enhance the accuracy of the $\dot{P}$ determination 8-times. The full list of the used data and special models of phase curves will be published elsewhere soon.

We give here only the results of the phase function fitted by the model of the quartic orthogonal polynomial, which is formally identical with that used for V901 Ori (see Eq. (29) and (30), and the references around). The parameters and coefficients $\beta_{jk}$ are listed in Table 1.

We found the following model parameters: $\tilde{M}_{0}=2\,454\,296.5980(2)$, the mean period $\tilde{P}=1.190\,836\,98(6),\tilde{P}^{\prime}=3.08(3)\times 10^{-9}=0.0981(9)$ s yr*-1*, and $\tilde{P}^{\prime\prime}=3(1)\times 10^{-13}$ d*-1*. The points in the O-C diagram are excellently fit by a simple parabola (see Fig. 9); the influence of a cubic term, if any, is obviously negligible. The change of the period is almost linear so that it could also be caused by mechanisms of the power law category. $\sigma$ Ori E is the hottest of the stars with varying periods, so its stellar winds could be dense enough to explain its moderate rotational breaking.

Up to 2010, it seemed that all rotationally unstable mCP stars belong to the hotter and thus more massive and younger mCP stars (Mikulášek et al., 2011b). Naturally, the observed affinity may only be the result of the wrong conclusion drawn from a very limited sample of such stars. Therefore we directed our analyses also to short-periodic moderately cool mCP stars, also called as SrCrRE or SiCr stars exhibiting light variations with a large amplitude. All studied moderately cool mCP stars turned out to be stable with one exception described below.

4.5 BS Cir – a moderately cool surprise

BS Cir = HD 125630 = HIP 70346 is a southern mCP star of the A2pSiCr type that is a representative of moderately cool magnetic chemically peculiar stars that display rather strong light variations in Strömgren index $c_{1}=(u-v)-(v-b)$ indicating large changes in the height of the Balmer jump. Combining available kinematic, photometric and spectroscopic data of the star, Mikulášek et al. (2014) derived the following astrophysical parameters: $T_{\mathrm{eff}}=8800\pm 500$ K, $L=41.7\pm 1.4\,\mathrm{L}_{\odot},\ M=2.32\pm 0.14\,\mathrm{M}_{\odot},\ \mathrm{age}=510^{+90}_{-150}$ Myr. A moderate global magnetic field with bipolar strength $B_{\rm{p}}$ of several kG is present (Kochukhov & Bagnulo, 2006; Hubrig et al., 2006).

BS Cir was observed in 1975-6 by Vogt & Faundez (1979) in uvby and then in 1980 by Manfroid & Renson (1980, 1983); Mathys & Manfroid (1985) have revealed the star to be photometrically variable with a period $P=2\fd 205\pm 0\fd 004$ (apparently not taking into account the photometry of Vogt & Faundez, 1979). A data set of Manfroid & Renson (1983) was later reanalysed and the formerly found period was confirmed by Mathys & Manfroid (1985). Catalano & Leone (1993) then added Manfroid & Renson data to their 56 precise $\mathit{uvby}$ measurements taken in 1991 and ‘improved’ the period to $P^{\prime}=2\fd 20552(6)$. As they again did not consider Vogt & Faundez data, they could be mistaken in the total number of cycles between 1980 and 1991 ($\Delta t\simeq 11\times 365\simeq 4015$ d) by one $(\Delta k=-1)$. The incorrect period is also given in the list of Dubath et al. (2011) containing periods and types of variable stars in the Hipparcos survey where BS Cir was treated as an eclipsing binary with the orbital period of 1\fd1020.

The first actual BS Cir period of 2\fd2043, derived from the analysis of all the above mentioned data sets, replenished by Hipparcos, ASAS, Pi of the Sky observations and our own observations in the South Africa Astrophysical Observatory by various photometric instruments in 2011-4, was briefly reported in the review paper of Mikulášek et al. (2014). The ephemeris of BS Cir, based on 14 488 photometric measurements in eleven data sets that cover, more or less evenly, a time interval of 38 years, also contained the square term. It showed undoubtedly that the instantaneous period of the star is rising with the rate $\dot{P}=5.6(4)\times 10^{-9}=0.181(13)\,\mathrm{s\,yr}^{-1}$.

Mikulášek et al. (2009) argued that periodograms of BS Cir, based on all the photometric material, admit the only mean period of 2\fd2042 (see Fig. 11). The nature of the star is then more extensively dealt with in the little monograph of Mikulášek et al. (2015a), based on the rich photometry and also the spectroscopy. However, the period analysis was based solely on the more reliable photometric data. Light curves obtained in uvbyUBVRIHp photometric colors were successfully fitted with the model expressed by Eq. (19) assuming two differently colored, centrally symmetric photometric spots centered at the phases $\varphi_{02}=0.0000(5)$ and $\varphi_{02}=0.4802(11)$, with half-widths of $d_{1}=0.136(2),\ d_{2}=0.117(3)$ (see Fig. 10).

The different spectral energy distribution in the spots shows that at least two diverse mechanisms of the redistribution of the radiative energy are active there. According to Eq. (23) and (4.1), we found coefficients for the quartic orthogonal polynomial $\beta_{jk}$, listed in Table 1. Then, we derived the quartic orthogonal ephemeris for the phase function $\vartheta(\vartheta_{0})$, formally identical with the same relations for $\sigma$ Ori E (Eq. (29)). The origin and the period of the linear approximation are: $\tilde{M}_{0}=2\,453\,943.3871(12)$, $\tilde{P}=2\fd 204\,284\,86(7)$. The mean lengthening of the observed period of $\dot{P}=5.4(4)\times 10^{-9}=0.170(13)$ s yr*-1* is well established, opposing to higher terms that can be neglected.

Currently, BS Cir belongs to the best and longest monitored mCP stars. However, we are not able to mine more information from the period analysis, because there are no signs of any elongation from the simple model assuming a steady rise of the period. Nevertheless, some information we can be gain from possible secular changes in the light curves themselves. Mikulášek et al. (2015a) have found small long-term trends in the shapes of the light curve (the amplitude of both spots is decreasing) while the phase distance between the centers of the spots remains constant. Such moderate changes (if real) could bear evidence of a slow free precession of the stellar body.

5 Comparison of mCP stars having unstable periods.

Discussion

The relative number of stars showing rotational period changes on the time scale of decades among single upper main sequence stars was theoretically quite unpredicted, and until now the true cause (or causes) of such peculiar behavior is not known. The main aim of the present analysis of more than 60 thousand mostly good photometric measurements of four magnetic CP stars, suspected of unstable rotation, was to prove those suspicions and find the properties of this phenomenon which could lead us to the solution of that puzzle.

The analysis confirmed without any doubts that all four stars mentioned above, namely CU Vir, V901 Ori, BS Cir, and $\sigma$ Ori E, display changes in their observed period of variability with a certainty of more than 13 $\sigma$. As we have shown in the case of V901 Ori (for details, see Sect. 4.3), the phase curves and their relative phase shifts remain constant during several decades, which indicates that the outer layers of mCP stars are extraordinarily stable including the magnetic geometry and the locations of photometric and spectroscopic features. This stability is very likely caused by the strong global magnetic field frozen into the atmospheric plasma that fastens the surface of the star, so it behaves like a solid body. Thus, the monitoring of the period changes tells us something about rotation of this outer, directly observable part of mCP stars. The found gradual variations in the angular velocity of several mCP stars could, but need not, inform us about the rotation of the inner part of the star. Nevertheless, the finding of the physical mechanism causing the observed rotational variations represents an interesting astrophysical challenge both for observers and theoreticians.

The studied subgroup of mCP stars consists of very disparate members with a minimum of common properties. CU Vir is a unique main-sequence radiopulsar, V901 Ori is one of the hottest mCP stars exhibiting a strong and extraordinarily complex magnetic field, $\sigma$ Ori E is a hybrid between a Be star and a He-strong star with the very strong magnetized stellar wind, and the last star is an old, moderate cool mCP star. It seems that the only shared common property is the fact that their periods are shorter than those of the majority of magnetic chemically peculiar stars (Mikulášek et al., 2009). Unfortunately, this might only be due to the selection effect following from relations in Eq. (31), according to which uncertainties in the determination of the period and period changes are proportional to the square of the period ($\sim P^{2}$).

The relevance of the individual parameters of the expansion can be judged by a comparison of their numerical values and the estimates of their uncertainties, $\delta P_{0},\ \delta\dot{P}_{0},\ \delta{\ddot{P}_{0}}$

[TABLE]

where $q$ is an auxiliary quantity, $\sigma$ is the mean scatter of data used for the ana-lysis, $n$ is the number of measurements, $A$ is the mean amplitude of variations, and $\Omega$ is the time interval covered by observations666Applying formula (31) to V901 Ori, one of the best photometrically monitored variables among mCP stars with parameters: $P=1.539$ d, $\sigma=0.0048$ mag, $A=0.040$ mag, $\Omega=40$ yr, $n=4500$, the uncertainties of the period $P$ and its first and second derivatives $\dot{P}$ and $\ddot{P}$, can be estimated as $\delta P=4.5\times 10^{-7}$ d $=0.04$ s; $\delta\dot{P}=2.4\times 10^{-10}=0.76$ s cen*-1*, $\delta\ddot{P}\simeq 2\times 10^{-13}$ d*-1*.. The relations show that the best chance of revealing period changes exists in continuously monitored stars with well defined phase curves and periods as short as possible.

Unfortunately, we do not have reliable statistics on the percentage of rotationally unstable stars among various types of mCP stars. This complicates our further speculation on the nature of the observed period variations. Nevertheless, it would be inspiring to discuss possible explanations and expectations more carefully.

In the following, we shall discuss the particular results of the period analyses, many of which are presented in Table 1. The main source of the phase information is the photometry which is most plentiful and relatively accurate.

The best-observed star, judging according to the ratio of the total number of observations $n$ to the square of the mean uncertainty $\sigma$, $n/\sigma^{2}$, is apparently $\sigma$ Ori E, which is mostly due to the unprecedentedly good and numerous observations from the MOST satellite (Townsend et al., 2013). Very well observed is also CU Vir thanks to long-time photometric monitoring of the star by Diana Pyper and Saul Adelman (Pyper et al., 2013, and many others).

The only star showing verifiably sinusoidal long-term changes in its period is CU Vir which has been monitored during the whole cycle of about 66 years. For expressing a satisfactory polynomial model of the phase function (see Sect. 4.1), a quintic parabola is needed, in the case of CU Vir, a cubic parabola for V901 Ori and a simple square parabola for the other two stars. Higher terms seem to be nearly negligible. Unfortunately, we cannot dig up more information on possible cyclicity in all the stars, excluding CU Vir. It is of particular interest that the mean rate of the period is rising in all the stars, although the rotation of both CU Vir and V901 Ori is currently accelerating.

Some physical information, excluding precession as the explanation of period changes in CU Vir, V901 Ori, and $\sigma$ Ori E is the amplitude of phase change $A_{\varphi}$. In the case of CU Vir we can read the amplitude $A_{\varphi}$ directly from its O-C diagram – 0.62 of the period (see Fig. 5); for the other stars we can estimate at least its lower limit according to a simple relation: $A_{\varphi}>\tilde{P}^{\prime}/2n_{\rm{rev}}^{2}$, where $n_{\rm{rev}}$ is the number of revolutions monitored. The limit for V901 Ori is 0.55 of the revolution, 0.26 in $\sigma$ Ori E and only 0.11 for BS Cir. This value is allowed by the model of free precession (see appendix in Mikulášek et al., 2008). However, the purely quadratic course of the modulation of the phase function indicates that the true phase amplitude is rather larger.

The nature of the period variations of the discussed stars is not known up to now, and we cannot exclude that their reasons may be different. The present analy-sis supports the idea (more broadly investigated in Mikulášek et al., 2017) that the tiny, nearly linear period changes of BS Cir can be due to precession of the magnetically distorted star, while the nearly sinusoidal period oscillations of at least CU Vir might be interpreted as a consequence of internal waves disseminating within the magnetic rotating stars. Period changes in hot mCP stars with stellar wind escaping from their extended magnetosphere can be explained by angular momentum loss ($\sigma$ Ori E), modulated by the gradual reconfiguration of the magnetic field firmly connected with the surface layers (V901 Ori).

6 Conclusions

We have presented a novel and commonly applicable method of period analysis based on the simultaneous exploitation of all available observational data containing phase information. This phenomenological method can monitor gradual changes in the observed instantaneous period very efficiently and reliably.

We have presented up to date results of the period monitoring of V901 Ori, CU Vir, $\sigma$ Ori E, and BS Cir, known to be mCP stars changing their observed periods and have discussed the physics of this unusual behavior. To compare the period behavior of those stars, we treated their data with an orthogonal polynomial model, which we specifically develop for this purpose.

We have confirmed period variations in all investigated stars and shown that they reflect real changes in the angular velocity of outer layers of the stars fastened by their global magnetic fields.

However, the nature of the observed rotational instabilities has remained elusive up to now. The discussed group of mCP stars is inhomogeneous to such extent that each of the stars may experience a different cause for its period variations.

Acknowledgements.

The research was supported by the project GA ČR 16-01116S. The author thanks to Stefan Hümmerich and Miroslav Jagelka for their careful reading of the manuscript and improving its language. The author of the review paper is also indebted to his close collaborators and all who provided him with data on investigated mCP stars, particularly to S. A. Adelman, K. Bernhardt, M. Chrastina, S. de Villiers, M. Drozd, G. W. Henry, S. Hubrig, S. Hümmerich, M. Jagelka, J. Janík, O. Kochukhov, J. Krtička, D. O. Kudryavtsev, P. Kurfürst, R. Kuschnig, J. Landstreet, J. Liška, T.  Lüftinger, M. Netopil, M. Oksala, E. Paunzen, T. Pribulla, M. Prvák, D. Pyper, T. Rivinius, A.  Reiners, I. I.  Romanyuk, T. Ryabchikova, E. Semenko, M. Shultz, D. Shulyak, N. Sokolov, K. Stepień, G. Szász, R.H.D. Townsend, C.  Trigilio, M. Vaňko, G. A. Wade, M. Zejda, P. Zielinski, J. Zverko, P. Zvěřina, J. Žižňovský, and many others.