The cross-matches between the B-type star sample and the Gaia EDR3 and 2MASS catalogs were done using TOPCAT2 (*Tool for OPerations on Catalogs And Tables* [54,55]). The search radius was set at 0.2 arcsec to avoid duplicated sources. We also selected stars with accurate photometry. This meant they satisfied the following criteria: (a) the magnitude *G* < 16 mag, (b) mean errors in *Bp* and *Rp* little than 0.01 mag for the Gaia EDR3, and (c) *JHKs* photometry with quality flags ("Qfl" or "ph\_qual") equal "AAA" for the 2MASS. In this way, we guaranteed that the *σmi* contributions at *σQ* (in Equation (4)) were as small as possible.

To calculate the reddening-free *Q* values (Equation (2)) we used the combination of all color indices: (*Bp* − *Rp*), (*Bp* − *J*), (*Bp* − *H*), (*Bp* − *Ks*), (*Rp* − *J*), (*Rp* − *H*), (*Rp* − *Ks*), (*J* − *H*), (*J* − *Ks*), and (*H* − *Ks*), which led to a total of *C*(10, 2) = 10! 2!(10−2)! = 45 *Q* parameters. The color-excess ratios were obtained from the relative absorption coefficients, *rma* = *Ama*/*A*V, for each of the selected five bands, which were estimated using the "extinction" code Fast Interstellar Dust Extinction Laws in Python3 [56]. With this tool, we estimated the *rma* values for different wavelengths and different extinction laws for each of the three available models [57–59], which also permitted the use of different values for *R*V. These values are given in Table 2. This code also calculated the extinction for Fitzpatrick and Massa [60]'s model but only for *R*V = 3.1. Therefore, this last model was not adopted in this work. It is important to stress that the dependence of *Bp* and *Rp* relative absorption coefficients with the *T*eff are still unknown. However, as the filters were narrower than the one of the G band, we would expect a lower effect.

**Table 2.** Values of relative absorption coefficients for the Gaia EDR3 and 2MASS filters

In order to estimate typical *Q* errors for normal B-type stars, we calculated ∆*Q* and *σQ* (Equations (3) and (4), respectively). On the one hand, we estimated |∆*Q*| where the *Q*0 values were calculated using Fitzpatrick's model for *R*V = 2.1, 4.1, and 5.1 while the *Q* values corresponded to *R*V = 3.1. On the other hand, we calculated *σQ* for each reddening law (see Table 2) using *R*V = 3.1 and *σE* = 0. Both |∆*Q*| and *σQ* values were calculated for the 45 combinations of color-excess ratios corresponding to our entire sample of normal B-type stars. These values are shown in Figure 3. The results showed that the mean values were: |∆*Q*| ∼ 0.05 and *σQ* ∼ 0.15. Therefore, any star would show anomalous reddening if its |*Q*| value was larger than 3|∆*Q*| ∼ 0.15 of the expected one for its spectral type. By analogy, we could define a *Q* range of 3 *σQ* ∼ 0.45 that would contain all the normal B-type stars. This criterion could be used to identify peculiar stars since they would be located outside the defined *Q* range.

**Figure 3.** Frequency histogram of error distributions for different reddening laws. The *σQ* (left panel) and |∆*Q*| (right panel) values were calculated for our entire sample of "normal" B-type stars using the 45 reddening-free *Q* parameters. To calculate *σQ* we used *R*V = 3.1 for the different reddening laws, while to calculate |∆*Q*| we used *R*V = 2.1, 4.1, and 5.1. See details in the text.

Figure 4 shows some examples of the behavior of the parameter *Q* with *T*eff and their respective color–color indices. These diagrams were selected considering the color indices frequently used in the literature. To study this behavior, we plotted our B-type star sample, together with the main sequence (MS) taken from the Modern Mean Dwarf Stellar Color and Effective Temperature Sequence4 (much of the content of this table, but not all, was incorporated into the Table 5 of [61]). The *Q* values for both the B star sample (red points) and the MS (solid gray line) were calculated using the model created by Fitzpatrick [59] with *R*V = 3.1. To show the effect of anomalous and normal reddening on color–color diagrams, the MS was reddened considering *A*V = 1 mag, adopting the model done by Fitzpatrick [59] for the different *R*V coefficients listed in Table 2, while in the *Q*-*T*eff diagrams the *Q* values for MS were calculated also using the Fitzpatrick's model with *R*V = 2.1, 4.1, and 5.1 (solid and dashed lines, see the caption to Figure 4). Then, considering that *σQ* = 0.15 (as set out above), we could define the following limits in the *Q*-*T*eff diagrams and CCDs (shown in brown dashed lines):

$$-0.10 \le Q\_{J\text{KR}pI} \le 0.05\tag{6}$$

$$-0.07 \le Q\_{HKBpRp} \le 0.08 \tag{7}$$

$$-0.11 \le Q\_{JHK} \le 0.04\tag{8}$$

that encloses almost all the B-type stars of the sample (around 86%, 93%, and 64% of a total of 163 stars for each of the limits defined above).

**Figure 4.** Behavior of *QJKRp J*, *QHKBpRp* and *QJHK* parameters with *T*eff (left panels) and their corresponding color–color diagrams (right panels). The red points represent the loci of the "normal" B-type star of the sample [45–50]. The main sequence band (solid gray line, [61]) was reddened considering *A*V = 1 mag adopting the model done by Fitzpatrick [59] for different *R*V coefficients: *R*V = 2.1 (pink dotted line), *R*V = 3.1 (solid blue line), *R*V = 4.1 (orange dashed line), *R*V = 5.1 (green dashed line). The effect of the reddening vector, i.e., *A*V = 2 mag, for different ISM laws, is indicated in the top-left part of the plots. The brown dashed lines enclosed the cloud of points within *σQ* < 0.15.

#### *3.2. The Q–Color and Q–Q Diagrams for Normal B-Type Stars*

Generally, the effective temperature of the B-type stars is a priori an unknown quantity and is often derived from spectroscopic data. Therefore, an alternative way to describe the behavior of the *Q* parameters with *T*eff, or the spectral types, is to use the color indices. Thus, the (*Bp* − *Rp*) color index is useful, for example, to classify stars according to their *T*eff. Jordi et al. [43] remarked that (*Bp* − *Rp*) plays the same role as (*V* − *I*C), thus the relationship between (*Bp* − *Rp*) and *T*eff could be almost equivalent to the well-known relation of *T*eff = *f*(*V* − *I*C). On the other hand, Mucciarelli et al. [62] also found that (*Bp* − *Ks*)0 and (*G* − *Ks*)0 were the best choice to derive precise and accurate *T*eff values among late-type stars.

Then, to explore the behavior of *Q* parameters with the color indices, we present some QCD in Figure 5. The *Q* parameters for the sample of B-type stars and the MS were constructed as explained in Section 3.1. In these plots, the MS was also reddened. These diagrams allowed early-type stars (red cloud of points) to be separated from late types on the MS. In particular, in the *QBp JRpK* vs. (*J* − *Ks*) and *QBpRpHK* vs. (*Bp* − *Rp*) diagrams (left panels), the B-type stars were located in a well-delimited region. This group was located at

$$-0.11 < Q\_{Bp/RpK} < 0.20 \qquad \text{and} \qquad -0.20 < (f - \text{Ks}) < 0.14 \tag{9}$$

$$-0.60 < Q\_{BpRpHK} < 0.50 \qquad \text{and} \qquad -0.30 < (Bp - Rp) < 0.40 \tag{10}$$

where the upper limits for the color indices correspond to an F0 V star. Another relevant characteristic of the diagrams of Figures 4 and 5 is that the differences between the MS and anomalous reddening ones seemed to be very small, at least for the selected B-type star group.

Our sample of B-type stars in the *QBpK JH* vs. (*Bp* − *Ks*) diagram (top-right panel of Figure 5) was located in the range −0.17 < *QBpK JH* < 1.0, where the lower limit corresponded to the position of an F0 V-type star. This region could be compared with that defined by Poggio et al. [63] in the (*J* − *H*) vs. (*G*DR2 − *Ks*). These authors made a preliminary selection of OB candidates that satisfied both (*J* − *H*) < 0.14 (*G*DR2 − *Ks*) +0.02 and (*J* − *Ks*) < 0.23 (*G*DR2 − *Ks*). The first condition could be expressed in terms of *Q* as *QJHG*DR2*K* = (*J* − *H*) − 0.14(*G*DR2 − *Ks*) < 0.02. Then, by multiplying both sides by (−0.14) −1 , we obtained *QG*DR2*K JH* = (*G*DR2 − *Ks*) − 0.07(*J* − *H*) > −0.14 which was similar to the region defined for *QBpK JH*.

As was already mentioned in Section 2, the reddening vector in the QCD has a horizontal direction (Figure 2), for this reason, the reddened B-type stars were overplotted to A and F-type MS stars. This effect was more evident in the *QJHK* vs. (*H* − *Ks*) diagram (bottom-right panel of Figure 5).

Finally, Figure 6 shows other examples of QQDs. Here, we find that in the *QJHK* vs. *QBpRpK* (top-left panel) and *QRp JK* vs. *QBpRpHK* (bottom-left panel) diagrams, the loci of the B-type stars were restricted to a small region around the coordinate point (0, 0):

$$-0.18 < Q\_{\rm JHK} < 0.01 \qquad \text{and} \qquad -0.12 < Q\_{\rm BpRpK} < 0.18 \tag{11}$$

$$-0.05 < Q\_{Rp/K} < 0.24 \qquad \text{and} \qquad -0.6 < Q\_{BpRpHK} < 0.7 \tag{12}$$

These relations could be used as possible criteria for selecting B-type stars. In particular, the first diagram is very sensitive to the ISM reddening law, and it could also be helpful for detecting stars with anomalous reddening.

On the other hand, the combination of parameters *QJHK* vs. *QBpRpHK* and *QBpKHK* vs. *QBpRpHK* (right panels) were not convenient at all, since they mixed late- and early-type stars. Notwithstanding, it drew attention to the fact that the B-type stars were located in a very narrow line that was almost perpendicular to the reddening direction.

**Figure 5.** Examples of the behavior of the *Q* parameter with the color indices. The red points represent the loci of the "normal" B-type star of the sample [45–50]. The main sequence band (solid gray line) was reddened considering *A*V = 1 mag for different *R*V coefficients: *R*V = 2.1 (pink dotted line), *R*V = 3.1 (solid blue line), *R*V = 4.1 (orange dashed line), and *R*V = 5.1 (green dashed line). See the text for details.

**Figure 6.** Examples of selected *Q*-*Q* diagrams using different *R*V absorption coefficients. The red points represent the loci of the "normal" B-type stars of the sample [45–50]. The main sequence band was calculated using different *R*V coefficients: *R*V = 2.1 (pink dotted line), *R*V = 3.1 (gray solid line), *R*V = 4.1 (orange dashed line), and *R*V = 5.1 (green dashed line). See the text for details.

#### **4. B-Type Stars with Emission Lines**

In the previous section, using several Gaia and 2MASS color–color, *Q*–color, and *Q*–*Q* diagrams, we delimited regions where normal B-type stars locate. In this section, we analyzed the location in the color–color, *Q*–color, and *Q*–*Q* diagrams of different kinds of emission-line B-type stars with respect to the normal ones. Our goal was to be able to find criteria that would allow us to separate among the diverse groups of stars with emission lines.

To analyze the position of the emission-line B-type stars in the different diagrams, we constructed our sample using stars classified by Lamers et al. [7], Kraus [12], Miroshnichenko [14], Cidale et al. [20], Oksala et al. [22], Zorec et al. [45,46], Aidelman et al. [47,48,49], Cochetti et al. [50], Clark et al. [64] and Guzmán-Díaz et al. [65]. In this way, in addition to the 163 normal B-type stars already studied, our sample was composed of 84 classical Be stars (CBe), 204 HAeBe and HAeB[e] stars, 33 LBVs, and 23 B[e]SGs. We searched for Gaia and 2MASS photometric data, and, finally, we obtained a sample of 156 emission-line B-type stars. It was striking that only 34 HAeBe stars, among the 204 selected stars, had photometric data in the Gaia EDR3 and 2MASS catalogs.

#### *4.1. Color–Color Diagrams for Emission-Line B-Type Stars*

The Gaia and 2MASS color–color diagrams are shown in Figure 7. Again, the solid gray line is the main sequence. The reddening vector (*A*v = 2 mag) is indicated in the upper left. The cloud of red and cyan points are, respectively, the normal B and CBe stars. Blue symbols are for HAeBes (and a few HAeB[e]s), orange triangles for B[e]SGs, and green squares for LBVs. Symbols surrounded by black open circles correspond to stars with accurate photometry. These stars satisfied the stated criteria (a), (b), and (c) in Section 3.1. We distinguished two groups of CBe stars in all the plots. The first group behaved as normal B-type stars do. These stars seemed to have small or no infrared excess. The second group encompassed stars with a large color excess, sharing similar photometric characteristics to some LBV stars. These CBe stars presented interstellar reddening and intrinsic color excesses (because the reddening vector points toward a slightly different direction). The large color excess observed in the last group of stars is discussed in Section 5. Particularly interesting are the (*Rp* − *Ks*) vs. (*Bp* − *Rp*) and (*J* − *Ks*) vs. (*Bp* − *Ks*) diagrams (top-left and top-right panels in Figure 7, respectively). These are very useful for searching for CBe stars with large infrared excesses because the normal B stars distribute along a narrow and well-defined band.

The best studied NIR color–color diagram in the literature is the (*J* − *H*) vs. (*H* − *Ks*) (see Figure 7, middle-left panel). Even when the stars are affected by reddening, this NIR two-color diagram permits separate LBVs from B[e]SG stars [22,23]. However, there is still significant contamination between LBVs and HAeBe stars. This diagram also allows distinguishing, although some mix is observed; the HAeBe stars from both the CBe stars [66] and the B[e]SGs [23]. For illustrative purposes, two straight brown dashed lines ((*J* − *H*) = 1.93(*H* − *Ks*) − 1.11 and (*J* − *H*) = 1.93(*H* − *Ks*) − 0.55) were plotted to indicate the regions mostly occupied by each group of stars: LBVs, HAeBes and B[e]SGs. Similarly, a (*J* − *H*) vs. (*Bp* − *Ks*) diagram (Figure 7, bottom-right panel) shows some mixture between young and evolved emission-line stars.

A relevant result was the behavior of the two-color diagram (*Bp* − *Rp*) vs. (*H* − *Ks*) (Figure 7, bottom-left panel), which separated groups of stars (but without recognizing their evolutionary stages) in regions or bands of different color excesses. This behavior was not caused by the reddening vector, which pointed almost perpendicularly to the horizontal shift shown for the different groups of stars. Moreover, this color–color diagram also separated the LBVs from the MS late-type stars. Similarly, the (*Rp* − *J*) vs. (*J* − *Ks*) diagram discriminated, but to lesser degree, LBVs from the other groups. Although, it seemed to gather groups of stars with different circumstellar reddening (Figure 7, middleright panel).

#### *4.2. Q–Color and Q–Q Diagrams for Emission-Line B-Type Stars*

As shown in Section 4.1, it was very difficult to delimit regions for HAeBe and B[e]SG stars using two-color diagrams. Particularly, in all these diagrams, these two types of stars were, to more or less degree, mixed. The same happened with the population of HAeBe and LBV stars. Part of this degeneracy could be solved using the *Q*-color and *Q*-*Q* diagrams.

Figure 8 presents several examples of QCD. The top-left panel, *QBpRpHK* vs. (*Bp* − *Rp*) located the LBVs over the line cut *QBpRpHK* = −4.5 (brown dashed line). Particularly, the *QJHK* vs. (*J* − *Ks*) diagram (top-right panel) separated the different groups according to the amount of IR emission. The greater the amount of hot dust, the higher the (*H* − *Ks*) color index, which implied a more negative value for *QJHKs* = (*J* − *H*) − *E*(*J*−*H*) *E*(*H*−*Ks*) (*H* − *Ks*). In addition, the dashed brown lines at *QJHK* = −1.11 and *QJHK* = −0.55 represent, respectively, the straight lines plotted on the (*J* − *H*) vs. (*H* − *Ks*) diagram. According to this, there were three groups of HAeBe stars, one of them located between the main sequence and the line cut *QJHK* = −0.55. It shared its properties with the LBV group of stars. A second group was located within a region with −1.11 < *QJHK* < −0.55, while the last group had properties similar to the B[e]SGs.

**Figure 7.** Selected Gaia and 2MASS color–color diagrams for B-type stars with emission lines. The main sequence was taken from Pecaut and Mamajek [61]. It was reddened using the extinction law calculated by Fitzpatrick [59], applying *A*V = 1 mag and different *R*V values. Symbols represent the normal B-type stars (in red), the CBe stars (in cyan), the HAeBes/HAeB[e]s (in blue), the LBVs (in green), and the B[e]SGs (in orange). Black open circles point to stars with accurate photometry (see Section 3.1). The reddening vector *A*V = 2 mag is indicated in the upper-left side of the plots.

The middle-left panel, *QJKHK* vs. (*H* − *Ks*), better separated HAeBes from B[e]SGs. There, the HAeBe stars were more reddened than the B[e]SGs. It is also interesting to note that in both *QJHK* vs. (*J* − *Ks*) and *QJKHK* vs. (*H* − *Ks*) diagrams, it was possible to distinguish a group of five strongly reddened LBV stars as they were shifted horizontally to the right.

The most interesting diagram was *QBp JHK* vs. (*Bp* − *Ks*) (middle-right panel), which enabled identification of LBV, B[e]SG, and HAeBe stars from each other. Almost all LBVs had *QBp JHK* > −7. The B[e]SG stars were also well-separated from HAeBe by the straight brown dashed line (*QBp JHK* = −4.33(*Bp* − *Ks*) − 0.57) drawn in the figure. However, in a *QRp JK* vs. (*J* − *Ks*) diagram (bottom-left panel), both CBe and B[e]SG stars seemed to lie on the same straight line (*QRp JK* = −1.99(*J* − *Ks*) + 0.08). Instead, the LBVs had values of *QRp JK* > −1.7. The *QBpK JK* vs. (*Rp* − *Ks*) diagram showed similar behavior.

**Figure 8.** *Q*-color diagrams of emission-line B-type stars. The main sequence was reddened applying *A*V = 1 mag and different *R*V values (see Section 3 for details). Symbols represent the normal B-type stars (in red), the CBe stars (in cyan), the HAeBes/HAeB[e]s (in blue), the LBVs (in green), and the B[e]SGs (in orange). Black open circles point to stars that have accurate photometry (see Section 3.1). The brown dashed lines delimit regions of stars that share similar properties.

Figure 9 illustrates several examples of QQD. LBVs and B[e]SGs which are clearly apart from each other. In most of these diagrams, the HAeBe stars are also separate from the B[e]SGs. For instance, on a *QJHK* vs. *QBpRpHK* diagram (top-left panel), most of the LBVs lie on the right of the vertical line at *QBpRpHK* = −4.5 and above the horizontal at *QJHK* = −0.55, thereby complementing information given in Figure 8 (top panels). In a diagram *QJHK* vs. *QBpK JK* (middle-left panel), the LBVs are inside a box delimited by −4.8 < *QBpK JK* < −0.25, and −0.55 < *QJHK* < 0.11. In addition, it is also possible to separate HAeBes from B[e]SGs through the line *QJHK* = 0.1 *QBpK JK* − 0.35. Finally, the *QJHK* vs. *QBpRpK* diagram (bottom-left panel) shows similar behavior.

Concerning the panels on the right side, they show a greater mix between the group of stars than the ones on the left side, e.g., in the *QRp JHK* vs. *QBpK JK* diagram (middle-right panel), the LBVs have values of *QRp JHK* > −2.4, and the HAeBes and B[e]SGs can be separated with the line *QRp JHK* = 0.48 *QBpK JK* − 1.4.

**Figure 9.** Example of *Q*–*Q* diagrams for B stars with emission lines. The main sequence was reddened applying *A*V = 1 mag and different *R*V values (see Section 3 for details). Symbols represent the normal B-type stars (in red), the CBe stars (in cyan), the HAeBes/HAeB[e]s (in blue), the LBVs (in green), and the B[e]SGs (in orange). Black open circles point to stars that have accurate photometry (see Section 3.1). The brown dashed lines delimit regions of stars that share similar properties.

#### **5. Discussion**

Traditionally, the positions of the stars in different photometric diagrams, or even in the HR diagram, have been adopted as selection criteria to classify stars according to their intrinsic properties, memberships in clusters or associations, variability types, etc. In this work, we focused on studying various color–color diagrams, taking advantage of the Gaia and 2MASS photometric surveys as powerful analysis tools to select and classify B-type stars with and without emission lines. Notably, one would expect that the NIR radiation excess emitted from gaseous and dust CEs of early-type stars would show different properties according to their chemical composition (gas-to-dust ratio), the absorption and heating mechanisms of particles, dust temperature, and the size of dust grains. As each diagram has its own limitations, mainly because the colors are affected by ISM, it was important to adopt multiple criteria. For this purpose, we also analyzed different reddening-free Q parameters, since the method presents high confidence. For instance, Mohr-Smith et al. [67] found that the extension of the *Q* method to the Sloan filter system had a confidence of 97%. The study was based on follow-up spectroscopy of a sample of OB stars selected in a previous study using the *Q* method [68].

In general, from the two-color diagrams of Figure 7, the HAeBe stars were very difficult to separate from B[e]SGs or LBVs. The best option was the (*J* − *H*) vs. (*H* − *Ks*) diagram that enabled distinguishing HAeBe from B[e]SGs [23], and LVBs from B[e]SGs [22], delimiting particular regions through straight lines (middle-left panel). However, there was still a substantial mixture between the HAeBe stars and the LBVs group (according to [23]). This preliminary classification could be complemented by the information derived from a (*Bp* − *Rp*) vs. (*H* − *Ks*) diagram that permitted, in addition, the identification of hot-HAeBe and B[e]SG stars with strong intrinsic IR emissions and no reddening interstellar excess. These stars were located within a region defined by 0 . (*Bp* − *Rp*) . 1 and 0.5 . (*H* − *Ks*) . 1.25. All this information could also be retrieved in the *QJHK* vs. *QBpRpHK* diagram, where the straight lines, plotted on the *JHKs* two-color diagram, were transformed into horizontal lines at *QJHK* = −0.55 and *QJHK* = −1.11. The box in the (*Bp* − *Rp*) vs. (*H* − *Ks*) two-color diagram resulted in the range −12.4 . *QBpRpHK* . −3.96. Furthermore, in this QQD the three groups of stars, that is the HAeBes, LBVs, and B[e]SGs, were well separated from each other. According to Oksala et al. [22], the B[e]SGs were surrounded by larger amounts of hot dust than the LBVs, whereas the LBVs hosted huge amounts of cool dust at far distances.

We also recognized three kinds of HAeBe stars that presented different amounts of NIR radiation excess. In all CCDs, QCDs, and QQDs, we observed a group of HAeBes located in the same region as LBVs and some CBe stars. Figure 9 (top-left panel) shows that this group of stars had *QJHK* > −0.55. The other two groups had −1.11 < *QJHK* < −0.55, and *QJHK* < −1.11, respectively. The latter corresponded to HAeBe stars with significant IR excess. Our result agreed with that found by Chen et al. [19] in a (*J* − *H*) vs. (*H* − *Ks*) diagram. These authors noted that the majority of HAeBe stars were located in the region of 0 < (*H* − *Ks*) < 2 and −0.2 < (*J* − *H*) < 2.0, which means −4.08 < *QJHK* < 2.  Finally, they found sources that had very large infrared excesses with (*H* − *K*) > 1.5. Sources with less excess IR radiation probably originated from free ionized gas emissions. Instead, those sources with large excesses of IR radiation were likely indicative of the influence of a bright nebulosity related to low temperatures or mainly thermal emissions from the surrounding dust [4].

A somewhat similar result was achieved for the CBe stars. Our sample of CBe stars encompassed bright Be field stars and Galactic cluster members. Among them, we recognized two types of CBe stars. One group, with small near-IR excess, resembled the normal B-type stars. The second group consisted of a few CBe stars (and a small group of stars classified as normal B-type stars) with a moderate IR excess, sharing properties with the LBVs. The latter could be better identified in the QCDs of Figure 8, particularly, in the

*QBpRpHK* vs. (*Bp* − *Rp*) (top-left panel) and *QBp JHK* vs. (*Bp* − *Ks*) (middle-right panel) diagrams. Zhang et al. [66] explained that a small radiation excess could be accounted for by hot circumstellar plasma, but a larger IR excess would require thermal emission from circumstellar dust. Lee and Chen [69] proposed that the dust grains were produced by condensation in the expanding stellar atmosphere, and, therefore, the dust particles should be very small in size to reprocess starlight efficiently. In contrast, the HAeBe stars showed a very large infrared excess which arose from relatively large grains as part of surplus star-forming materials. However, Rivinius et al. [5] asserted that dust had never been found in the envelopes of Be stars and that it was precisely the presence of dust that made it possible to distinguish a CBe from a HAeBe. As a consequence, stars with large color excesses are likely misclassified and would not be CBe.

#### **6. Conclusions**

This work discussed the properties of different color–color diagrams and equivalent Johnson's reddening-free *Q* parameters built using a combination of the Gaia and 2MASS photometry. Although the *Q* values were free from the ISM reddening effect, they still showed a dependence on the reddening law in the sky region. Nevertheless, using appropriate combinations of color indices, it was possible to minimize this effect among the O, B, and A stars. Furthermore, by analyzing the different diagrams, we were able to establish criteria to identify stars of spectral type B (Equations (8), (10) and (11)), as well as to identify early-type stars with emission lines of different nature. Main selection criteria are summarize in Table 3.

**Table 3.** Summary of classification criteria for normal B-type stars, LBVs, and B[e]SGs.

In addition to the well-known (*J* − *H*) vs. (*H* − *Ks*) color–color diagram, which is an excellent dust tracer of circumstellar envelopes, we also found that the (*Bp* − *Rp*) vs. (*H* − *Ks*) diagram could be used to identify emission-line B-type stars corrected by interstellar reddening, but with strong intrinsic IR emission. Particularly interesting were the (*Rp* − *Ks*) vs. (*Bp* − *Rp*) and (*J* − *Ks*) vs. (*Bp* − *Ks*) diagrams, as they might be helpful to search for CBe stars with large infrared excesses, because the normal B-type stars distribute along a narrow and well-defined band.

The main advantage of the QCDs over the color–color ones was that the B-type stars did not mix with the M-type stars, as shown, particularly in the *QBpRpHK* vs. (*Bp* − *Rp*), and *QBp JHK* vs. (*Bp* − *Ks*) diagrams. Moreover, early-type stars that have anomalous reddening, due to the presence of circumstellar envelopes, can be easily recognized. This is because the effect of interstellar extinction causes stars of the same spectral types to shift to the right of the diagram. At the same time, the presence of the circumstellar envelope makes the *Q* value more negative (when using an IR color index). We suggest using different selection criteria, such as those that arise from the *QJKHK* vs. (*H* − *Ks*) and *QBp JHK* vs. (*Bp* − *Ks*) diagrams, to search for, and classify, B-type stars with emission lines. These diagrams also allow the LBV stars to be well-recognized from HAeBes and B[e]SGs. In addition, the region of the HAeBe stars can be separated from that of the B[e]SG stars by a straight line. However, there was still significant contamination between HAeBes and LBVs, that it would be desirable to get rid of when making source selections.

As was already stressed, QQDs are free from interstellar extinction, although they are sensitive to anomalous reddening. However, the presence of a circumstellar envelope produces a shift of the *Q*-values in a direction almost perpendicular to the effect of the anomalous reddening. Then, the QQDs are powerful tools to identify stars with intrinsic IR emission.

Using the *Q* parameters allowed us to recognize two classes of CBe stars. One group resembled the normal B-type stars with small near-IR excess, and the other group had a few objects showing moderate to large IR excess. This second group might show dusty envelopes formed as a result of the evolution and would not be CBe stars. We also detected three kinds of HAeBe stars. The first group had small IR radiation excess. Stars in this group shared properties with CBe stars and LBVs. A second group showed a moderate contribution from dust, and a third group showed a significant IR excess comparable with the properties of B[e]SGs. All these groups need further study to understand their dust grain properties and formation processes.

Thus, based on the location of the stars in different CCDs, QCDs, and QQDs, we were able to state various selection criteria that allowed separating, as much as possible, the diverse groups of B-type stars with a moderate and large IR excess (e.g., CBes, HAeBes, LBVs, and B[e]SGs). Both the diagrams and stated criteria can be very useful tools for automated designs of machine learning and optimal search algorithms.

**Author Contributions:** Conceptualization, Y.A. and L.S.C.; methodology, Y.A.; formal analysis, Y.A.; investigation, Y.A.; writing—original draft preparation, Y.A. and L.S.C.; writing—review and editing, L.S.C.; visualization, Y.A.; supervision, L.S.C.; funding acquisition, Y.A. and L.S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by CONICET (PIP 1337) and the Universidad Nacional de La Plata (Programa de Incentivos 11/G160), Argentina. In addition, this project has received funding from the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant Agreement No. 823734.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This publication made use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work presented results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular, the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia (the last accessed date was 1 February 2023). The Gaia archive website is https://archives.esac.esa.int/gaia (the last accessed date was 1 February 2023).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Notes**

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Synthetic Light Curve Design for Pulsating Binary Stars to Compare the Efficiency in the Detection of Periodicities**

**Aldana Alberici Adam 1,\* ,† , Gunther F. Avila Marín 2 , Alejandra Christen 2 and Lydia Sonia Cidale 1,3**

	- 2 Instituto de Estadística, Universidad de Valparaíso, Valparaíso 2340000, Chile;
	-<EMAIL_ADDRESS>(G.F.A.M<EMAIL_ADDRESS>(A.C.)
	- 3 Facultad de Ciencias Astronómicas y Geofísicas, La Plata B1900FWA, Argentina
	- **\*** Correspondence<EMAIL_ADDRESS>	- † Current address: Observatorio Astronómico, Paseo del Bosque, La Plata B1900FWA, Argentina.

**Abstract:** B supergiant stars pulsate in regular and quasi-regular oscillations resulting in intricate light variations that might conceal their binary nature. To discuss possible observational bias in a light curve, we performed a simulation design of a binary star affected by sinusoidal functions emulating pulsation phenomena. The Period04 tool and the WaveletComp package of R were used for this purpose. Thirty-two models were analysed based on a combination of two values on each of the *k* = 6 variables, such as multiple pulsations, the amplitude of the pulsation, the pulsation frequency, the beating phenomenon, the light-time effect, and regular or quasi-regular periods. These synthetic models, unlike others, consider an ARMA (1, 1) statistical noise, irregular sampling, and a gap of about 4 days. Comparing Morlet wavelet with Fourier methods, we observed that the orbital period and its harmonics were well detected in most cases. Although the Fourier method provided more accurate period detection, the wavelet analysis found it more times. Periods seen with the wavelet method have a shift due to the slightly irregular time scale used. The pulsation period hitting rate depends on the wave amplitude and frequency with respect to eclipse depth and orbital period. None of the methods was able to distinguish accurate periods leading to a beating phenomenon when they were longer than the orbital period, resulting, in both cases, in an intermediate value. When the beating period was shorter, the Fourier analysis found it in all cases except for unsolved quasi-regular periods. Overall, the Morlet wavelet analysis performance was lower than the Fourier analysis. Considering the strengths and disadvantages found in these methods, we recommend using at least two diagnosis tools for a detailed time series data analysis to obtain confident results. Moreover, a fine-tuning of trial periods by applying phase diagrams would be helpful for recovering accurate values. The combined analysis could reduce observational bias in searching binaries using photometric techniques.

**Keywords:** methods: statistical methods: numerical; binaries: eclipsing; stars: oscillations (including pulsations)

#### **1. Introduction**

Observational surveys of massive stars of spectral type O show that most are forming binary systems [1–3]. Studies carried out in the region of 30 Dor also show that the binarity frequencies among dwarf and giant B-type stars agree with those of the O-type stars [4,5]. However, in contrast with the high binary incidence among dwarf massive stars, binary systems with B supergiant components are very rare. Barbá et al. [6] reported a noticeable sharp drop in the number of spectroscopic binaries between the supergiants of spectral types O9.7 and B0. These authors attributed this peculiarity to observational bias but suggested also discussing the action of possible evolutive scenarios such as binary mergers, binary disruptions, etc., to explain this effect.

**Citation:** Alberici Adam, A.; Avila Marín, G.F.; Christen, A.; Cidale, L.S. Synthetic Light Curve Design for Pulsating Binary Stars to Compare the Efficiency in the Detection of Periodicities. *Galaxies* **2023**, *11*, 69. https://doi.org/10.3390/ galaxies11030069

Academic Editor: Oleg Malkov

Received: 25 March 2023 Revised: 19 May 2023 Accepted: 22 May 2023 Published: 31 May 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

On the other hand, massive stars pulsate in regular and quasi-regular (strange modes) oscillations (cf. [7]). They may show *β* Cephei type modes (low-order p and g modes) and Slowly Pulsating B (SPB) type variability (high-order g modes) [8,9]. Strange modes might also be excited with long periods, of the order of hours or days [10–13]. The pulsation mechanism is, in general, well understood (e.g., the *κ* mechanism, *e* mechanism, or strange-mode instabilities) but for other groups, such as the *γ* Dor stars, roAp stars, and S Dor stars, no consensus has yet been reached [14]. These pulsations are revealed by photometric variations that can be identified by asteroseismological methods, while in the stellar spectrum, they generate variability in the broadening of the line profiles [15,16]. However, when pulsating stars are found in binary systems, their luminosity variations are modulated by the orbital motion [17]. This is the light-time effect, equivalent to a periodic Doppler shift of the pulsation frequency. Therefore, these effects generate complex patterns in the light curves that could make it difficult to detect or estimate periods or even conceal a binary detection.

The importance behind the detection of pulsating stars in massive binary systems lies in their impact on the study of the mass discrepancy problem, the theory of stellar interiors, and the verification of stellar evolution models [18], as they allow the masses of each component to be measured accurately and independently.

The search for periodicities is one of the most basic tasks in time series data analysis. It is a crucial feature for classifying variable stars and deriving stellar parameters. It is usually solved using an estimation function called a periodogram. Various types of periodograms are used in practice; for example, the Lomb–Scargle periodogram [19,20] is a known algorithm for detecting periodic signals in unevenly and regularly sampled time series. When periodic or quasi-periodic fluctuations arise in an interrupted or transient event [21], the wavelet transform is well suited for detecting changes in the parameters of signals due to its characteristic of focusing on a limited time interval of the data [22]. Nevertheless, the wavelet analysis is appropriate for time series with regular sampling. Today this problem has been partially addressed since satellite missions (Kepler, TESS, among others) provide photometric data of pulsating stars and binary systems with almost fixed cadences of a few seconds or minutes. Hence, using techniques for regularly sampled data is pertinent in this context [21,23–25].

On the other hand, it is well-known that power spectra can often be misleading. The performance of many different methods for period detection and estimation in light curves has been tested in the literature, for example by Graham et al. [26]. Those authors examined a huge number of sources from three datasets (CRTS, MACHO, ASAS) and ran the data through eleven different algorithms such as the Fourier transform, Lomb– Scargle, phase dispersion minimisation, among others. The study of different observed light curves showed that, at best, period-finding algorithms could retrieve the period of a regularly periodic object with a reasonable degree of accuracy in only about 50% of cases. However, the situation is much worse for stellar objects with semi-periodic, quasi-periodic or multi-periodic variability; typically, only around 10–20% of the cases are the periods successfully recovered, evidencing the low level of accuracy for most of the methods, as shown in Figure 19 given by Graham et al. [26]. Although the authors performed an excellent comparative study among many methods for period searching, it still has the disadvantage that the authentic periods of the signal could be miss-detected because, even when they are accurate, they might be inexact. Based on this weakness, it is necessary to explore simulation designs for which the true values of the periods are known but compare the performance of a method not analysed in the above work, such as wavelet analysis, with the revisited Fourier analysis. This is the primary goal of the work.

In addition, one of the usual drawbacks when analysing observed data is the noise level of the signal, which generates consequently noisy periodograms. Thus, specific peaks in a periodogram might be spurious and not due to the presence of any real periodic phenomena [27,28]. Other drawbacks can also arise, for example, when the sample data are not evenly spaced in time. This situation is very frequent in the case of astronomical

observations, where the data are out of phase in time due to various natural phenomena. It is well known that these non-equispaced data produce the so-called alias [29], i.e., false periods, where it is difficult to discern which peak of the periodogram is the real one and which is its alias. Unfortunately, simple and rigorous methods are unavailable to solve the alias discrimination task. In practice, we need to have prior information about the phenomenon under investigation to discuss the number and reliability of the possible periods. For example, for binary systems, the half period is often the most significant peak in a periodogram [26].

In summary, for a thorough analysis of the light curves of massive stars, it is necessary to have (i) a valid mathematical model to describe the time series, (ii) a period analysis tool that supports the detection of variable frequencies, and (iii) a homogeneous time distribution of the data to ensure a stable solution. Although the Fourier transformation provides a confident signal frequency spectrum, we often lose information about when this frequency happens. An alternative method is the wavelet analysis [30–32] that allows searching for local periods in a given window of time to obtain both frequency and time resolution.

In this work, we proposed to use a simulation design developed to study how different periodic or quasi-periodic simulated phenomena are reproduced in a wavelet analysis scalogram and a Fourier periodogram. For comparison purposes, we generated multiperiodic synthetic light curves (with pre-fixed periods) for a binary system with a pulsating companion. These synthetic models, unlike others, consider a statistical noise, a gap, and a slightly irregular sampling. In addition, phenomena such as beating and light time effects were also emulated. We obtained surprising and interesting results that will be useful in future studies of the periodicities of light curves and which have possible implications for the low incidence rate of binary systems among B supergiants.

Section 2 introduces the simulation design by detailing the variables and their features. Section 3 briefly describes the time series analysis methods used. Section 4 shows the results obtained from the wavelet and Fourier analyses. Section 5 discusses our results in the context of the expected objectives and future work. Finally, Appendix A shows the complete frequency list for each model and all the periodograms resulting from the Fourier analysis. Appendix B shows the list of periods detected for the synthetic light curves with the wavelet analysis and the power spectrum for each one.

#### **2. Simulation Design**

In statistics, the area of experimental designs is concerned with finding the design that best helps to determine and study the factors that influence a given variable. There are many different types of designs: the factorial, fractional factorial, block, Latin square, and nested designs, among others [33,34]. Factorial designs are those in which experiments corresponding to all possible combinations of the "levels" or values of the factors or measurements of interest are carried out. The 2 *k* designs are factorial designs with *k* factors or variables for which only two levels or values are considered for each. These designs with few values for each factor are often used when it is desired to determine which factors are affecting a certain process or phenomenon. Extreme values of quantitative variables (such as temperature or concentration) or extreme categories of qualitative variables (such as the presence-absence of an attribute, among other possibilities) are considered for the levels. This type of design makes it possible to determine which factors or variables really influence the change or variation of the process or phenomenon. They are used in the first stage of the study when many variables are known about a process, and it is desired to determine its most influential factors. Once the relevant factors have been determined, it is possible to move on to a more exhaustive design, which considers more levels or values for each factor or even a more complex design. Each power-of-2 design considers 2 *k* different runs, which can be a considerable amount depending on the number *k*. Because it is not always feasible to perform 2 *k* experiments or runs, one can consider fractional factorial designs (e.g., 2 *k*−*p* ), which consist of running only a fraction of the total amount of runs given by 2 *k*/2 *p* . Therefore, it is important to optimally select this fraction. For this purpose,

the fraction of the design is chosen in a way that allows obtaining a complete factorial design when *p* factors are discarded, for instance, because they are considered negligible after running the experiment. In this case, the number of factors *p* to be removed agrees with the power *p* in the denominator of 2 *k*/2 *p* . This methodology is based on the belief that if one considers a large number of variables, only a smaller number is relevant, and it does not produce significant variations in the phenomenon. In a later step, the remaining executions could be carried out to obtain a complete 2*k* design.

One of the most important uses of fractional factorial experiments is in factorscreening1 [33,34]. Furthermore, one of the benefits of using a statistical methodology is to be able to have hypothesis tests that will give us ways to quantify the significance of the findings. In this work, a fractional factorial design was used to establish the models to be simulated and underlined through Fourier and wavelet analysis. The design is detailed in the following section.

#### *Designs of Synthetic Light Curves*

In this stage, instead of performing a real experiment, we considered using numerical simulations of light curves to know in advance the result that we should find and thus be able to make a diagnosis of the goodness of each of the methods, Fourier analysis and wavelet analysis, to be used to detect periodicity. In the context of this research, we performed the simulations of synthetic light curves of a pulsating star in an eclipsing binary system. To build them, we used the photometric data of a known binary, HD 19,356, obtained by the TESS (Transiting Exoplanet Survey Satellite) mission [35] with a temporal resolution of 120 s and a gap of about 4 days. Then, to obtain the synthetic light curve, we modelled the data with the software PHOEBE (PHysics Of Eclipsing BinariEs) [36,37], adopting the orbital period of 2.87 days corresponding to the selected binary system [38]. The resulting synthetic curve, referred to as model No. 0, is shown in Figure 1. This model will be the basis for the complete simulation design.

**Figure 1.** Synthetic light curve of an eclipsing binary affected by ARMA (1, 1) noise.

As the time series of many physical processes show stochastic and auto-correlated properties [39,40], we contaminated the synthetic curves with an ARMA (1, 1) noise [41]. Simulations of light curves with ARMA noise have already been used by authors such as Caceres et al. [40,42] to study different phenomena, such as exoplanetary transits. According to these authors, autoregressive modelling is more effective after minimising systematic variations due to instrumental or atmospheric conditions. Alberici [37] performed numerical simulations of light curves with three different types of stochastic processes to model noise: two stationary, a Gaussian process (white noise), an ARMA (1, 1) process (dependent of the past), and a non-stationary GARCH (1, 1) process (also dependent of the past). This author independently tested the Fourier and wavelet analysis and found no difference among the diverse types of noise models when they were analyzed with the same techniques.

Therefore, according to Alberici [37] and Caceres et al. [40,42], we opted for using an ARMA (1, 1) noise model, which has some dependence on the prior values, for the random fluctuation of the star luminosity [40]. The ARMA processes were used in this work to model the synthetic light curves, considered without systematic effects. The magnitude of the variability of the ARMA process was determined at a low value for simplicity, and it is left for future work to perform a study with higher variability of the noise process (i.e., where the variability of the luminosity could affect the detection of the binary star eclipse or pulsations, for example).

Therefore, a simulation design was planned with six variables and two categories (factors) each, defined as follows:

#### 1. Beating phenomenon (*X*1).

Beating is the periodic variation in amplitude at a certain wave point due to the superposition of two waves having slightly different frequencies *fi* and *fj* , producing a time disturbance between constructive and destructive interference. For simplicity, we considered two waves of equal amplitude, *A*, and wave number, *k*, but with slightly different angular frequencies:

$$y\_i(\mathbf{x}, t) = A \cos(k\mathbf{x} - \omega\_i t), \qquad y\_j(\mathbf{x}, t) = A \cos(k\mathbf{x} - \omega\_j t). \tag{1}$$

The resulting wave is the sum of *yi* + *yj* , which gives

$$y(\mathbf{x},t) = 2A\cos\left(\frac{\omega\_i - \omega\_j}{2}t\right)\cos\left(k\mathbf{x} - \frac{\omega\_i + \omega\_j}{2}t\right),\tag{2}$$

where

$$
\omega = 2\pi f.\tag{3}
$$

Simulations with and without beating phenomena were considered.

2. Period of the pulsation (*X*2).

In this variable, we explored the effects on the light curve considering regular pulsation cycles with a period longer or shorter than the binary's orbital period.

3. Amplitude of the pulsation with respect to the depth of the primary eclipse (*X*3).

Two categories were considered for this variable. The first is that the pulsation amplitude is 20% of the depth of the primary eclipse, which we call a large amplitude. The second is that the pulsation amplitude is 4% of the depth of the primary eclipse, referred to as a small amplitude. The values of the amplitudes of the pulsations were chosen in this way because they agree with those observed in TESS light curves of some pulsating binary stars [43].

4. Number of pulsation periods (*X*4).

Scenarios where the light curve is affected by simultaneously one or two pulsation periods were considered.

5. Light-time effect (*X*5).

This scenario considers a star in a binary system that pulsates sinusoidally with a single frequency. Its luminosity varies with time as a consequence of this pulsation. However, the orbital motion of the star leads to a periodic variation of the distance between the observer and the star, so the phase of the observed luminosity variation also varies with the orbital period [17]. The light-time effect leads to frequency multiplets in the Fourier spectra for pulsating binary stars. To model this effect in our synthetic light curves, we relied on the simplest case of a pulsating binary star with a circular orbital motion. Furthermore, to emulate the splitting of the pulsation period, we proceeded to affect the curve with three sinusoidal pulsations: the intrinsic pulsation frequency of the star and two additional oscillation frequencies due to the light-time effect, where the spacing between each of them

and the intrinsic frequency of the pulsation is equal to the orbital frequency (for more details see [17]).

Let *Po* be the orbital period of the star and *Pp* the period of the pulsation, then the multiple components of the luminosity variations will have the following frequencies:

$$f\_{+} = \frac{1}{P\_o} + \frac{1}{P\_p}, \qquad f\_{-} = \left| \frac{1}{P\_o} - \frac{1}{P\_p} \right|. \tag{4}$$

Furthermore, the amplitudes and phases of the multiple components can be used to derive meaningful information about the orbital elements and the radial velocity curve. Scenarios with and without the light-time effect were considered in this work.

6. Quasi-regular pulsations (*X*6).

Let us now consider a set of harmonic oscillators. Periodic motion means that there must be some time *τ* at which all the different oscillators have completed an integer number of cycles, i.e.:

$$\begin{aligned} \omega\_0 \tau &= 2\pi n\_0, \\ \omega\_1 \tau &= 2\pi n\_1, \\ \omega\_2 \tau &= 2\pi n\_2, \\ &\dots \end{aligned} \tag{5}$$

This means that all *nk* are integers; therefore, the quotient *ωi*/*ωj* must be a rational number. If, by contrast, it is irrational, any oscillator *i* or *j* might not complete an integer number of cycles in the time *τ*, giving rise to quasi-periodic pulsations (QPP). This kind of behaviour is a common feature in the solar atmosphere and flare stars where the flaring energy is released by magnetic reconnection (cf. [44]). Quasi-periodic pulsations are also observed in the *α* Cyg variables and are attributed to strange-mode oscillations [45,46]. A star oscillating with some periodic phenomenon would appear in a power spectrum as a peak at precisely that frequency (Dirac-*δ* function). Otherwise, it will appear as a broader (sometimes Lorentz-shaped) peak.

Simulations with and without quasi-regular periods were considered.

Since we have defined *Xi* variables (with *i* = 1, · · · , 6) and each has two levels or categories, the total design is 2 6 = 64 models. Since the design is optimal, we can consider only a fraction of the total design without losing information about the variables involved. Thus, in this work, we analysed 2 6/2 1 = 32 models. The design parameter model is defined in Table 1.

**Table 1.** Simulation design features. The table lists the phenomena, properties, and periods fixed to simulate the pulsations in each model. The light-time effect and the beating phenomenon generate periods that are labeled with the indices "*t*" and "*b*", respectively. In all cases, these effects are present only in the lowest pulsation period.

**Table 1.** *Cont.*

a The pulsation amplitude is with respect to the depth of the primary eclipse. c e: Euler number.

#### **3. Time Series Analysis Methods**

Signals were differentiated into stationary and non-stationary. The first ones were localised in the time since their frequencies did not vary. The spectral analysis of this type of signal was carried out using the Fourier transform, which allows its decomposition into infinite sinusoidal terms, transforming the signal from the time base to the frequency base and back again. When switching from the frequency domain to the time domain, valuable information is lost, which, due to the stationary nature of the signal, is irrelevant. However, in non-stationary signals, i.e., those whose frequencies vary with time, the loss of information becomes very relevant, and it is difficult to determine when a frequency change occurs [47]. This is why wavelet analysis is presented as a tool to obtain a detailed decomposition and reconstruction of signals with abrupt changes, using a multi-resolution analysis with windows of a variable length adapted to the change in frequency of the signal. This is the main property of the wavelet, as it minimises the computational cost of the analysis by not needing infinite terms. Moreover, it improves the detail of the detection by allowing the use of long time intervals in those segments where greater precision is required at low frequencies and smaller intervals where information is required at high frequencies. Unlike Fourier, where the basis functions are sines and cosines of infinite length, in the wavelet analysis, the basis functions are localised functions in frequency and time [48]. Thus, they are a suitable method for studying quasi-regular, transient and discontinuous phenomena, such as light curves of binary stars.

The software used to carry out data analysis are detailed below.

#### *3.1. The Period04 Tool*

The Period04 tool [49] is a hybrid program, written in Java and C++, dedicated to the statistical analysis of large astronomical time series which contain gaps or missing data over a significant time interval. The program allows the extraction of individual frequencies from the multi-periodic content of the time series through an analysis based on the Discrete Fourier Transform (DFT). The light curves were fitted using the formula:

$$\chi(t) = Z + \sum A\_i \sin[2\pi(f\_i t + \phi\_i)],\tag{6}$$

where *Z* is the zero-point, *Ai* is the amplitude, *fi* is the frequency, and *φi* is the corresponding phase.

#### *3.2. The WaveletComp Package*

WaveletComp is an *R* package [50–52] for continuous time series analysis based on the Morlet wavelet family [32]. This family of wavelets leads to a wavelet transform of the time series. The transform can be separated into real and imaginary parts, thus providing information about the local amplitude and instantaneous phase of any periodic process over time.

The Morlet "mother" wavelet [32], in the version implemented by WaveletComp, is:

$$
\phi(t) = \pi^{-1/4} e^{i\omega t} e^{-t^2/2}.\tag{7}
$$

In this frame, the angular frequency *ω* is set to 6 [51], therefore the relation between the wavelet scale, *a* obtained with wavelet analysis, and the "equivalent Fourier period" [53] is given by:

$$\frac{1}{f} = \frac{4\pi a}{\omega + \sqrt{2 + \omega^2}},\tag{8}$$

where

$$\frac{1}{f} \sim 1.033 \, a \qquad \text{if} \qquad \quad \omega = 6.1$$

This correction is given directly by the WaveletComp package.

The Morlet wavelet transform of a time series *f*(*t*) is defined as the convolution of the series with a set of "daughter wavelets" generated by the mother wavelet through translations in time, *τ*, and dilatations in frequency, *a*, being:

$$\mathcal{W}(\tau, a) = \int\_{-\infty}^{\infty} \mathbf{x}\_t \frac{1}{\sqrt{a}} \phi^\*(\frac{t-\tau}{a}),\tag{9}$$

where *φ* ∗ denotes the complex conjugate function.

The local amplitude of any periodic component of the series can be expressed as:

$$\text{Ampl}(\tau, a) = \frac{1}{a} |W(\tau, a)|. \tag{10}$$

Then, the square of the amplitude can be interpreted as the energy density of the wave in the time-frequency domain and is called the *wavelet power spectrum*:

$$\text{Power}(\mathbf{r}, a) = \frac{1}{a^2} |\mathcal{W}(\mathbf{r}, a)|^2. \tag{11}$$

This power spectrum is represented by a scalogram, see Figure 2, where a color code, from blue (low power) to red (high power), is used to indicate the power spectrum range. Maximum probability values are represented with black lines. White lines delimit probability regions with a 0.1 significance level against a white noise null. The WaveletComp

package returns by default a scalogram whose axes are period vs index (ordinate and abscissa, respectively), where index refers to the amount of data in the time series and is directly related to the time axis of the series. More information can be obtained from the average wavelet power over time, which will depend on how many cycles a certain period completes in the series and not only on the magnitude of the wavelet power. The red curve represents values with a significance level less than or equal to 0.05, all of which differ from white noise. The values with the highest average powers indicate the most relevant periods for the signal. A practical example of wavelet analysis for multi-periodic time series with periods 30, 38, and 80 is shown in Figure 2. Periods 30 and 38 generate a beating effect in the series due to their proximity. All periods with higher power are detected with solid black lines in the scalogram, and the significance decreases towards nearby periods. The color code on the right-hand sidebar indicates a higher significance for the period 30 and decreases until the period 80 (see central panel). This is due to the number of cycles each period completes in the time range of the series. In the right panel, on the average wavelet power plot, periods 30 and 38 are more difficult to identify individually, unlike period 80, although they show a striking feature in the scalogram over the entire timeline identifying the presence of the beating. The results show the mixing effect produced in power obtained with the wavelet analysis for nearby periods. In agreement with the power spectrum, the maximum average wavelet power is around period 30 and decreases for longer periods.

**Figure 2.** Example of applying wavelet analysis for a multi-periodic time series. Left panel: sinusoidal series with periods 30, 38, and 80. Central panel: Scalogram or wavelet power spectrum of the series. The top in the color code indicates a higher significance for periods. Periods with a higher power are detected with solid black lines. Right panel: average wavelet power of the series in logarithm scale. In this diagram, the beating periods 30 and 38 are more difficult to identify individually, unlike period 80, but they show a striking feature in the scalogram. Although all the periods represented in red have a significance level less than or equal to 0.05, the most relevant ones, about 30, 38, and 80, are the three highest average wavelet power. They were obtained using the WaveletComp Monte Carlo sampling-based test.

Estimating the spectral power using the continuous wavelet transform of a nonperiodic time series can generate false spectral structures at the beginning and the end of the wavelet spectrum [54]. This is a known problem in the periodogram analysis caused by replacing the integrals over an infinite continuous signal by finite sums [55]. The "zero padding" method is used to solve this problem, which completes with zero values on both edges of the time series after a sample mean is removed. Then, the wavelet spectrum is estimated [53,54]. One of the drawbacks of applying this method is that it produces discontinuities at the edges of the wavelet spectrum. The region of the wavelet spectrum where padding effects become important is known as the cone of influence (COI) [53]. Therefore, the portions of the wavelet spectrum between the time axis and the COI should be considered unreliable because statistical significance tests in that region cannot be trusted [53,56]. Additionally, wavelet analysis performed with the WaveletComp package can only produce accurate results if the time series is equispaced.

#### *3.3. Period Detection Criteria*

To identify the periods detected with both tools in the synthetic light curves, the following selection criteria were established:

#### **4. Results**

Each synthetic light curve generated with the design model of Table 1 is shown in Figure 3. These curves were analysed using the wavelet and Fourier analysis. To give a detailed description of the results obtained and how they are interpreted, we selected a group of four models with specific characteristics in common and compared them with Model No. 0 (shown in Figure 1). For the remaining 28 models, we will comment on the general results obtained, following the same analysis criteria. The results for the complete set of models are in Appendices A and B.

Table 2 lists the independent periods detected for Models No. 0, 1, 5, 17, and 21 with each method and their corresponding errors. Among these five models, the magnitude of the pulsation amplitude, the presence or absence of the beating phenomenon, and the values of the pulsation periods vary. The respective wavelet power spectrum, average wavelet power, and periodograms are illustrated in Figures 4 and 5.

**Figure 3.** Synthetic light curves obtained from the simulation design. The model number is placed according to the list in Table 1.

**Table 2.** Comparison of Models No. 0, 1, 5, 17 and 21 with similar properties. The independent periods detected with wavelet (w) and Fourier (F) analysis are provided and their corresponding errors. For the periods obtained with the wavelet analysis, the weighted average with respect to the amount of data for each part of the curve is listed.

For Model No. 0, the wavelet and Fourier analysis recognises the binary's orbital period as a fundamental frequency and finds its harmonics. The wavelet scalogram for both parts of the curve shows a structure (seen as "faces") at the position of the minima of the light curve, also reproducing the eclipses of the binary (compare the timeline and structures of Figures 1 and 4). However, the period detected by the wavelet tool is systematically slightly shorter than the original orbital period. This shift could be related to the fact that the time scale, adopted from a TESS light curve observation, is not completely equispaced.

When examining the wavelet power diagrams for Model No. 1, we noticed some differences from Model No. 0, particularly in the features that identify the eclipses. For example, it loses detail in finding harmonics with periods 1.43 and 0.95 days, maybe due to the size of pulsation amplitude (20%). Even so, it finds the binary 2.87 days period with an error of 0.09, although in the second part of the curve and with much lower power than in Model No. 0. Additionally, a period of about 1.15 days is found, and its trace shows a sinusoidal modulation in the scalogram. Although it could be a spurious value, we associated this period as an independent value, likely with a harmonic of 8.96 or 10 days (8 × 1.15 = 9.2), where the value 8.96 days were artificially introduced to obtain the beating. However, this shift could be due to three near pulsation periods: 3e, 8.96, and 10 days. Furthermore, a period of 7.852 days is also present with high power and could be related to the 3e days period. Even though the wavelet analysis does not find those exact values, the software likely detects the presence of significant periods in that range.

Model No. 5 is affected by the same pulsation periods as Model No. 1 but with a smaller amplitude. Therefore, the model's scalogram and power diagram show a lower significance, but the same periods and scalogram features are maintained. As shown in Figure 2, the wavelet analysis performed with the WaveletComp package presents complications in differentiating nearby periods, particularly cases with a beating, such as Models No. 1 and 5, resulting in an intermediate value.

**Figure 4.** Results of the wavelet analysis for Model No. 0, 1, 5, 17, and 21. From left to right, the wavelet power spectrum and the average wavelet power for the second part of each synthetic light curve.

**Figure 5.** Periodograms obtained for the models No. 0, 1, 5, 17, and 21.

When analysing Model No. 17, we observed that the structures present in the scalogram are very similar to those of Model No. 1, both with a pulsation amplitude of 20% with respect to the depth of the primary eclipse. Since they share two periods, the 7.9 days period and the harmonic of 1.15 days, they could be related to the periods of 3e and 10 days.

Finally, the scalogram and average wavelet power of Model No. 21 is very similar to that found in Model No. 5 since they share two pulsations with amplitudes of 4%. However, in this case, it detects the binary period and another with a 7.92 days period, but with lower power than in Model No. 5.

The periodograms resulting from the Fourier analysis show remarkable similarities between the models affected by low-amplitude pulsations regardless of the periods involved. This is particularly seen in models No. 5 and 21. At the same time, slight variations are observed with respect to model No. 0. The pattern changes markedly for models No. 17 and 21, as seen in Figure 5). A similar situation occurs with models No. 6–8.

Quasi-periodic oscillations decrease the power of detection of the Fourier Transform method, tending not to detect them. The Period04 tool accurately sees the fundamental frequency of the binary orbital motion (*P*orb = 2.87 ± 0.01 days) and its harmonics in all cases. However, it observes two pulsation periods with a significant error close to *P* = 10 days, detected as 9.26 days and 12.9 days in Model No. 1 and 9.47 days in Model No. 5. The 9.26 days (or 9.47 days) period could be likely related to the beating phenomenon with *P* = 8.96 or the pulsation with *P* = 10 days, or a mix of both. This is supported because, for example, the period 9.26 days also appears in Model No. 2 (see Appendix A) that has regular oscillations (with *P* = 7 days and *P* = 10 days, and a beating of 7.7 days), but in this case, the period 7 days is inexactly seen as *P* = 6.657 days (see also Model No. 6). Similarly to the wavelet analysis, a period of 1.1 days (*f* = 0.9036 days−1 ) is present in Models No. 1 and 17. Models with the same set of parameters but low amplitude oscillations (i.e., Models No. 5 and 21) behave similarly. Model No. 21, similar to Model No. 17, detects an accurate period for *P* = 10 days.

Generally, the orbital period is well-detected in most models with Fourier analysis; only a few detect half of the period (i.e., Models No. 4, 11, 20, and 27). These models present regular oscillations with short periods or have a quasi-regular period √ 2. The latter is similar to the half-orbital period. Models with quasi-regular pulsations show more dispersion of values, mainly when the time light effect is included. The Fourier analysis finds all periods in 14 models (Models No. 10, 12, 14, 16–19, 21–23, 26, 28, 30 and 32), 9 of which have low pulse amplitude (4%).

On the other hand, the wavelet analysis found the binary's period satisfactorily, although it underestimated the value with an average of 2.81 days (of 31 models), which suggests a period shift of 0.06 days. Concerning the rest of the periods, this method found an equal or smaller number of periods than the Fourier analysis. Although the mean value has a systematic error, its variance is similar to that obtained with the Fourier method. It detected a few periods, mainly in models with low pulsation amplitude and where the period values coincide with a harmonic of the binary's orbital period. However, it found all periods in four models (Models No. 17, 18, 26 and 28). These have a high pulsation amplitude (20%) and absence of beating phenomenon.

A comparison of the methods is shown in Figure 6. The boxplots of the hit ratios of all the prefixed periods obtained by wavelet and Fourier analysis are in the left and right panels, respectively. We see that the median percentages (thick black line inside the box) are much lower for the wavelet analysis (<0.4) than for the Fourier analysis (about 0.8).

Figure 7 presents a boxplot comparison for both techniques, of the hit ratio segmented by the categories of five of the variables included in the simulation design: beating phenomenon (*X*1), period of the pulsation (*X*2), amplitude of the pulsation (*X*3), number of pulsations periods (*X*4), quasi-regular pulsations (*X*6), calculated with a wavelet and Fourier analysis, respectively. Since we were dealing with the half-model design, the result does not provide information about one variable.

Regarding the results obtained for the wavelet analysis, we can observe that the variable that influences its performance the most is the pulsation amplitude. At low amplitude, it is difficult to detect the intrinsic period of the model. The presence of the beating phenomenon and the fact that the period of the pulsation is quasi-regular also have some influence. The hit ratio increases when the pulsation period is greater than the orbital period and when there is more than one pulsation (see Figure 7). The analysis of the variance of the percentage of hits of a quasi-periodic oscillation with a large pulsation amplitude (i.e., 20%) finds the period only the 55% of the time. While for a regular oscillation, it detects the period 85% of the time. Instead, if the pulsation amplitude is low, it finds the period about 30% of the time in both cases. Quasi-regular or regular period detection is also affected by their interaction with the orbital period, depending if they are shorter or longer than it. The regular period detection rate varies from 60% (for shorter periods) to 50% (for longer periods) and for quasi-regular periods from 65% to 30%.

Unlike wavelet analysis, Fourier analysis is unaffected by the pulsation amplitude. The proportion of detected periods (an approximate value) increases if the period pulsation is longer than the orbital one, but they are more accurate when it is shorter (see Figure 7). The presence of beating and quasi-period decreases the detection capability a little. In both cases, the medians (thick black line in the box) is about 75%, which changes to 100% without the phenomenon. From an analysis of the variance of the hit ratio, we obtain that if the pulsation period is shorter than the orbital period, the presence of the beat phenomenon improves period detection (70% vs 80%). Otherwise, if the pulsation period is longer, the presence of this effect worsens it (95% vs 75%). Quasi-regular or regular period detection is also affected by their interaction with the orbital period, depending if they are shorter or longer than it. The period detection rate for regular oscillations varies from 100% (for shorter periods) to 80% (for longer periods) and for quasi-regular periods from 50% to 85%. For pulsation periods longer than orbital periods, the performance is similar.

**Figure 6.** Comparison between the proportion of correctness of the original period obtained by wavelet (cyan) and Fourier (grey) analysis. The thick solid black line represents the median of the data, and the red symbol stands for an outlier or extreme value.

**Figure 7.** Boxplots of hit ratio of the original periods with the wavelet (cyan) and Fourier (grey) analysis for each variable involved. The thick solid black line represents the median of the data, and the red symbol stands for an outlier or extreme value. The hit rate is on the *y*-axis and the two categories (indicated by −1 and 1) for each variable are shown on *x*-axis: beating phenomenon (absence/presence), period of the pulsation (minor/major), amplitude of the pulsation (4%, 20%), number of pulsations periods (1/2), and quasi-regular pulsations (absence/presence).

#### **5. Discussion and Conclusions**

It is widely known that power spectra can be misleading and frequently lead to inaccurate conclusions. For this reason, many works present a comparative analysis of popular algorithms applicable to period searching. Sometimes the study of one particular light curve with various numerical methods resulted in a high dispersion or undetected periods as reported by Graham et al. [26] and Distefano et al. [59].

The objective of the current work was to evaluate the response of two methods, based on the Morlet wavelet function and continuous Fourier transform, in the production of bias, alias, or shifts in period detection using a simulation design with pre-fixed and known periods. For this purpose, we performed a simulation design of synthetic light curves of a pulsating star in an eclipsing binary system with an orbital period of 2.87 days. It is a 2 *k*−*p* -statistical design with *k* = 6 variables, in which the beating phenomenon and the light-time effect were considered. Of the total design method, only one-half was executed. The 32 time series models were analysed using the WaveletComp and Period04 tools. Although many time series analysis tools are available, the application of such a large simulation design to more than two techniques would be unattainable.

The obtained synthetic light curves, shown in Figure 3, present similarities with the light curves of pulsating binary stars obtained with the TESS mission (see Figure 8) in shape and cadence. Among them, we can highlight the similarity between HD 152,248 with model No. 9, TIC 80,042,405 with model No. 27, and TIC 97,467,902 with models with low amplitudes with respect to the depth of the primary eclipse (models No. 5–8, 13–16, 21–24, 29–32).

**Figure 8.** Light curves obtained by TESS of pulsating binary stars.

Comparing both applied methods, we observed that the orbital period and its harmonics were well detected in most models. Even though the Fourier method provided more accurate period detection, the wavelet analysis found it more times. When estimating the orbital period, the WaveletComp package trends to provide a low period value or period shift of approximately 0.06 days in all the models, but this has a similar dispersion to the orbital period obtained from the Fourier analysis. In agreement with this finding, Foster [22] also reported through numerical simulations that for a simple sinusoid, the maximum frequency obtained with wavelet analysis has an offset with respect to the signal frequency, which can be compensated. This author also commented that under irregular sampling, the detection of the period by the wavelet analysis varies with time, making it difficult to give an exact value of the same, especially with the increase in missing data in the sampling, which can reach up to 20% at high frequencies or find spurious periods at low frequencies. Although we considered "almost" a regular sampling, it may be that some of the still sampling irregularity present in the timeline scale selected from TESS observations has affected the performance of the WaveletComp package.

As shown in the previous section, pulsation amplitude is a very powerful variable in period detection. In the case of the wavelet analysis, more satisfactory results are found recovering the periods of the pulsations with large amplitudes. On the other hand, the Fourier analysis did not show difficulties in recognising the periods when the amplitude is low, but they were detected with significantly lower power than models with large amplitudes.

Neither tool was able to distinguish accurately the periods leading to a beating phenomenon when they were longer than the orbital period, resulting in both cases in an intermediate value. However, when the beating was formed by periodic pulsations shorter than the orbital period, the Fourier analysis found it in all cases, but quasi-regular periods remained unresolved. Concerning the beating period √ 2, it is necessary to remember that this period coincided (unfortunately, but possible in a real scenario) with the first harmonic of the selected orbital period of 1.43 days. Therefore, when a value of 1.4 days was detected, it was considered the harmonic of the orbital period. If the values set in the simulation design were unknown, it would be impossible to distinguish between them. The same holds for the 0.94 days period. By analogy, wavelet analysis does not detect these periods either.

Finally, the wavelet power spectrum revealed in all models a particular structure located at the position of the primary eclipses of the binary in the synthetic light curve, even showing the instant of time in which they occur. Then, if a similar structure is found in a real data set, it could indicate binarity. In addition, as mentioned before, several harmonics of the orbital period were detected by both tools, which is also an indicator of binarity.

In general, the Fourier analysis by means Period04 tool had a higher proportion of periods correctly and accurately detected than the analysis performed with the WaveletComp package. However, the discrepancy could be attributed to the less information available in the latter when performing the period analysis since splitting the light curve into two parts was necessary. Despite the high cadence, the TESS data outside the large gap also present a slightly irregular sampling, which affects the analysis. Furthermore, this package uses a somewhat coarse resolution level, preventing it from accurately determining the periods encountered. Foster [22] also found some kind of asymmetry in period detections, with more values at low frequencies and suppression at high frequencies.

Considering the period hit rate analysed using synthetic light curves for eclipsing binaries, a bias in detecting some pulsation periods (beating phenomenon and light time effects) in multi-periodic light curves would be expected. The hit rate depends on the amplitude and frequency of the pulsation with respect to eclipse depth and orbital period. The situation would be more unfavourable for purely elliptical double stars, partial eclipses, or systems influenced by the gravitational attraction of a third body that might behave similarly to a sinusoidal pulsation. In some cases, only a light modulation could likely be acting because the binary orbit is not aligned with the observer's direction. In these cases, the period hit rate would be lower, masking the binary detection, for example, among B supergiants.

Although the recovery of all pre-fixed periods was not feasible using WaveletComp and Period04, valuable insights were obtained regarding how each technique responds to the various modelled situations. For instance, we observed that WaveletComp can only detect an intermediate value when confronted with two close periods due to the poor spectral resolution of this package. Additionally, we have noted that Period04 fails to recognise quasi-regular periods, but has good accuracy in detecting high-frequency oscillations. Furthermore, we found with both tools structures that indicate binarity (in the periodograms and scalograms), information that will be very useful when performing these analyses on the light curves of binary star candidates, among other findings.

Considering the advantages and disadvantages of the techniques analysed, we recommend using at least two diagnosis tools to conduct a detailed analysis of time series data to obtain confident results. Additionally, a second criterion should be considered to test the accuracy of the found frequencies or periods since the methods might find a harmonic, inaccurate period, or spurious value. Particularly, the wavelet method might predict only a harmonic of the true period. Since approximated values were detected with at least one method, a fine-tuning of trial periods around a certain percentage of its value, by applying phase diagrams, would help to recover accurate values.

Evaluating the behaviour of other tools based on wavelets that do not have this weakness (e.g., Weighted wavelet Z-transform (WWZ) [22]) is one of the goals to follow. This method also enables an improvement in the resolution of the period searches, allowing a high degree of accuracy. In addition, considering even more realistic models or observed light curves that present a gap and time-varying periods would allow us to take full advantage of the main feature of wavelet analysis.

**Author Contributions:** Formal analysis, A.A.A., G.F.A.M., A.C. and L.S.C.; Investigation, A.A.A., G.F.A.M., A.C. and L.S.C.; Methodology, A.A.A., G.F.A.M., A.C. and L.S.C.; Writing—original draft, A.A.A.; Writing—review & editing, G.F.A.M., A.C. and L.S.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project has received funding from the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant Agreement No. 823734. AA and LC thanks financial support from CONICET (PIP 1337) and the Universidad Nacional de La Plata (Programa de Incentivos 11/G160), Argentina. AC and GA acknowledges support from Centro de Estudios Atmosféricos y Astroestadística (CEAAS), Universidad de Valparaíso.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

TESS Transiting Exoplanet Survey Satellite PHOEBE PHysics Of Eclipsing BinariEs

#### **Appendix A. Fourier Analysis Results**

For each synthetic light curve, the frequencies with a signal/noise ratio ≥ 5 and their errors, are listed in Table A1. The fixed listed periods of the design models are listed in column 2. It does not include the period of the binary system (2.87 days) present in all the models, except for model No. 0, where is the only period present. Columns 3 and 4 list the frequency and the corresponding period in days. Column 5 indicates the independent frequencies found and the linear relationship among them. The respective periodograms are shown in Figures A1 and A2.

**Table A1.** *Cont.*

**Table A1.** *Cont.*

**Table A1.** *Cont.*

**Figure A1.** *Cont.*

**Figure A1.** Periodograms obtained from the Fourier analysis for models No. 0 to No. 23.

**Figure A2.** Periodograms obtained from the Fourier analysis for models No. 24 to No. 32.

#### **Appendix B. Wavelet Analysis Results**

The periods detected for each synthetic light curve with the wavelet analysis, and their errors, are listed in Table A2. The list of set periods in the simulation design is listed in column 2. It does not include the period of the binary system (2.87 days) present in all the models, except for model No. 0, where is the only period present. Columns 3 and 4 lists the periods present in the first and second parts of the curve. Column 5 indicates the independent periods found and the linear relationship among them. The wavelet power spectrum and the average wavelet power are illustrated in Figures A3–A8.

**Table A2.** Wavelet analysis results reporting the significant periods and harmonics. The detected independent periods are highlighted in blue. The symbol e is the Euler number.

**Table A2.** *Cont.*

**Table A2.** *Cont.*

**Table A2.** *Cont.*

**Figure A3.** Results of the wavelet analysis for models No. 0 to No. 5. From left to right the wavelet power spectrum and the average wavelet power for the first and second part of each synthetic curve.

**Figure A4.** Results of the wavelet analysis for models No. 6 to No. 11. From left to right the wavelet power spectrum and the average wavelet power for the first and second part of each synthetic curve.

**Figure A5.** Results of the wavelet analysis for models No. 12 to No. 17. From left to right the wavelet power spectrum and the average wavelet power for the first and second part of each synthetic curve.

**Figure A6.** Results of the wavelet analysis for models No. 18 to No. 23. From left to right the wavelet power spectrum and the average wavelet power for the first and second part of each synthetic curve.

**Figure A7.** Results of the wavelet analysis for models No. 24 to No. 29. From left to right the wavelet power spectrum and the average wavelet power for the first and second part of each synthetic curve.

**Figure A8.** Results of the wavelet analysis for models No. 30 to No. 32. From left to right the wavelet power spectrum and the average wavelet power for the first and second part of each synthetic curve.

#### **Note**

1 The search for the most important factors that may be varied in a simulated system.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **A Mini Atlas of H-Band Spectra of Southern Symbiotic Stars**

**Paula Esther Marchiano 1,2,\* , María Laura Arias 1,2,† , Michaela Kraus 3 , Michalis Kourniotis 3 , Andrea Fabiana Torres 1,2,† , Lydia Sonia Cidale 1,2,† and Marcelo Borges Fernandes 4**

**Abstract:** Symbiotic stars are interacting binary systems composed of an evolved star (generally a late-type red giant) and a degenerate or dwarf companion in orbit close enough for mass transfer to occur. Understanding the status of the late-type star is important for developing binary models for the symbiotic systems as it affects the transfer of matter needed to activate the hot component. Infrared observations have been very useful in probing the nature of late-type stars in symbiotic systems. This work presents a set of symbiotic stars observed with SOAR/OSIRIS (*R*∼3000) in the H-band. We aimed to search for possible molecular circumstellar emission, to characterize the cool companion in these systems, and to confront the new findings with those obtained from the previous K-band classifications. We detected molecular emission from just one object, BI Cru, which displays the second-overtone CO-bands. To fit the observed photospheric CO absorption bands, we used the MARCS atmosphere models. We present our results as a mini atlas of symbiotic stars in the near-infrared region to facilitate the comparison among different observed symbiotic systems.

**Keywords:** techniques: spectroscopic; stars: symbiotic, binaries; stars: infrared

#### **1. Introduction**

Symbiotic stars display a combined spectrum that shows high excitation emission lines superposed on the absorption lines of two stellar components: a late-type giant and a hot compact star. This characteristic led to defining symbiotic stars as strongly interacting binaries consisting of a cool giant of spectral type M (sometimes K- or G-type giant) as a donor, which transfers material to an accretor compact star, most commonly a hot white dwarf (or, in some systems, a neutron star). Such a binary system is surrounded by a circumstellar nebula enriched by the material lost from both components [1,2]. In addition, there are other contributions to the spectral energy distribution (SED) of a symbiotic star, such as emission from heated dust and jets or collimated bipolar outflows. Therefore, the entire SED of symbiotic systems is a superposition of all these contributions, which dominate the emission spectrum in different wavelength regions. With additional processes such as mass loss, accretion, and ionization, the spectrum of each symbiotic object is unique, and as a group, they appear rather heterogeneous. In this sense, Reference [3] proposed a physical definition where "*a symbiotic system is a binary in which a red giant transfers enough material to a compact companion to produce an observable signal at any wavelength*".

Regarding the near-infrared region, symbiotic stars have been divided into two types [4]. If they display an infrared (IR) excess emission that results from dust, they are called D-type (dust). In these systems, the cool giant is a very evolved Mira variable surrounded by warm dust. If they display an IR continuum emission from the cool companion, the symbiotics belong to the S-type (stellar), in which the cool star is a regular red giant,

**Citation:** Marchiano, P.E.; Arias, M.L.; Kraus, M.; Kourniotis, M.; Torres, A.F., Cidale, L.S.; Borges Fernandes, M. A Mini Atlas of H-Band Spectra of Southern Symbiotic Stars. *Galaxies* **2023**, *11*, 80. https://doi.org/ 10.3390/galaxies11040080

Academic Editor: Oleg Malkov

Received: 29 April 2023 Revised: 4 June 2023 Accepted: 12 June 2023 Published: 22 June 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

often filling its Roche lobe. The orbital periods of S-type symbiotics are about 2–3 years, while those of the D-types are at least an order of magnitude longer [5]. There is a third type of symbiotic system introduced by Allen et al. [6]: D0 -type systems, which are a sub-type of the classical D-types. In these symbiotics, a significantly hotter star of spectral types F and G constitutes the cool component. Furthermore, like D-type symbiotics, D0 -type systems also exhibit IR dust emission. However, their SEDs reveal distinct characteristics, with a nearly flat profile indicative of a cool companion and a dusty shell. This is in contrast to D- and S-types, which exhibit intensity peaks at wavelengths between 2 and 4 µm and between 0.8 and 1.7 µm, respectively [7].

As previously mentioned, the accretion object for most symbiotic stars is a white dwarf. However, in a small number of systems, it is believed to be a neutron star. These symbiotic X-ray binaries (SyXBs) represent a special class of symbiotic binaries that have garnered significant interest among astronomers. One notable example is IGR J16194-2810, which was observed in the low/hard state using the Suzaku X-ray satellite. This particular SyXB has been classified as a low-mass X-ray binary [8,9] since the system is composed of an M-type giant and probably a neutron star [10].

Symbiotic systems show a complex emission spectrum with plenty of lines of different ionization degrees and forbidden transitions and also show a different phenomenology, including stellar pulsations (semi-regular and Mira variables), orbital variations, and slow changes due to dust. These systems are interesting objects to study the formation and evolution of gaseous and dusty circumstellar or circumbinary environments, providing clues to investigate other related objects, such as those showing the B[e] phenomenon. Symbiotic stars may also help to solve one of the timeliest problems in modern astrophysics, the missing progenitors of Type Ia supernovae (SNe Ia).

Observations in the IR are particularly important for those systems in which the cold companion is hidden in the optical spectral range. For example, the water vapor and carbon monoxide absorption bands observed in the 1–4 µm spectral range allow us to describe the characteristics of the late-type stars in these systems. It is interesting to note that, in systems with extremely dense environments, the CO-bands might turn into emissions. So far, BI Cru has been known for a long time to show intense CO-bands in emission [11]. Furthermore, H 1-25 (Hen 2-251) and RX Pup were reported by Schmidt and Mikołajewska [12] as symbiotic stars whose spectra show evidence of weak CO emission instead of absorption.

In a previous work, the medium-resolution K-band spectra of a sample of eight symbiotic objects were presented [13]. These observations were used to characterize the cool companions based on their CO first-overtone absorption bands along with measured equivalent widths of Na I and Br*γ*. Now, in order to obtain a more complete picture of each source, we present the H-band spectra of the same objects and follow a different approach. We derived the effective temperature of the late-type components using synthetic spectra computed with the MARCS atmosphere model' code [14]. However, considering the significant deviation of effective temperatures in red giant stars from temperatures derived using blackbody models, we propose to implement the correction recommended by Akras et al. [7], thereby obtaining what we will refer to as an equivalent temperature *T*E. These results were used to estimate a more solid spectral type determination of the red giant.

The paper is structured as follows: Section 2 provides a description of the observations and data reduction process. Section 3 presents the fitting with the MARCS model spectra, the determination of the equivalent temperature *TE*, the estimation of the spectral type for each symbiotic object, along with a concise overview of each object. In Section 4, we analyze the advantages and disadvantages of the spectral classification in different spectral bands. Finally, Section 5 ends the paper with our conclusions.

#### **2. Observations and Data Reduction**

On 8 June 2014, a selection of nine southern symbiotic stars was observed with the Ohio State InfraRed Imager/Spectrometer (OSIRIS), coupled to the 4.1 m telescope of the Southern Observatory for Astrophysical Research (SOAR, Cerro Pachón, Chile) under the proposal ID: SO2014A-009. The selected instrumental configuration consisted of a singleorder long-slit spectrometer and the camera f/7 to obtain medium-resolution (*R*∼3000) spectra in the H- and K-bands. The spectral coverage reached *λλ* 1.50–1.77 µm and *λλ* 2.0–2.4 µm for the H- and K-bands, respectively. The data of eight symbiotic stars obtained in the K-band were published by Marchiano et al. [13].

During all observations, several offset patterns following an ABBA cycle were taken, and the AB pairs were subtracted to remove the sky emission. Each spectrum was flatfielded, telluric-corrected, and wavelength-calibrated. The reduction process was carried out with the IRAF 1 software package.

Relevant information on the observed targets is given in Table 1. There, we list the star names, the stellar coordinates, the H and K magnitudes, and the IR classification type (D or S) [7,15], and in the last six columns, we show the spectral type of the cool component determined from near-IR spectra (using wavelengths between 7000 and 10,000 Å), the effective temperature *Tcool*, the references where we obtained these data, the spectral type of the cool component derived from the CO-bands at 2.3 µm by Marchiano et al. [13], the spectral type of the cool component obtained in this work, and the equivalent temperature *TE*, which we will explain in the next section.

**Table 1.** Stellar sample data with literature classifications.

References: 1. Schulte-Ladbeck ([16]), 2. Mürset and Schmidt ([17]), 3. Gałan et al. ([18]), 4. Zamanov et al. ([19]), 5. Pujol et al. ([20]).

#### **3. Results**

The H-band spectral observation of each symbiotic system listed in Table 1 is presented in Figure 1. This figure shows the complete set of spectra, which were normalized and vertically shifted to ease inspection and comparison. The spectra are first displayed according to the presence of CO-bands in emission and then according to the strength of the atomic absorption lines.

**Figure 1.** Comparison between the H-band spectra of our sample stars observed with SOAR/OSIRIS (black) and the best synthetic spectra computed with the MARCS model atmospheres (in different colors). Brackett series are indicated in gray, emission/absorption of CO in red, Mg I in green, and Si I in light blue. Ca I lines (1.6150 µm, 1.6157 µm, and 1.6197 µm) are indicated with a pink band and Al I lines (1.67504 µm, 1.67505 µm, 1.67189 µm and 1.67633 µm) with a yellow band. The stellar parameters of each model are also provided.

Our observations showed that all lines in BI Cru were in emission, and the remaining objects all showed CO-bands in absorption. More accurately, carbon monoxide molecules were observed in strong absorption in all stars, confirming the typical atmosphere of cool giants [21].

As a special case, the same as what was seen in the K-band [13], was RS Oph, which showed a weaker CO-band intensity than expected for the spectral type of the cool component. This weakening might be a consequence of the presence of dust with a temperature of ∼1000 K. Such a warm, dusty envelope veils the stellar photospheres and contributes to the total SED with emission of the absorbed and scattered stellar light at wavelengths longer than 2 µm [22].

Brackett series emission was significant in the spectra of BI Cru and Hen 3-1341. In addition to the lines of the hydrogen series, we also identified lines from Mg I, Ca I, Si I, and Al I. The Mg I line at *λ* 1.7109 µm was present in eight of the nine stars of our sample in absorption and in emission only in BI Cru. Note that the Si I line at *λ* 1.5888 µm was blended with a line of the Brackett series; it is important to take this into account in future observations with higher resolution to better distinguish the contributions of each line.

We fit the observations with the MARCS model spectra [14] to obtain the effective temperature (*T*eff). A grid was computed for M giant stars, covering the temperature range 2500 K ≤ *T*eff ≤ 4600 K, with a step of 100 K. Models were generated with a power resolution of *R* = 3000, solar metallicity [M/H] = 0, *v* sin*i* = 8 km s−1 , and 30 km s−1 and log *g* = 3.0, 2.5, and 2.0. For temperatures lower than 3500 K, the microturbulence was taken with a value of 4 km s−1 and the macroturbulence equal to 8 km s−1 . For a hotter grid (temperatures higher than 3500 K), lower values for microturbulence (1 km s−1 ) and macroturbulence (3 km s−1 ) were considered, based on typical values adopted by studies of red giants and symbiotic stars [23,24]. On the one hand, it is important to note that the projected rotational velocity *v* sin*i* = 8 km s−1 , or 30 km s−1 , taken as a fixed parameter, does not represent the majority of late-type giants in a symbiotic binary system ([25–27]). On the other hand, due to the instrument's spectral resolution used, the minimum detectable velocity would be approximately 100 km s−1 . Therefore, this particular instrument would not resolve any velocity change below 100 km s−1 .

In order to achieve a more accurate fit, we took into account not only the intensity of the CO absorptions, but also the relative relationships between these lines and other absorptions observed in each spectrum. The estimated values of *T*eff and log *g* of each symbiotic star are shown on the right side of Figure 1. The error of *T*eff was estimated to be around 100–200 K. Temperatures *T*E between parentheses refer to the correction made to the effective temperature of red giant stars with *T*eff < 4000 K according to the polynomial functions applied in the work of Akras et al. [7] and van Belle et al. [28]. We used this correction to the *T*eff values derived from the MARCS model to enable their comparison or equivalence with temperatures calculated via the blackbody approximation. Once this correction is obtained, which we refer to as the equivalent temperature *T*E for convenience, we applied the empirical relationships for red giants from the work of van Belle et al. [28] to determine the corresponding spectral type.

A brief summary of each individual symbiotic system is provided below with a description of the emission and absorption characteristics observed in this sample. It is worth clarifying that all systems belong to the catalog of symbiotic stars published by Belczynski et al. [15], which was used for the selection of the nine observed objects. Furthermore, there exist more recent catalogs of symbiotic stars, such as those compiled by Akras et al. [7] and Merc et al. [29].

**BI Cru (= Hen 3-782 = WRAY 15-967):**

BI Cru is a D-type symbiotic object. It is one out of three currently known symbiotic stars in which CO molecular emission was detected in the near-infrared at 2.3 µm [30], which seems to be stable over long time intervals. From the modeling of the CO-band emission detected in a high-resolution (*R*∼45,000) K-band spectrum, Marchiano et al. [31] derived the kinematics and the physical properties of the molecular circumstellar medium. This CO emission is associated with a dense disk, which would form after one or more episodes of intense mass loss from the red giant.

McCollum et al. [32] reported an IR shell of BI Cru more than five-times larger in arc size than the star's optical lobe. The temperature of this IR dust emission associated with our object was estimated to be 1300 K [33,34].

#### **RS Oph (= HD 162214 = MWC 414):**

RS Ophiuchi is a symbiotic recurrent nova consisting of a white dwarf with a mass near the Chandrasekhar limit that orbits inside the outer wind of a red giant. The system has had numerous outbursts recorded in 1898, 1933, 1958, 1967, 1985, and 2006, and more recently, in August 2021, it underwent its seventh optical eruption. Brandi et al. [35] redetermined the spectroscopic orbit of this star based on the optical spectra over the decade of 1998–2008.

The spectral type of the red giant and its temperature were estimated by Ribeiro et al. [36] and Pavlenko et al. [37]. In Table 1, we mark the reference of Zamanov et al. [19] because these authors applied a more precise method to classify the star. However, the spectral types obtained both by Marchiano et al. [13] and through this work (see Table 1) indicated an earlier spectral type than Zamanov et al. [19].

#### **Hen 3-1341 (= V2523 Oph = SS73 75):**

Hen 3-1341 is one of about ten symbiotics that shows hints of jets. According to Stute et al. [38], its optical and ultraviolet spectra show strong emissions in the Balmer continuum from N V, [Fe VII], and He II and the band at 6830 Å due to Raman scattering by neutral hydrogen. Furthermore, the object's infrared colors are appropriate for a cool giant (Spectral Type M4) without circumstellar dust. In quiescence, Hen 3-1341 resembles Z And, the prototype of symbiotic stars [39].

So far, no H-band observations have been reported. Our K-band observations [13] revealed an earlier spectral type for the red giant of this symbiotic star (see Table 1). From our fitting of the spectral features in the H-band, we found *T*eff = 3700 K and *T*E = 4170 K if we applied the temperature scale correction [7,28]. With this correction, we obtained the same spectral type found in the K-band. It is important to note that CO-bands in the spectrum of this star (see Figure 1) could be contaminated with emission lines. Furthermore, Hen 3-1341 presents emission lines that resemble those of BI Cru. New observations of this star in this band will be necessary to better confine its temperature.

#### **CL Sco (= Hen 3-1286):**

Infrared radial velocities were used by Fekel et al. [40] to compute the orbital parameters of this star. They obtained physical parameters adopting a period *P* = 625 days for a circular solution, which agrees with Kenyon and Webbink [41]. Besides, the orbital solution for both components was presented by Montané et al. [42].

In the near-IR, Gałan et al. [18] measured the photospheric chemical abundances of this object from high-resolution (*R*∼50,000) spectra. They presented the observed Kand H-band spectra of CL Sco together with synthetic spectra calculated using their own abundance estimates. According to our fitting in the H-band, the obtained *T*E value would indicate a later spectral type than that obtained by Marchiano et al. [13] from the analysis of the K-band; however, this coincides with that obtained by Gałan et al. [18].

#### **SY Mus (= HD 100336 = Hen 3-667):**

Dumm et al. [43] obtained the orbital parameters of the M star in this eclipsing symbiotic system and observed an asymmetry in the UV continuum light curve, which they explained as being caused by an asymmetric distribution of the wind from the cool component.

SY Mus was observed at high resolution (*R*∼50,000) in the HA-, K-, and Kr-bands using the Phoenix cryogenic echelle spectrometer on the Gemini South telescope byMikołajewska et al. [23]. All the spectra cover a narrow spectral range of ∼100 Å. Spectrum synthesis employing standard local thermal equilibrium analysis and atmosphere models was used to perform an analysis of the photospheric abundance of CNO for the red giant component. They obtained an isotopic ratio of 12C/ 13C ∼ 6 − 10, which indicates that the giant has experienced the first dredge-up.

In this work, we found that SY Mus has an effective temperature *T*eff = 2600 K. Applying the temperature scale correction, we obtained *T*E = 3490 K. This result is close to that published by Mürset and Schmid [17] and Gałan et al. [18], who indicated that the cold companion has a spectral type M4.5. However, this differs from what was obtained by Marchiano et al. [13], who found an earlier spectral type (see Table 1).

#### **RT Cru (= HV 1245):**

RT Crucis was classified as a symbiotic star by Cieslinski et al. [44], and the same authors reported the presence of rapid variations in its brightness. Gromadzki et al. [45] using optical and IR observations, detected two periodicities in the light curves: *Po* = 325 ± 9 days and *Pp* = 63 ± 1 days, which corresponded to the orbital and pulsation periods of the red giant, respectively.

RT Cru also attracted attention due to the discovery of hard X-ray emission detected in 2003–2004 with INTEGRAL/IBIS, which promoted the study of how accumulated material advances through a disk and reaches the surface of the white dwarf [46].

Using the first-overtone band of CO in the K-band, Marchiano et al. [13] found that the spectral type of the red giant companion was earlier than the one published by Pujol et al. [20]. A result in agreement with Marchiano et al. [13] was obtained in this work using synthetic spectra computed by MARCS model atmospheres in the H-band, as can be seen with the best fit in Figure 1 and in Table 1.

#### **V347 Nor (= Hen 2-147 = WRAY 16-208):**

This is a symbiotic Mira star. It has a pulsation period of 373 days and lacks the obscuration events typical of many other D-type symbiotics [47]. Its extended nebula was first modeled as an expanding ring inclined with respect to the plane of the sky by Corradi et al. [48]. Utilizing Doppler shift measurements from VLT integral-field spectroscopy in combination with the HST images, Santander-García et al. [47] demonstrated that the intrinsic geometry of the nebula is indeed that of a circular, knotty ring of ionized gas. This ring was found to be inclined by 68 degrees with respect to the line of sight and to expand with a velocity of ∼90 km s−1 .

We present here the first H-band spectrum for this symbiotic system. Figure 1 also shows that the best-fitting spectrum gave a spectral type not only earlier than the one found by Mürset and Schmid [17], but also earlier than the one found by Marchiano et al. [13] (see Table 1).

#### **KX TrA (= Hen 3-1242):**

KX TrA is a high-excitation S-type symbiotic binary. Its optical and near-IR spectra were studied by Ferrer et al. [49] and Marchiano et al. [50]. With observations taken between 1995 and 2007, they determined an orbital solution to the absorption lines of the red giant through the radial velocity curve. Since 2004, the system has entered into a period of activity.

Gałan et al. [18,51] presented high-resolution (*R*∼50,000) near-IR spectra for a sample of symbiotic systems. They employed stellar atmosphere models using a standard local thermal equilibrium analysis and derived the chemical abundance of several systems including KX TrA. Adopting the spectral type shown in the seventh column in Table 1, their analysis revealed a slightly subsolar metallicity ([Fe/H] ∼ −0.3) for KX TrA. Based on our analysis (see Figure 1 and Table 1), the spectral type was slightly later than the spectral type reported by Marchiano et al. [13] and was earlier than the one adopted by Gałan et al. [51] cited above.

#### **V694 Mon (= MWC 560):**

The most-spectacular features in the optical spectrum of V694 Mon are the broad absorption lines, most prominently at the Balmer transitions. V694 Mon is a symbiotic star in which the accretion disk drives a powerful high-velocity jet during outbursts, producing broad, blue-shifted, variable absorption lines from atomic transitions that extend up to thousands of km s−1 , not only in the optical range, but also in the IR, near-ultraviolet, and far-ultraviolet [52,53].

By fitting the observations of the donor star, we derived for the first time its effective temperature of *T*eff = 3000 K (see Figure 1). We also obtained, through the black body approximation, the *T*E value, which allowed us to estimate the spectral type of the cool companion (see Table 1).

#### **4. Discussion**

Analysis of the near-infrared range in symbiotic systems is key to their study since it has been used to classify these in S- and D-type stars based on their division into color–color diagrams, (J-H) versus (H-K), and (H-K) vs. (K-L) [54]. In addition, the techniques used in this region of the electromagnetic spectrum have proven to be the most-suitable for finding the spectral type of the late-type component in these binary systems. Moreover, the near-infrared is typically not contaminated by the nebula and the emission from the hot component, which both strongly affect the spectra at optical wavelengths.

Keenan and Hynek [55] introduced the spectral classification for the cool component of symbiotic binaries using the red TiO-bands, which increase in strength with decreasing temperature. Sharpless [56] solved the spectral classification using the characteristics of the most-important molecular bands between 7500 and 8900 Å. His temperature classification was also based on the growth of the TiO- and VO-bands, but he found that the CN-bands at 7916, 7941, 7878, and 8068 Å are useful luminosity discriminators in this spectral range as well. According to Schulte-Ladbeck [16], numerous observations of symbiotic systems are subject to contamination by emission lines from the nebula. This represents a source of error because these lines weaken the absorption of molecular bands in relation to the continuum. The longer wavelength bands are expected to be less affected since the IR brightness of the red giant increases rapidly in that range.

The K-band has been widely employed in most spectroscopic studies of cool or obscured objects, due to their higher brightness in this band compared to the J- or Hbands. Moreover, for the spectral classification of cool stars, several atomic features of Mg I, Ca I, and Na I are utilized, along with the band heads of the first-overtone CO-bands, which dominate the K-band spectra. Nevertheless, observations have revealed that stars of spectral types K3-M5 also exhibit significant temperature and luminosity-sensitive features in the H-band, including the lines of Mg I, Al I, and OH, and the band heads of the CO second overtone bands [57].

The circumstellar dust in symbiotic systems (more noticeable in D-type than S-type stars) often causes a continuum emission excess in the IR, which can reach down to wavelengths of around 2 µm. If this is the case, the excess emission can significantly hinder the extraction of precise information about the photosphere of the cool component. In such systems, shorter-wavelength spectra are required to identify and characterize the underlying star.

After considering the pros and cons of performing spectral classification in different near-infrared spectral bands, we propose that the H-band can provide valuable complementary information that can help to better constrain the parameters in symbiotic systems, in particular the temperature of the red giants in dusty symbiotic systems whose K-band spectra might be contaminated by the dust emission. It is also important to clarify that the differentiation of the spectral types obtained in the H- and K-band for each symbiotic object is based on observations that have been taken with the same instrument and within the same observing run. The studies themselves, however, employed different methods. The spectral classification found from the K-band spectra is based on measurements of equivalent line widths [13], unlike in this study, where we used a more robust method. Thus, an estimate of the similarities or differences between the spectral types obtained in each band is not entirely accurate if a different methodology is used, but we are confident that the presented analysis for the H-band spectra of symbiotic objects provides more decent results regarding the classification of their cool components.

#### **5. Conclusions**

To summarize, the IR spectra of symbiotic stars are a very useful tool to perform spectral classification of the cool components of these systems [16]. They are also excellent sources of information on the physical conditions of the disks from which molecular emission (if detected) originates. On the other hand, for many years, different studies, dedicated to the near-IR spectra of symbiotic stars have utilized low-resolution data, and in general, many of those studies were restricted to a small sample of stars, some examples being Schulte-Ladbeck [16], Kenyon and Gallagher [58], Schild et al. [59]. In the last two decades, new instruments have given access to the IR spectral range and have made it possible to obtain high-quality spectra with the necessary spectral resolution to improve the understanding of these complex objects.

With the aim to enlarge the sample of symbiotic objects with near-IR classifications, we presented the medium-resolution H-band spectra of nine symbiotic systems obtained with the OSIRIS spectrograph. We summarize our results in two points:

**Author Contributions:** Conceptualization, P.E.M., M.L.A., M.K. (Michaela Kraus) and L.S.C.; methodology, P.E.M., M.L.A, M.K. (Michaela Kraus), A.F.T. and L.S.C.; software, P.E.M., M.K. (Michalis Kourniotis); investigation P.E.M., M.L.A, M.K. (Michaela Kraus), M.K. (Michalis Kourniotis), A.F.T., L.S.C. and M.B.F.; resources, P.E.M., M.L.A, M.K. (Michaela Kraus), M.K. (Michalis Kourniotis), A.F.T., L.S.C. and M.B.F.; data curation, M.L.A.; writing—original draft preparation, P.E.M.; writing—review and editing, P.E.M., M.L.A, M.K. (Michaela Kraus), M.K. (Michalis Kourniotis), A.F.T., L.S.C. and M.B.F.; funding acquisition, A.F.T., M.L.A, M.K. (Michaela Kraus) and M.B.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** A.F.T. and M.L.A. acknowledge financial support from the Universidad Nacional de La Plata (Programa de Incentivos 11/G160) and CONICET (PIP 1337), Argentina. M.K. (Michaela Kraus) acknowledges financial support from the Czech Science Foundation (GACR, Grant Number ˇ 20-00150S). The Astronomical Institute of the Czech Academy of Sciences Ondˇrejov is supported by the project RVO:67985815. This project has received funding from the European Union's Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowska-Curie Grant Agreement No. 823734. M.B.F. acknowledges financial support from the National Council for Scientific and Technological Development (CNPq) Brazil (grant number: 307711/2022-6)

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data underlying this article will be shared upon reasonable request to the corresponding author.

**Acknowledgments:** This research made use of the NASA Astrophysics Data System (ADS) and of the SIMBAD database, operated at CDS, Strasbourg, France. This paper is based on observations obtained at the Southern Astrophysical Research (SOAR) telescope, which is a joint project of the Ministério da Ciência, Tecnologia, e Inovação (MCTI) da República Federativa do Brasil, the U.S. National Optical Astronomy Observatory (NOAO), the University of North Carolina at Chapel Hill (UNC), and Michigan State University (MSU) under Program ID SO2014A-009.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notes**

1 IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

#### **References**

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

### *Article* **Revealing the Binarity of HD 36030—One of the Hottest Flare Stars**

**Olga Maryeva 1,\* , Péter Németh 1,2 and Sergey Karpov 3,4**

**Abstract:** The Kepler and TESS space missions significantly expanded our knowledge of what types of stars display flaring activity by recording a vast amount of super-flares from solar-like stars, as well as detecting flares from hotter stars of A-F spectral types. Currently, we know that flaring occurs in the stars as hot as B-type ones. However, the structures of atmospheres of hot B-A stars crucially differ from the ones of late types, and thus the occurrence of flaring in B-A type stars requires some extension of our theoretical views of flare formation and therefore a detailed study of individual objects. Here we present the results of our spectral and photometric study of HD 36030, which is a B9 V star with flares detected by the *TESS* satellite. The spectra we acquired suggest that the star is in a binary system with a low-mass secondary component, but the light curve lacks any signs of periodic variability related to orbital motion or surface magnetic fields. Because of that, we argue that the flares originate due to magnetic interaction between the components of the system.

**Keywords:** stars flares; stars activity; stars binaries spectroscopic; early-type stars; stars individual; HD 36030

**Citation:** Maryeva, O.; Németh, P.; Karpov, S. Revealing the Binarity of HD 36030—One of the Hottest Flare Stars. *Galaxies* **2023**, *11*, 55. https://doi.org/10.3390/ galaxies11020055

Academic Editors: Lydia Sonia Cidale, Michaela Kraus and María Laura Arias

Received: 14 March 2023 Revised: 1 April 2023 Accepted: 10 April 2023 Published: 12 April 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Stars classified as B-type span a very wide range of temperatures (from 25,000 to 11,000 K) and masses (from 50 M for B0 Ia0 down to 3.5 M for B9 V). They usually display the variability related to changes of the photospheric structure, e.g., pulsations of different varieties [1], such as: slowly pulsating B-type stars (SPB) [2,3], variables of the *β* Cephei type [4], variables of the *α* Cygni type [5], or classical Be (CBe) stars [6]. Less frequently, B-type stars show variability related to extended atmospheres and outflows, usually observed in peculiar B-type stars, such as B[e] [7,8] or luminous blue variables (LBVs) [9,10]. Quite unexpectedly, some B stars have also been shown to manifest very energetic stellar flares (e.g., [11] and references therein)—unpredictable dramatic increases in brightness for a few dozen minutes.

Such flaring activity is a well-known effect that is often seen both on the Sun and on socalled flare stars. Energetic flaring events are attributed to the release of the energy stored in coronal magnetic fields in magnetic reconnection events that accelerate the particles (electrons and ions) to high energies, thus heating dense regions of the stellar atmosphere and generating the flare-like emission in a wide range of frequencies, from radio to gammarays [12,13].

Although the majority of flare stars were found among objects of late spectral types, primarily cool K-M type red dwarfs [14,15], such outbursts were also detected in some hot B–A type stars, both in the optical range [16–18] and in X-rays [19,20]. In cool stars (F5 and later), the flares are linked with strong magnetic fields produced by a dynamo mechanism working in sufficiently deep outer convection zones [21]. On the other hand, the interiors of hot early-type stars are mostly radiative, with no convective motions to power the dynamo, and therefore magnetic field production is not expected in them, and other explanations of the eruptive activity are needed.

Pedersen et al. [22] suggested that all flaring hot stars are parts of binary systems and that the flares occur on unresolved cool components. However, Švanda and Karlický [23] and Bai and Esamdin [24] found that the frequency distribution of flare energies for A-type stars is steeper (with more flares having large energies) than the one for stars of later spectral types. It may suggest that the nature of flares on hot stars is indeed different, and they are not occurring on their colder companions. Balona [17,18] studied the variability of flaring B-A type stars and argued that rotational modulation seen in them confirms the presence of strong surface magnetic fields. One more possible mechanism for the flares may be related to magnetic interaction between the components in a close binary system, as was suggested by Yanagida et al. [20] for the pre-main sequence stars HD 47777 and HD 261902.

Currently, more than one hundred hot flaring B-A type stars are known from observations of the *Kepler* and *TESS* missions [11,17,18]. Van Doorsselaere et al. [25] demonstrated that about 2.45% of all B-A stars display detectable flares. However, the nature of their flares is still elusive. Thus, the study of these stars, and their possible binarity, is of utmost importance.

This paper is devoted to the spectral and photometric study of HD 360301 , classified as B9 V [26] and belonging to the galactic open cluster ASCC 21 (or [KPR2005] 21 [27,28]). For the first time, HD 36030 was mentioned as a flare star in Maryeva et al. [11], where the authors searched for flare stars among confirmed members of galactic open clusters using high-cadence photometry from the *TESS* mission. Maryeva et al. [11] found eight B-A type flaring objects; among them, HD 36030 was the hottest. Independently, Balona [18] and Yang et al. [29] also selected HD 36030 as a flare star while searching for flares in *TESS* data. Despite being included in a couple dozens of catalogs, and being detected as an X-ray source by the Swift X-ray Telescope (XRT) [30]), to date HD 36030 lacks any specific studies devoted to it.

In Section 2, we present new spectral data, their analysis, and the results of atmospheric modeling using the TLUSTY code. In Section 3, we perform the analysis of photometric data on different time scales, and in Section 4, we discuss the results, the parameters of the second component, and the nature of flares. Finally, short conclusions are given in Section 5.

#### **2. Spectroscopy of HD 36030**