#  The investigation on Magnetic Activity of HD 134319 based on TESS
Photometry and Ground-based Spectroscopy

Fukun Xu${}^{\href https://orcid.org/0000-0002-8618-3551}$,1,2 Shenghong
Gu1,2,3111shenghonggu@ynao.ac.cn(SG) and Panogiotis Ioannidis4
1Yunnan Observatories, Chinese Academy of Sciences, Kunming 650216, China
2Key Laboratory for the Structure and Evolution of Celestial Objects, Chinese
Academy of Sciences, Kunming 650216, China
3School of Astronomy and Space Science, University of Chinese Academy of
Sciences, Beijing 101408, China
4Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029
Hamburg, Germany
Email: xufukun@ynao.ac.cn(FX); shenghonggu@ynao.ac.cn(SG)

(Accepted xxx. Received xxx; in original form xxx)

###### Abstract

We present analysis of the starspot properties and chromospheric activity on
HD 134319 using high precision photometry by TESS in sectors 14–16 (T1), 21–23
(T2) and high-resolution spectroscopy by OHP/ELODIE and Keck/HIRES during
years 1995–2013. We applied a two-spot model with GLS determined period
$P=4.436391\pm 0.00137$ day to model chunks sliding over TESS light curve, and
measured relative equivalent widths of Ca II H and K, H$\beta$ and H$\alpha$
emissions by subtracting overall spectrum from individual spectrum. It was
found that a two-spot configuration, i.e. a primary, slowly evolving and long-
lasting spot (P) plus a secondary and rapidly evolving spot S, was capable of
explaining the data, although the actual starspot distribution is unable to
derived from collected data. Despite the spot radius-latitude degeneracy
revealed in the best-fit solutions, a sudden alternation between P and S radii
followed by gradual decrease of S in T1 and the decrease of both P and S from
T1 to T2 were significant, corresponding to the evolution of magnetic
activity. Besides, S revealed rotation and oscillatory longitude migration
synchronized to P in T1, but held much larger migration than P in T2. This
might indicate the evolution of internal magnetic configuration. Chromospheric
activity indicators were found to tightly correlated with each other and
revealed rotational modulation as well as a long-term decrease of emissions,
implying the existence and evolution of magnetic acitivity on HD 134319.

###### keywords:

starspots – stars: chromospheres – stars: magnetic fields – stars: rotation

††pubyear: 2022††pagerange:  The investigation on Magnetic Activity of HD
134319 based on TESS Photometry and Ground-based Spectroscopy– The
investigation on Magnetic Activity of HD 134319 based on TESS Photometry and
Ground-based Spectroscopy

##  1 Introduction

Stellar magnetic activity arising from the enhancement of magnetic flux
generated at the tachocline, penetrated the outer convection zone and then
exceeded the surface of late-type stars manifests itself in abundant
phenomena, such as starspots, plages and flares, causing series of distortions
detectable in both photometric and spectroscopic observations. Periodic or
quasi-periodic variations in photometric measurements of cool stars are widely
attributed to the rotational modulation of starspots. This makes starspot a
good tracer of stellar rotation and detector of differential rotation which
plays a critical role in generating and maintaining the solar-like magnetic
activity through an $\alpha\Omega$ dynamo (Işık et al., 2011). Meanwhile,
emissions of chromospheric activity indicators such as the Ca II H & K
resonance and Balmer lines are generally used as diagnoses of magnetic
activity in stellar outer layers and also show modulation under stellar
rotation (Vida et al., 2015).

In analogy to the case of the Sun, the magnetic activity on active stars is
believed to reveal both spatially inhomogeneous distribution and temporally
evolution. The former creates inhomogeneous configuration of the magnetic
field accompanying by the emergence of photospheric starspot and chromospheric
excitation (Hempelmann et al., 2016; Balona et al., 2019), while the latter
exhibits as evolution of starspot in terms of emerging, migrating, decreasing
and vanishing during its life as well as the variation of chromospheric
proxies in multiple timescales (Mittag et al., 2019). By long-term monitoring
on photometric and spectroscopic variations, one can derive the starspot
distribution, evolution and its connection with the chromospheric activity,
which were reported on several targets and believed to be widely existed on
late-type stars (e.g. Vida et al., 2015; Morris et al., 2018; Xu et al.,
2021). However such a study is usually non-trivial due to the limit in
acquiring high precision data of target in distance.

In recent years, space-borne telescopes such as MOST (Walker et al., 2003),
CoRoT (Auvergne et al., 2009), Kepler (Koch et al., 2010) and most recent
Transiting Exoplanet Survey Satellite (TESS, Ricker et al., 2015) produced
bulk of long-term photometric data. Since its initiation, TESS was designed to
cover $>80\%$ of the whole sky and deliver light curves (LCs) spanning from
$27$ days to years for millions of stars, providing us the opportunity to
analyse the starspot configuration and its evolution on stellar photosphere
with high-precision (Reinhold et al., 2013; Namekata et al., 2020).

The star HD 134319 (TIC 202426247, GJ 577, Gl 577, HIP 73869, with magnitude
$V=8.41$ mag, $G=8.23$, and color $B-V=+0.677$), is a young ($\sim
0.625\,\mathrm{Gyr}$, Montes et al., 2001; McCarthy et al., 2001; Decin et
al., 2003) active G5 main-sequence star of BY Dra type (Soderblom, 1985) which
can be taken as a young solar analogue with a shorter period of about
$4.448\,\mathrm{day}$ (Messina et al., 1998). Locating in the north semi-
sphere, was observed in six sectors by TESS (14–16 and 21–23) spanning over
$300$ days, revealing significant active regions surviving over 90 days as
well as evidences of evolution in both short and long timescales. HD 134319’s
stellar parameters were greatly summarized in
SIMBAD111http://cdsweb.u-strasbg.fr/ and
ExoFOP222https://exofop.ipac.caltech.edu/tess/target.php?id=202426247 and some
of them are listed in table 1. The age of HD 134319 was firstly deduced as
$0.625\,\mathrm{Gyr}$ by classifying it as a member of the Hyades supercluster
in terms of space motion by Montes et al. (2001) and McCarthy et al. (2001),
while other values from high-resolution echelle spectra were determined by
Valenti & Fischer (2005) ($6.3\,\mathrm{Gyr}$), Takeda et al. (2007)
($2.36\,\mathrm{Gyr}$) and lastly Isaacson & Fischer (2010)
($0.06\,\mathrm{Gyr}$). Besides, a young nearby companion binary system, with
GAIA magnitude of $G=14.25$, and separation of $5.39\pm 0.02\,\arcsec$ from HD
134319, was detected by infrared imaging (McCarthy et al., 2001; Mugrauer et
al., 2004) and adaptive optics (Lowrance et al., 2003) (see section 2.2). The
radial velocity (RV) of HD 314319 was firstly determined as
$V_{\mathrm{R}}=-3.8\,\mathrm{km/s}$ by Wilson & Joy (1950), and most recently
derived from OHP/ELODIE data with value $V_{\mathrm{R}}=-6.362\pm
0.011\,\mathrm{km/s}$, which was included in catalogue of RV standard stars by
GAIA (Soubiran et al., 2018). Its long-term stability from 2002 to 2013 was
reported by Butler et al. (2017) in their RV exoplanet survey using precision
RV measurements by Keck/HIRES.

The chromospheric activity of HD 134319 was recognized by several studies
(Soderblom, 1985; Duncan et al., 1991; Wright et al., 2004; López-Santiago et
al., 2010; Butler et al., 2017; Morris et al., 2019), by measuring its Ca II H
& K emissions in terms of S-index and/or its corresponding derivative $\lg
R_{\mathrm{HK}}^{\prime}$, giving typical median values of $S\sim 0.42$ and/or
$\lg R_{\mathrm{HK}}^{\prime}\sim-4.12$, which indicates HD 134319 to be
likely an ultra-active star (Boro Saikia et al., 2018). While the spot
coverage on HD 134319 was reported to be very small through estimation of TiO
molecular band absorption by Morris et al. (2019) in their investigation of
the connection between activity level and spot coverage, they noted HD 134319
as an outlier with high activity but small spot coverage.

The starspots on HD 134319 from photometric analysis were firstly proposed by
Messina & Guinan (1998) who clarified the increase of peak-to-peak amplitude
of LC towards decreasing wavelengths to the presence of of spots. Soon Messina
et al. (1998) used rotational modulation to predict long-lasting active
longitudes and spot covering fraction of at least $fs\geq 0.16$ for this BY
Dra type variable by employing light curve inversion to multi-band (u, v, b
and y, Strömgren, 1966) photometric data spanning from 1991 to 1995. Later
Messina et al. (2001) inverstigated the rotation-activity connection from
photometry, including HD 134319.

In a word, the youth and high level of chromospheric and photospheric magnetic
activity on HD 134319 make it a good proxy for the young Sun not far after its
arrival at the zero age main sequence. In this paper we present a detailed
analysis of starspot configuration, its evolution and chromospheric activity
on HD 134319, using high-precision LC newly observed by TESS and collected
spectroscopic data. In section 2, we introduce both the photometric and
spectroscopic observations. We describe our approach for sliding spot
modelling and the method for measuring relative equivalent widths from
spectroscopic data in section 3. We then give the results and respective
discussions in section 4, and finally the conclusions are summarized in
section 5.

##  2 Data and REDUCTION

Table 1: Stellar parameters of HD 134319 in literatures.

Parameter | Value
---|---
RA, DEC | 15h05m49.90423, +64°02’49.9415"${}^{\ref{ft:SIMBAD}}$
$V$ mag | V=8.41${}^{\ref{ft:SIMBAD}}$
Radius ($R_{\sun}$) | 0.94a, 0.92984998b, 0.92985${}^{\ref{ft:exofop}}$
Mass ($M_{\sun}$) | 1.01${}^{\ref{ft:exofop}}$
$T_{\text{eff}}$ (K) | 5668.04${}^{\ref{ft:exofop}}$ , 5635.87a, 5662c, 5636d
Metallicity [Fe/H] | 0.03${}^{\ref{ft:SIMBAD},\ref{ft:exofop}}$
log$g$ | 4.50556${}^{\ref{ft:exofop}}$
$v~{}\text{sin}i$ (km/s) | 10.6c, 17.89e, 11.39$\pm$0.06f, 10.9d

a(Gaia Collaboration, 2018), b(Stassun et al., 2018), c(Valenti & Fischer,
2005), d(Morris et al., 2019), e(White et al., 2007), e(López-Santiago et al.,
2010).

###  2.1 Spectroscopy

The spectroscopic data for HD 134319 collected online consist of two parts.
One is the observations publicly accessible and available for download from
Keck Observatory Archive333https://koa.ipac.caltech.edu/, which were carried
out during years 1999 - 2013 by the HIRES spectrometer on the Keck I-10m
telescope (Keck/HIRES, Vogt et al., 1994). The spectral resolution is about
$R=67000$. The other is the data from ELODIE Archive444http://atlas.obs-
hp.fr/elodie/ which was observed by ELODIE spectrograph on the 1.93m telescope
of Observatoire de Haute-Provence (OHP/ELODIE, Baranne et al., 1996; Moultaka
et al., 2004). The spectral resolution is $42000$. Totally 26 spectra of
OHP/ELODIE observed during years 1995 - 1997 and 22 spectra of Keck/HIRES
observed during years 1999 - 2013 were collected.

Reduction of the raw data from Keck/HIRES was carried out using standard tasks
in the IRAF package555IRAF is distributed by the National Optical Astronomy
Observatories, which is operated by the Association of Universities for
Research in Astronomy (AURA), Inc, under cooperative agreement with the
National Science Foundation., which basically included flat-field division,
scattered light correction and wavelength calibration. The bias and dark
corrections were omitted due to their absence in some spectra, and that flat
corrections were done using flats with wide orders (there existed another type
of flats with narrow order in a fraction of observations).

OHP/ELODIE provides two types of products: the combined one dimensional
spectrum (s1D) which was uniquely resampled in wavelength and given in
instrumental relative flux covering $4000$ to $6800\,\mathring{A}$, and the
two dimensional image (s2D) which contains the extracted and deblazed spectrum
in 67 orders.

According to wavelength coverage, we chose the spectral portions around Ca II
H & K, H$\beta$ and H$\alpha$ lines. For Keck/HIERS, H$\alpha$ is available
since 18 August 2004. For OHP/ELODIE s2D data, H$\beta$ portion was selected
from order 32 (H$\beta$ in order 31 was omitted due to bad pixels H$\alpha$
portion was selected from order 64, and Ca II H & K (in orders 3 and 2,
respectively) were omitted due to low signal to noise ratio (SNR).

###  2.2 GJ 577 B/C

A young nearby proper motion companion of HD 134319, with separation of
$5.39\pm 0.02\,\arcsec$ to the west and GAIA magnitude of $G\sim 14.25$, was
detected by infrared imaging (McCarthy et al., 2001; Mugrauer et al., 2004)
and adaptive optics (Lowrance et al., 2003). Mugrauer et al. (2004) determined
its spectral type as $M4.5$, mass as $0.16-0.20\,M_{\sun}$, age as $>100$ Myr
and found its strong H$\alpha$ emission and deep TiO and VO molecular
absorption bands. Lowrance et al. (2003) further recognized it as a binary
system (GJ 577 B/C) with separation $0.082\,\arcsec$, both of which lie on the
stellar/substellar boundary with spectral type between $M5$V and $M6$V. This
scenario of binary system was also noted by Burgasser et al. (2005) and Martin
et al. (2017) and confirmed by proper motion anomaly (Kervella et al., 2019).

The separation of $5.39\,\arcsec$ is apparently far enough in spectroscopy on
Keck/HIRES666https://www2.keck.hawaii.edu/inst/hires/slitres.html and
OHP/ELODIE777http://www.obs-hp.fr/www/guide/elodie/elodie.html and thus
unlikely to contaminate the primary’s spectra. However the companion is near
enough in TESS photometry (spatial resolution of
$21\,\arcsec/\mathrm{px}$888https://heasarc.gsfc.nasa.gov/docs/tess/observing-
technical.html) to be crowded in the same pixel as the primary. For simplicity
we assume that it will not put detectable distortion, but at most a small
offset about $0.40\%$, on resultant TESS LC because it is much fainter
($G=14.25$) than the primary ($G=8.23$) as measured by GAIA (Gaia
Collaboration et al., 2018).

###  2.3 TESS Light Curve

The TESS mission (Ricker et al., 2015) was designed to observe an extremely
wide sky area in sectors measuring $24\times 96\degr$, utilizing four wide-
field cameras aligned in a mosaic, extending from near the ecliptic equator to
beyond the ecliptic pole. The initial two-year observation strategy of TESS
includes 26 sectors, each of which is observed in two highly elliptical
$13.7\,\mathrm{day}$ orbits around the Earth, spanning about $27$ days. TESS
provides simple aperture photometry (SAP) and pre-search data conditioning
(PDC) processed LC generated from two-minute cadence data, processed by NASA’s
TESS Science Processing Operations Center (SPOC) pipeline (Jenkins et al.,
2016) which is an successor of Kepler’s
pipeline999https://heasarc.gsfc.nasa.gov/docs/tess/, as well as the calibrated
target pixel file (TPF, which contains the pixel level photometric time
series). All of them are publicly accessible and available for download from
Mikulski Archive for Space Telescopes (MAST) at STSci data
base101010https://mast.stsci.edu/ to ensure simultaneous usage of final
products and check of the data quality in pixel level.

Considering that the PDC process was primarily designed for planet hunting but
not for the study of stellar variabilities which usually show non-unique and
much broader periodic durations and variations than planetary transits, and
could be ignored or misinterpreted by the pipelines, one is generally
recommanded to check the systematic correction by PDC and redo it when
necessary (Still & Barclay, 2012; Vinícius et al., 2017). We used calibrated
TPFs to check whether the pipeline chosen apertures, which are different
between sectors, can overcome the possible contamination from the crowed
neighbours and temporal variation of photons distribution, i.e., whether the
PDC processed LC is suitable for stellar activity investigation.

The determination of optimal aperture in wide field-of-view photometry like
TESS (Ricker et al., 2015) is non-trivial due to the crowdedness of objects
versus wide spread image of stars on CCD. Crowdedness contributes to extra
flux on target and should be avoided by employing smaller aperture, which,
however, causes wastage of the target flux and leads to lower SNR. Estimations
of flux fraction and crowding matrix were done by PDC to recover the real
brightness of target in Kepler and later TESS (Batalha et al., 2010). However
the automatically predetermined aperture and its correction in pipeline might
be unsuitable for individual case, as widely revealed (for example, the large
jumps between adjacent sectors) on many targets (e.g. Özavcı et al., 2018; Xu
et al., 2021), and should be checked with caution, especially in investigation
of stellar variation.

HD 134319 was observed by TESS in sectors 14–16 and 21–23 during its second
year’s routine. Besides GJ 577 B/C, there are other three neighbours
(1620024559829845504, 1620023597757174144 and 1620024692973605760) matched
with Gaia DR2 (Gaia Collaboration et al., 2016, 2018) around HD 134319
recorded on the CCD image of TPFs. The last two are far enough from HD 134319
and properly excluded by PDC. The first one ($\text{Gmag}=13.9574$, fainter
than HD 134319 by about two orders of magnitude) was excluded outside the PDC
apertures in sectors 14–16, but not in sectors 21–23 due to the variation of
flux spread. Considering that no recognizable difference in profile of whole
LC between with and without this neighbour was found, we concluded that it
dose not contribute detectable contamination on HD 134319’s LC.

Moreover, we assessed the validity of PDC determined apertures by comparing
the LC profiles derived from different apertures. One can manually define a
smaller TPF masks, i.e. photometric apertures, by setting a larger value of
parameter threshold, which defines the threshold in choosing the pixels to be
integrated in Lightcurve, and vice versa. We found that the profile of
integrated LC was stable with masks both slightly smaller and larger than the
PDC, indicating that the spread function of HD 134319 was stable over time and
the apertures determined by PDC were suitable for our study.

Consequently, we used the PDC corrected LC, i.e. PDCSAP long cadence data, in
our analysis. There was no empirical evidence of improper correction by PDC
for all sectors except two parts with improper offsets. One is the first half
of sector 22 (1899 – 1927 after $\mathrm{BTJD}-2457000$) and was subtracted by
an offset of $400$ ADU, the other is sector 23 and was added by an offset of
$1000$ ADU. This kind of offset, to remove the discontinuities between the
improper segment and its surrounding LCs, was estimated qualitatively due to
the too short time length of segment to adopt a semi-quantitative method (e.g.
Özavcı et al., 2018). The possible under- or over-estimation of the offset
attributes to a common deviation of LC, which results in either spot size
variation or additional polar spot but has not much to do with the spot
evolution. At last the LC was normalized to unity by dividing its median over
the whole time span (figure 1).

Figure 1: Resultant light curve.

##  3 Model and analysis

###  3.1 Modelling the light curve

Periodic or quasi-periodic variation in photometric LC is generally attributed
to the starspot rotating into and out of view from Earth (Valio et al., 2017).
As a manifestation of internal magnetic flux on the stellar surface, starspot
generates notable distortion on LC which makes it as good indicator of the
internal magnetic structure and tracer in measuring the surface rotation.
However, reconstructing two-dimensional distribution of the starspots from
one-dimensional disk-integrated photometric time series is always not easy due
to the weak constraints from the observation itself (Lanza et al., 2016) as
well as the large parameter space and high degeneracy between parameters of
spot model. One should to reduce the number of parameters as small as possible
in spot modelling and treat the results with caution.

####  3.1.1 Rotational period of HD 134319

In literatures (table 2), the rotation period of HD 134319 was reported as
$4.448\,\mathrm{day}$ by Messina & Guinan (1998) and $4.43\,\mathrm{day}$ by
Wright et al. (2011) from photometric observations, and estimated as
$5.0\,\mathrm{day}$ by Wright et al. (2004) and $3\,\mathrm{day}$ by Isaacson
& Fischer (2010) from statistical determination of chromospheric activity
fluctuations. With the favour of long-term precise photometry by TESS, we can
obtain an accurate measurement of the period revealed by the rotational
modulation. We used the generalized Lomb-Scargle periodogram (GLS, Zechmeister
& Kürster, 2009) to measure the period of time series in sectors 14–16 and
21–22 of TESS (sector 23 was temporarily excluded due to its empirically
questionable LC) as $P=4.436391\pm 0.00137\,\mathrm{day}$, which is slightly
larger than the one of Wright et al. (2011) and smaller than the one of
Messina et al. (1998).

Table 2: Rotational periods of HD 134319 measured by different authors.

Reference | Period (days) | Memo
---|---|---
Messina et al. (1998)a | $4.448\pm 0.005$ | b
Wright et al. (2004) | $5.0$ | cf. $R_{\mathrm{HK}}^{\prime}$
Isaacson & Fischer (2010) | $3$ | cf. $R_{\mathrm{HK}}^{\prime}$
Wright et al. (2011) | $4.43$ | c
This work | $4.436391\pm 0.00137$ | d

aAlso in Messina & Guinan (1998); Messina et al. (2001), bphotometric Scargle-
Press period (Scargle, 1982), cphotometric period by FEPS (Meyer et al.,
2006), this value was cited by Butler et al. (2017); Mittag, M. et al. (2018);
Morris et al. (2019), dphotometric period by GLS (Zechmeister & Kürster,
2009).

In figure 2 we show the phase-folded LC of HD 134319 using this rotation
period. Two dips in brightness per rotation with different depths are visibile
throughout the observation age. The primary dip ("P") is visible over the
whole observation season, starting from $\text{Phase}\sim 0.7$ at onset and
shifting gradually to $\text{Phase}\sim 0.5$ in the end. The secondary dip
("S"), emerging at $\text{Phase}\sim 0.3$ at onset, is recognizable before
time $1763.5$, then becomes ambiguous after that due to its small amplitude.
On the other hand, notable decrease in brightness loss of both "P" and "S" can
be found with increasing time, revealing notable evolution of spot
configuration.

Figure 2: Phase-folded, median smoothed LCs for HD 134319 with GLS estimated
rotation period $P=4.436391\,\mathrm{day}$, showing a larger, primary "P" plus
a smaller, secondary "S" condensed regions under starspot modulations and
their notable evolution over time. The vertical position for each curve
corresponds to the start time of chunk. Each time window is scaled to larger
amplitude for easier view by subtracting $1$ then multiplying it by $800$,
i.e. $f^{\prime}=800(f-1)$. The long-time gap from $1763.5$ to $1869.5$ is
marked by breaks in vertical axes.

####  3.1.2 GEMC_LCM

As an analytical model towards photometric observation, light curve modelling
(LCM, Budding, 1977; Dorren, 1987) simulates the LC with assumption that the
drop in light intensity is caused by a small number of circular spots
occupying on the stellar photosphere and thus is capable of reducing the
parameter space to some extent. LCM consists of parameters relative to the
star (the unspotted intensity $U$, the linear limb-darkening effect $u$ and
the inclination $i$) and individual circular spot (the spot-to-photosphere
intensity ratio $\kappa_{w}$, the spot latitude $\beta$, longitude $\lambda$,
radius $\gamma$ and its rotational modulation period $P$). Practically stellar
parameters $i$, $U$ and $u$ are constants and parameter $\kappa_{w}$ can be
approximately constants and common to all spots. Taking into account the
rotation, the periodical variation of longitude can be better defined as a
function of of time $t$ by the rotational period $P$ and initial longitude
$\lambda(0)$ as $\lambda(t)=2\pi t/P-\lambda(0)$ for convenience (Xu et al.,
2021).

An efficient program, called "GEMC_LCM" (Xu et al., 2021), which was designed
to inherit both the superior global optimization power of genetic algorithm
and the high efficiency on parameter space exploration of Markov Chain Monte
Carlo algorithm, was employed for spot modelling using LCM of Budding (1977).
GEMC_LCM is able to run parallelly on mulit-core CPU with the favour of
OpenMP.

####  3.1.3 Parameter degeneracies

Besides the large parameter space, the parameter degeneracies in spot
modelling are another kind of problem in photometry modeling, which block the
algorithm converging to the true solution. Practically it is difficult to
overcome this problem in spot modelling in the absence of independent
constraints. However it is reasonable to divide such degeneracies into two
types according to their sources, i.e., they are inherent or produced. The
former comes from the model itself, while the latter comes from measurement
related effects which can be suppressed by high precision photometry.

An inherent degeneracy between stellar inclination and spot latitudes was
identified from numerical simulation by Walkowicz et al. (2013), who also
showed the possibility of failure in detecting the differential rotation when
LC is dominated by one large active region. Actually, when it comes to the
case of LCM with single spot, one can easily find that this degeneracy can be
expressed analytically as

$I_{c}(i,\beta)=I_{c}\left(90\degr-|\beta|,\mathrm{Sign}\left[\beta\right](90\degr-i)\right)$
(1)

Where $I_{c}$ is the simulated LC, $\mathrm{Sign}\left[\beta\right]$ means the
sign of $\beta$, and specially $\mathrm{Sign}\left[0\right]=\pm 1$. In other
words, the spot at latitude $\beta$ with inclination $i$ creates definitely
the same LC as spot at latitude $\mathrm{Sign}\left[\beta\right](90\degr-i)$
with inclination $90\degr-|\beta|$.

A produced degeneracy between spot latitude and radius due to the low
constraint of photometry versus noise and contamination was noted as an "ill-
posed" issue (Lanza et al., 2016). Fox example, a larger spot at higher
latitude distorts the LC with similar amplitude as a smaller spot at lower
latitude by means of projection to light of sight. However, it is possible to
minimize such a defect to some extent with the mercy of high precision
photometry as they lead to subtle but different profiles in principle
(Walkowicz et al., 2013). Besides, the contamination coming from the residual
systematic errors in long-term observations such as Kepler and TESS can mixed
with stellar variations in multiple timescales. Thus it is usually recommanded
to do the systematic error corrections, circumstances alter cases, in stellar
variation investigations (Still & Barclay, 2012; Vinícius et al., 2017). We
have checked the LC as described in section 2.3.

To derive the properties of spot configuration independent on parameter
degeneracies, it is feasible to employ some assumptions a priori. A predefined
stellar inclination is capable of breaking its degeneracy with spot latitude.
The stellar inclination can be calculated from measured $v\sin i$, stellar
radius $R$ and rotation period $P$ by

$i=\arcsin\left(v\sin i*P/{2\pi R}\right)$ (2)

With the measured period $P=4.436391\,\mathrm{day}$ as derived in section
3.1.1, the stellar radius in table 1 implies a maximal $v\sin i\sim
10.7\,\mathrm{km/s}$. Considering that the projected rotational velocity
$v\sin i$ was estimated in the range between $10.6\,\mathrm{km/s}$ and
$17.89\,\mathrm{km/s}$ (table 1), HD 134319 is preferred to have a high
inclination. As a comparison, we adopted fixed inclination of $i=75\degr$ and
$i=90\degr$ in modelling.

####  3.1.4 Surface differential rotation

The differential rotation was thought to play a critical role in generating
and maintaining the stellar magnetic field. The surface differential rotation
(SDR) on the Sun was observed in relative motion of sunspots and can be
expressed by a quadratic law

$P(\beta)=P_{\mathrm{eq}}/(1-\alpha\sin^{2}\beta)$ (3)

Where $P(\beta)$ is the stellar rotation period at latitude $\beta$,
$P_{\mathrm{eq}}$ is the reference period at equator, and $\alpha$ is relative
rate representing the strength of SDR, which was measured as $\alpha=0.2$ on
the Sun. A similar law of SDR is generally assumed on stars other the Sun by
analogy (Henry et al., 1995). By this definition a zero $\alpha$ means a rigid
rotation, a positive $\alpha$ represents a Solar-like SDR while a negative
$\alpha$ represents an anti-Solar SDR.

####  3.1.5 Spot modelling

The whole LC was split into chunks with $5$ days duration each, slightly
longer than one rotation period, and advanced by $1$ day forward, providing a
total of $129$ such chunks for spot modelling. In such a short timescale, the
spot positions are not expected to evolve, while their radii are allowed to
change linearly indicated by the apparent variation of LC within adjoint
rotations.

The modelling was done under hypothesis of a stellar sphere with uniform
surface brightness occupied by circular dark spots with uniform temperature.
The pre-determined parameters were chosen as follows. The limb darkening
effect was assumed to be a linear law with coefficient $u=0.5119$ from the
table of Claret (2018) for a star with $T_{\mathrm{eff}}\sim 5600\,K$ and
$\log{\it{g}}\sim 4.5$. The uniform spot-to-photosphere intensity ratio was
fixed to $\kappa_{w}=0.22$ empirically derived from statistics (figure 7 of
Berdyugina (2005)). The stellar surface brightness was also fixed as the
maximal value of LC, $U\sim 1.02$, which is artificial to some extent because
the target is always occupied by spots. At last, by comparison, we employed
two values of stellar inclination $i=75\degr$ and $i=90\degr$.

The spot configuration in each chunk was derived under a two-spot model. The
number of spots was determined by two reasons. One is that there are two
minima in almost all chunks, probably inferring a distribution of two
concentrated active regions on the surface. The other comes from our tests in
checking the global optimizing capability of models with two and three spots.
It was found that three (or more) spot model provided slightly better fiting
as expected, while resulted more likely in families of solutions with equally
good fits which reduces its reliability, comparing with two-spot model.
Therefore a two-spot model was sufficient to obtain acceptable fits. Note
that, strictly speaking, the real number of spots (or active regions) might or
might not eqaul to the number of minima in LC (Jeffers & Keller, 2009; Basri &
Shah, 2020), but we would rather to fit the LC with a simple spot-model for
efficiency and convergence of optimization, and if the spot distributions on
the star can be described by such a two-spot model it would be the result.

The determination of rotational modulation period corresponding to spot is
non-trivial. Usually it can be determined with high precision due to the
strong constraint of long-term photometric time series on longitude, despite
exceptions, e.g. when one spot dominates the LC (Walkowicz et al., 2013).
However, such a constraint becomes weak in our case when we do spot modelling
for single chunk, i.e. in pretty short timescale. Instead we adopted fixed
period in modelling after several numerical experiments, and thus the derived
longitude might vary due to two effects, i.e. the evolution of spot longitude,
and the difference between the chosen and actual rotational modulation period
of spot.

As a conclusion, we employed a two-spot model to inverse the longitudes,
latitudes, linearly evolving radii of spots (i.e. eight free parameters in
total) for each chunk, using fixed period $P=4.436391\,\mathrm{day}$. For the
sake of higher probability of globally optimized solution, each chunk was
fitted for many runs, and the results were then divided into groups
corresponding to different spot configurations. Finally we elected the group
with minimal combined residuals $\sum_{j}(\chi^{2}_{j})$ summing over all
chunks ($j$) as the final result.

###  3.2 Relative variations of chromospheric activity

Core emissions in strong optical spectral lines, such as the Ca II H & K and
H$\alpha$ lines, are proven proxies for magnetic flux on the Sun (Eberhard &
Schwarzschild, 1913). Their strengths as well as variations were used for
monitoring stellar activity levels and detecting long-term activity cycles
similar to the solar $11\,\mathrm{year}$’s cycle. For example, the S-index,
developed by the Mount Wilson Observatory HK Project (Duncan et al., 1991), is
a measure of the emission in the Ca II H & K line cores in the lower to middle
chromosphere due to magnetic heating, and was widely used in measuring the
stellar rotation period as well as the long-term activity cycles (e.g.
Hempelmann et al., 2016; Butler et al., 2017; Mittag, M. et al., 2018; Mittag
et al., 2019).

####  3.2.1 Approach

According to the wavelength coverage of spectra, we adopted the Ca II H & K,
H$\beta$ and H$\alpha$ lines to measure chromospheric activity. The spectral
subtraction technique (Barden, 1985; Montes et al., 1995) has been widely used
in measuring their emissions with the favour of comparison star or synthesized
spectrum (e.g. Montes et al., 1995; Gu et al., 2002; Cao et al., 2019).
However, noting the difficulty in finding the continuum around the Ca II H & K
lines (Shkolnik et al., 2005) while the requirement of precise measurement on
temporal variations of chromospheric activity indicators in our case, we
alternatively employed a scheme based on Xu et al. (2021), which is capable of
deriving relative variation of equivalent widths of chromospheric activity
indicators in series of spectra with high precision.

The scheme consists of following steps. $\mathbf{1}$. Delete outliers in each
spectrum contaminated by cosmic rays, telluric lines, CCD bad pixels, etc.
$\mathbf{2}$. Normalize the spectrum with highest SNR to unity as reference
spectrum $F_{R}$. $\mathbf{3}$. Fit other individual spectra
$F^{\prime}(\lambda^{\prime})$ by a product of $F_{R}(\lambda)$ and a
polynomial $P(\mathbf{a};\lambda-\lambda_{0}$), i.e.
$F^{\prime}(\lambda)=P(\mathbf{a};\lambda-\lambda_{0})F_{R}(\lambda)$ where
$\lambda_{0}$ is the reference wavelength, and then apply a shift
$\Delta\lambda$ due to RV, i.e. $\lambda^{\prime}=\lambda+\Delta\lambda$. Each
individual spectrum is inversely transformed to unity $F_{U}$ using derived
parameters $\mathbf{a}$ and $\Delta\lambda$. $\mathbf{4}$. The overall
spectrum $F_{M}$ is obtained as the median over a subset of spectra with high
SNR. Then the residual spectrum is $f=F_{U}-F_{M}$. $\mathbf{5}$.
Correspondingly the relative equivalent width is $\Delta W=W-W_{M}=-\int
fd\lambda$. By fitting $f$ with a Gaussian function,
$g(\lambda)=A\exp\left[-(\lambda-\lambda_{0})^{2}/(2\sigma^{2})\right]$, we
have $\Delta W=-\sqrt{2\pi}A\sigma$. Note that minus value means emission.

Figure 3 shows two example spectra observed by HIRES/Keck at 20060416 and
20110615 (table 4), when all of Ca II K, H, H$\beta$ and H$\alpha$ emisson
lines exhibited almost the most significant negative and positive residuals,
respectively. The relative variation was small compared to the overall spectra
but can be extracted with reliability by above scheme.

Figure 3: Examples of spectra observed at 20060416 (red line) and 20110615
(cyan line) reduced by our scheme. Each panel shows spectra $F^{\prime}$
derived by step $\mathbf{3}$ with the overall spectrum $F_{M}$ overlaid by
black dashed line, along with their respective residual spectra
$f=F^{\prime}-F_{M}$ derived by step $\mathbf{5}$.

By using subsets of Keck/HIRES spectra, Wright et al. (2004), Isaacson &
Fischer (2010) and Butler et al. (2017) measured S-index of HD 134319
independently. Comparison between our measurements and their results is shown
in figure 4, overlaid with respective linear fittings. Note that Wright et al.
(2004) employed a scheme similar to us by using the highest SNR observation as
a template in measuring the "sensitive differential S-values". Our result fits
well with the ones of Wright et al. (2004) and Isaacson & Fischer (2010),
indicating the reliability of our scheme in measuring the relative variation
of chromospheric activity indicators.

Figure 4: Comparisons between our measurements, $\Delta W_{\text{K}}+\Delta
W_{\text{H}}$, and previous measurements of S-index from Wright et al. (2004)
(W04), Isaacson & Fischer (2010) (I10) and Butler et al. (2017) (B17). Linear
fittings to the three cases are overlaid. The horizontal error bar represents
the typical uncertainty of our measurement.

####  3.2.2 Residual spectra

The parameters in above scheme were determined by numerical experiments and
listed in table 3. In step $\mathbf{2}$, The reference spectrum was the one
observed at Jan 27, 2008 (record 41 in table 4) by Keck/HIRES. In step
$\mathbf{3}$, the background portions, whose length should be short enough for
a good fitting by a low-order polynomial while long enough to include enough
absorption lines to measure RV shift, were chosen as a compromise and the gap
centred on respective emission core was temporally included for background
fitting. In step $\mathbf{4}$, the subset used in calculating the overall
spectrum $F_{M}$ consists of records 31–37, 39–41, 43, and 46–48 (table 4)
with high SNRs and common to all indicators. In step $\mathbf{5}$, the width
of Gaussian function, in accordance with $\sigma$, was estimated by plotting
residual spectra $f$ together and fixed as $\sigma=\mathrm{FWHM}/2.35482$ for
each index. The fitting interval within which Gaussian fitting was done was
used to estimate the fitting errors.

Table 3: Spectral fitting parameters. Index | $\lambda_{0}$ | Parameters in step $\mathbf{3}$ (Polynomial fit) | Parameters in step $\mathbf{5}$ (Gaussian fit)
---|---|---|---
| ($\mathring{A}$) | Background portions ($-\lambda_{0}$) | Order | Fitting interval ($-\lambda_{0}$) | $\sigma$
Ca II K | $3933.52$ |  $-7$ – $-0.75$ and $0.75$ – $7$ | $2$ |  $-7$ – $7$ | $0.140$
Ca II H | $3968.32$ |  $-7$ – $-0.75$ and $0.75$ – $7$ | $2$ |  $-7$ – $7$ | $0.140$
H$\beta$ | $4861.17$ |  $-8$ – $-2$ and $2$ – $8$ | $3$ |  $-8$ – $8$ | $0.425$
H$\alpha$ | $6562.55$ |  $-13$ – $-3$ and $3$ – $13$ | $3$ |  $-13$ – $13$ | $0.828$

Normalizations of Ca II H and K lines to unity through continuum division are
usually questionable due to the difficulty in finding their continua, an
alternative way was employed by fitting a straight line to the edges of
spectral portion (Shkolnik et al., 2005), $7\,\mathring{A}$ wide, centred on
emission cores, respectively. Figure 5 shows the residual spectra $f$ and
their respective overall spectra $F_{M}$ from Keck/HIRES. Emissions of
resonance from $F_{M}$ and variations of chromospheric activity from $f$ are
notable, indicating a high level of magnetic activity.

Figure 5: The residual spectra $f$ (up panel) and the overall spectrum $F_{M}$
(lower panel) ($f=F_{U}-F_{M}$) for Ca II K (left) and H (right). The typical
width of Gaussian fit (corresponding to $\sigma$) in measuring relative
equivalent width $\Delta W$ was over-plotted by black lines.

H$\alpha$ and H$\beta$ can be normalized to unity by general continuum
division with reliability. Moreover, with the high quality of H$\beta$ and
H$\alpha$ lines by Keck/HIRES, we recognized an improper trend between fitted
residuals and spectral absorption depth, i.e. the background portions in
residual $f$, which might result from improper reduction, e.g. scattered light
residuals, reflected by high SNR data. An offset is capable of fitting out
such a trend mathematically, i.e. $F^{\prime}(\lambda)\equiv
o+P(\mathbf{a};\lambda-\lambda_{0})F_{R}(\lambda)$, and was included in
processing H$\beta$ and H$\alpha$ of Keck/HIRES spectra. Figure 6 shows
example of H$\beta$ with the most notable difference between with and without
offset observed at Jan 14, 2003.

Figure 6: Comparison between residuals $f$ with and without an offset $o$
included in fitting $F^{\prime}(\lambda^{\prime})$ by
$P(\mathbf{a};\lambda-\lambda_{0})F_{R}(\lambda)+o$ versus $1-F_{M}$: the
example of H$\beta$ was observed by Keck/HIRES on Jan 14, 2003. The overlaid
larger triangles and circles, horizontally misaligned from each other, were
sliding average of respective smaller symbols, the error bar represents the
standard deviation of corresponding interval. The overall spectrum is plotted
in sub-window with background portions in table 3 overlaid with black lines.

Figure 7 shows the residual spectra $f$ and respective overall spectrum
$F_{M}$ of H$\beta$ and H$\alpha$ lines, from which we can also see variable
emissions with different amplitudes.

Figure 7: Same as figure 5, but for H$\beta$ (left) and H$\alpha$ (right).

##  4 Results and Discussion

It’s unfeasible to decide the real number of starspots from disk integrated
photometric observation. The number of dips in LC corresponds only to the
lower limit of starspots number, because LC can always be fitted by model with
as many starspots as one likes (Jeffers & Keller, 2009; Basri & Shah, 2020).
One has to adopt a priori compromise to reduce the free parameter space and
increase the possibility of convergence in optimization so as to obtain a
reasonable solution in fitting. For example, LCM simulates the LC by an
analytically model with a few circle spots, while light curve inversion (LCI)
usually employs either the maximum entropy or Tikhonov criterion for a unique
and stable solution (e.g. Lanza et al. (2006)). Thus, if the photosphere of
HD134319 is covered by or equivalent to two dominant starspots (or active
regions) versus surroundings with non-spot or stable distributed starspots
contributing no rotational modulation, our result should reveal the
distribution and relative evolution of starspots, or else it would rather be
the description of hemispheric asymmetry.

###  4.1 Results comparison between different inclinations

Using the best-fit and manually checked solutions from many independent runs
for each chunk, we were able to increase the probability of optimization to
fall into global minima in limited iterations, and thus to trace the sizes and
locations of starspots over the entire span of TESS data with reliability.

Practically, each chunk was modelled independently by $50$ runs for
inclination $i=75\degr$ and $40$ runs for inclination $i=90\degr$. The
solutions could be divided into a few groups characterized by special spot
configuration. In most cases it was easy to elect the globally optimal
solution by choosing minimal residual $\chi^{2}$, however outliers existed
where local minima had equally good or even better fit than the predicted one.
To overcome this problem we additionally assumed that spot configuration among
adjacent chunks should be similar to each other as they shared about $80\%$
data points.

By comparison, figure 8 displays the spot longitudes, latitudes and radii
corresponding with $i=75\degr$ and $i=90\degr$ respectively. Under assumption
of linear evolution, we took the value of spot radius at mid-epoch of each
chunk in plot and discussion hereafter.

Figure 8: Temporal variations of spot longitudes (top), latitudes (middle) and
radii (bottom) for inclination $i=75\degr$ (green crosses and red diamonds)
and $i=90\degr$ (blue pluses and orange circles), respectively. The radius of
each spot plotted, hereafter, takes its value at mid-epoch of respective chunk
under the assumption of linearly evolving size. "P" represents the larger spot
(primary) and "S" represents the smaller spot (secondary) in each case. The
sizes of symbols are proportional to the spot radii.

The spot longitudes were derived with high stability for both cases as
expected due to both the strong constraint of photometry on rotational phase
and the orthogonality between stellar inclination and longitudes. Two active
longitudes can be easily recognized i.e. the primary one (named "P") centring
around longitude $\lambda\sim 220\degr$ and the secondary one (named "S")
varying around longitude $\lambda\sim 130--50\degr$.

Two features can be found from the distribution of spot latitudes. First, the
case of inclination $i=75\degr$ generally shows spots at lower latitudes than
$i=90\degr$, which implies a mathematical deviation of configuration
corresponding to input inclination, i.e., an evidence of parameter degeneracy.
Second, the result shows similar relative variations of spot latitudes between
input inclinations, which indicates a possible constraint of high accuracy by
TESS photometry on spot latitudes versus noise, to some extent.

A degeneracy between spot radius and latitude is recognizable in both cases.
Figure 9 shows the distribution of spot radius versus its latitude for "P" and
"S" in different time ranges for inclination $i=90\degr$, and their comparison
to simulations for a supposed spot at latitude $\beta=70\degr$ initially with
series radii $\gamma=13,20,23$ and $28\degr$ through setting the spot’s
latitude ranging from $40$ to $85\degr$. The highly coincidence between
results and simulations within time range 1869–1995 indicates probable
existence of degeneracy, while the apparent deviation within time range
1683–1764 from such a coincidence means reliable variation of either spot
radius or latitude.

Figure 9: Spot radius–latitude ($\gamma-\beta$) correlation from best-fit
solutions, different symbols represent spots "P" and "S" at different time
ranges. The simulations for an artificial spot at latitude $\beta=70\degr$
initially with series radii $\gamma=13,20,23~{}\text{and}~{}28\degr$ through
varying its latitude value from $\beta=40\degr$ to $85\degr$ are overlaid with
dash-dot lines.

Figure 10 shows one fitting example in each sector for cases of inclination
$i=75\degr$ and $i=90\degr$. Both cases reveal that good fits have typical
residuals around $0.0001-0.0002$, implying that the spot configuration on HD
134319 can be reasonably described by two condensed active regions. As a
whole, fitting of $i=75\degr$ resulted in slightly larger $\chi^{2}$ than
$i=90\degr$ in most chunks, which indicates an extremely high inclination of
HD 134319, coinciding with the prediction from spectroscopy. Hereafter spot
modelling result with inclination $i=90\degr$ is employed in discussion.

Figure 10: Illustration of modeling on inclinations $i=75\degr$ (blue dashed
line) and $i=90\degr$ (red line) for TESS LC (cyan dots), one chunk for each
sector, with epoch of each chunk indicated. The horizontal axis of each panel
was modified to $0-5$ by subtracting its epoch. Residuals are over-plotted at
value of $0.96$.

###  4.2 Starspot Evolution and Differential Rotation

####  4.2.1 Configuration of two active regions

In figure 11 we reproduce the continuous phase versus time evolution map of
flux, with spot longitudes from the best-fit solution and their linear
fittings overlaid. Two active regions characterized by their amplitudes of dip
in brightness, i.e. the primary "P" located at phase about $0.71$ and the
secondary "S" located at phase of about $0.34$, separating by $140\degr$ in
longitude, are revealed to be located at high latitudes varying between $50$
and $80$ degrees (figure 8). The solutions also indicate the difference
between spot sizes, and the smaller one evolves more rapidly than the bigger
one, despite the possible degeneracy between latitude and radius (figure 9).

Figure 11: Continuous phased LC map of the entire TESS long cadence data set
with period $P=4.436391\,\mathrm{day}$, with the spot longitudes from our two
spot model as well as their linear fittings overlaid. The pixel shade, from
dark to light, indicates the flux in each (time/phase) chunk, and is
interpolated and plotted as contour map. Vertical white gaps correspond to
times with no TESS data, and the long-time gap of observation are marked as
breaks in horizontal axes. The longitude corresponding to the primary active
region, "P", and its linear fitting are overlaid with orange circles and red
line, while the longitude corresponding to the secondary active region, "S",
and its linear fittings to the first and last half data separately, are
overlaid with blue diamonds, blue dashed and dark dash-dot lines,
respectively. The symbols are scaled proportional to the spot radii.

In literature, the photometric study of HD 134319 was presented by Messina et
al. (1998) and Messina & Guinan (1998), who analysed HD 134319 LC from 1991 to
1995 and reported a two spot configuration existed over $5$ years, despite the
large gap between sections. Such a configuration was widely found on other
stars with different rotation periods. Rapid rotator AB Dor ($P\sim
0.51479\,\mathrm{day}$, K0V, $70\sim 100$ Myrs) was reported to maintain two
long-lasting longitudes over tens of years from long-term photometry (e.g.
Berdyugina & Järvinen, 2005; Ioannidis & Schmitt, 2020) and evolving SDR from
Dopple imaging (e.g. Donati & Collier Cameron, 1997; Petit et al., 2002;
Collier Cameron & Donati, 2002; Jeffers et al., 2007) as well as ultra-active
phenomena like superflares (Schmitt et al., 2019). The rapidly rotating M4
type star GJ 1243 ($P\sim 0.59\,\mathrm{day}$, mass $\sim 0.24\,M_{\sun}$) was
found to be dominated by one stable and another evolving active longitudes
during Kepler and TESS epoch (Davenport et al., 2020). Detailed studies on FK
Com ($P\sim 2.4\,\mathrm{day}$, G4III) revealed two active regions with
apparently different sizes and their "flip-flop", phase-jumps, spot emerges,
drifts and so on (e.g. Korhonen et al., 2001; Oláh et al., 2006; Hackman et
al., 2013). On BY Dra type star LQ Hydrae ($P\sim 1.6\,\mathrm{day}$,
$\text{age}\sim 60\,\mathrm{Myr}$, K1V), long-term photometry revealed
dominance of two active longitudes over 20 years so that "the first evidence
of flip-flop on single dwarf" and multiple activity cycles, in both non-
axisymmetric and axisymmetric (as in the Sun) modes, in time-length of years
were reported (Berdyugina et al., 2002; Lehtinen et al., 2012). Another
analogy, Kepler-17 (with slower rotation, $P\sim 12.01\,\mathrm{day}$, G2V)
also showed two active longitudes, one of them survived over Kepler’s
observation period, and an evidence of cycle (Lanza et al., 2019). Thus one
can reasonably expect that such a two spot configuration exists widely on
stars, and we can suggest an internal non-axisymmetric geometry for HD 134319,
possessed by almost constant mean field configuration which was believed to
indicate a small differential rotation (Messina & Guinan, 1998).

####  4.2.2 Rotational modulation periods of starpots

As shown in figures 8 and 11, the spot configuration on HD 134319 revealed
long-lasting features as well as significant evolution in their locations and
sizes. For easier discussion, we divided the observation time into two
durations, i.e. "T1", from 1683 to 1764 (sectors 14–16) and the other "T2",
from 1869 to 1955 (sectors 21–23). The primary region "P" tended to be a long-
lasting feature undergoing slow evolution from "T1" to "T2", while the
secondary region "S" migrated not too much in "T1", like "P", but exhibit a
much rapid evolution in "T2".

The rotational modulation period of individual spot can be measured from the
phased LC (Davenport et al., 2015) by fitting its longitude with a linear
function

$P_{i}=P_{0}/(1-m_{i}P_{0})$ (4)

where $P_{0}$ is the phase folding period, the subscript $i$ represents the
index of spot and $m_{i}$ is the slope. Note by this definition a negative
slope yields a smaller value than $P_{0}$.

The fittings are overlaid in figure 11 with three lines. Fitting of "P" gave
an average rotational modulation period of $P_{\mathrm{P}}=4.4236654\pm
0.0003350\,\mathrm{day}$, starting at phase about $\lambda_{\mathrm{P}}=0.71$
at epoch 1683 (red line). Fitting of "S" in "T1" gave
$P_{\mathrm{S1}}=4.4226344\pm 0.0016632\,\mathrm{day}$ and
$\lambda_{\mathrm{S1}}=0.34$ at epoch 1683 (blue dashed line). And fitting of
"S" in "T2" gave $P_{\mathrm{S2}}=4.3580376\pm 0.0064554\,\mathrm{day}$ and
$\lambda_{\mathrm{S2}}=0.19$ at epoch 1689 (dark dash-dot line).
$P_{\mathrm{P}}$ and $P_{\mathrm{S1}}$ were estimated with high confidence due
to their stabilities, however $P_{\mathrm{S2}}$ was determined with larger
uncertainty due to its rapid evolution. Note that these values are smaller
than the phase folding period, i.e. $P_{0}=4.436391\,\mathrm{day}$, derived by
GLS determination of the whole LC.

####  4.2.3 Starspot evolution

The deviation of spot longitude from its rotational modulation period
represents the spot longitudinal migration from average location, as shown in
figure 12. In "T1", the two spots exhibited highly synchronized oscillatory
variations with amplitude of about $15\degr$ and period of about $40$ days.
This phenomenon can also be seen in figure 2 where the two dips in LC varied
periodically with time. In "T2", "P" remained the same level of evolution like
before, however "S" underwent a much sudden evolution than "T1", i.e. it had
an obviously small average rotational period than "P" and exhibited a much
larger longitudinal oscillation with amplitude of about $40\degr$ and period
of $30-40$ days, which broke the previous synchronization and seemed to be a
new distribution of the spot feature. This indicates that "P" survived over
time spanning from "T1" to "T2", while "S" survived over "T1" but might be a
new one in "T2".

Figure 12: Subtraction of spot longitude ($\lambda$) by its linear fitting,
$\lambda-$(linear fitting), i.e. spot longitudinal migration from its average
location, versus time. The symbols are scaled proportional to the spot radii,
as in figure 11.

On the other hand, although correlated with recovered latitudes as shown in
figure 9, the spot sizes were found to evolve with time in terms of relative
variations. Firstly, compared with a supposed spot located at average latitude
$\beta=70\degr$, "P" has a maximal radius of about $28\degr$ in "T1" and
decreases to about $20\degr$ in "T2", while "S" has an average radius of about
$23\degr$ in "T1" and a smaller radius of only $13\degr$ in "T2". So the spot
configuration was likely to be dominated by one larger plus a smaller active
region, and both of them decreased to smaller scales from "T1" to "T2", in the
sense of projected area. Secondly, as shown in figure 8, "S" radius in "T1"
increased rapidly around epoch $1692$ from $26\degr$ to $31\degr$ and then
decreased gradually to small value of about $12\degr$ at the end of "T1",
which was in similar size as in "T2", while "P" underwent an opposite
variation which had a rapid decrease at epoch $1692$ from $29\degr$ to
$22\degr$ followed by a gradually recovery to large size. For the case of
"T2", in contrary, the two spots showed synchronized variation of radii over
time.

Considering the spot radius-latitude correlation, the variations in "T2" might
result from the parameter degeneracy due to weak constraints of photometry,
however one can still recognize two evidences of spot evolution in "T1" with
reliability. One is the simultaneous but opposite variations between "P" and
"S" around $1692$, which indicates either a switch of the activity level from
one spot to another or the latitudinal migration to opposite directions. This
phenomenon was also reported and analysed in detail on HQ Hydrae by Lehtinen
et al. (2012) who discriminated the "flip-flop" from more commonly "switch"
phenomenon. The other is the decrease of "S" since $1720$.

Figure 13 presents the long-term evolution of the LC mean brightness $M$ and
peak-to-peck amplitude $A$, which were usually employed in estimating the
stellar activity variations due to spots (Berdyugina et al., 2002; Lehtinen et
al., 2012; Ioannidis & Schmitt, 2020). The variation in $M$ could be explained
by activity difference, or spot size variations, between "T1" and "T2". While
the decrease in $A$ attributes mainly to the decay of spot "P" to small size
from "T1" to "T2", supporting the above analysis.

Figure 13: The long term variation of the mean brightness $M$ (top) and peak-
to-peak amplitude $A$ (bottom) of HD 134319 measured from figure 2.

Considering the coherence of the spot longitudinal migration and radial
evolution, the primary active region "P" is likely to be the same one which
migrated and oscillated from phase $0.71$ and lived longer than the
observation season ($280$ days), while the secondary active region "S"
exhibited different types of evolution and thus is likely to represent
different spot features at least between "T1" and "T2". This inference is in
accordance with the scenario that spot lifetime is proportional to its size
and large spot can survive for many years (see reviews of Berdyugina (2005)
and Strassmeier (2009)).

####  4.2.4 Differential rotation

Despite the weak constraint on spot latitude, the multiple rotational
modulation periods in abundant photometric time series were widely taken as
indicator of SDR on many stars other than the Sun. By assuming the latitudinal
occupation of spots towards both equator and pole as much as possible, the
seasonal period variation from extreme long-term photometry was applied to
estimate the lower-limit of SDR (e.g. Henry et al., 1995; Messina & Guinan,
2003; Balona & Abedigamba, 2016). At the mercy of precise photometric data,
there exist studies attempting to inverse spot rotational modulation period
and latitude simultaneously and thus directly estimating the SDR (e.g. Croll
et al., 2006; Fröhlich et al., 2012; Lanza et al., 2016). However, studies
aiming at examining the reliability of SDR detection based on photometry alone
revealed the probability of misleading results especially in evolving stages
(Aigrain et al., 2015; Basri & Shah, 2020).

The difference of spot rotational modulation periods between "P" and "S" in
"T1" is quite small, i.e. $\Delta P=0.0010310\pm 0.0016966\,\mathrm{day}$,
corresponding to a lower-limit of SDR as $\alpha=0.000233$. Meanwhile, the
difference between "P" and "S" in "T2" yields $\Delta P=0.0656278\pm
0.0064641\,\mathrm{day}$, corresponding to a lower-limit of SDR as
$\alpha=0.0148$, which is apparently smaller than prediction from statistics
(Reinhold et al., 2013; Balona & Abedigamba, 2016) and dynamo studies (Küker &
Rüdiger, 2008; Kitchatinov & Olemskoy, 2012). Such a small difference is
probably due to the nearby latitudinal distribution of spots during the
observation season (or, spot pair is located on opposite hemispheres, which is
indistinguishable in the case of inclination $i=90\degr$). As in figure 8, the
best-fitted spot latitudes were close to each other, i.e. between
$\beta=50\degr$ and $85\degr$ and varied temporally not too much, supporting
above situation.

Due to the weak constraint of photometric observation on latitudinal
information and the possible migration of starspots, the estimation of SDR is
pretty difficult. The lower limit of SDR shear is reliable only when the
magnetic configuration can be described by such a two-spot model and all
starspots are centered in certain longitudes, however both of which cannot be
confirmed from collected data, and thus the above discussion should be treated
with caution.

###  4.3 Chromospheric activity in short and long timescales

The relative equivalent widths ($\Delta W$s) of Ca II H and K, H$\beta$ and
H$\alpha$ lines are listed in table 4 and plotted in figure 14. Observations
by OHP/ELODIE within about one year have good phase coverage, the measurements
yield typical median uncertainties of $0.036088$ for H$\beta$ and $0.015585$
for H$\alpha$, except two extremely low SNR records 23 and 24. Observations by
Keck/HIRES have long time baseline over 14 years but are sparsely sampled, the
measurements yield typical median uncertainties of
$0.025996,~{}0.018091,~{}0.006445,~{}\text{and}~{}0.005927$ for Ca II H and K,
H$\beta$ and H$\alpha$, respectively.

Table 4: The relative equivalent widths measured from the residual spectra, which represent the difference between individual spectrum and the overall spectrum for each chromospheric activity indicator (see section 3.2). The last colume is the instrument, "E" represents OHP/ELODIE and "H" represents Keck/HIRES. No. | Date | Epoch | $\Delta W_{\text{K}}$ | $\Delta W_{\text{H}}$ | $\Delta W_{\text{H}\beta}$ | $\Delta W_{\text{H}\alpha}$ | Inst.
---|---|---|---|---|---|---|---
| yyyymmdd | MJD-50000 | ($\mathring{A}$) | ($\mathring{A}$) | ($\mathring{A}$) | ($\mathring{A}$) |
1 | 19951106 | 27.758241 | … | … | -0.031373$\pm$0.040562 | -0.018857$\pm$0.015511 | E
2 | 19960503 | 206.975949 | … | … | -0.029332$\pm$0.027837 | -0.033498$\pm$0.013834 | E
3 | 19960503 | 206.989294 | … | … | -0.030835$\pm$0.026398 | -0.030377$\pm$0.013863 | E
4 | 19960504 | 208.035718 | … | … | -0.027873$\pm$0.021146 | -0.029040$\pm$0.010667 | E
5 | 19960505 | 209.994572 | … | … | -0.015346$\pm$0.020592 | -0.032421$\pm$0.010319 | E
6 | 19960506 | 210.940266 | … | … | -0.017242$\pm$0.026576 | -0.013584$\pm$0.010882 | E
7 | 19960508 | 212.998855 | … | … | -0.040752$\pm$0.032760 | -0.021943$\pm$0.016166 | E
8 | 19960628 | 262.953877 | … | … | -0.018798$\pm$0.036594 | -0.001461$\pm$0.016825 | E
9 | 19960628 | 262.962755 | … | … | -0.025168$\pm$0.039243 | -0.004372$\pm$0.017009 | E
10 | 19960629 | 263.881146 | … | … | -0.035705$\pm$0.035581 | -0.022479$\pm$0.016555 | E
11 | 19960629 | 263.890013 | … | … | -0.039269$\pm$0.037212 | -0.027851$\pm$0.015313 | E
12 | 19960630 | 264.888473 | … | … | -0.030101$\pm$0.033388 | -0.038060$\pm$0.015731 | E
13 | 19960630 | 264.897384 | … | … | -0.018138$\pm$0.038456 | -0.034214$\pm$0.016794 | E
14 | 19960701 | 265.879850 | … | … | -0.011007$\pm$0.031238 | -0.006143$\pm$0.014964 | E
15 | 19960701 | 265.892223 | … | … | -0.020166$\pm$0.038616 | -0.014393$\pm$0.017049 | E
16 | 19960702 | 266.930995 | … | … | -0.051228$\pm$0.070634 | -0.042379$\pm$0.024720 | E
17 | 19960702 | 266.939896 | … | … | -0.021975$\pm$0.074569 | -0.004292$\pm$0.023377 | E
18 | 19960703 | 267.924133 | … | … | -0.006526$\pm$0.037805 | 0.009826$\pm$0.015257 | E
19 | 19960703 | 267.932987 | … | … | -0.023558$\pm$0.042209 | 0.010922$\pm$0.015659 | E
20 | 19960828 | 323.834537 | … | … | -0.024840$\pm$0.051999 | -0.022079$\pm$0.017712 | E
21 | 19960829 | 324.823485 | … | … | -0.026197$\pm$0.029949 | -0.025716$\pm$0.012245 | E
22 | 19960901 | 327.855637 | … | … | -0.031695$\pm$0.033460 | -0.015077$\pm$0.013060 | E
23 | 19961229 | 447.207779 | … | … | -0.197406$\pm$0.339053 | 0.181427$\pm$0.068684 | E
24 | 19970125 | 474.192547 | … | … | -0.128411$\pm$0.360773 | 0.169545$\pm$0.141722 | E
25 | 19970128 | 477.205753 | … | … | -0.019608$\pm$0.031867 | 0.001652$\pm$0.019156 | E
26 | 19970129 | 478.191632 | … | … | -0.022774$\pm$0.019649 | -0.019219$\pm$0.009338 | E
27 | 19990424 | 1292.561160 | -0.084215$\pm$0.069323 | -0.034360$\pm$0.031120 | -0.035413$\pm$0.007033 | … | H
28 | 19990425 | 1293.526583 | -0.001237$\pm$0.068204 | -0.013785$\pm$0.054094 | -0.026791$\pm$0.011333 | … | H
29 | 19990519 | 1317.430287 | 0.008699$\pm$0.045056 | 0.000250$\pm$0.036867 | 0.001960$\pm$0.009238 | … | H
30 | 19990808 | 1398.261886 | 0.023564$\pm$0.053873 | 0.000142$\pm$0.041517 | 0.007955$\pm$0.010172 | … | H
31 | 20030114 | 2653.679294 | 0.054104$\pm$0.027859 | 0.031481$\pm$0.017816 | 0.022237$\pm$0.005083 | … | H
32 | 20030314 | 2712.478795 | -0.025440$\pm$0.028476 | -0.012997$\pm$0.018487 | -0.014811$\pm$0.005179 | … | H
33 | 20030616 | 2806.324191 | 0.004558$\pm$0.017696 | 0.005198$\pm$0.010904 | 0.000763$\pm$0.003301 | … | H
34 | 20030616 | 2806.375226 | 0.020166$\pm$0.024592 | 0.007891$\pm$0.017010 | 0.010546$\pm$0.006181 | … | H
35 | 20030713 | 2833.330383 | 0.004602$\pm$0.016301 | 0.001800$\pm$0.010753 | 0.005014$\pm$0.004157 | … | H
36 | 20030728 | 2848.266842 | 0.043567$\pm$0.016419 | 0.022013$\pm$0.013506 | 0.019795$\pm$0.004684 | … | H
37 | 20040625 | 3181.446408 | -0.091720$\pm$0.017171 | -0.058185$\pm$0.011094 | -0.025069$\pm$0.004413 | … | H
38 | 20040710 | 3196.314068 | 0.002408$\pm$0.053364 | -0.010890$\pm$0.043795 | 0.006985$\pm$0.009848 | … | H
39 | 20050225 | 3426.617754 | -0.023167$\pm$0.022615 | -0.017766$\pm$0.015543 | -0.001114$\pm$0.005827 | -0.003765$\pm$0.005488 | H
40 | 20060416 | 3841.393636 | -0.038543$\pm$0.022988 | -0.027235$\pm$0.017408 | -0.021964$\pm$0.006786 | -0.028862$\pm$0.005404 | H
41 | 20080127 | 4492.683551 | 0.006933$\pm$0.019474 | 0.005565$\pm$0.013863 | 0.009798$\pm$0.005215 | 0.030829$\pm$0.004608 | H
42 | 20090604 | 4986.408533 | 0.033238$\pm$0.038755 | 0.020317$\pm$0.027897 | 0.012719$\pm$0.010399 | 0.055268$\pm$0.008783 | H
43 | 20100204 | 5231.671857 | -0.034793$\pm$0.026827 | -0.022060$\pm$0.017126 | -0.009158$\pm$0.008149 | -0.020032$\pm$0.005889 | H
44 | 20110615 | 5727.448702 | 0.050378$\pm$0.044145 | 0.023823$\pm$0.030052 | 0.028416$\pm$0.012267 | 0.094395$\pm$0.010165 | H
45 | 20110615 | 5727.458635 | 0.058932$\pm$0.042421 | 0.033285$\pm$0.027970 | 0.018301$\pm$0.010821 | 0.063508$\pm$0.008704 | H
46 | 20120305 | 5991.668097 | -0.012412$\pm$0.025164 | -0.006934$\pm$0.018366 | 0.002887$\pm$0.006191 | 0.002646$\pm$0.006686 | H
47 | 20120601 | 6079.331947 | 0.011411$\pm$0.019098 | 0.003193$\pm$0.013857 | 0.008064$\pm$0.005063 | 0.048057$\pm$0.005655 | H
48 | 20130709 | 6482.252854 | -0.043213$\pm$0.024948 | -0.020351$\pm$0.022389 | -0.003531$\pm$0.006699 | -0.010131$\pm$0.005964 | H
Figure 14: The relative equivalent widths versus time as table 4. The linear
fitting on H$\beta$ is overlaid with dashed line, and the time durations used
for analysing rotational modulation (figure 16) are marked with dark lines.
Hereafter, the axis corresponding to equivalent widths is inverted for easier
view.

####  4.3.1 Long-term evolution

As shown in figure 14, during observations in near 20 years, the chromospheric
activity indicators show variations exceed $0.1\mathring{A}$. Evolutions in
both short and long timescale are discriminable by comparing their relative
variations of equivalent widths in different timescales. The short-term
evolution within one or a few rotations can be attributed to the rotational
modulation as revealed by OHP/ELODIE observation with good phase coverage,
which will be discussed in more detail in section 4.3.3. While the notable
variation in Keck/HIRES observing season shows the maximal scale, indicating
evolution in long timescale. Linear fitting of H$\beta$ (overlaid by green
dashed line) reveals a trend of its average $\Delta W$ from $-0.012$ to
$0.011\mathring{A}$ with time, implying a decrease of activity level in long
timescale.

####  4.3.2 Correlations between indicators

Correlations between chromospheric activity indicators were widely studied on
active stars (e.g. Strassmeier et al., 1990; Montes et al., 1995; Cincunegui
et al., 2007; Martínez-Arnáiz et al., 2011; Scandariato et al., 2017). The
excess emissions of Ca II H and K ($3968$ and $3933\,\mathring{A}$) resonance
were found to be tightly correlated with each other and thus were practically
measured together, for example, the widely used S-index (Duncan et al., 1991),
in estimating the stellar activity. Two other well studied chromospheric
proxies, H$\alpha$ ($6363\,\mathring{A}$) and H$\beta$ ($4861\,\mathring{A}$)
lines, which can be feasibly measured with high SNR on even pretty cool stars
where the Ca II resonance becomes progressively faint, resulted from the two
most probable transitions of electrons between energy levels of Hydrogen which
is the most abundant element inside the Sun and other stars, were also found
to be correlated with each other. However, correlation between H$\alpha$ and
Ca II resonance was more complicated, because H$\alpha$ reveals not only a
non-monotonically variation with increasing heating rate (Linsky, 2017) but
also the dependence on stellar physical parameters such as pressure and
temperature (Cincunegui et al., 2007).

Figure 15 shows the correlation between $\Delta W$s of Ca II K and H, H$\beta$
and H$\alpha$ lines. Linear fitting to each relation was done for estimating
the relative difference between indicators and overlaid in this figure by
dashed line. Note by the linear fitting a positive slope yields a positive
correlation. We can find that the relative variation of any indicator is
strongly correlated with the others, despite the quantitatively different
uncertainty levels.

Figure 15: Relative equivalent widths of Ca II K versus H (top), Ca II K
versus H$\beta$ (middle) and H$\alpha$ versus H$\beta$ (bottom), the
uncertainties are plotted by horizontal and vertical errorbars, respectively.
Weighted linear fittings are overlaid with dashed lines, the ratio of
$W_{\text{K}}/W_{\text{H}}\sim 1.698487$ is overlaid with dark dots.

Due to the lack of comparing star as template, it is difficult to estimate the
absolute value of H$\beta$ and H$\alpha$ emissions, while we can estimate
equivalent widths of Ca II H and K lines due to their narrow emission
profiles. Measurements gave $W_{\text{H}}=-0.342627\,\mathring{A}$ and
$W_{\text{K}}=-0.581948\,\mathring{A}$ for Ca II H and K lines, respectively,
by directly integrating their overall spectra over $1.1\,\mathring{A}$ width
portions centred on the emission cores, and yielded a ratio of
$W_{\text{K}}/W_{\text{H}}=1.698487$ which is close to the slope of linear
fitting (top panel in figure 15). This infers that the relative variations of
Ca II H and K lines are proportional to their respective absolute emission
strengths.

The larger variation of H$\alpha$ than H$\beta$ in bottom panel of figure 15
could be attributed to its lower energy required in electron transitions in
principle, while the more commonly mentioned correlation between Ca II and
H$\alpha$ lines was found to depend on stellar parameters such as the stellar
effective temperature, metallicity and pressure etc. (e.g. Cincunegui et al.,
2007; Walkowicz & Hawley, 2009; Scandariato et al., 2017). By analogy to the
study on the Sun by Gebbie & Steinitz (1974), who proposed that the formation
of Ca II is direct collision dominated due to turbulent velocities while the
formation of H$\alpha$ is photoionization dominate due to radiation which
decreases in late-type stars, one could reasonably expect different
correlations between Ca II and H$\alpha$ on different type of stars. Thus the
ratio between Ca II and H$\alpha$ (simply as $(\Delta W_{\mathrm{K}}+\Delta
W_{\mathrm{H}})/\Delta W_{\mathrm{H}\alpha}\sim 2.238$) is larger than one,
might indicate that the chromospheric excitation on HD 134319 is dominated by
collision due to turbulent velocities.

####  4.3.3 Rotational modulation of indicators

In short timescale, variations of chromospheric proxies under rotational
modulation were widely reported on stars (e.g. Wright et al., 2004; Isaacson &
Fischer, 2010; Flores Soriano & Strassmeier, 2017). Such variations can also
be found on HD 134319. Figure 16 shows the phase folded $\Delta W$s of Ca II H
and K, H$\beta$ and H$\alpha$ within chosen time durations no longer than half
a year. Among the five durations, (a) and (b) have good phase coverage in
about one rotation cycle, while (c) – (f) are phase–folded for data collected
in months. With the good phase coverage of OHP/ELODIE observations, two peas
around phases $0.2$ and $0.7$ are recognizable from panels (a) – (b) during
time from May to July of 1996, corresponding to two active regions
characterized by enhancement of chromospheric emissions. Later, for the
precise data observed by Keck/HIRES in year 1999 and 2003, at least one active
region around phase $0.4$ can be recognized, despite the poor phase coverage.

Figure 16: Phase–folded relative equivalent widths ($\Delta W$s) of Ca II K
and H, H$\beta$ and H$\alpha$ lines in five chosen time durations shorter than
half a year (marked in figure 14). Peaks in the phased curves represent the
phase of enhanced chromospheric emissions under rotational modulation.
Moreover, the y-band LC taken from figure 3 of Messina et al. (1998) is re-
folded by newly measured rotation period and plotted in (a) and (b) with dark
dots, overlaid by smoothed dotted lines, and is inverted so that the possible
spatial connection between chromospheric emissions and photospheric spots is
manifested as the positive correlation between inverted $\Delta W$ and LC.

The long-term behaviour of the chromospheric emissions under rotational
modulation implies that HD 134319 was dominated by two active regions during
years from 1996 to 1997, and the pea around phase $0.2-0.4$ might be due to a
long-lived active region lasting over years from 1996 to 2003.

####  4.3.4 Connection between chromospheric activity and photospheric spots

The possible connection between chromospheric activity and photometric
variation would indicate a configuration of the magnetic field spreading over
photosphere and chromosphere. A clear anti-correlation between them on LQ
Hydrae was reported and interpreted as the spatial connection between
photospheric dark spots and chromospheric plages (Cao & Gu, 2014; Flores
Soriano & Strassmeier, 2017), the similar phenomenon was also reported on FK
Comae (Vida et al., 2015).

It is difficult to analyse this relation exactly on HD 134319 due to the lack
of overlap between photometric and spectroscopic observations analyzed in this
paper. In figure 16 (a) and (b), we over-plot the phased LC taken from figure
3 of Messina et al. (1998) observed around 1995.38, which is the closest to
the spectroscopic data in time, seperated by more than one year. A positive
correlation between $\Delta W$ and inverted photometric LC can be found but at
low significance. Thus no reliable conclusions can be concluded with collected
data.

##  5 Summary

In this paper, we present the analysis of the starspot configuration,
evolution and chromospheric activity on HD 134319, by employing the high
precision photometry by space-based TESS telescope in sectors 14–16 (epoch
1683–1764, called "T1") and 21–23 (epoch 1869–1955, "T2") and the
spectroscopic data observed by OHP/ELODIE and Keck/HIRES from year 1995 to
2013.

We firstly measured the rotation period of the star from TESS LC as
$P=4.436391\pm 0.00137$ days which is used in subsequent analysis, using GLS
method (Zechmeister & Kürster, 2009). Then, the starspot configurations at
series of epochs were derived by splitting the whole TESS LC into 129 chunks
in $5$ days length and then fitting them separately with a two-spot model.
Besides, based on analysis of the high-resolution spectroscopic data, we
derived the relative equivalent widths ($\Delta W$s) of the Ca II H and K,
H$\beta$ and H$\alpha$ lines to investigate its chromospheric activity.

We focus on magnetic characteristics on HD 134319 using a simple two-spot
model. One should keep in mind that the actual number of starspots (or active
regions) could not be decided with collected data and might much larger than
the number of dips in LC. However our model revealed reasonable good fits, and
if the star really has only two spots or can be described by such a model, it
would be the result. Our main results are summarized as follows:

(I) As revealed by analysis of TESS LC data, a two-spot configuration, i.e., a
long-lasting primary spot "P" plus a secondary spot "S", was capable of
explaining the LC variation on HD 134319 during observation. Furthermore, the
primary spot was likely to survive at least longer than TESS observation
duration, i.e. 300 days, and might survive over years.

(II) Linear fittings on spot longitudes give the average rotational modulation
period of spots as: $P_{\mathrm{P}}=4.4236654\pm 0.0003350\,\mathrm{day}$
emerging at phase $\lambda_{\mathrm{P}}=0.71$ at epoch 1683,
$P_{\mathrm{S1}}=4.4226344\pm 0.0016632\,\mathrm{day}$ emerging at phase
$\lambda_{\mathrm{S1}}=0.34$ at epoch 1683 and $P_{\mathrm{S2}}=4.3580376\pm
0.0064554\,\mathrm{day}$ emerging at phase $\lambda_{\mathrm{S2}}=0.19$ at
epoch 1689. This corresponds to a low-limit of SDR, if starspots underwent no
migration at long-time scale, as $\alpha=0.0148$ which is weaker than
prediction, which might due to a nearby distribution of spot latitudes.

(III) Reliable evolutions of spot radii can be derived despite the radius-
latitude degeneracy (figures 8 and 9). A sudden increase of "S" from radius
$\gamma=26\degr$ to $31\degr$ and simultaneous sudden decrease of "P" from
$\gamma=29\degr$ to $22\degr$ around epoch $1692$ indicate an exchange of
activity strength between spots. Since epoch $1692$, "S" radius underwent a
gradual decrease to about $12\degr$ at the end of "T1", similar to its radius
in "T2". Besides, decrease of spot radii from epoch $1683-1764$ to $1869-1955$
was notable (figure 9). Comparing with a reference spot at latitude
$\beta=70\degr$, "P" had radius $\gamma\sim 28\degr$ and "S" had radius
$\gamma\sim 23$ in epoch $1683-1764$, while "P" had radius $\gamma\sim
20\degr$ and "S" had radius $\gamma\sim 13\degr$ in epoch $1869-1955$. This
indicates the decrease of magnetic activity from epoch $1683-1764$ to
$1869-1955$.

(IV) The spots also exhibited longitudinal migrations. "P" evolved slowly and
migrated in oscillation around its average longitude with amplitude of about
$15\degr$ and period of about $40$ days, "S" in "T1" showed migration tightly
synchronized with "P", while "S" in "T2" showed much larger oscillation with
amplitude of about $40\degr$ and period of $30-40$ days (figures 11 and 12).

(V) Relative variations of chromospheric activity indicators from 1995 to 2013
revealed both short-term rotational modulation (figure 16) and long-term
decrease of activity (figure 14), impling the existence and evolution of
magnetic activity, and thus the distortion of LC is likely due to starspots on
HD 134319.

## Acknowledgments

This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France (http://cdsweb.u-strasbg.fr/), the data from the European
Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed
by the Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium), the Keck Observatory
Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA
Exoplanet Science Institute (NExScI), under contract with the National
Aeronautics and Space Administration and the ELODIE archive at Observatoire de
Haute-Provence (OHP, http://atlas.obs-hp.fr/elodie/). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. TESS photometric data
presented in this paper were obtained from the Mikulsky Archive for Space
Telescopes (MAST). STScI is operated by the Association of Universities for
Research in Astronomy, Inc., under NASA contract NAS5-26555. This paper
includes data collected by the TESS mission. Funding for the TESS mission is
provided by the NASA’s Science Mission Directorate. We made use of Lightkurve,
a Python package for Kepler and TESS data analysis (Lightkurve Collaboration
et al., 2018) launched in Jupyter Notebook environment
(https://github.com/jupyter/notebook/), built on top of libraries including
NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020) and Matplotlib
(Hunter, 2007) and relative to astropy (Astropy Collaboration et al., 2018),
astroquery (Ginsburg et al., 2019), celerite (Foreman-Mackey et al., 2017) and
tesscut (Brasseur et al., 2019). NumPy, Matplotlib, SciPy and brokenaxes
(https://github.com/jerrylikerice/brokenaxes) were used in preparing figures
of this paper. We also thank Dr. Longcheng Gui${}^{\href
https://orcid.org/0000-0001-7766-4815}$ for his great help in computing
resources. This work is supported by National Natural Science Foundation of
China through grants Nos. 10373023, 10773027, 11333006, U1531121, and
11903074. We acknowledge the science research grant from the China Manned
Space Project with NO. CMS-CSST-2021-B07. The joint research project between
Yunnan Observatories and Hamburg Observatory is funded by Sino-German Center
for Research Promotion (GZ1419).

## Data Availability

The data underlying this article are available in the article. Part of the
data underlying this paper are in the public domain and available in: the Gaia
Archive (https://gea.esac.esa.int/archive/), the Keck Observatory Archive
(KOA; https://koa.ipac.caltech.edu/), the ELODIE Archive (http://atlas.obs-
hp.fr/elodie/) and Mikulski Archive for Space Telescopes (MAST;
https://archive.stsci.edu).

## References

  * Aigrain et al. (2015) Aigrain S., et al., 2015, MNRAS, 450, 3211
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, aj, 156, 123
  * Auvergne et al. (2009) Auvergne M., et al., 2009, A&A, 506, 411
  * Balona & Abedigamba (2016) Balona L. A., Abedigamba O. P., 2016, MNRAS, 461, 497
  * Balona et al. (2019) Balona L. A., et al., 2019, MNRAS, 485, 3457
  * Baranne et al. (1996) Baranne A., et al., 1996, A&AS, 119, 373
  * Barden (1985) Barden S. C., 1985, ApJ, 295, 162
  * Basri & Shah (2020) Basri G., Shah R., 2020, ApJ, 901, 14
  * Batalha et al. (2010) Batalha N. M., et al., 2010, ApJ, 713, L109
  * Berdyugina (2005) Berdyugina S. V., 2005, Living Reviews in Solar Physics, 2, 8
  * Berdyugina & Järvinen (2005) Berdyugina S. V., Järvinen S. P., 2005, Astronomische Nachrichten, 326, 283
  * Berdyugina et al. (2002) Berdyugina S. V., Pelt J., Tuominen I., 2002, A&A, 394, 505
  * Boro Saikia et al. (2018) Boro Saikia S., et al., 2018, A&A, 616, A108
  * Brasseur et al. (2019) Brasseur C. E., Phillip C., Fleming S. W., Mullally S. E., White R. L., 2019, Astrocut: Tools for creating cutouts of TESS images (ascl:1905.007)
  * Budding (1977) Budding E., 1977, Ap&SS, 48, 207
  * Burgasser et al. (2005) Burgasser A. J., Kirkpatrick J. D., Lowrance P. J., 2005, AJ, 129, 2849
  * Butler et al. (2017) Butler R. P., et al., 2017, AJ, 153, 208
  * Cao & Gu (2014) Cao D.-t., Gu S.-h., 2014, AJ, 147, 38
  * Cao et al. (2019) Cao D., et al., 2019, MNRAS, 482, 988
  * Cincunegui et al. (2007) Cincunegui C., Díaz R. F., Mauas P. J. D., 2007, A&A, 469, 309
  * Claret (2018) Claret A., 2018, A&A, 618, A20
  * Collier Cameron & Donati (2002) Collier Cameron A., Donati J. F., 2002, MNRAS, 329, L23
  * Croll et al. (2006) Croll B., et al., 2006, ApJ, 648, 607
  * Davenport et al. (2015) Davenport J. R. A., Hebb L., Hawley S. L., 2015, ApJ, 806, 212
  * Davenport et al. (2020) Davenport J. R. A., Mendoza G. T., Hawley S. L., 2020, AJ, 160, 36
  * Decin et al. (2003) Decin G., Dominik C., Waters L. B. F. M., Waelkens C., 2003, ApJ, 598, 636
  * Donati & Collier Cameron (1997) Donati J. F., Collier Cameron A., 1997, MNRAS, 291, 1
  * Dorren (1987) Dorren J. D., 1987, ApJ, 320, 756
  * Duncan et al. (1991) Duncan D. K., et al., 1991, ApJS, 76, 383
  * Eberhard & Schwarzschild (1913) Eberhard G., Schwarzschild K., 1913, ApJ, 38, 292
  * Flores Soriano & Strassmeier (2017) Flores Soriano M., Strassmeier K. G., 2017, A&A, 597, A101
  * Foreman-Mackey et al. (2017) Foreman-Mackey D., Agol E., Angus R., Ambikasaran S., 2017, ArXiv
  * Fröhlich et al. (2012) Fröhlich H. E., Frasca A., Catanzaro G., Bonanno A., Corsaro E., Molenda-Żakowicz J., Klutsch A., Montes D., 2012, A&A, 543, A146
  * Gaia Collaboration (2018) Gaia Collaboration 2018, VizieR Online Data Catalog, p. I/345
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A1
  * Gebbie & Steinitz (1974) Gebbie K. B., Steinitz R., 1974, in Athay R. G., ed., Vol. 56, Chromospheric Fine Structure. p. 55
  * Ginsburg et al. (2019) Ginsburg A., et al., 2019, AJ, 157, 98
  * Gu et al. (2002) Gu S.-h., Tan H.-s., Shan H.-g., Zhang F.-h., 2002, A&A, 388, 889
  * Hackman et al. (2013) Hackman T., et al., 2013, A&A, 553, A40
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hempelmann et al. (2016) Hempelmann A., Mittag M., Gonzalez-Perez J. N., Schmitt J. H. M. M., Schröder K. P., Rauw G., 2016, A&A, 586, A14
  * Henry et al. (1995) Henry G. W., Eaton J. A., Hamer J., Hall D. S., 1995, ApJS, 97, 513
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Işık et al. (2011) Işık E., Schmitt D., Schüssler M., 2011, A&A, 528, A135
  * Ioannidis & Schmitt (2020) Ioannidis P., Schmitt J. H. M. M., 2020, A&A, 644, A26
  * Isaacson & Fischer (2010) Isaacson H., Fischer D., 2010, ApJ, 725, 875
  * Jeffers & Keller (2009) Jeffers S. V., Keller C. U., 2009, in Stempels E., ed., American Institute of Physics Conference Series Vol. 1094, 15th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun. pp 664–667, doi:10.1063/1.3099201
  * Jeffers et al. (2007) Jeffers S. V., Donati J.-F., Collier Cameron A., 2007, MNRAS, 375, 567
  * Jenkins et al. (2016) Jenkins J. M., et al., 2016, in Chiozzi G., Guzman J. C., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9913, Software and Cyberinfrastructure for Astronomy IV. p. 99133E, doi:10.1117/12.2233418
  * Kervella et al. (2019) Kervella P., Arenou F., Mignard F., Thévenin F., 2019, A&A, 623, A72
  * Kitchatinov & Olemskoy (2012) Kitchatinov L. L., Olemskoy S. V., 2012, MNRAS, 423, 3344
  * Koch et al. (2010) Koch D. G., et al., 2010, ApJ, 713, L79
  * Korhonen et al. (2001) Korhonen H., et al., 2001, A&A, 374, 1049
  * Küker & Rüdiger (2008) Küker M., Rüdiger G., 2008, in Journal of Physics Conference Series. p. 012029, doi:10.1088/1742-6596/118/1/012029
  * Lanza et al. (2006) Lanza A. F., Piluso N., Rodonò M., Messina S., Cutispoto G., 2006, A&A, 455, 595
  * Lanza et al. (2016) Lanza A. F., Flaccomio E., Messina S., Micela G., Pagano I., Leto G., 2016, A&A, 592, A140
  * Lanza et al. (2019) Lanza A. F., Netto Y., Bonomo A. S., Parviainen H., Valio A., Aigrain S., 2019, A&A, 626, A38
  * Lehtinen et al. (2012) Lehtinen J., Jetsu L., Hackman T., Kajatkari P., Henry G. W., 2012, A&A, 542, A38
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration et al., 2018, Lightkurve: Kepler and TESS time series analysis in Python (ascl:1812.013)
  * Linsky (2017) Linsky J. L., 2017, ARA&A, 55, 159
  * López-Santiago et al. (2010) López-Santiago J., Montes D., Gálvez-Ortiz M. C., Crespo-Chacón I., Martínez-Arnáiz R. M., Fernández-Figueroa M. J., de Castro E., Cornide M., 2010, A&A, 514, A97
  * Lowrance et al. (2003) Lowrance P. J., Nicmos Environments Of Nearby Stars Team STIS 8176 Team 2003, in Martín E., ed., Vol. 211, Brown Dwarfs. p. 295 (arXiv:astro-ph/0208072)
  * Martin et al. (2017) Martin E. C., et al., 2017, ApJ, 838, 73
  * Martínez-Arnáiz et al. (2011) Martínez-Arnáiz R., López-Santiago J., Crespo-Chacón I., Montes D., 2011, MNRAS, 414, 2629
  * McCarthy et al. (2001) McCarthy C., Zuckerman B., Becklin E. E., 2001, AJ, 121, 3259
  * Messina & Guinan (1998) Messina S., Guinan E. F., 1998, Information Bulletin on Variable Stars, 4576, 1
  * Messina & Guinan (2003) Messina S., Guinan E. F., 2003, A&A, 409, 1017
  * Messina et al. (1998) Messina S., Guinan E. F., Lanza A. F., 1998, Ap&SS, 260, 493
  * Messina et al. (2001) Messina S., Rodonò M., Guinan E. F., 2001, A&A, 366, 215
  * Meyer et al. (2006) Meyer M. R., et al., 2006, PASP, 118, 1690
  * Mittag, M. et al. (2018) Mittag, M. Schmitt, J. H. M. M. Schröder, K.-P. 2018, A&A, 618, A48
  * Mittag et al. (2019) Mittag M., Schmitt J. H. M. M., Hempelmann A., Schröder K. P., 2019, A&A, 621, A136
  * Montes et al. (1995) Montes D., Fernandez-Figueroa M. J., de Castro E., Cornide M., 1995, A&A, 294, 165
  * Montes et al. (2001) Montes D., López-Santiago J., Gálvez M. C., Fernández-Figueroa M. J., De Castro E., Cornide M., 2001, MNRAS, 328, 45
  * Morris et al. (2018) Morris B. M., Curtis J. L., Douglas S. T., Hawley S. L., Agüeros M. A., Bobra M. G., Agol E., 2018, The Astronomical Journal, 156, 203
  * Morris et al. (2019) Morris B. M., Curtis J. L., Sakari C., Hawley S. L., Agol E., 2019, The Astronomical Journal, 158, 101
  * Moultaka et al. (2004) Moultaka J., Ilovaisky S. A., Prugniel P., Soubiran C., 2004, PASP, 116, 693
  * Mugrauer et al. (2004) Mugrauer M., et al., 2004, A&A, 417, 1031
  * Namekata et al. (2020) Namekata K., et al., 2020, ApJ, 891, 103
  * Oláh et al. (2006) Oláh K., Korhonen H., Kővári Z., Forgács-Dajka E., Strassmeier K. G., 2006, A&A, 452, 303
  * Özavcı et al. (2018) Özavcı I., Şenavcı H. V., Isık E., Hussain G. A. J., O’Neal D., Yılmaz M., Selam S. O., 2018, MNRAS, 474, 5534
  * Petit et al. (2002) Petit P., Donati J.-F., Collier Cameron A., 2002, MNRAS, 334, 374
  * Reinhold et al. (2013) Reinhold T., Reiners A., Basri G., 2013, A&A, 560, A4
  * Ricker et al. (2015) Ricker G. R., et al., 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Scandariato et al. (2017) Scandariato G., et al., 2017, A&A, 598, A28
  * Scargle (1982) Scargle J. D., 1982, ApJ, 263, 835
  * Schmitt et al. (2019) Schmitt J. H. M. M., Ioannidis P., Robrade J., Czesla S., Schneider P. C., 2019, A&A, 628, A79
  * Shkolnik et al. (2005) Shkolnik E., Walker G. A. H., Bohlender D. A., Gu P.-G., Kürster M., 2005, ApJ, 622, 1075
  * Soderblom (1985) Soderblom D. R., 1985, AJ, 90, 2103
  * Soubiran et al. (2018) Soubiran C., et al., 2018, A&A, 616, A7
  * Stassun et al. (2018) Stassun K. G., et al., 2018, AJ, 156, 102
  * Still & Barclay (2012) Still M., Barclay T., 2012, PyKE: Reduction and analysis of Kepler Simple Aperture Photometry data, Astrophysics Source Code Library (ascl:1208.004)
  * Strassmeier (2009) Strassmeier K. G., 2009, A&ARv, 17, 251
  * Strassmeier et al. (1990) Strassmeier K. G., Fekel F. C., Bopp B. W., Dempsey R. C., Henry G. W., 1990, ApJS, 72, 191
  * Strömgren (1966) Strömgren B., 1966, ARA&A, 4, 433
  * Takeda et al. (2007) Takeda G., Ford E. B., Sills A., Rasio F. A., Fischer D. A., Valenti J. A., 2007, ApJS, 168, 297
  * Valenti & Fischer (2005) Valenti J. A., Fischer D. A., 2005, ApJS, 159, 141
  * Valio et al. (2017) Valio A., Estrela R., Netto Y., Bravo J. P., de Medeiros J. R., 2017, ApJ, 835, 294
  * Vida et al. (2015) Vida K., Korhonen H., Ilyin I. V., Oláh K., Andersen M. I., Hackman T., 2015, A&A, 580, A64
  * Vinícius et al. (2017) Vinícius Z., Barentsen G., Gully-Santiago M., Cody A. M., Hedges C., Still M., Barclay T., 2017, KeplerGO/PyKE, doi:10.5281/zenodo.835583, https://doi.org/10.5281/zenodo.835583
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogt et al. (1994) Vogt S. S., et al., 1994, in Crawford D. L., Craine E. R., eds, Proc. SPIEVol. 2198, Instrumentation in Astronomy VIII. p. 362, doi:10.1117/12.176725
  * Walker et al. (2003) Walker G., et al., 2003, PASP, 115, 1023
  * Walkowicz & Hawley (2009) Walkowicz L. M., Hawley S. L., 2009, AJ, 137, 3297
  * Walkowicz et al. (2013) Walkowicz L. M., Basri G., Valenti J. A., 2013, ApJS, 205, 17
  * White et al. (2007) White R. J., Gabor J. M., Hillenbrand L. A., 2007, AJ, 133, 2524
  * Wilson & Joy (1950) Wilson R. E., Joy A. H., 1950, ApJ, 111, 221
  * Wright et al. (2004) Wright J. T., Marcy G. W., Butler R. P., Vogt S. S., 2004, ApJS, 152, 261
  * Wright et al. (2011) Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011, ApJ, 743, 48
  * Xu et al. (2021) Xu F., Gu S., Ioannidis P., 2021, MNRAS, 501, 1878
  * Zechmeister & Kürster (2009) Zechmeister M., Kürster M., 2009, A&A, 496, 577