TYC 3340-2437-1: A Quadruple System with a Massive Star

Hierarchical massive quadruple systems are ideal laboratories for examining the theories of star formation, dynamical evolution, and stellar evolution. The successive mergers of hierarchical quadruple systems might explain the mass gap between neutron stars and black holes. Looking for light curves of O-type binaries identified by LAMOST, we find a (2+2) quadruple system: TYC 3340-2437-1, located in the stellar bow-shock nebula (SBN). It has a probability of over 99.99% being a quadruple system derived from the surface density of the vicinity stars. Its inner orbital periods are 3.390602(89) days and 2.4378(16) days, respectively, and the total mass is about (11.47 + 5.79) + (5.2 + 2.02) = 24.48 M ⊙. The line-of-sight inclinations of the inner binaries, B1 and B2, are 55.°94 and 78.°2, respectively, indicating that they are not coplanar. Based on observations spanning 34 months and the significance of the astrometric excess noise (D > 2) in Gaia Data Release 3 data, we guess that its outer orbital period might be a few years. If it were true, the quadruple system might form through the disk fragmentation mechanism with outer eccentric greater than zero. This eccentricity could be the cause of both the arc-like feature of the SBN and the noncoplanarity of the inner orbit. The outer orbital period and outer eccentric could be determined with the release of future epoch astrometric data of Gaia.


INTRODUCTION
Massive stars are progenitors of black holes (BHs) and neutron stars (NSs).Despite more and more BHs and NSs being found by the observations of LIGO, Virgo, Xray telescopes, and radio telescopes, there still is a mass gap between NSs and BHs in the range of about 2 − 5 M ⊙ (e.g.Farr et al. 2011).It is ambiguous whether the gap is caused by observation bias or as a result of massive stellar evolution (Abbott et al. 2019;Fryer et al. 2012).In general, the evolution of a massive OB star Corresponding author: Jiao Li; Chao Liu lijiao@bao.ac.cn; liuchao@nao.cas.cn with a mass of 8-25 M ⊙ is thought to possibly experience a hydrogen-rich type IIP supernovae and leave behind a neutron star (Smartt 2009;Vink 2021).However, most of the massive stars are in binaries, and over 21% of the massive stars probably lead to a merger, which could deeply influence their evolution route and fate (Sana et al. 2012).A massive wide binary with a period longer than 10 years would also produce a triple-ring structure similar to that of SN 1987A (Morris & Podsiadlowski 2007).
If massive stars were members of quadruple systems, their evolution would be more exotic.The successive mergers of hierarchical quadruple systems might explain the mass gap between NSs and BHs (Hamers et al. 2021).In the 3+1 quadruple system, the mass gap probably arises if two NSs in the innermost orbit merge due to secular evolution forming a mass gap object, and becoming a triple system.The resulting inner binary may successively merge due to continued secular evolution, leading to the formation of a BH outside the mass gap (Safarzadeh et al. 2020).In the 2+2 quadruple system, the merger remnant of an inner binary probably receives a recoil kick from the anisotropic emission of gravity waves, triggering a dynamical interaction where another merge would happen (for details, see Fragione et al. 2020;Hamers et al. 2021).They are ideal laboratories to study stellar formation and evolution in a dynamically complex environment where the constituent binary stars can be close enough to interact with each other.It can also be easy to imagine that, at a specific condition for 2+2 quadruple, the less massive component binary might achieve a certain velocity to be a runaway binary if the more massive component merged and exploded as a supernova (Gao et al. 2019).
Currently, there are three main categories of formation theories for multiple systems: turbulent (core and filament) fragmentation, disk fragmentation, or dynamic interaction.The turbulent fragmentation model presumes that multiple systems generate from overdensities that develop and collapse within molecular filaments, producing initially widely separated systems (> 500 au).The disk fragmentation model posits that multiple systems arise from a massive accretion disk, forming initially smaller separations (< 500 au), and similar orientations.The dynamic interaction mode can establish gravitational bonds among nearly independently single stars at large distances through energy and angular momentum exchange with the surrounding gas cloud or individual circumstellar disks.(Offner et al. 2023 and references therein).So, the semi-major axis of the outer orbit in quadruple systems and the orientations of their inner binaries could serve as possible indicators for discerning the formation mechanisms of multiple systems, particularly in distinguishing between the turbulent and disk fragmentation models (Lorenzo et al. 2017).Moreover, the dynamic interaction mode possesses the capability to rearrange the hierarchy and multiplicity of systems formed via turbulent or disk fragmentation.
Despite the fraction of quadruple system being about 38% among O-type stars (Moe & Di Stefano 2017), only a handful of massive quadruples are resolved with orbital parameters (e.g.Lorenzo et al. 2017;Southworth 2022;Nazé et al. 2022).Detection and characterization of a quadruple can be challenging.It is very difficult to detect eclipse light curves, ellipsoidal light curves, or radial velocity variations for both inner binaries of a quadruple, especially, for the quadruple that inner binaries have distinct orbital inclinations.
We discover a quadruple system by cross-matching Otype stars found by the Large Sky Area Multi-Object Spectroscopic Telescope survey (LAMOST: Liu et al. 2019;Li 2021) with the light curve of Transiting Exoplanet Survey Satellite (TESS: Ricker et al. 2015): TYC 3340-2437-1.The system has two inner binaries.One is composed of O+B type stars, and the other is composed of B+A type stars.Its distance is about 2.31 kpc (see Section 5).The variation of the radial velocity of LAMOST Medium Resolution Survey ( LAMOST-MRS: Liu et al. 2020;Zhang et al. 2021) is larger than 200 kms −1 .TYC 3340-2437-1 has been identified as G151.9227-00.5651,located within the stellar bow-shock nebula (SBN), by Kobulnicky et al. (2016).SBN is an arcuate structure created by the interaction between the stellar wind and the surrounding interstellar medium (ISM).TYC 3340-2437-1 would be an exotic object to study the evolution of massive stars.

TESS LIGHT CURVE
TYC 3340-2437-1 was observed by TESS in Sector 19 with a cadence of around 30 minutes.The light curve (hereafter, we denote LC as the TESS light curve of the object) has been extracted by using the MIT Quick Look Pipeline (Huang et al. 2020), as shown in panel (a) of Figure 1.
A Lomb-Scargle periodogram of LC shows a domination peak at 1.69759(59) day (detail see Section A.1).It is similar to an ellipsoidal light curve by folding LC into twice this peak period of 3.39518 days.In addition, we have observed ten narrow valleys in LC, which are represented as blue dots with blue triangles in Figure 1 (a).These valleys appear to have arisen due to the eclipsing of a binary system.We fit each of the narrow valleys using the Gaussian function and obtain the conjunction points of the eclipsing binary, shown as the blue triangles of Figure 1 (a).Then, fitting the times of the conjunction points with a linear function, where i = [0, 1, 2...10], t i is the ith conjunction time.We obtain a period P 2 = 2.43781(16) days and a transit time T c2 = 2458817.49056(75)BJD for the eclipsing binary.
After folding LC with P = 2.43781 days and T c = 2458817.49056BJD (see Figure 1 b), and excluding the data (the blue dots in Figure 1 b) in phases of (−0.04 < ϕ < 0.04) or (0.46 < ϕ < 0.54), a new light curve (LC 1 ) without eclipses is obtained (see Figure 1 c).LC 1 is folded into a period of P 1 = 3.390497 days, which  (Huang et al. 2020); the blue triangles are the minimum points of the narrow "valley"; the blue dots are the points in the phase of (−0.04 < ϕ < 0.04) or (0.46 < ϕ < 0.54) when the LC is folded with a period of 2.43781 days.Panel (b), the phase is folded with a period of 2.43781 days and to the transit time of 2458817.49056BJD.Panel (c), the LC1 shows the removal of the blue points in the Panel (a) is derived by fitting the radial velocity curve (see Section 3.1).The resulting light curve exhibits a clear sinelike feature, which is indicative of the ellipsoidal effect (as shown in the middle panel of Figure 2 b).By dividing LC by a model ellipsoidal light curve (or a smoothed light curve of LC 1 ), the eclipsing light curve (LC 2 ) can be obtained (see Figure 2 c).
It is highly unlikely that the periods of LC 1 and LC 2 originate from a triple system since building a stable triple system with these two periods is extremely difficult (Georgakarakos 2008;Harrington 1975).Instead, they may be contributed to a binary system in which rotation and binary orbit are asynchronous.In this case, the periods of the ellipsoidal effect and radial velocities should be 3.390497 and 2.43781 days, respectively.Alternatively, they may have arisen from two isolated binary systems or a 2+2 quadruple system.This issue can be figured out by the radial velocity curve of the object.In the following, we will see that LC 1 is generated from the ellipsoidal effect of an inner binary with O+B type stars (referred to as binary B 1 hereafter) of a quadruple system, while LC 2 is induced by the other inner binary consisting of B+A type stars (referred to as binary B 2 hereafter).We define their orbital periods as P 1 and P 2 .The schematic configuration of TYC 3340-2437-1 is shown in Figure 3.

SPECTRA AND RADIAL VELOCITY CURVE
The spectra of TYC 3340-2437-1 have been observed through LAMOST low-resolution survey (LAMOST-LRS) and the binarity and exotic star project of LAMOST-MRS (Li et al. 2023).LAMOST is a reflecting Schmidt telescope with an effective aperture of about 4 m, located at the Xinglong Station of the National Astronomical Observatories of the Chinese Academy of Sciences (Cui et al. 2012;Zhao et al. 2012).However, the radial velocities measured by the spectra were inadequate to derive a reliable period when we found the star interesting, and could not cover all phases of the orbital period P 1 (as shown by the star and triangle markers in the top panel of Figure 2 b ).To address this issue, follow-up spectra were observed by using the BFOSC E9+G10 instrument of the Xinglong 2.16-m telescope (Zhao et al. 2018) and the YFOSC E9+G10 instrument of the Lijiang 2.4-m telescope (observation and data reduction see Section A.2).

Radial Velocity of Star B 1a
The radial velocities (see Table A.2) are measured using cross-correlation function (CCF) with OSTAR2002 grid of TLUSTY (Lanz & Hubeny 2003), implemented through the package of laspec (Zhang et al. 2020(Zhang et al. , 2021)).We subsequently employed a Markov Chain Monte Carlo (MCMC) approach to determine the orbital parameters of the binary with these radial velocities.The posterior distribution is where v n and σ vn are the radial velocities and errors;   (Huang et al. 2020); the red lines are light curve calculated by ellc (Maxted 2016).Panel (b) shows that the ellipsoidal light curve of binary B1 (LC1) and corresponding radial velocity curve, both of which are folded by a period of 3.390497 days.In the top panel, blue solid and dash lines are Kepler orbital radial velocities of binary B1a and binary B 1b , respectively; the markers of the open square, dot, and star are radial velocities measured from spectra of BFOSC, YFOSC, and LAMOST-LRS (R ∼ 1800), respectively.The open and solid triangles stand for radial velocities of star B 1b to star B1a, which are measured from spectra of LAMOST-MRS (R ∼ 7500).The red triangles stand for the phases where the spectra are used to measure the radius ratio of star B 1b to star B1a.In the middle panel, the gray dots donate the light curve removed the data near the eclipses of binary B2.The bottom panel shows the folded residual of the TESS LC.Panel (c) shows the relative eclipse light curve of binary B2 (LC2) and folded by a period of 2.43781 days, which is obtained by dividing the TESS LC into the model light curve of binary B1.The remarkable deviations at the phase of ∼ 0.4 of LC2 might be attributed to the de-trending process of the TESS LC, which can be seen at the edges of Panel (a).Note: the residuals of Panel (b) and (c) are the same residual of Panel (a) but folded with different periods.
ing errors such as errors of the instrument, wavelength calibration, and even radial velocity model.
Our analysis yielded a period of 3.390497(30) days, which is consistent with the period (3.39518 days) derived from the periodogram of the light curve (see Section 2).The period of 3.390497(30) is adopted as P 1 because it is more precise, thanks to the radial velocities collected over a longer observation time.The other pa-rameters can be found in Table 1.Upon folding the radial velocity curve and LC 1 with the period of 3.390497 days, we find that the folded LC 1 exhibits two peaks, with the peak near phase 0.25 being slightly higher than the one near phase -0.25 (see Figure 2 b), due to the Doppler boosting effect.These provide evidence that LC 1 arises from the ellipsoidal effect of a binary with an orbital period of 3.390497 days.In other words, LC 1 and LC 2 are generated by two distinct binary systems.
If the components of binary B 1 were bright enough and the spectrum resolution was greater than 2000, we would expect to observe double peaks in the CCF profiles when the spectra were observed at an orbital phase close to 0.25 or -0.25, given that the semi-amplitude of the radial velocity K 1a is larger than 100 km s −1 .However, upon examining the CCF profiles of all spectra, we only found a single peak.Therefore, we use the LAMOST-LRS spectrum to measure the atmosphere parameters of star B 1a and consider the contamination of the spectrum by the other three stars (details see Section A.3), obtaining: T eff1a = 30279 ± 1600 K, log g 1a = 3.74 ± 0.25 dex, Z 1a ∼ 0.35 Z ⊙ , and v sin i 1a = 111±42 kms −1 .The exemplary spectroscopic fit is shown in Figure 4

Radial Velocity of Star B 1b and Radius Ratio
We then carefully examined the LAMOST-MRS spectra taken at the orbital phases of -0.187 and 0.245 and discovered very faint lines that were moving in antiphase with main absorption lines (such as He I 4922, Hα, and He I 6678, as shown by the blue lines in Figure 5), which means the faint lines arise from the other component star B 1b .Furthermore, these lines displayed a higher velocity shift, indicating that star B 1b has less mass than star B 1a .
The star B 1b should be a B-type star because the HeI and Hα adsorption lines of B 1b can be seen, as shown in Figure 5.The uncertainties of atmospheric parameters measured from LAMOST-MRS for B-type stars are very large.Additionally, the depths of these absorption lines in star B 1b are comparable to the noise level, further complicating the determination of its atmospheric parameters.Fortunately, its radial velocity can still be determined using the LAMOST-MRS spectra.
Its lower mass compared to star B 1a implies that it is likely to be a main sequence star.Assuming specific atmosphere parameters, we could derive its radial velocity and radius ratio (r 1 = R 1b /R 1a ).For the sake of simplicity, we set the atmosphere parameters of B 1a as fixed values: T eff1a = 30279 K, log g 1a = 3.74 dex, Z 1a = 0.35 Z ⊙ .Since the radius ratio depends on atmosphere parameters of B 1b , we drive the radial velocity of B 1b and radius ratio by exploring a grid of T eff1b and log g 1b , where T eff 1b ranges from 19000 K to 29000 K with a step of 2000, and log g 1b ranges from 3.8 dex to 4.5 dex with a step of 0.1, Z 1b = 0.35 Z ⊙ .
Therefore, the spectrum of binary B 1 depends on radial velocities (v), the projected rotation velocities (v sin i) and radii (R) of its two components.It can be as a function of these parameters: where f 1a and f 1b are synthetic spectra of stars B 1a and B 1b specific atomashpere paramters, respectively; d is the distance of TYC 3340-2437-1.We use an MCMC approach to fit the normalized spectra of LAMOST-MRS in wavelength ranges of [6520, 6600] and [6650, 6700] with the normalized broadened synthetic spectrum with a resolution of 7500 (the posterior probability distribution is similar to Eq. A1).This enables us to derive the radial velocities of star B 1b and the radius ratio of star B 1b to B 1a across the grid of T eff1b and log g 1b .The radial velocities of B 1b cannot be measured when the spectra were observed near phase 0, as the top panel of      A.2.Moreover, the radial velocities were nearly independent of T eff1b and log g 1b .
In order to obtain an accurate radius ratio, we first selected 12 spectra that were observed near phases -0.25 and 0.25 (their locations in the folded radial velocity curve are indicated by the red triangles in the top panel of Figure 2 b).We then calculate the median values of the radius ratios of each point in T eff1b and log g 1b grid, with their errors estimated by the bootstrap method using the package of Scipy.The result is listed in Table 2.The radius ratio was found to decrease with an increase of log g 1b .Furthermore, the radius ratio was observed to decrease with an increase of T eff1b , except when log g 1b was equal to 4.4 dex and 4.5 dex, as shown in Figure 6.

Mass Ratio of Star B 1b to B 1a
Since we have measured the radial velocities of B 1a and B 1b at various time points, the mass ratio q 1 = m 1b m −1 1a can be derived.We fit the radial velocities using the MCMC approach.The posterior distribution is where v yn is the measured radial velocity of star B 1a or We have determined K 1a = 111.3± 20.2 km s −1 , K 1b = 201.3± 1.8 km s −1 , and calculated q = 0.55 ± 0.10.The value of K 1a is consistent with the result of 105.4 ± 1.3 km s −1 obtained in Section 3.1 within the error range, although it is slightly higher.

EVIDENCE OF QUADRUPLICITY
Since the pixel size of TESS is 21 ′′ , the TESS light curve might be contaminated by its vicinity stars.So we extract a light curve from a single pixel from the full-frame image (FFI) centered at TYC 3340-2437-1 by Tesscut (Brasseur et al. 2019), the ellipsoidal and eclipse features still exist.
It's worth noting that we initially presumed TYC 3340-2437-1 to be a binary system (B 1 ), and the eclipsing light curve originated from a nearby eclipsing binary B 2 .We inferred the minimum luminosity of binary B 2 by comparing the eclipse depth of LC 2 with the G RP = 10.12 mag magnitude of TYC 3340-2437-1 (for more details, refer to Section A.4).
To validate this assumption, we searched for light curves of stars within a 60 ′′ angular radius in the Zwicky Transient Facility (ZTF: Masci et al. 2019;Bellm et al. 2019) Data Release (DR9).Specifically, we searched for eclipsing binaries with a period of P 2 = 2.43781 days within the magnitude range G RP ∈ [10.12, 15.12] mag.However, none was found.TYC 3340-2437-1, therefore, might be a quadruple system composed of binary B 1 and binary B 2 .
Moreover, binary B 2 should be fainter than star B 1b , since no signal of binary B 2 can be detected from the LAMOST-MRS spectra (see Figure 5).We assume that the distribution of nearby stars around TYC 3340-2437-1 is uniform in terms of solid angle.Counting from Gaia DR3, there are 1695 stars located in the radius of 0.5 • centered at TYC 3340-2437-1 with G RP from 10.12 to 15.12 mag, so the steradian density is ρ = 7.08 ± 0.17 × 10 6 sr −1 .The probability of optical double P ∼ ρπR 2 ang = 8.4 ± 0.2 × 10 −5 adopted effective angular resolution of Gaia R ang ∼ 0.4" = 1.94 × 10 −6 rad 1 .In other words, TYC 3340-2437-1 has more than 99.99% probability to be a 2+2 quadruple system.This is reinforced by the fact that the parameters of astrometric excess noise ϵ = 0.1 mas and significance of the noise D = 19.68 in Gaia DR3, because ϵ > 0 mas and D > 2 represent an additional intrinsic scatter term in the astrometric solution (Lindegren et al. 2012).In this case, it might be caused by orbital motion between binary B 1 and binary B 2 .

Characterizing the System
Combining the binary mass function of star B 1 and 1 https://www.cosmos.esa.int/web/gaia/dr2Since the variables on the right-hand side of Eq. ( 6) were obtained from LAMOST-LRS spectra and the radial velocity curve of binary B 1 (for details, see Section 3), and R 2 1a sin 3 i 1a remains constant.The characteristics of an ellipsoidal light curve are predominantly influenced by the mass ratio, orbit inclination, and the ratios of star radii to the semi-major axis (for details, see Section 3.2.1 of Prša 2018).For TYC 3340-2437-1, LC 1 is primarily a result of the ellipsoidal variations of star B 1a , given that R 1b < 0.5R 1a (see Figure 6) and q 1 ∼ 0.5.It depends mainly on q 1 , i 1 , R 1a , and the smi-major axis A 1 of binary B 1 .Notably, combining indicates that A 1 is solely determined by R 1a in this particular case.Therefore, LC 1 is expressed as a function of i 1a and R 1a .
Considering that the luminosity ratio of binary B 2 to B 1 is less than 0.02 (see Section A.4), it can be ignored that LC 1 is influenced by the light coming from binary B 2 .With the two variables and the two functions, we can constrain both i 1 and R 1a by a simultaneous fitting of the ellipsoidal curve and the radial velocity curve.This allows us to calculate the masses of stars B 1a and B 1b .
Based on the atmosphere parameters (T eff1a = 30249 K, log g 1a = 3.74 dex, and Z 1a = 0.35 Z ⊙ ) of star B 1a , it appears to be a main sequence star.Additionally, we do not detect any signal of star B 2a and star B 2b from our current spectra, indicating that their luminosities are less than that of star B 1b .So, we assume that star B 1b , B 2a , and B 2b are also main sequence stars.The luminosity of star B 1b can be obtained according to the relation between the luminosity and mass of the main sequence star, given a specific mass.We use isochrones of PAdova and Trieste Stellar Evolution Code (PARSEC) stellar model (Bressan et al. 2012;Chen et al. 2014) to interpolate the luminosity from the mass of star B 1b .The masses of star B 2a and star B 2b can then be interpolated based on their luminosities, which can be estimated by the relative eclipse depths of binary B 2 (see Section A.4 and Figure 2 (c) for clarity).
Consequently, the semi-major axis of binary B 2 can be calculated with its masses and period.Additionally, the eclipse widths of LC 2 are influenced by the radii of star B 2a and star B 2b (for details, refer to Section 3 of Prša 2018).Furthermore, their luminosities would increase as the inclination of binary B 2 (i 2 ) decreases.For example, given a luminosity ratio of star B 2b to star B 2a , if i 2 = 90°, the primary eclipse of the binary would be the deepest, and the least luminosity of the star is needed to reach the relative depth.By simultaneously fitting the eclipse widths of LC 2 with their masses and orbital period, we can determine their radii and inclination of binary B 2 .
We assumed that the other three stars have the same metallicity as star B 1a and fixed orbital periods of binary B 1 (3.390497 days) and binary B 2 (2.43781 days), the transit time of binary B 2 (2458817.49056BJD), and the effective temperature (30279 K) of star B 1a to simplify the analysis.Then an MCMC approach is performed to fit the TESS LC and radial velocity curves of binary B 1 with these assumptions.The posterior distribution can be obtained using these priors; for details, see Section A.5.The best-fit model is shown in Figure 2 (a).We obtain that the masses of TYC 3340-2437-1 are (11.47+0.28  −0.26 + 5.79 +0.13 −0.12 ) + (5.20 We noted that T eff2a , log g 2a and R 2a are similar to the corresponding parameters of star B 1b .However, despite these similarities, we did not detect any discernible signal in the LAMOST-MRS spectrum corresponding to B 2a .This lack of detection could be attributed to two reasons.Firstly, B 2a doesn't exhibit a noticeable anti-phase movement with respect to star B 1a , making it challenging to distinguish its spectral lines since the depths of the absorption lines of B 2a are comparable to (or lower than) the noise level in our observations.It is also plausible that the luminosity of B 2a is overestimated, although this cannot be conclusively determined with our spectra.
From the bottom panel of Figure 2 (a), we can see that there are some residuals with absolute values exceeding three times their errors, see bottom panel of Figure 2 (a).When we fold the residuals into orbital periods of binary B 1 and binary B 2 , we do not observe any period signals, as shown in the bottom panels of Figure 2 (b)  and (c).This indicates that these larger residuals are not generated by the orbital motions of binary B 1 or binary B 2 .There are several possible reasons for this, such as 1) the de-trending process of the light curve; 2) pulsations of the four components; 3) the orbital motion between binary B 1 and B 2 ; 4) underestimated errors of the light curve.The small errors in the derived parameters (see Figure A.6) may suggest the light curve errors were underestimated.However, we cannot figure out the exact cause of the large residual with the current observation data.

DISTANCE FROM SED FITTING
The SED fitting was performed following a strategy similar to the one used in the Python package Speedyfit.We prepared grids of fluxes passing through filters of Gaia, APASS, 2MASS, and allWISE using SEDs of TLUSTY OSTAR2002 (Lanz & Hubeny 2003), TLUSTY BSTAR2006 (Lanz & Hubeny 2007), and Kurucz (Castelli & Kurucz 2003).
Fixing the effective temperatures (30279, 19299, 18023 and 8140 K) and surface gravities (3.717, 4.271, 4.256 and 3.955 dex) of the four components of TYC 3340-2437-1 derived from the light curve and radial velocity curve fitting, we interpolated fluxes for stars B 1a , B 1b , B 2a , and B 2b from the flux grids of TLUSTY OSTAR2002, TLUSTY BSTAR2006, TLUSTY BSTAR2006, and Kurucz, respectively.We assumed that their radii followed Gaussian priors with means of 7.77, 2.92, 2.83, and 2.02, and standard deviations of 0.09, 0.12, 0.14, and 0.57 R ⊙ .We used an MCMC approach to fit the observed fluxes of Gaia EDR3 G BP -, G-, and G RP -bands, APASS B-, V -, G-, R-, and I-bands, 2MASS J-, H-and K s -bands, and ALLWISE W 1 -and W 2 -bands.The best-fit model is shown in Figure 7.
Our analysis yielded a distance of 2.31 ± 0.03 kpc and a color excess of E(B −V ) = 1.198±0.002mag for TYC 3340-2437-1.This distance is less than the one (d Gaia = 2.75 ± 0.14 kpc) derived from Gaia DR3 parallax.This discrepancy indicates that the astrometric solution of TYC 3340-2437-1 cannot be adequately explained by the single-star model again.
The angle distance between star B 1a and star B 1b (star B 2a and star B 2b ) is about 49.4 µas (29.7 µas).This angle distance is too small to have a significant impact on the astrometric solution for Gaia.Therefore, the solution of Gaia parallax might be influenced by the orbital motion between binary B 1 and binary B 2 .A.2) and obtain v Na = −18.5 +6.9 −11.8 km s −1 .Since we know the Galactic coordinates and the distance (d = 2.31 ± 0.03 kpc derived from SED fit) of the star, we calculate the radial velocity of the clouds where the star was born v ev ∼ −18.6 km s −1 with respect to the local standard of rest (LSR) by assuming that the Sun's motion with respect to the LSR is (11.1, 12.24, 7.25) km s −1 and the circular speed is 238 km s −1 (Schönrich et al. 2010;Schönrich 2012).It suggests that the comoving CO cloud of the star has a radial velocity similar to the ISM.
We download the 12 CO(J = 1 − 0) data around TYC 3340-2437-1 from the Milky Way Imaging Scroll Painting (MWISP) project (Su et al. 2019) and integrate emission of 12 CO(J = 1 − 0) from -30.3 to -11.7 km s −1 , which are the 16th and 84th percentile of the radial velocities of the Na D1 and Na D2 absorption lines.The radial velocity is also consistent with the systematic radial velocity of binary B 1 of -14.4 km s −1 .As shown in Figure 8, TYC 3340-2437-1 is located in a region with Table 3. Stellar, orbital, and astrometric parameters of TYC 3340-2437-1.The stellar and orbital parameters are derived using an MCMC approach to fit the radial velocity curve of binary B1 and the light curve of the quadruple system.Astrometric parameters are from Gaia EDR3 (Gaia Collaboration et al. 2022), with a zero-point correction applied to the parallax (Lindegren et al. 2021).Quantities shown are effective temperature T eff , surface gravity log g, stellar mass m, radius R, Roche lobe radius RL, luminosity log L, metallicity Z, the semi-amplitude of radial velocity K, the orbital period P , the zero-point of the ephemeris Tc, the systemic velocity of the binary vγ, the orbital inclination i, the smi-major axis A, the orbital eccentricity e, the mass ratio q, the longitude of periastron ω, the total mass of the quadruple system m total , parallax ϖ, the Renormalised Unit Weight Error RUWE, the astrometric excess noise ϵ (astrometric excess noise), the significance of astrometric excess noise D (astrometricexcess noise sig), the ICRS coordinates (α, δ), the Galactic coordinates (l, b), proper motions in the right ascension µα, and declination µ δ , the magnitude of Gaia (G, GBP, GRP).
12 CO integrated intensity of ∼ 20 K km s −1 .The cloud might be the reservoir of star formation to generate the quadruple system.

DISCUSSION
If it were true that the Gaia DR3 astrometric solution of TYC 3340-2437-1 was affected by the orbital motion between binaries B 1 and B 2 , then the outer orbital period P out would be comparable to the coverage time of the astrometric time series.Considering that the Gaia DR3 catalog was created from raw data collected by the Gaia instruments during the first 34 months (Gaia Collaboration et al. 2022), the outer period could be a few years or a decade.
Assuming P out is 5 (or 10) years and it is a circular orbital, the outer semi-major axis a out ∼ 1826(2899) R ⊙ between binary B 1 and B 2 .Following the Equation (12) of Georgakarakos (2008), we obtain the largest stable criteria (a out /a in ) c ∼ 7, assuming binary B 1 as the inner binary, B 2 as the third star with a mass of (m 2a + m 2b ), and the third star exhibiting prograde motion.Therefore, the quadruple system would be stable due to a out /a in ∼ 74 > 7.
The orbital velocities of binary B 1 and B 2 would be 14.9 (or 11.8) and 35.7 (or 28.3) km s −1 .If the outer orbit had an eccentricity e out > 0, the orbital velocity of binary B 1 at the periapsis would be greater than 14.9 km s −1 and even greater than 20 km s −1 .This raises the probability that the arc-like feature of the bow-shock nebula may result from the interaction between the surrounding ISM and stellar wind expelled when binary B 1 near periapsis.
In such a scenario, one would expect the bow-shock nebula to exhibit discrete features in the direction of the arc-like structure, as these shocks could form at different times when binary B 1 approaches periapsis.However, no such discrete features are evident in the WISE image (see insert panel of Figure 8).This absence could be attributed to limitations in the spatial resolution of the WISE telescope.Alternatively, the bow-shock nebula might have formed from the unitary motion of the quadruple system, similar to the bow-shock nebula of a runaway star.It is also plausible that the bow-shock nebula is an "in situ" feature, supported by interactions between TYC 3340-2437-1 and an outflow of hot gas from a nearby star-forming region or the H II region (the bow-shock theory can be found from Kobulnicky et al. ( 2016) and the references therein).
If P out were 5 (or 10) years, the semi-major axis would be about 8.5 (or 13) au.This is significantly less than 500 au, suggesting that TYC 3340-2437-1 might have formed through disk fragmentation (Offner et al. 2023).However, it's noteworthy that the orbital inclinations of binaries B 1 and B 2 are 55.9 +0.55 −0.49 and 78.20 +3.16 −0.96 degrees, respectively.These inclinations indicate that the orbital orientations of the inner binaries are not coplanar, which contrasts with the idea that disk fragmentation mechanisms typically result in multiple systems with initially similar orientations.
However, if the eccentricity of the outer orbit were greater than zero, an interesting possibility arises.The Kozai-Lidov mechanism can induce large-amplitude oscillations of eccentricities and inclinations for highly inclined triple systems (Naoz et al. 2013 and references therein).In a similar context, it's conceivable that binaries B 1 and B 2 were initially formed with similar orientations, but with non-zero eccentric of outer obit, which would lead to changes in their orientations over time through secular evolution of the system.This issue might be figured out when the epoch positions of Gaia are released in the future.

CONCLUSION
Cross-matching TESS LC with the O-type stars found by LAMOST, we find a hierarchical quadruple system TYC 3340-2437-1.It is excluded to be an optical double star and has a probability of more than 99.99% to be a (2+2) quadruple system by analyzing the data of Gaia DR3 and ZTF.Its inner orbital periods are 3.390497(30) days and 2.43781( 16) days.Assuming the four stars are main sequence stars, the total mass of the quadruple system is estimated to be about 24.48 = (11.47+5.79)+(5.2 + 2.02) M ⊙ by a simultaneous fitting of both light curve and radial velocity curve.
The line-of-sight inclinations of the inner binaries are estimated i 1 = 55.94 +0.55 −0.49 and i 2 = 78.20 +3.16 −0.96 .TYC 3340-2437-1 is a quadruple system in SBN and surrounded by cloud with 12 CO integrated intensity of ∼ 20 K km s −1 which could be the reservoir for forming TYC 3340-2437-1.The outer orbital period of TYC 3340-2437-1 might be a few years or a decade.The quadruple system might formed through the disk fragmentation mechanism.The arc-like feature of SBN and non-coplanar of the inner binary might be due to an eccentric outer orbit.The outer orbit period and outer eccentric could potentially be determined when the epoch positions of Gaia are released in the future, which can also be used to constrain the mechanism of star formation.

ACKNOWLEDGMENTS
The authors thank Hao Tian, Man I Lam, Jie Lin, Sheng-Hong Gu, Fang Yang, Er-Lin Qiao, Jian-Rong Shi, Hai-Liang Chen, Zheng-Wei Liu, Xiang-Cun Meng, Yan Gao, Bo Wang, Shu Wang, Xiao Zhou, Zhi-Cun Liu, Yan-Jun Guo and Jian-Ping Xiong for valuable discussion; and thank Jin-Ming Bai, Yu-Feng Fan, Xiao-Guang Yu, Jian-Guo Wang, Kai-Xing Lu, Ju-Jia Zhang, Jie Zheng, Hong Wu for the suggestion and assistance of observation.We acknowledge the support of the staff of the Lijiang 2.4m telescope and the Xinglong 2.16m telescope.We thank the anonymous reviewer for the valuable suggestions and comments.
This work is supported by the National Key R&D Program of China grant (Nos. 2021YFA1600401, 2021YFA1600400, and 2019YFA0405501.L.C.), the Natural Science Foundation of China (Nos. 12090040, 12090043, 11873016, 11873054, 12303068 Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by CASs.Funding for the project has been provided by the National Development and Reform Commission.LAMOST is operated and managed by the National Astronomical Observatories, CASs.The LAMOST FELLOWSHIP is supported by Special Funding for Advanced Users, budgeted and admin-istrated by Center for Astronomical Mega-Science, CASs.This work is supported by Cultivation Project for LAMOST Scientific Payoff and Research Achievement of CAMS-CAS and Special research assistant program of CASs.
This paper includes data collected by the TESS mission.Funding for the TESS mission is provided by the NASA Explorer Program.This work presents 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 Multi Lateral Agreement (MLA).The Gaia mission website is https://www.cosmos.esa.int/gaia.The Gaia archive website is https://archives.esac.esa.int/gaia.This research made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multi-line survey in 12CO/13CO/C18O along the northern galactic plane with PMO-13.7mtelescope.The authors thank all the members of the MWISP working group, particularly the staff members at the PMO-13.7mtelescope, for their long-term support.MWISP was sponsored by National Key R&D Program of China with grant 2017YFA0402701 and CAS Key Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047.

A.2. Observation and data reduction
TYC 3340-2437-1 spectra were collected from three observation facilities, including LAMOST, the Xinglong 2.16-m telescope at the Xinglong station of the National Astronomical Observatories, CAS, and the Lijiang 2.4-m telescope at the Lijiang station of the Yunnan Observatories, CAS.We collected spectra observed in both the low (R ∼ 1800) and medium (R ∼ 7500) resolution surveys from the LAMOST.Additional spectra were observed using the BFOSC E9+G10 instrument with a 1.6 ′′ on the Xinglong 2.16-m telescope (R ∼ 2000), and the YFOSC E9+G10 instrument with a 0.58 ′′ slit on the Lijiang 2.4-m telescope with R ∼ 2500.The BFOSC spectra were reduced and extracted by using the developed Python package BFOSC.The YFOSC spectra were reduced using iraf.Standard reduction procedures, such as bias subtraction, flat fielding correction, spectral order extraction, and wavelength calibration using the Fe/Ar lamp for BFOSC observations (He/Ne lamp for YFOSC) were included to reduce the data.Spectral orders were combined using the Python package pyrafspec.Flux measurements contaminated by skylight were  Table A1.Basic information of each observation facility.R is the spectrum resolution.N obs is observation number.Wavelength stands for the effective wavelength coverage.Note: the effective wavelength coverage of YFOSC spectrum is less than that of BFOSC because the absence of emission lines of HeNe lamp in the wavelength range 3900-5200 Å.
excluded beyond the wavelength range larger than 7200 Å.In Table A1, we list the telescope, instrument, adopted slit, spectral resolution, selected calibration lamp, wavelength coverage, the number of observations, and the location of the telescope.

A.3. Atmosphere Parameters of Star B 1a
Despite the spectral resolution of LAMOST-LRS being R ∼ 1800, for early-type stars, the uncertainties of atmosphere parameters obtained by LAMOST-LRS spectra are better than that of LAMOST-MRS spectra (R ∼ 7500) since the distinctness optical feature of early-star concentrate on the wavelength range of [3700, 5000] Å (the wavelength coverage shown in Table A1).The wavelength of LAMOST-LRS spectrum range from [3700, 9000] Å, while LAMOST-MRS wavelength only is in the intervals of [4900,5350] Å (blue arm) and [6300, 6800] Å (red arm).Moreover, the signalto-noise ratio (SNR) of BFOSC spectra are less than the SNR of LAMOST-LRS spectrum.Therefore, we used LAMOST-LRS spectrum, which has been corrected for its wavelength to the rest frame, to measure the atmosphere parameters.We perform a Markov chain Monte Carlo approach (MCMC) to determine the parameters: effective temperature T eff , surface gravity log g, metallicity Z, and projected rotation velocity v sin i. is, where f i and σ fi represent flux and flux error of each pixel of LAMOST-LRS spectrum, respectively.Pixels within the wavelength range of (3760 < λ < 3900) or (4000 < λ < 4400) or (4460 < λ < 5000) were used to estimate the atmosphere parameters (as the residual spectrum shown in the bottom panel of Figure 4).f syn denote the synthetic spectrum interpolated using the Stellar Label Machine (SLAM) (Zhang et al. 2020;Guo et al. 2021) with a TLUSTY spectra grid.SLAM is a data-driven method based on support vector regression, which can be used to interpolate spectra with a trained model.The model was trained by spectra of a subgrid of TLUSTY OSTAR2002.We have selected a subset: 27500 ≤ T eff ≤ 37500 K with a step of 2500, 3.25 ≤ log g ≤ 4.5 dex with a step of 0.25, four chemical compositions of 0.2, 0.5, 1 and 2 Z ⊙ , and microturbulence velocity ξ t = 10 km s −1 , from the TLUSTY OSTAR2002 grid (Lanz & Hubeny 2003).Each spectrum within the subset grid was broadened with projected rotational velocities 0 ≤ v sin i ≤ 200 km s −1 with a step of 10 (using the broadening kernel described by Eq. 18.14 of Gray ( 2021) ).The limb-darkening coefficient was fixed ϵ = 0.2.Subsequently, the broadened spectra were degraded to an instrument resolution of 1800 and then resampled with a wavelength range of 3700 ≤ λ < 6700 Å, with logarithmic steps of 0.0001 based on a base 10 logarithm.This corresponds to the dispersion per pixel of LAMOST-LRS spectra.
We obtain T eff = 30279 ± 230 K, log g = 3.74 ± 0.03 dex, Z = 0.35 ± 0.02 Z ⊙ , and v sin i = 111 ± 8 kms −1 .Their errors are relatively small (see Figure A.2). To assess this, we applied the Cramér-Rao bound (Cramér 1946;Rao 1945) for the trained SLAM model.The theoretically achievable uncertainties of T eff , log g, Z and v sin i were calculated to be 210 K, 0.027 dex, 0.019 Z ⊙ and 6.5 km s −1 , respectively.They are comparable to the errors of the MCMC approach.Therefore, the low errors are reasonable, since only precision measurements were taken into account.

A.3.1. Contamination of the other stars
The spectrum of star B 1a is contaminated by the other three stars.To simulate the LAMOST-LRS spectrum of the quadruple system, we assigned the following parameters: for star B 1a , T eff = 30000 K, log g = 3.7 dex, Z = 0.35 Z ⊙ , v sin i = 0 km s −1 , and a radius of 7.8 R ⊙ ; for the other three stars, T eff = 19000 K, log g = 4.3 dex, Z = 0.35 Z ⊙ , v sin i = 0 km s −1 , and a radius of 3 R ⊙ .The radial velocities of stars B 1a and B 1b were set to 9 and -18 km s −1 , respectively.The semi-amplitudes of the radial velocities of stars B 2a and B 2b were 83 and 215 km s −1 , respectively.
We varied the center mass radial velocity of binary B 2 relative to binary B 1 from -20 to 20 with a step of 10 km s −1 and the orbital phase of binary B 2 from 0 to 0.9 with a step of 0.1.For each center mass and orbital phase point of binary B 2 , we generated 100 spectra and randomly resampled flux to ensure that each pixel of the mock spectra had the same SNR as the LAMOST-LRS spectrum using a Gaussian distribution.Subsequently, the atmospheric parameters of the mock spectra were predicted using SLAM.
We observed biases of about 200 K, 0.18 dex, -0.05 Z ⊙ , -5 km s −1 in T eff , log g, Z and v sin i respectively, for star B 1a (see back line of Figure A.4).When considering only the contamination of the spectrum by star B 1a , the biases in T eff , log g, Z and v sin i were about -72 K, 0.04 dex, -0.03 Z ⊙ and -5 km s −1 (see blue line of Figure A.4).We also observed that the biases are almost invariable for all orbital phases and center mass velocities of binary B 2 .
For a more comprehensive evaluation that considers both precision and accuracy, we have accepted the uncertainties reported by Guo et al. (2021) for σ(T eff ) = 1600 K, σ(log g) = 0.25 dex, and σ(v sin i) = 42 km s −1 .It should be noted that the resulting estimated metallicity could potentially be unrealistic due to the weakness of metallicity lines in the optical wavelength range, especially given the high temperature of the spectrum.

A.4. Vicinity Eclipse Binary?
Assuming that binary B 2 is a vicinity eclipse binary with an inclination angle of 90°, the primary eclipse relative depth of B 2 is magnitude of binary B 2 is 15.12 mag based on the G RP = 10.12 mag of TYC 3340-2437-1.Crossing match Gaia DR3 with an angular radius of 60 ′′ at the center of TYC 3340-2437-1, we find most of the sources are fainter than 17 mag, expect two sources have magnitude G RP = 14.68 (38 ′′ away; α = 62.7854177°, δ = 50.6993595°)and G RP = 15.22 (54 ′′ away; α = 62.7887683°, δ = 50.6935972°).We obtain their light curves from ZTF DR9 and then fold them into the period of 2.43781 days (see Figure A.5).The folded light curves are almost uniformly distributed in the period phase, indicating that they do not have the orbital period of 2.43781 days.Binary B 2 , therefore, is impossible to be a vicinity eclipse binary.
A.5. Posterior distribution of fitting light curve and radial velocity curve In the process of fitting, we fixed the orbital periods of binary B 1 and B 2 , the eccentricity and transit time of B 2 , and the effective temperature of star B 1a (P 1 = 3.390497 days, P 2 = 2.43781 days, e 2 = 0, T c2 = 2458817.49056BJD, T eff1a = 30279 K).The posterior distribution of MCMC fit by simultaneously using LC and radial velocity curves of binary B 1 .
where θ = (T c1 , √ e 1 cos ω 1a , √ e 1 sin ω 1a , q 1 , K 1 , v γ1 , s, R 1a , r 1 , log g 1a , H 1b , i 2 , L 2a , R 2a , L 2b , R 2b ).T c is the time of transit; e is eccentricity; ω is the longitude of periastron; q is the mass ratio (q 1 = m 1b m −1 1a ); K is radial velocity semi-amplitude; v γ is systemic velocity; R is star radius; r 1 is radius ratio (r 1 = R 1b R −1 1a ); g 1a is the gravity of star B 1a ; H 1b represents the albedo effect of star B 1b , which influences the minima of ellipsoidal light curve.As shown in the middle panel of Figure 2 (b), the minimum at the phase of 0 is less than the minimum at the phase of 0.5 (or -0.5).i 2 is orbital inclination of binary B 2 ; L is luminosity; s is an additional "jitter" variance to description system errors of radial velocities or radial velocity model (In this case, v 1a and v 1b might be influenced by the binary B 2 ).L(x|θ) is the likelihood function of x, x are the TESS light curve, radial velocities of star B 1a and star B 1b .where θ l = (T c1 , √ e 1 cos ω 1a , √ e 1 sin ω 1a , q 1 , K 1 , R 1a , r 1 , log g 1a , H 1b , i 2 , L 2a , R 2a , r sb2 , R 2b ); l n is the observed light curve flux at time t n .σ ln is error of

Figure 1 .
Figure 1.TESS light curve.Panel (a) shows the LC of SAP FLUX extracted by the MIT Quick Look Pipeline(Huang et al. 2020); the blue triangles are the minimum points of the narrow "valley"; the blue dots are the points in the phase of (−0.04 < ϕ < 0.04) or (0.46 < ϕ < 0.54) when the LC is folded with a period of 2.43781 days.Panel (b), the phase is folded with a period of 2.43781 days and to the transit time of 2458817.49056BJD.Panel (c), the LC1 shows the removal of the blue points in the Panel (a)

Figure 2 .
Figure2.Panel (a) shows the TESS light curve (LC, dark points) of SAP FLUX extracted by using the MIT Quik Look Pipeline(Huang et al. 2020); the red lines are light curve calculated by ellc(Maxted 2016).Panel (b) shows that the ellipsoidal light curve of binary B1 (LC1) and corresponding radial velocity curve, both of which are folded by a period of 3.390497 days.In the top panel, blue solid and dash lines are Kepler orbital radial velocities of binary B1a and binary B 1b , respectively; the markers of the open square, dot, and star are radial velocities measured from spectra of BFOSC, YFOSC, and LAMOST-LRS (R ∼ 1800), respectively.The open and solid triangles stand for radial velocities of star B 1b to star B1a, which are measured from spectra of LAMOST-MRS (R ∼ 7500).The red triangles stand for the phases where the spectra are used to measure the radius ratio of star B 1b to star B1a.In the middle panel, the gray dots donate the light curve removed the data near the eclipses of binary B2.The bottom panel shows the folded residual of the TESS LC.Panel (c) shows the relative eclipse light curve of binary B2 (LC2) and folded by a period of 2.43781 days, which is obtained by dividing the TESS LC into the model light curve of binary B1.The remarkable deviations at the phase of ∼ 0.4 of LC2 might be attributed to the de-trending process of the TESS LC, which can be seen at the edges of Panel (a).Note: the residuals of Panel (b) and (c) are the same residual of Panel (a) but folded with different periods.

Figure 4 .
Figure4.The LAMOST-LRS spectrum.Top, the black line is observed spectra, the red line is the best-fit model of the MCMC approach; Bottom, the residual between the model and the observation spectrum, and the gaps in the spectrum are excluded when we determine the atmosphere parameters.The gaps are contaminated by the absorption lines of interstellar material (ISM): Ca II K, Ca II H, and diffuse interstellar bands (4430 Å).

Figure 5 .
Figure5.LAMOST-MRS spectra observed at the orbital phases of 0.24 and -0.1868 (P = 3.390497 days) of binary B1; the three panels are absorption lines of He I 4922, Hα and He I 6678 from left to right; the black, red and blue thin lines are LAMOST-MRS spectra, synthetic spectra of B1a, and B 1b ; the gray thick lines are co-added synthetic lines of B1a, and B 1b .

Figure 2
Figure2(c) shows that the number of solid triangles is five less than the number of open triangles.The radial velocities can be found in TableA.2.Moreover, the radial velocities were nearly independent of T eff1b and log g 1b .In order to obtain an accurate radius ratio, we first selected 12 spectra that were observed near phases -0.25

Figure 6 .
Figure6.The radius ratio of star B 1b to star B1a (r1 = R 1b /R1a) varies with T eff1b and log g 1b .The vertical lines stand for the minimum and maximum errors.

Figure 7 .
Figure 7. Comparson of synthetic and observed photometry.Left axis is the spectral energy distribution; the open dots represent the fluxes from Gaia EDR3 GBP-, G-and GRP-bands, APASS B-, V -, G-, R-, I-bands, 2MASS J-, H-and Ks-bands, and ALLWISE W1-, W2-bands, which are downloaded from CDS portal; the black curve stands the best-fitting SED of TYC 3340-2437-1 summed by SEDs of stars B1a, B 1b , B2a, and B 2b .Right axis shows the transmission curve normalized at the maximum.The transmission curves are labeled as the same colors of open dots.

Figure 8 .
Figure 8. Integrated emission of 12 CO(J = 1−0) of MWISP in the interval of [-30.3, -11.7] kms −1 , which is the 68% confidence interval of the radial velocity measured by Na D1 and Na D2 absorption lines.The blue square is the pixel size of MWISP image.Insert panel is WISE image colored by bands of W1 (blue), W2 (green), and W4 (red).It is the same size as the black quadrangle of the main panel.The black curves (arc) of the main and insert panels have the same radius of ∼ 24.6" (0.33 pc), which go through the intensity peak of WISE W4 and indicate the position of the SBN.
.B.Z. ), and the Science Research Grants from the China Manned Space Project (NOs.CMS-CSST-2021-A10, CMS-CSST-2021-A08).XC acknowledges the National Science Fund for Distinguished Young Scholars (No. 12125303 ).CJ acknowledges the Strategic Pioneer Program on Space Science, Chinese Academy of Sciences (CASs) through grant NOs.XDA15052100, and the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB0550200.
, Python, Numpy (van der Walt et al. 2011; Harris et al. 2020), Scipy (Virtanen et al. 2020; Jones et al. 2001), PyAstronomy (Czesla et al. 2019), Speedyfit.APPENDIX A. APPENDIX INFORMATION A.1.Lomb-Scargle periodogram A Lomb-Scargle periodogram of the TESS LC shows a peak.So we use the Gaussian function to fit the peak and obtain the peak period is 1.69759(59) day (see Figure A.1).The peak period and its error are the mean and square root of the estimated covariance of the mean of the Gaussian function.This is performed by a Python function curve fit.
Figure A.1.Periodogram of the TESS light curve.The black and red lines are the Lomb-Scargle periodogram and the fitted Gaussian curve, respectively.

Figure
Figure A.3.The LAMOST-MRS spectrum.The spectrum was observed at the orbital phase of 0.028665 of binary B1.The right and left panels are blue and red arm spectra of LAMOST-MRS.The black line is interpolated by using the best-fitting atmosphere parameters of the LAMOST-LRS spectrum.The left bottom panel contains blank wavelengths where the DIBs are located.

Figure A. 4 .Figure
Figure A.4.The atmosphere parameters biases of star B1a due to the contamination of spectrum by the components stars.The blue and black lines represent the bias caused by only star B 1b and the other three stars, respectively.

Figure A. 7 .
Figure A.7.The model light curves of best-fit parameters but for different color excess E(B − V ) = 0 and 1.2.The residual is calculated as the difference between the normalized fluxes F E(B−V )=0 − F E(B−V )=1.2 .

Table 2 .
The grid of radius ratio of star B 1b to star B1a (r1 = R 1b /R1a) varied with T eff1b and log g 1b .

Table A2 .
The posterior distribution Radial velocities, RV1 is radial velocity of stat B1a, which is measured by CCF.RV2 is the radial velocity of star B 1b , which is measured by MCMC approach (details see Section 3.2).RVNaD1 and RVNaD2 are RVs of ISM absorption lines Na D1 (5890 Å) and Na D2 (5896 Å), which are measured by the Gaussian fit.