Optical and Near-Infrared Monitoring of Gamma-ray Binaries Hosting Be Stars

Optical and near-infrared observations are compiled for the three gamma-ray binaries hosting Be stars: PSR B1259-63, LSI+61 303, and HESS J0632+057. The emissions from the Be disk are considered to vary according to the changes in its structure, some of which are caused by interactions with the compact object (e.g., tidal forces). Due to the high eccentricity and large orbit of these systems, the interactions -- and, hence the resultant observables -- depend on the orbital phase. To explore such variations, multi-band photometry and linear polarization were monitored for the three considered systems, using two 1.5 m-class telescopes: IRSF at the South African Astronomical Observatory and Kanata at the Higashi-Hiroshima Observatory.


Introduction
Gamma-ray binaries are the group of X-ray binaries emitting strong, high-energy (0.1-100 GeV) and very high-energy (>100 GeV) gamma-rays [1]. About 10 systems have been classified as gamma-ray binaries in the Galaxy and Large Magellanic Cloud. The optical/nearinfrared companion stars of the compact objects are either an O star or a Be star-a B-type star exhibiting emission-line profiles in Balmer lines from its circumstellar disk (Be disk) [2]. Companion massive stars dominate the emissions from the binaries in the optical and near-infrared bands. In a gamma-ray binary, the relativistic particles from the compact object interact with the photon field and the low-density fast outflow (wind), emitting UV radiation from the massive star. In the case that the compact object is the rotation-powered pulsar, when the pulsar wind collides with the stellar winds and/or the Be disk, particle acceleration is thought to occur at the shock region, where non-thermal emissions increase: this is the case considered in a pulsar-wind model. On the other hand, in the system where mass is transferred from the stellar wind and/or the Be disk to the compact object, sufficient mass transfer leads to jet formation, such that particle acceleration occurs (microquasar model). In either model, the stellar wind and/or the Be disk play important roles in the mass/plasma distribution, photon field, and other structures determining binary interactions.
The most fascinating fact about gamma-ray binaries is that despite many discussions (both observational and theoretical), the nature of the compact object (and, hence, the binary) has not been well established in most of the systems, except for those where a pulse has been detected (e.g., PSR B1259−63 [3] and LMC P3 [4]); however, recently, possible pulse detection has been reported for LS 5039 [5] and LSI+61 303 [6]. Another interesting issue is why only a few High-Mass X-ray binaries (HMXBs) have shown such high-energy emission, arXiv:2109.00975v1 [astro-ph.SR] 2 Sep 2021 among hundreds of X-ray binaries. Furthermore, due to the diversity of the orbital period (i.e., binary separation) and the eccentricity, gamma-ray binaries provide a good experimental environment to study the mechanisms of such high-energy emissions. As HMXBs generally have a high eccentric orbit, the separation between the compact object and the massive star changes within each orbital cycle. The varying positions of the stars and the separation change the contact surface between the outflows and visibility from the observer, which causes an orbital modulation in the emissions. For instance, light curves in X-and gammarays of LS 5039 reach their maximum and minimum around conjunctions. In PSR B1259−63, outbursts occur from radio to gamma-rays around the periastron, where the pulsar approaches the Be disk.
To search for variabilities in the optical and the near-infrared bands, we monitored five gamma-ray binaries from 2010 to 2018. The optical and near-infrared emissions of the gamma-ray binaries were dominated by those from the OB stars and the circumstellar matter. Variations in these bands and their correlation with other wavelength bands, therefore, are expected to provide key evidence of binary interactions, as well as stellar activities contributing to the high-energy emissions. We also monitored several HMXBs and compared their observational features with those of the gamma-ray binaries, to investigate whether any features are unique to the gamma-ray binaries.
In this paper, we report our observations of three gamma-ray binaries (PSR B1259−63, LSI+61 303, and HESS J0632+057) hosting a Be star. If the companion is of Be type, its equatorial decretion disk adds complexity to the circumstellar environment. The Be disk consists of a high-density, partially ionized plasma, from which the optical and near-infrared emission lines and the IR excess arise. Be disks in binary systems are considered to be subject to such processes as precession, warping, tidal deformation, and truncation under the influence of the compact object [2]. The disk may affect the pressure balance of the winds according to its density and geometrical parameters; for instance, the pulsar wind may be suppressed when it approaches the high-density region of the Be disk [7]. On the other hand, in a system with an O star, mass transfer and shock originate from the stellar winds. As the stellar wind is symmetric, compared to the Be disk, the orientation between the two objects and the line of sight are more related to orbital variations. The observations of the gamma-ray binaries hosting O stars, and HMXBs will be reported in another paper, which is currently in preparation.
PSR B1259−63 is a 48 msec period rotation-powered radio pulsar [3] in a highly eccentric (e = 0.87) and wide (P orb = 1236.72 days) binary orbit with a Be star, LS 2883 [8]. The pulsed emission becomes undetectable around the periastron, and the pulsar is assumed to pass behind the Be circumstellar disk during the pulse eclipse, crossing the disk twice in each orbit: before and after periastron. The system is the first binary to be detected with TeV gamma-rays. The radio, X-ray, GeV, and TeV flares around the periastron have been observed for multiple orbital cycles [9]. The double-peak structure of the radio, X-ray, and possibly TeV light curves are interpreted to be associated. However, the GeV flare shows a puzzling feature; it starts after the post-periastron passage, and has no counterpart in other multi-wavelengths. Discussion of the GeV emission mechanism has led to the suggestion that the infrared excess from the Be star may work as the inverse Compton target of the high-energy electrons, contributing to the optical depth in the GeV emissions [10]. LSI+61 303 has the smallest orbit among the three systems considered in this work; P orb = 26.496 days and e = 0.537 [11]. The nature of this system has been discussed by many studies, in the frame of both the microquasar model and the pulsar-wind model. The system is well known with both orbital modulation and super-orbital modulation (∼1700 days) from radio to gamma-rays. Zamanov et al. (2014) studied the correlation between Hα lineprofile parameters, V-magnitude, and Fermi-LAT flux. They found anti-correlation between V-magnitude and Fermi-LAT flux, and correlation between Equivalent Width (EW) and Fermi-LAT flux [12]. Line-profile variability in Hα indicates perturbations in the Be disk within one orbit, and a bright, high-speed (300-443 km s −1 ) component in the red side of the double-peaked profile before the apastron [13]. More recently, orbital variability in its optical photometry and polarization has been studied [14]. The polarization results showed a periodicity of 13.244 days-half of the orbital period. In the B-, Vand R-band light curves orbital modulation with the peak shift was observed, suggesting super-orbital modulation.
HESS J0632+057 is the most puzzling system of the three considered systems, due to the uncertainty of its orbital parameters. Since its discovery [15], its orbital period has been refined by many authors, ranging from 313 to 321 days, although the proposed periods are consistent within the error. Based on the radial velocity of the Be star, two sets of orbital parameters have been proposed. Both show very eccentric orbit, but periastron phase are different; e = 0.83 and φ periastron = 0.967 [16], or e = 0.64 and φ periastron = 0.663 [17]. Here φ = 0 is the the starting point of the Swift observation (JD 2454857.5) [15], which is used by most of the previous studies. In either orbit, two outbursts occur in one orbital cycle away from the periastron: A primary outburst at the phase φ∼0.3-0.4, with a strong peak, and a secondary outburst at φ∼0. 6-0.9. The X-ray and gamma-ray flux decrease to a minimum ("dip") between the two outbursts: φ∼0. 5. There have been few reports regarding orbital modulation in the optical and near-infrared bands to date, but polarization variation around the periastron has been reported [18]. Line-profile variabilities have been shown in the Balmer lines, suggesting Be disk perturbation at φ∼0.5-1, with the timescale of ∼60 and ∼150 days [19,20]. 1.10-day periodic variation has been detected through TESS photometry (6000-10,000 Å), caused by the stellar rotation [21].
In Sections 2 and 3, we describe the observations and the results, respectively. We discuss our results in Section 4, and give our concluding remarks in Section 5.

Observations
The objective of this work is to search for variabilities in the optical and near-infrared bands in the gamma-ray binaries, which is related to the existence of the compact object. Such variabilities are often expected to have orbital modulations. The considered targets host Be stars, but their orbital period and the X/gamma-ray activities differ widely. We aim to compare variabilities in the optical and near-infrared bands with those in X-rays and gammarays to see if the variations are related to binary interactions. Searching for commonality and difference between the three targets will provide hints to the relationship between interactions and the orbital parameters.
We monitored the three gamma-ray binaries using two telescopes, in South Africa and Japan, from 2010 and 2018. Among the three systems reported in this work, HESS J0632+057 was visible from both telescopes, but the other systems were visible from either telescope, due to their high declination. The observation log is summarized in Table 1. We also revisited our spectroscopic observations of HESS J0632+057 taken with high-dispersion spectrographs (HIDES and ESPaDOnS) [17,19] from 2013 to 2018, to compare the photometry and linear polarization with line-profile variabilities. We derived EW of Hα and Hβ. Please note that part of EW(Hα) has already been published [19].

IRSF
The near-infrared monitoring observations of PSR B1259−63 and HESS J0632+057 were performed from 2010 to 2018 using the IRSF (InfraRed Survey Facility) 1.4 m telescope located at the Sutherland station of the South African Astronomical Observatory, using a three-channel infrared camera: SIRIUS (Simultaneous InfraRed Imager for Unbiased Survey) [22]. SIRIUS offers J, H, and Ks bands simultaneously, with a field of view of 7 .7 × 7 .7.
For polarimetric observations, the single-beam polarimeter, SIRPOL (SIRius POLarimetry mode) [23] was installed upstream of the camera. The typical accuracy of the polarization degree δ P was about 0.3%, depending on the stability of the sky. We observed PSR B1259−63 around the 2010 and 2014 periastron passages for about a month in each orbital cycle, with some observations in the quiescent phase taken as reference. For HESS J0632+057, we first focused on the phase zero of the Swift observation [15], which had been considered to be near the periastron [16] . After 2015, we attempted to obtain a larger coverage of the orbital phase.
The primary data reduction was carried out using the standard pipeline software for the SIRIUS and SIRPOL. Aperture photometry was followed using the IRAF APPHOT task, with the aperture size adjusted to the FWHM of stars in the respective image. Differential photometry with multiple reference stars in the FoV was calculated, and the light curves were deduced with reference to a quiescent observation. The standard error of the zero point of instrumental magnitude was adopted for calculation of the error for each observation, including the uncertainty due to the instability of the reference stars, as well as statistical errors. For HESS J0632+057, the target was brightest by a magnitude of about 3, compared to the stars in the FoV, which made it difficult to set an appropriate exposure time for the observation of both the target and the reference stars, as well as to reduce the errors of the light curves. To handle this issue, we took defocused images in the latter period of the observations. For the polarimetric observations, the average of each night was calculated. It is noted that the number of the dataset differs night by night, so the size of the error bar (standard deviation) differs accordingly. The data on the night when only one or two sets were taken have no error bar in the plotted figures. It is also noted that polarization was not calculated for the data in 2018.

Kanata
Multiple-band photometric observations of LSI+61 303 and HESS J0632+057 were carried out from 2012 to 2016, using the 1.5 m telescope Kanata equipped with HOWPol and HONIR at Higashi-Hiroshima Observatory in Japan. These systems were monitored along with other HMXBs, including LS 5039, to explore orbital variabilities as well as to conduct follow-up comparisons with X-ray and gamma-ray activities. For this purpose, the targets were observed every few days, without focusing on a certain phase.
HOWPol (Hiroshima One-shot Wide-field Polarimeter) [24] is the optical instrument capable of imaging, polarimetry, and low-resolution spectroscopy. Imaging mode has a 15' diameter of FoV. Among the five filters, VR c I c filters were used to monitor the two systems. HOWPol was also used to monitor R-band linear polarization of these systems.
HONIR (Hiroshima Optical and Near-InfraRed camera) [25] is also capable of imaging, polarimetry and low-resolution spectroscopy, covering the optical and near-infrared bands. The imaging mode has a 10' square of FoV. HONIR has two arms 1 , each with selectable filters. Optical and near-infrared data were captured simultaneously. VRJK filters were used mainly, but the H filter was also used on several nights at the beginning of the monitoring period.
Data reduction of the Kanata data was carried out using IRAF packages for bias/dark subtraction and flat fielding. The same reference stars are used to measure the magnitude of LSI+61 303 using aperture photometry. For HESS J0632+057, there were few proper comparison stars which were close to the target and bright enough in both the optical and near-infrared bands. Therefore, different sets of the comparison stars were used between the two instruments. Similar to IRSF/SIRIUS photometry, differential magnitude from a given time was calculated, to compare the observations (see Appendix A.). For polarization data, two or three sets of images were taken, while rotating the half-wavelength plate (0 • , 22.5 • , 45 • and 67.5 • ) and the average of these sets was calculated. The parameters Q and U, followed by Table 1. Summary of the observations. Column 1 is the target name. Columns 2 and 3 list the used instrument and filters, respectively. Columns 4 and 5 list the observational period (in terms of dates and MJD) and the number of the nights during it.

Results
The observed light curves and polarization showed variabilities or indications of variability. In the following sections, the observed features of the individual systems are described.

PSR B1259−63
In the near-infrared band, PSR B1259−63 showed the ∼0.1-mag brightening around periastron [26]. The light curves in the J-, H-, and Ks-bands were almost identical in the 2010 and 2014 orbital cycles. The brightening started no later than 10 days before the periastron and reached maximum brightness at 12-17 days after periastron. In 2014, the brightness was nearly back to the level in the quiescent phase at 70 days after periastron. A difference in the temporal feature was observed between the light curves of J-, Hand Ks-bands; thus, a characteristic track appeared on the near-infrared color-magnitude diagram. The brightness in the Ks-band increased more rapidly than in the H-band, making the infrared color (H−Ks) redder at first. However, the flux increase of the H-band gradually caught up with that of the Ks-band, which turned blue. The infrared color became redder again as the flare decayed. The time lag between the bands indicated that the variation in the Be circumstellar disk first occurs in an outer region. Figure 1 shows the linear polarization in the JHKs bands, taken with IRSF/SIRPOL around the 2010 periastron. The polarization degree in the near-infrared bands did not change significantly around the periastron, although it may have increased by <1% in the J and H bands at the time of brightening. Please note that the polarization degree observed 198 and 301 days before periastron was the same as that at periastron, within the error margins observed. The average polarization degree in the J, H, and Ks bands were 1.4 ± 0.6%, 1.5 ± 0.7% and 1.3 ± 0.8%, respectively. Johnston et al. (1996) have reported that the linear polarization degree in the radio bands was very high (several tens of %), but depolarization was observed around the periastron, when the pulsar was thought to be behind the Be disk [27].

LSI+61 303
Figures 2 and 3 show the photometric variation of LSI+61 303 in the optical and nearinfrared bands. Orbital modulation is clearly seen in the V and R bands by HOWPol; the brightness increases from around periastron to apastron, then decreases from apastron to periastron. The binned light curve (shown in red) showed the amplitude of the variation is ∼0.05 mag, which is consistent with the previous study [14]. Although the phase coverage is not sufficient, the same sinusoidal feature is seen in the I band light curve. HONIR data showed larger scattering than HOWPol data. In particular, high scattering in the near-infrared data made it difficult to discuss the associated variations. This may be caused by stellar activities on shorter/larger time scales, but the lack of data made it impossible to reach any conclusion. This sinusoidal variation of HOWPol light curve agreed with previous studies ( [12,14] for instance). Both HOWPol and HONIR light curves suggest a short, ∼0.1-0.2-mag brightening a few days after periastron (φ ∼ 0.3-0.35) in the optical bands.  Figure 4 shows the polarimetric variation in the R and J bands. The polarization degree (∼1%) is very small, comparable to the order of the error of the instruments. No orbital modulation was suggested. Periodic variability has been reported to have a period of 13.244 days-half of the orbital period [14]. The average polarization degree and angle were 1.2 ± 0.2% and 144.47.2 degree in the R band, and 0.6 ± 0.2% and 151.7 ± 7.3 degree in the J band. The result for the R band was consistent with that of previous studies; 1.25 ± 0.04% and 135.4 ± 1.0 degree [14].       Figure 5 displays the H-band magnitude folded with three orbital periods proposed for HESS J0632+057 with which the orbital parameters has been calculated: 321 days [15], 316.8 days [28], and 313 days [17]. As described in Section 1, the orbital parameters of HESS J0632+057 still have large uncertainties. We have selected the orbital periods which represent the range of proposed period to check the effect of its uncertainty. The flux decreases around φ∼0.3-0.4, where the primary X-ray and gamma-ray outbursts occur, though the variation is not sufficiently significant considering the error bars and the scarce data points. Our data do not offer enough sensitivity to determine the orbital period, but, as seen in the figures, the several-day difference in the orbital period led to a difference in the phase distribution of the observed data, at φ of 0.8-1.3 in particular, when the data over several orbital cycles are combined. Hereafter, we use P orb = 313 days which was deduced by the radial velocity of the observations with a wider coverage of the orbit to discuss the orbital variability. Figure 6 shows the J-, H-, and Ks-band photometry data, compared with the Equivalent Width (EW) of Hα and Hβ lines. A decrease of magnitude is also seen in the J-band. The EW of Hα did not show a clear orbital modulation. The absolute value of EW of Hβ, on the other hand, showed a hint of modulation with two minima in one orbital cycle; φ∼0.3 and φ∼0. 9. The phase at its minimum is different from those of the J and H flux (φ∼0.3-0.4). Figure 7 compares the color change (∆J − ∆H) with respect to the variation in the J band (∆J). The color-magnitude diagram suggests that when the flux decreases at φ∼0.3-0.4, the system becomes redder. In contrast to PSR B1259−63 [26], it seems that no color transition was observed while the magnitude was increasing and decreasing.   [17,19]. From the top to bottom, EW(Hα), EW(Hβ), J, H, and Ks magnitude are plotted. The orbital period (313 days) to fold the data, and the phase at periastron (φ periastron = 0.663) are cited from Moritani et al. [17]. See Figure 5 for the description of the gray and striped area.    Figure 11 shows the result of linear polarization in the Rand J-bands observed with Kanata. The averaged value of the polarization degree observed with IRSF/SIRPOL was 2.7 ± 0.5% (J), 2.1 ± 0.8% (H) and 2.3 ± 0.5% (Ks), and that with Kanata/HOWPol+HONIR is 2.6 ± 0.4% (J) and 3.8 ± 0.2% (R). The polarization degree in the J-band was the same between IRSF/SIRPOL in 2011-2011 and Kanata/HONIR in 2014-2016. Polarization in the near-infrared bands was smaller than that in the R band, but the trend appeared similar between them; the polarization degree in both of the Jand R-bands decreased and increased again around MJD 57,400. Over several years of the observation period, no significant change was observed in linear polarization. This agreed with the fact that the Hα profiles did not change their intensity drastically within the observational periods [16,17,19].

HESS J0632+057
The feature of orbital modulation is slightly seen in the R band; the polarization degree increased by 0.5% when the phase changed from 0.5 to 0.8 and then decreased, although the amplitude was comparable to measurement error. The polarization angle did not change clearly. Yudin et al. (2017) have reported multicolor polarimetry in the optical band at different orbital phases [18]. This work shows a result consistent with theirs: the polarization degree and angle in the R-band were 3.83 ± 0.15%, 171.3 ± 3 degree (JD 2,457,369) and 3.74 ± 0.15%, 167.9 ± 3 degree (JD 2,457,369). Additionally, note that intra-night variation in the polarization degree has also been suggested, with the amplitude of ∼0.5% (in the V band) [29].   Figure 9. The same as Figure 3, but for HESS J0632+057 with Kanata/HONIR. See Figure 5 for the shaded and striped area.  R J Figure 11. The same as Figure 4, but for HESS J0632+057 with Kanata/HOWPol (until MJD 57,042, R-band) and Kanata/HONIR (from MJD 57,330, R-band, and from MJD 57,139, J-band). See Figure 5 for the shaded and striped area.

Discussions
Multi-year monitoring of the gamma-ray binaries in the optical and near-infrared bands was considered to cover enough orbital cycles and phases to enable comparison of the orbital variability among these systems. The observed light curves showed a variety of features which differed from system to system, as well as from wavelength band to band. Table 2 summarizes the variations in the three gamma-ray binaries. vs X/gamma-rays Outbursts at periastron 1 [26] Anti-correlation with the primary outburst?
In the gamma-ray binaries, the fluxes in the optical and near-infrared bands were dominated by that of the massive star. Be disks exhibit emission-line components in Hydrogen (H I), Helium (He I), and other metallic lines (e.g., Fe II, Si II) in their spectra, whose profiles represent a Keplerian rotation disk. Free-free and free-bound emissions from the Be disk can lead to excess in the near-infrared continuum. Disk absorption affects the observed spectral energy distribution depending on viewing angle of the system [2].
The Be disk-isolated or in a binary system-changes its structure over various time scales, causing variations in the observed magnitude, color, spectra, and polarization. When the equatorial slow outflow from the photosphere stops, the Be disk dissipates and the emission-line components, as well as the IR excess, disappear. Such variability is on a timescale from years to decades, much longer than the orbital period of the gamma-ray binaries, and likely occurs independently of the interaction with the compact object; and, hence, independently of the orbital phase. In addition to the variability due to stellar activities, the existence of the compact object gives extra sources of perturbation to the Be disk in a binary system, through the tidal force and/or the plasma wind. Due to the large, eccentric orbit, such a perturbation is expected to occur at a dedicated phase, thus resulting in orbital modulation.
In spectroscopy, line-profile variability indicates the structural changes of the Be disk such as density perturbation, precession, and warping. In several Be stars 'S-shaped' variation was observed; an emission component has moved blue-ward and red-ward in the line profile. Moreover, many Be stars have shown 'V/R oscillation'; relative intensity of the blue peak to the red peak of double-peaked profile has oscillated. These variabilities are thought to be caused by the one-armed oscillation [2,20,32,33]. Episodic change between the single-and double-peaked profiles in Hα line has been observed in the Be/X-ray binary 4U 0115+63, implying the precession of the Be disk [34]. A triple-peaked profile-indicating a precessing, warped Be disk-has been observed in the Be/X-ray binary A 0535+262 [35].
In photometry, disk growth causes brightening of the magnitude in the optical and near-infrared bands for most Be stars. The disk absorbs and scatters part of the stellar flux and if the disk is seen close to edge-on, its growth may result in a net dimming of the system. Since the radius of the emission region of the individual band differs [2], the color changes as the disk size changes. Indeed, EW(Hα) and color (J − Ks) are correlated, as both observables represent the disk size [36]. Our observation showed that the average brightness in the optical and near-infrared bands did not change over several years. This suggests that the overall size of the Be disk was stable over the observational period.
Linear polarization in the optical and near-infrared bands originates from the scattering of free electrons in the Be disk. The polarization degree is roughly anti-correlated with the opacity, and the polarization angle is thought to be perpendicular to the major axis of the Be disk [2]. Statistical analysis of 672 Be stars has shown that the intrinsic polarization degree is <1.5% [37]. Simulated polarization in the continuum band of the Be disk has predicted that the polarization parameters oscillate when the Be disk has precession and/or an azimuthally asymmetric structure [38]. It has also been shown through simulation that the polarization degree becomes saturated as the Be disk grows [39]. Spectropolarimetric observations of Be stars over decades have shown variation in V-band polarization as the Be disk dissipates and a new disk shows up [40]. Polarimetric observation of the Be/X-ray binary 4U 0115+63, covering a giant outburst, has shown that the polarization degree and angle changed significantly while the Be disk warped [41].

Orbital Modulation
At the periastron passage, PSR B1259−63 showed brightening of the near-infrared bands, indicating structural change of the Be disk due to tidal force from the pulsar [26]. van Soelen et al. (2016) has reported an increase of EW (Hα) around the periastron, suggesting expansion or elongation of the Be disk [42]. Furthermore, the double-peaked profile of He I line from their observation suggested a hint of oscillation in the peak separation and V/R ratio (i.e., the ratio of the peak intensity of the blue part to that of the red part), indicating excitation and propagation of the density wave in the Be disk by tidal force from the pulsar. These variations occurred at a similar time-about 10 days before periastron-and showed a peak around 15 days after periastron.This implies that the timescale of the tidal force impact on the Be disk of this system is around ∼25 days.
Similar to PSR B1259−63, LSI+61 303 also showed a possible brightening in the V and R bands at the post-periastron passage (Figures 2 and 3), while it was not seen in HESS J0632+057 at either proposed periastron (Figures 6, 8 and 9). However, S-shaped variation in the Balmer lines has been reported in HESS J0632+057 around φ∼0.5-1.0 [19,20]. If the periastron is at φ = 0.663 [17], the line-profile variability suggests excitation of the density wave. If this were the case, brightening may be observed if precise photometric monitoring were carried out for this system. Considering the timescale of the disk reconstruction seen in PSR B1259−63 (∼25 days for P orb = 1236.72 days, corresponding to the phase duration of ∼0.02), monitoring with a high cadence around periastron is needed to further discussion of this possibility.
In contrast to PSR B1259−63, HESS J0632+057 showed no turning point in color change (Figure 7; as the J-band brightness decreases, the J − H color becomes reddened). The tendency of color-magnitude variation may be similar to "the negative correlation" classified by Harmanec [43] for shell stars; the stronger H I line emission, the fainter the star. According to Harmanec [43], the edge-on Be stars should also grow redder as the star dims due to obscuration by a growing dense disk. The Balmer lines of HESS J0632+057 do not exhibit a shell profile [19] but the inclination angle i of 47 • -80 • has been estimated [16]. Our result may suggest a variation in size and/or density of the highly inclined disk to the line of sight. In spectroscopy, LSI+61 303 [12] and HESS J0632+057 (this work) did not show EW change at a dedicated phase as PSR B1259−63 did; EW variation looked rather sinusoidal within one orbital cycle.
Orbital variations in polarization suggests structure changes in the Be disk, or an additional source of polarization or depolarization at a certain phase. The lack of change in the polarization degree within 1 % in PSR B1259−63 despite the variability in photometry and spectroscopy at periastron (Figure 1), suggested that the Be disk is larger than 10 R * [39]. Although our monitoring did not detect polarimetric variation in LSI+61 303, another longterm monitoring study has shown such variation with a period of half of the orbital period [14]. A small change (∼0.5%) in polarization degree in HESS J0632+057 seemed to occur on a timescale on the order of several tens of days to one orbit (Figures 11), thus not being a short event at a dedicated phase (∼10 days), such as the periastron event of PSR B1259−63. These features suggest that the disk structure changes in this timescale; i.e., changing and recovering the structure while the compact object revolves around the orbit. Therefore, the tidal force is rather unlikely to be the origin of such polarimetric variation.

Long-Term Be Disk Activity
As described above, the Be disk itself changes structure due to stellar activities. The timescale of Be activity is longer than the orbital period and is thought to cause super-orbital variation, as seen in LSI+61 303. Our monitoring observation did not show a sign of superorbital modulation for the other systems, although the coverage of orbital cycles and/or accuracy of the measurement are not sufficient to determine this conclusively. Please note that the Hα emission line of HESS J0632+057 strengthened in 2018-after our monitoring period-thus indicating Be disk growth [44].

Comparison with High-Energy Emissions and Orbital Parameters
From the viewpoint of the relationship with X-ray and gamma-ray activities, the nearinfrared brightness in PSR B1259−63 increased, while high-energy emissions increased around periastron. On the contrary, LSI+61 303 showed anti-correlation between the optical brightness and high-energy emissions, while HESS J0632+057 showed anti-correlation at the primary outburst. The color changed monotonously as the brightness increased and decreased in HESS J0632+057 (Figure 7), in contrast to what is expected from reconstruction by the tidal force, as seen in PSR B1259−63 [26]. This feature suggests that the size of the Be disk had simply changed, due to either mass transfer or stripping by the pulsar wind.
As shown in Table 2, the orbit size was quite different among the three considered gamma-ray binaries. The Be disk is larger with a larger orbit, as the Be disk is truncated by the compact object. Here, we estimated the Be disk size in PSR B1259−63 using the relationship between EW and emission region of Hα ( [30] eq. 8), citing EW at the 2010 periastron [45]. We also estimated the Be disk size by the emitting region of Hβ of PSR B1259−63, from the peak separation of the double-peaked profile ( [31] ∼116 km s −1 ) assuming a Keplerian disk. The disk size of LSI+61 303 and HESS J0632+057 were cited from the previous study [30]. The observed variation was also different, which seems to be related to the differences in orbital parameter and nature of the system. The light curve showed a clear sinusoidal pattern in LSI+61 303, the binary with the smallest orbit, where the Be disk fills a large fraction of area of the orbit compared with the other systems. On the other hand, PSR B1259−63, where the Be disk covers only a small area of the orbit, showed variations only at periastron. HESS J0632+057 shows similar variations with both PSR B1259−63 and LSI+61 303. Compared to PSR B1259−63, HESS J0632+057 has ∼8 times smaller orbit, whereas the Be disk is ∼3-4 times smaller. On the other hand, compared to LSI+61 303 HESS J0632+057 has ∼40 times larger orbit, whereas the Be disk is ∼3-4 times larger. These features suggest that HESS J0632+057 is in the middle of PSR B1259−63 and LSI+61 303 from the viewpoint of the size if the Be disk in comparison to the orbit; although the alignment of the Be disk plane with respect to the orbital plane also affects the interactions.
Besides the Be disk size compared to the orbit, other parameters also have impacts on binary interaction, such as the orbit eccentricity, misalignment angle between the Be disk and orbit, and Be disk base density. The impact of the disk truncation is weaker if the Be disk is misaligned to the orbit of the companion star, which enables the disk to become larger.
Smoothed particle hydrodynamics (SPH) simulation of the density wave in the Be disk in binary systems has shown that spiral arms and the radial density distribution depends on the misalignment angle between the Be disk and the orbital plane, as well as the viscosity of the Be disk [46]. In a Be disk with higher misalignment angle and higher viscosity, the spiral arm becomes more tight and radial density distribution becomes flatter, without a clear cut-off.
We monitored several Be/X binaries, which have similar orbit size and eccentricity to the gamma-ray binaries but have accreting pulsars. Line-profile variabilities implying the density wave have also been observed in several Be/X-ray binaries. Comparison of orbital modulation between these systems may provide clues to discuss the relationships between these systems.

Concluding Remarks
Three gamma-ray binaries hosting a Be star were monitored, in terms of their optical and near-infrared brightness and linear polarization, from 2010 to 20108. The long-term light curves indicated the orbital modulations (or their sign) in all the systems, although the variations differed from system to system. All three binaries showed variation due to the tidal force from the compact object, although more precise monitoring with a higher cadence is required to determine more conclusive results. The diversity of the variation is possibly related to the difference in the size and base density of the Be disk, the disk size compared to the orbit, and the alignment of the Be disk plane with respect to the orbital plane. In this case, error in ∆m Kanata should be RMS of each measurement.

1
By design, HONIR has three arms, among which two arms were available during the observational period.