\section{Introduction}

GRO\,J2058+40 was discovered as a 198 s X-ray pulsar by BATSE on board the
Compton Gamma-Ray Observatory ({\it CGRO}) during a giant outburst in 1995
September-October \citep{wilson98}. This giant (type II) outburst was followed
by a series of weaker (type I) outbursts, whose intensity peaked at intervals of
about 55 days \citep{wilson98,wilson05}. Because BATSE detected the odd
outbursts, counting from the giant outburst, to be brighter than even outbursts,
it was proposed that \gro\ was undergoing periastron and apastron outburst in a
110 day orbit \citep{wilson98}. However, the all-sky monitor (ASM)
on board the {\it Rossi X-Ray Timing Explorer (RXTE)} did not see such a pattern.
Consequently, the orbital period was established to be 55 days \citep{wilson05}.

\gro\ was X-ray active from 1996 to 2002, but only detected at periastron
passages. Pulsations were not detected during the epoch 2003--2004 
\citep{wilson05}.
\gro\ entered a new period of enhanced X-ray activity in May 2008
\citep{krimm08}. Although the X-ray flux doubled in less than a week, the
luminosity was relatively low, $L_x=1.5\times 10^{36}$ \ergs\ assuming a
distance of 9 kpc \citep{wilson05,reig15}. This is more consistent with  the lower-luminosity type I outbursts. However, because the 2008 X-ray data were
limited to just one {\it Swift}/XRT pointed observation, it is likely that 
the peak of the outburst was missed.


Another giant X-ray outburst occurred in March 2019 with a peak luminosity in
the energy range 3--78 keV of $L_x=5.6 \times 10^{37}$ \ergs\ \citep{kabiraj20},
similar to the discovery luminosity back in 1995.  {\it NuSTAR} and {\it
AstroSat} observations performed during this outburst revealed the possible
presence of a cyclotron line at 10 keV together with some harmonics
\citep{molkov19,mukerjee20}. During {\it AstroSat} observations, a
quasi-periodic oscillation at 0.090 Hz was detected during the decay of the
outburst \citep{mukerjee20}. The {\it NuSTAR} observations only detected the
cyclotron line in a very narrow range of the spin phases of \gro\ and covering
only $\sim10$\% of the entire spin period \citep{molkov19}, but it did not seem
to be present in the pulse-average spectrum \citep{kabiraj20}.

The optical counterpart to GRO\,J2058+40 is an O9.5-B0e IV-V star
\citep{reig05a,wilson05}. The massive companion is a long-term variable system
with photometric and spectroscopic changes on timescales of years
\citep{kiziloglu07b,reig15,reig16}. The underlying B star also displays fast
optical multiperiodic variability that is  attributed to nonradial pulsations
\citep{kiziloglu07b,reig22}. We present the most complete and
detailed study of the optical variability of the BeXB \gro\ performed so far. 


\section{Observations and data analysis}

We obtained optical photometry, spectroscopy, and polarimetry of the optical
counterpart to \gro. The observations were made from the 1.3 m telescope at the
Skinakas observatory\footnote{https://skinakas.physics.uoc.gr/en/} (SKO) in
Crete (Greece). Additional spectra were retrieved from the data archive of the 
2.0 m Liverpool Telescope\footnote{https://telescope.livjm.ac.uk/}(LT).


\subsection{Photometry}

Photometric observations with the Johnson-Cousins $B$,
$V$, $R,$ and $I$ filters were made with the 1.3 m telescope of the
Skinakas Observatory. The telescope
was equipped with a  2048$\times$2048 ANDOR CCD with a 13.5 $\mu$m pixel
size. In this configuration, the  plate scale is 0.28$\arcsec$/pixel,
hence providing a field of view of $9.5 \times 9.5$ arcmin$^2$. Standard
stars from the Landolt list \citep{landolt09} were used for the
transformation equations.  The data were reduced in the
standard way using the IRAF tools for aperture photometry\footnote{A User's Guide to CCD
Reductions with IRAF, Philip Massey, February 1997. https://iraf.net/irafdocs} 
\citep{tody86}. After the
standardization process, we finally assigned an error to the calibrated
magnitudes of the target given by the rms of the residuals between the
cataloged and calculated magnitudes of the standard stars. The photometric
magnitudes are given in Table~\ref{photres} (see Appendix~\ref{app}).


\subsection{Polarimetry}

Polarimetric observations were made with the RoboPol photopolarimeter attached
to the focus of the 1.3 m telescope of the Skinakas Observatory.  In the
polarimetry configuration a plate scale of  0.43 $\arcsec$/pixel is achieved
with a 2048$\times$2048 ANDOR CCD with a 13.5 $\mu$m pixel size.  RoboPol is an
imaging photopolarimeter that simultaneously measures the Stokes parameters of
linear polarization  of all sources in the field of view
\citep{king14,ramaprakash19}. RoboPol splits the incident light in two beams,
each half incident on a half-wave retarder followed by a Wollaston prism.  The
fast axis of the half-wave retarder in front of the first prism is rotated by
$67.5^{\circ}$ with respect to the other retarder. Every point in the sky is
thereby projected to four points on the CCD with different polarization state. 
The photon counts in each spot, measured using aperture photometry, were used to
calculate the $U$ and $Q$ parameters of linear polarization.  To optimize the
instrument sensitivity, a mask was placed in the telescope focal plane. The
absence of moving parts allows RoboPol to compute all the Stokes parameters of
linear polarization in one shot. The intrumental polarization and polarization
angle zeropoint were controlled with regular measurements of polarimetric
standards as described in \citet{blinov21}. All observations were corrected for the
instrumental polarization, and uncertainties were propagated accordingly. The
results of the polarimetric observations are given in Table~\ref{polres}
(Appendix~\ref{app}).

\subsection{Spectroscopy}

The 1.3\,m telescope at SKO  was equipped with a
2048$\times$2048 (13.5 micron) pixel ANDOR IKON CCD and a 1302 l~mm$^{-1}$
grating, giving a nominal dispersion of $\sim$0.9 \AA/pixel.  Spectra of
comparison lamps were taken before and after each exposure in order to perform
the wavelength calibration and account for small variations produced by the
tension and flexture of the telescope at different zenith angles during the
night that could affect the wavelength calibration The 2.0 m LT is a fully
robotic telescope at the Roque de Los Muchachos on the Canary Island of La Palma
(Spain). We downloaded and analyzed all the spectra available at the LT data
archive for this source.  These spectra were obtained  with the Fibre-fed
RObotic Dual-beam Optical Spectrograph (FRODOSpec), which is  a dual-beam design
multipurpose integral-field input spectrograph that splits the beam before the
entrance to the individually optimized collimators. We used the
spectra taken with the red arm, which cover the wavelength range 5900--8000 \AA\
with a dispersion of 0.6 \AA/pixel. A Xenon arc exposure was obtained after the
science exposure. The log of the spectroscopic observations is given in  
Table~\ref{specres} (see Appendix~\ref{app}).

The \ha\ line is the prime indicator of the state of the disk. The structure and
size of the disk affect not only the strength of the \ha\ line, but also its
shape. We investigated the evolution of the strength and line
profile over about an 18 year period. We took the equivalent width of the \ha\
line (\ew) as a proxy of the strength of the line. To study the changes
experienced by the shape of the line, we fit the line profile with one or two
Lorentzian functions, depending on whether the line displayed a single-peak or a
double-peak profile. In the case of a split profile, the peak at shorter
wavelength is called the blue or violet peak, and it is denoted by $V$, and the
peak at the longer wavelength is called the red peak and is represented by $R$.

The three fitting parameters are the central wavelength ($\lambda_i$), the
full width at half maximum ($FWHM_i$), and the relative intensity of the peak
($I_i$). Here the subindex $i$ refers to either the blue ($V$) or red peak
($R$). When the line showed a double-peak profile, we obtained the
following quantities: $\Delta V=\lambda_{\rm R}-\lambda_{\rm V}$
and $V/R=\log(I_{\rm V}/I_{\rm R})$, where $I_{\rm V}$ and $I_{\rm R}$ are the
intensity of the blue and red peaks relative to the continuum. For single-peak
profiles, $\Delta V=0$ and $\log(V/R)$ was not defined.

The \ew\ was computed by means of a numerical integration method, the
trapezoidal rule, as implemented in the scipy python package. To ensure an
homogeneous processing, all the spectra were normalized with respect to the
local continuum, which  was rectified to unity by employing a polynomial fit.
Because the definition of the continuum is the main source of uncertainty in the
computation of the equivalent width, we averaged over different selections of the
continuum and different polynomial fits (different grade). A total of 24
measurements were obtained for each spectrum. The error in \ew\ is the standard
deviation of these measurements.


For the line fits, we used the {\tt deblend} task of the SPLAT package
\citep{skoda14} of the Starklink project \citep{currie14}. The fits require the
continuum to be at zero level. Therefore, the best polynomial fit of the
continuum was subtracted from the spectra prior to the Lorentzian fits.

\begin{figure} 
\resizebox{\hsize}{!}{\includegraphics{./fig1.eps}} 
\caption[]{Evolution of the optical observables. From top to bottom, we show the \ha\
equivalent width, continuum emission in the $R$ band, $B-V$ color index, and
polarization degree in the $R$ band (not corrected for ISM polarization).}
 \label{opt}
\end{figure}


\section{Results}





\subsection{Brightness variability}
\label{res}

Figure~\ref{opt} shows the long-term optical variability of \gro. In general, 
the correlation among all optical observables is good. As the strength of
the line emission increases, the continuum emission and the polarization 
increase.  
\gro\ went through two high and three low optical states during the course of
our observations. 

A closer look at Fig.~\ref{opt} reveals some interesting results. First, the
duration of the line emission in the low-optical states is significantly shorter
than the continuum emission, that is, the changes in \ha\ occur on faster timescales than in the continuum. While the continuum emission may remain stable for
years during the optical minimum, the \ha\ line emission reaches a minimum and
begins to increase after just a few weeks. Second, the continuum emission reaches
its peak by about $\sim400$ days before the \ew\ does. Third, the source
exhibited stronger line and continuum emission  during the 2019 high state than
during the 2009 high state: \ew$=-13.6$ \AA\ and $R=13.74$ mag compared to
\ew$=-10$ \AA\ and $R=13.92$ mag. 


The transition from low to high states and then back to low states occurred
smoothly over a period of several years. Overall, the rate of change in the \ew\ is
roughly similar during the two transitions from the low to the high states as
it is for the two declining phases. However, the changes in the optical
observables during the declining phases are in general
more abrupt that during the rising phases. To quantify this result, we performed a
linear fit to the \ew\ for the four different phases, as shown in
Fig.~\ref{ewfit}. The results are summarized in Table~\ref{linfit}. The slopes of
the fits are 1.5--2 times larger during the decay phase.


\begin{figure} 
\resizebox{\hsize}{!}{\includegraphics{./fig2.eps}} 
\caption[]{Evolution of the \ha\ equivalent width. The circles represent
the data, and the solid lines are the best linear fit to the data points of the
corresponding interval. The red arrows mark the occurrence of an X-ray
outburst. The dash style simply denotes a weaker outburst.}
 \label{ewfit}
\end{figure}


\begin{table*}
\caption{Results of the linear fits to the different parts of the long-term
evolution of the \ew. See Fig.~\ref{ewfit} to identify the various regions.}
\label{linfit}
\center
\begin{tabular}{lccccc}
\hline
\hline
Region  &MJD range      &slope                  &Ordinate       &Corr. Coef.  &$N_{\rm obs}$  \\
\hline
\hline
Rise 1 (red)    &53544--55103   &$0.00598\pm0.00006$    &$-320\pm3$     &0.99   &30   \\
Rise 2 (blue)   &56220--58751   &$0.00410\pm0.00004$    &$-230\pm2$     &0.95   &67   \\
Decay 2 (green) &55103--56220   &$-0.00870\pm0.00011$   &$491\pm6$      &-0.92  &33   \\
Decay 3 (brown) &58751--59894   &$-0.00938\pm0.00010$   &$560\pm5$      &-0.98  &31   \\
\hline
\end{tabular}
\end{table*}
\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig3.eps} 
\caption[]{Representative profiles of the \ha\ line. }
\label{lineprof}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\includegraphics[width=10cm]{./fig4.eps} 
\caption[]{Evolution of the \ha\ line parameters. {\em Top}: Separation between the blue and 
red peaks (left axis) and the \ha\ EW (right axis). 
{\em Bottom}: Evolution of the $V/R$ ratio.  }
\label{linevol}
\end{center}
\end{figure}


\subsection{\ha\ line profile variability}
\label{prof}

In contrast to what it is typically seen in other BeXBs, the shape of the \ha\ line
is rather symmetric.  While the strength of the line changes substantially over
time, as described in the previous section, the Be star in \gro\ always displays
a quite symmetric double-peak profile (Fig.~\ref{lineprof}). The ratio of the
peak intensities of the blue and red peaks is almost always $V/R\sim1$
(Fig.~\ref{linevol}, bottom panel), which it is rather unusual as most BeXBs
exhibit $V/R$ variability \citep{reig00,reig10b,yan16,alfonso17,zamanov20}. 

The
separation between the blue and red peaks $\Delta V$ also evolves with time
 and is anticorrelated with \ew\ (Fig.~\ref{linevol}, top panel). At the optical
high states, the peak separation approaches the spectral resolution of our data
($\sim200$ km s$^{-1}$), and  the two peaks merge into a single-peak profile in
some spectra. 



\subsection{Optical polarimetry}

The continuum polarization in Be stars is attributed to Thomson scattering of
unpolarized starlight in the disk
\citep{poeckert79,waters92,wood96,yudin01,halonen13a}.   In general, it is
difficult to compute the actual intrinsic polarization because most of the
observed polarization is expected to have an interstellar origin. For example,
using the relation between polarization and extinction $P_{\rm max,ISM} (\%) =
3.5 \times E(B-V)^{0.8}$ \citep{fosalba02}, we estimate that the maximum
contribution of the interstellar medium (ISM) to the measured optical
polarization toward \gro\ is $P_{\rm ISM}= 4.6$\%,  which is the observed
polarization, on average. However, the fact that  the polarization degree varies
by more than 1\% (from $\sim4$\% to $\sim5.2$\%) and that this variation
correlates with \ew\ (it increased as the \ew\ increased) supports the idea that
part of the observed polarization is intrinsic to the source and originates in
the disk. To estimate the contribution of the ISM, we observed three nearby
stars in the field of view of \gro. Table~\ref{polfsinf} lists the field stars
that we used, while Table~\ref{polfs} gives the polarimetric data for these
stars.

Figure~\ref{extdist} shows the change in extinction with distance and the
location of the field stars and the source. The data were taken from the  dust
maps by \citet{green19}. The figure indicates that at least three different
molecular clouds are in between the observer and the BeXB along the line of
sight. Clearly, as the distance increases, so does the polarization. The field
star fs2 is too close to Earth to represent a reliable measurement of the ISM
polarization at the source distance.  On the other hand, gfs2 displays  stronger
polarization than the source, which might indicate that the star is
intrinsically polarized and hence again only a poor representation of the ISM
polarization. We note that at very short and long distances, the computation of
the extinction is very uncertain due to the lack of sufficient number of stars.
If we use fs1 to estimate the polarization from the ISM (by vector subtraction,
the $q$ and $u$ Stokes parameters of the source and the field star), then we
obtain that the intrinsic polarization varied from 0.3\% to 1.3\%. The lowest
value of the polarization degree corresponds to the 2013 and 2022 observations,
when the \ew\ displayed minimum values.  The level of intrinsic polarization
that we measure agrees with the expected  polarization in an axisymmetric
circumstellar disk predicted by single-scattering plus attenuation models, which
is typically $\simless 2$\% \citep{waters92}. It also agrees with the measured polarization
in a sample of classical Be stars \citep{yudin01}.


\begin{table}
\caption{Selected field stars in the field of view of \gro\ used to derived the
ISM polarization. }
\label{polfsinf}
\begin{center}
\begin{tabular}{@~l@~c@~ccc}
\noalign{\smallskip}    \hline\noalign{\smallskip}
ID      &RA             &DEC           &Distance        &{\it Gaia} ID  \\
        &               &              &(pc)            &\\
\noalign{\smallskip}    \hline\noalign{\smallskip}
fs1     &314.72292      &41.77444       &$1894^{+67}_{-48}$     &2065652877361883648  \\
fs2     &314.69167      &41.76639       &$642^{+7}_{-6}$        &2065653598916381568  \\
gfs2    &314.77043      &41.76491       &$5787^{+1771}_{-1095}$ &2065652671203453696  \\
\noalign{\smallskip}    \hline\noalign{\smallskip}
\end{tabular}
\end{center}
\tablefoot{Distance from \citet{bailer-jones21}}
\end{table}

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig5.eps} 
\caption[]{Extinction as a function of distance (data from \citet{green19}).
The dots mark the assumed location of the stars on the extinction curve,
according to their distance. Distances are from \citet{bailer-jones21}. We also
indicate the polarization degree of each star (see Table~\ref{polfs}).}
\label{extdist}
\end{center}
\end{figure}

\subsection{Reddening and distance}

Owing to the presence of the circumstellar disk, the photometric magnitudes and
colors are affected by disk emission.  The overall effect is to make the star
appear redder, that is, to be seen under higher extinction than in the absence of
the disk. As a result, the estimate of the distance obtained from photometric
magnitudes that are thought to be contaminated by disk emission represents a lower
limit.  Some attempts have been made to quantify the contribution of the disk
to the photometric indices \citep{dachs88,fabregat90,riquelme12}. Here we take
advantage of the low optical states in \gro\ (Fig.~\ref{opt}). Because the
source of long-term optical variability is the disk, low optical states
correspond to instances in which the disk is almost absent. Therefore,  the low
optical states give us the opportunity to observe the underlying B star without
a major contribution from the disk. 

To estimate the distance through the distance-modulus relation, $V-M_V-A_V=5
\log(d)-5$, the amount of interstellar extinction $A_{\rm V}=R \times E(B-V)$ to
the source has to be determined.  The most direct method for estimating $E(B-V)$ is
to use the calibrated color of the star according to the spectral type. The
color excess is defined as $E(B-V)= (B-V)_{\rm obs}-(B-V)_0$, where $(B-V)_{\rm
obs}$ is the observed color and $(B-V)_0$ is the intrinsic color of the star.
The expected color for an O9.5-B0V star is $(B-V)_0=-0.29\pm0.02$. This value is
the average of the calibrations from
\citet{johnson66}, \citet{fitzgerald70}, \citet{gutierrez79}, \citet{wegner94},
and \citet{pecaut13}. The error is
the standard deviation of the five values.   We applied this method to the three
low optical states independently and obtained $8.6\pm0.9$ kpc, $7.8\pm0.8$ kpc,
and $8.0\pm0.9$ kpc for the 2005--2006, 2012--2013, and 2022--2023 low states.
We note, however, that the average \ha\ EW was lowest during
the 2005--2006 period at \ew=--0.9 \AA, compared to \ew=--3.3 \AA\ and \ew=--1.3 \AA\
in the other two low states, respectively. Hence we take $8.5\pm0.8$ kpc as the
final value for the distance. This value compares to $8.7^{+1.1}_{-0.9}$ kpc
from {\it Gaia} \citep{bailer-jones21}. The average color excess is estimated to be
$E(B-V)=1.42\pm0.02$.



\section{Discussion}

BeXBs display variability on all timescales and at all wavelengths. The fastest
variability is detected in the X-rays (on the order of seconds) and corresponds
to the rotation of the neutron star, which manifests as pulsations. In the
optical band, the fastest timescales (on the order of hours or days) are
attributed to changes in the stellar photosphere (pulsation and rotation).
Long-term variations (on the order of months to years) are related to structural
changes in the disk. These long-term variations are the subject of the present
work.

\subsection{Disk formation and dissipation}

Be star disks are known to form and dissipate on timescales of years. The
disk emission increases with wavelength and becomes particularly significant in
the infrared band.
\gro\ went through two high (2008-2009 and 2018-2019) and three low optical
states (2004-2005, 2012-2013, and 2022-2023) during the course of our
observations. The transition from  low to high states and then back to low
states occurred smoothly over a period of several years. We attribute these
long-term changes to the formation and dissipation of the circumstellar disk. 
As expected in the disks of Be stars, the amplitude of the variability in the
continuum increases as the effective wavelength increases: the difference in
magnitudes from maximum to minimum is $\Delta B=0.32$, $\Delta V=0.46$, $\Delta
R=0.56$, and $\Delta I=0.70$.


\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig6.eps} 
\caption[]{Polarization angle and polarization degree of \gro\ for the period
2013-2022. The red arrow marks the occurrence of a type II outburst. Filled brown circles denote all the observations, and the open blue circles
represent the weighted mean calculated using the observations over one year. }
\label{pol}
\end{center}
\end{figure}


\subsubsection{Line and continuum emission}

In Sect.~\ref{res} we pointed out that the line emission variability timescales are faster than the continuum emission. The $BVRI$ magnitudes can remain
at a mean value with very small fluctuations during the low optical states. For
example, in the period MJD52800--54000, the mean and standard deviation of 14
measurements of the $V$ magnitude were 14.94 and 0.04 mag, respectively. For MJD
55800--56900 (six measurements), these values were 14.90 and 0.02 mag,
respectively. In contrast, \ew\ does not show a stable behavior like this and varies
more strongly. For the same periods, the average and standard deviation in
\ew\ was --1.8 and --0.8 \AA\ and --2.7 and --2.1 \AA, respectively. 

Another interesting result is the fact that the peak in the \ew\ occurred about
400 days after the peak in the photometric magnitudes. Although the
observational gaps cannot be ignored, the excellent spectral coverage  clearly
indicates that the delay between the \ew\  and $R$-band maxima is real. These
delays imply that the continuum and line emission are anticorrelated for a certain
period of time (i.e., the time elapsed between the two peaks). This result has
been reported in many other BeXBs: Swift\,J0243.6+6124 \citep{liu22},
RX\,J0440.9+4431 \citep{yan16}, 1A\,0535+26 \citep{clark99,yan12b}, 4U\,0115+63
\citep{reig07b}, and SAX\,J2103.5+4545 \citep{camero14}. Interestingly, the
X-ray outbursts occur close to the \ew\ maxima.


This delay can be explained as due to the different locations in the disk in which 
the line and continuum emission in Be stars originate. The outer parts of the
disk contribute most to the emission of the \ha\ line, whereas the origin of the
optical continuum resides in the inner parts. According to the viscous decretion
disk model \citep[see][for a review]{rivinius13a}, about 90--95\% of the
$V$-band flux comes from inside $1.8-2.5\,R_*$ \citep{haubois12}, while \ha\
emission reaches this percentage above $\simmore 6-10\, R_*$ \citep{carciofi11}.
Therefore, the fact that  the optical brightness has already started to decrease
by the time the \ew\ reached the maximum indicates that the dissipation of the
disk began from the inner parts of the disk. Disk dissipation occurs when the
mass ejection mechanism ceases. As the angular momentum supply stops, the inner
parts of the disk are reaccreted onto the star. However, viscosity still
transports angular momentum outward. Thus the outer parts of the disk continue
to expand, while the inner parts move inward \citep{haubois12,carciofi12}. Since
the flux in the $V$ (or $R$) band is mainly produced in the inner parts of the
disk, we expect the visual flux to start decreasing before the \ha\ line flux,
which mainly originates in the outer parts. 

The long-term trend of the polarization degree also supports a disk
formation and dissipation cycle. Polarization in Be stars results from the
scattering of stellar radiation by electrons in the circumstellar disk.
The degree of the polarization provides information about the number of
scatterers. Therefore, we should expect a significant decrease in the
polarization degree during low optical states, as is indeed observed
(Fig~\ref{opt}). Since these low optical states represent instances when the
disk almost vanished, we can assume that the polarization degree roughly
corresponds with that of the ISM. With this assumption, we estimate that the ISM
polarization from the observer to \gro\ is $\sim4$\%.

Strictly speaking, none of low optical states corresponded to the complete loss
of the disk because the \ha\ line did not revert fully into absorption (see
Fig.~\ref{lineprof}).  A normal (i.e., non-\ha\ emitting) B0V star is expected to
have \ew$=+3.5$ \AA\  \citep{jaschek87}, while the lowest values of \ew\ during
the low optical states were around $-[0.3-1]$ \AA. Nevertheless, the low optical
states do represent highly debilitated disk phases. 


\subsubsection{Connection with X-rays}

As we mentioned above, the \ew\ during the 2018-2019 maximum reached higher
values than during the 2008-2009 high state. Since \ew\ provides a measure of
the size of the disk, we conclude that the Be star developed a larger disk
during the 2018-2019 event. The accretion of a large amount of matter onto the
neutron star led to the type II outburst observed by the {\it NuSTAR} and Astrosat space
missions. In principle, it might be argued that since the disk was smaller during the
2008-2009 optical peak, the X-ray emission would be lower, as observed.   On the
other hand, given the similarities and regularity of the low and high optical
states (Figs.~\ref{opt} and \ref{ewfit}), the absence of a giant X-ray outburst
prior to the 2008--2009 optical maximum is surprising. The {\it Swift}/BAT hard X-ray
transient monitor detected an increase in 15--50 keV flux by a factor of two in
less than a week. A pointed observation with {\it Swift}/XRT was made on 9 May 2008.
The 0.3--10 keV flux was estimated to be $1.5 \times 10^{-10}$ erg cm $^{-2}$
s$^{-1}$, which at a distance of 9 kpc corresponds to a luminosity of $L_x=1.4
\times 10^{36}$ \ergs\ \citep{krimm08}. Although this luminosity is one order of
magnitude lower than the 2019 event, we cannot be sure that it corresponded to
the peak of the X-ray event. The 2008--2009 is reminiscent of a type I outburst,
and the lack of data can be attributed to the weakness of the detection due to
the large distance to the source.

One possible explanation for the different luminosity of the two events could be
a different geometry in a misaligned and tilted disk. Current models that
explain X-ray outbursts in BeXBs invoke warped, tilted, and misaligned disks.
Because disks in BeXBs are truncated \citep{reig97a,reig16}, the only way for
large amount of material to be transferred to the neutron star is in highly
asymmetric configurations in which a misaligned disk becomes warped and eccentric.
Highly distorted disks result in enhanced mass accretion when the neutron star
moves across the warped part \citep{martin11,martin14a,okazaki13}. Evidence of these warped disks has been reported for 1A\,0535+262 \citep{moritani13} based
on a spectral line analysis and for 4U\,0115 \citep{reig18b} based on variations in the polarization
angle. Although we observe a small change in the polarization angle
during the 2019 X-ray outburst (Fig.~\ref{pol}), it is not very significant.
Likewise, the symmetry of the \ha\ line profile argues against a highly
distorted disk. Therefore, it remains to be seen how a type I outburst can lead
to the almost complete destruction of the circumstellar disk.

We also noted in Sect~\ref{res} that the dissipation of the disk is faster than
its formation. The rate of change of the \ha\ equivalent width is about
1.8 \AA/yr during disk formation and about 3.3 \AA/yr during
disk dissipation (see Fig.~\ref{ewfit} and Table~\ref{linfit}). This result
contrasts with what it is observed in classical Be stars, in which the
timescales for disk growth are shorter than the timescales for disk dissipation
\citep{haubois12}. A crucial difference is the presence of the neutron star in
BeXBs. As the disk expands and reaches periastron distance, a substantial amount
of matter will be accreted onto the neutron star. Therefore, the dissipation of
the disk in BeXBs may occur faster owing to the interaction between the neutron
star and the disk. 

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig7.eps} 
\caption[]{ ($(B-V)-V$ color magnitude diagram. Empty symbols correspond to
the disk formation phases, and filled symbols correspond to the dissipation
phases.}
\label{vbv}
\end{center}
\end{figure}


\subsection{Inclination angle}
\label{ang}


We showed in Sect.~\ref{prof} that the source displays a double-peak profile even
during the high optical state (see Fig~\ref{lineprof}). This result suggests that we
see the circumstellar disk at an intermediate inclination angle. Typically,
single-peak profiles are associated with high inclination angles (pole-on
stars) and shell profiles with high inclination angles (edge-on stars;
\citealt{slettebak79,hanuschik95,rivinius06,silaj14,sigut22}). An example of a low
inclination system is V\,0332+54, whose inclination angle is estimated to be
$i<20^{\circ}$ \citep{negueruela99,zhang05}, and a BeXB that exhibits shell
profiles is IGR\,J21343+4738 \citep{reig14a}.

The evolution of the photometric magnitudes and colors can also be used to
constrain the inclination angle. Figure~\ref{vbv} shows the color-magnitude
diagram of \gro. Empty symbols correspond to the disk formation phases, and
filled symbols correspond to the dissipation phases. The positive correlation
between the color index $(B-V)$ and the $V$ magnitude implies that the system is
viewed at an intermediate inclination angle, where the disk emission to the
optical colors and magnitudes contributes significantly. As the disk grows, the
system brightens ($V$ decreases), while the overall (disk plus star) emission
becomes redder ($(B-V)$ increases).  At very high inclination angles (equator-on
stars), the inner parts of the Be envelope partly block the stellar photosphere,
while the small projected area of the disk on the sky keeps the disk emission to
a minimum. Thus, stars viewed at very high inclination angles ($i > i_{\rm
crit}$) would show an inverse correlation \citep{harmanec83,harmanec00}.   The
value of the critical inclination angle is not known, but a rough estimate based
on available data suggests $i_{\rm crit} \sim 60-70^{\circ}$
\citep{hanuschik96a,sigut13}.
In summary, the correlation between the photometric colors and magnitudes and
the stable double peak profile of the \ha\ line in \gro\ suggest an inclination
angle in the range $30^{\circ}-60^{\circ}$.


\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig8.eps} 
\caption[]{Peak separation as a function of \ew.  }
\label{deltav-ew}
\end{center}
\end{figure}


\subsection{Disk size and kinematics}
\label{velaw}




It was already proposed by \citet{struve31} that disks in Be stars are supported by rotation. The main evidence for this comes from the
correlation between the projected rotational velocity and the width and shape of
the emission line \citep{andrillat83,dachs86b,hanuschik88,hanuschik89,dachs92}.

In Sect~\ref{prof} we showed that the peak
separation $\Delta V$ and \ew\ are anticorrelated. This is the expected behavior in rotation-supported disks, in which the rotational velocity of the gas particles in the disk
varies as the inverse of the radius, $v_{\phi}\propto r^{-j}$  (e.g., for
Keplerian disks $j=1/2$).  In these disks, the peak separation roughly
corresponds to twice the projected rotational velocity of the particles in the
disk $\Delta V\approx 2 v_{\phi} \sin i$. As the radius of the disk increases
(i.e., as \ew\ increases), the rotational velocity decreases, and so does the peak
separation.  The index $j$ can be estimated from the relation between peak
separation and EW \citep{hanuschik88},

\begin{equation}
log(\Delta V/2 \sin i) = a \times \log(EW(H\alpha))+b
,\end{equation}

\noindent where the slope $a$ is related to the $j$ index as $a=j/2$. Figure \ref{deltav-ew} shows
this relation for \gro. When we restrict the fitting range to EWs between $1-8$ \AA, then we find $j=0.58\pm0.04$, which is not far from the
expected value of 0.5 for a Keplerian disk. Lower and higher values of the \ew\ were removed from the fit because below $\log(EW(H\alpha))\sim
0$ the disk is too small and has probably not reached a stable configuration,
while above 8 \AA, the $\Delta V$ saturates, which might be a spectral
resolution effect.

After showing that the particles in the disk of \gro\ follow a Keplerian
law, we can estimate the disk size. The disk velocity adopts the form 
\citep[see. e.g.,][]{huang72} 

\begin{equation}
\label{vphi}
v_\phi=v_0\left(\frac{r}{R_*}\right)^{-1/2}
,\end{equation}

\noindent where $R_*$ is the equatorial stellar radius, and $v_0$ is a initial
value of the disk velocity close to the stellar surface. In the limiting case,
$v_0$ would be the stellar rotational velocity. For simplicity, we
assumed $v_0=v_*$. By inverting eq.~(\ref{vphi}) and since 
$v_{\phi} \approx (\Delta V/2 \sin i),$

\begin{equation}
\frac{R_d}{R_*} = \left(\frac{2 v_* \sin i}{\Delta V}\right)^2 
,\end{equation}

\noindent where $r=R_d$ is the radius of the \ha\ emission region.
Figure~\ref{radew} shows the size of the circumstellar disk as a function of
time and \ew. As mentioned above, the leveling off at $\sim 50\,~{\rm R}_\odot$ marks the
minimum distance between the peaks that our spectrograph is capable of
discerning.

The Be star projected rotational velocity in \gro\ has been estimated by
\citet{kiziloglu07b} to be in the range $240-310$ km s$^{-1}$. Here we assumed $v_*
\sin i = 275$ km s$^{-1}$ , which corresponds to half the maximum peak separation.
We note, however, that \ha\ double-peak separation exceeding $2 v_* \sin i$ has
been measured in Be stars during phases of very weak \ha\ emission
\citep{dachs92}.

The largest disk radius is estimated to be $R_{\rm disk,max}\sim5-6 \, R_*$ or
$60\, ~{\rm R}_\odot$, assuming a stellar radius of an O9.5-B0e IV-V star of $10\, ~{\rm R}_\odot$
\citep{martins05}.  From Kepler's third law, we can estimate the orbital
separation $a\sim 160 \, ~{\rm R}_\odot$ for $P_{\rm orb}=55$ days and typical values of
the mass $M_*=18\, ~{\rm M}_\odot$. If we now assume that in order to produce an X-ray
outburst, the neutron star must interact with the disk, that is, the disk radius
must be similar to the periastron distance ($a(1-e)$), then the eccentricity of
the system should be $e\sim 0.6$.

\begin{figure}
\begin{center}
\includegraphics[width=8cm]{./fig9.eps} 
\caption[]{Disk radius as a function of \ew. }
\label{radew}
\end{center}
\end{figure}


\section{Conclusion}

We have performed the most detailed study of the long-term optical variability
of the BeXB \gro.  The system displays correlated variability on timescales of
years in all the optical observables. The source goes through optical high and
low states. The high optical state is characterized by large \ew\, bright
continuum emission, and significant polarization. We identify this state with a
well-developed decretion disk. The smooth increases and decreases of the optical
emission are interpreted as phases of formation and dissipation of the
circumstellar disk. The optical maxima coincide with the occurrence of X-ray
outbursts, after which all optical indicators decreased, leading to the
dissipation of the disk.  However, a disk-loss episode never occurred.  The
almost permanent double-peak line profile indicates that the systems is seen at
an intermediate inclination angle, whereas the absence of $V/R$ variability is
attributed to a relatively small and stable disk. We did not find evidence for a
warped disk during the 2019 type II X-ray outburst.

 


\begin{acknowledgements}

Skinakas Observatory is run by the University of Crete and the Foundation for
Research and Technology-Hellas.   The Liverpool Telescope is owned and operated
by the Astrophysics Research Institute of Liverpool John Moores University. The
Starlink software is currently supported by the East Asian Observatory. IRAF is
distributed by the National Optical Astronomy Observatories, which is operated
by the Association of Universities for Research in Astronomy, Inc. (AURA) under
cooperative agreement with the National Science Foundation. This research has
made use of the SIMBAD database, operated at CDS, Strasbourg, France and of
NASA's Astrophysics Data System operated by the Smithsonian Astrophysical
Observatory.

\end{acknowledgements}


\bibliographystyle{aa}
