On the stability and pulsation in models of B[e] star MWC 137

B[e] type stars are characterised by strong emission lines, photometric $\&$ spectroscopic variabilities and unsteady mass-loss rates. MWC 137 is a galactic B[e] type star situated in the constellation Orion. Recent photometric observation of MWC 137 by TESS has revealed variabilities with a dominant period of 1.9 d. The origin of this variability is not known but suspected to be from stellar pulsation. To understand the nature and origin of this variability, we have constructed three different set of models of MWC 137 and performed non-adiabatic linear stability analysis. Several low order modes are found to be unstable in which models having mass in the range of 31 to 34 M$_{\odot}$ and 43 to 46 M$_{\odot}$ have period close to 1.9 d. The evolution of instabilities in the non-linear regime for model having solar chemical composition and mass of 45 M$_{\odot}$ leads to finite amplitude pulsation with a period of 1.9 d. Therefore in the present study we confirm that this variability in MWC 137 is due to pulsation. Evolutionary tracks passing through the location of MWC 137 in the HR diagram indicate that the star is either in post main sequence evolutionary phase or about to enter in this evolutionary phase.


INTRODUCTION
Massive stars significantly affect the chemical enrichment and dynamics of galaxies (Nomoto et al. 2013;Goswami et al. 2022).MWC 137 is a massive galactic B[e] -type star which is surrounded by a optical nebula .Sharpless (1959) first classified this nebula as an H II region and named it as Sh 2-266.For MWC 137, Ciatti & Mammano (1975) have found strong H I, faint Fe II and He I lines.In addition to that, emissions of O I and Ca II were also observed in the same study.Ciatti & Mammano (1975) has suggested that the nebula around MWC 137 is caused by the expansion of ejected material from the star.Ulrich & Wood (1981) has reported the presence of He I recombination lines ( 5876 and 10830) in MWC 137 and linked it with chromospheric emissions.This star has been classified as Herbig Ae/Be star by Hillenbrand et al. (1992), Berrilli et al. (1992) and Thé et al. (1994).From the VLA/Australia Telescope survey, Skinner et al. (1993) have reported MWC 137 as a possible non-thermal radio source with a significant decrease in radio flux on the time scale of ≤ 5 months.Miroshnichenko (1994) has taken 19 UBVRĲHK and 12 UBVRI observation for MWC 137 using a 1-m telescope at Assy Observatory.The author has found variability in all the photometric band and classified this star as a B0 V type, having a mass -loss of 1.1 × 10 −6 M ⊙ yr −1 .A quasi-periodic variabil-★ E-mail: sugyannitr4142@gmail.com † E-mail: yadavap@nitrkl.ac.in ity of 4.07 days is reported for MWC 137 by Bergner et al. (1994).
From narrow band H image and high resolution spectroscopy of MWC 137 and its nebula, Esteban & Fernandez (1998) concluded that it is not a Herbig Ae/Be star rather a B[e] type supergiant with luminosity L = 2.36 × 10 5 L ⊙ .The same authors have suggested that the nebula around MWC 137 (S 266) is likely to be produced from the interaction of stellar winds and interstellar medium or unprocessed ejected material.Henning et al. (1998) have discovered an extended bipolar millimeter continuum source in the vicinity of central star MWC 137.Fuente et al. (2003) have reported the presence of circumstellar disk in three B0 stars including MWC 137.These authors have derived an upper mass limit of 0.007 M ⊙ for the disk around this star.Marston & McCollum (2006) and Marston & McCollum (2008) have noted the presence of a diffuse bipolar structure in the ring of nebula around MWC 137.Spectral lines of rings around the central star MWC 137 show that the ring is expanding with velocity of 10-15 km s −1 indicating a dynamical age of 10 5 yr (Esteban & Fernandez 1998).Muratore et al. (2015) have reported that the central star MWC 137 is surrounded by cool and dense ring of CO gas.The derived temperature, density and rotational velocity of the ring are horizontal 1900 (± 100) K, 3 (±1) × 10 21 cm −1 and 84 (±2) km s −1 , respectively.Muratore et al. (2015) have argued that this star is not a Herbig Be star as the gas present in the ring is enriched with isotope 13 C.
The ratio of 12 C and 13 C has become a diagnostic tool to distinguish young, pre-main sequence stars and evolved stars (see e.g., Kraus 2009;Liermann et al. 2010;Muratore et al. 2010;Oksala et al. 2013; Kraus 2019).The isotope ratio 12 C / 13 C will decrease during the evolution of stars from the initial unprocessed interstellar value of 90 to less than 5 (Kraus 2019).The derived value 25 ± 2 of the isotope ratio by Muratore et al. (2015) is appropriate for proto-planetary nebula, main-sequence or supergiant evolutionary phase.Using VLT MUSE IFU data, Mehner et al. (2016) have discovered a collimated outflow from MWC 137.The emerging jet from the star has many knots and the electron density of 10 3 cm −3 .These authors have concluded that the MWC 137 is in post main sequence evolutionary phase and proposed several possible scenarios for the presence of jets and disk including binary system and stellar merger.Kraus et al. (2017) have performed multi-wavelength study of this star and its surrounding nebula.The authors also have noted that northern part of the nebula is mostly blue shifted while the southern part is primarily red shifted.In the same study, it has also been concluded that the large scale nebula may not be created and shaped only by MWC 137.Kraus et al. (2021) have noted a substantial variation in the electron density (ranging from 0 to 800 cm −3 ) across the nebula.The overall morphology of the nebula is found to be unchanged within 18 yr.Using the time series photometric data of Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), Kraus et al. (2021) have reported a stable variability of 1.93 d for MWC 137.The cause of this variability is unknown.The period of 1.93 d is too short to be originated by rotation or orbital motion due to any companion.Therefore it is speculated to be linked with pulsation in MWC 137.In several Btype massive stars, photometric and spectroscopic variabilities have been observed (see e.g.Frost 1902;Guthnick 1913) and found to be caused by stellar pulsation (e.g., Yadav et al. 2021;Yadav & Glatzel 2016).Since the luminosity to mass ratio for MWC 137 is greater than 10 3 in solar unit therefore dynamical instabilities can be expected in models of this star.In models of massive O and B type stars, it has been found that these instabilities may lead to pulsation, surface eruptions and re-arrangement of stellar structure (see e.g., Glatzel et al. 1999;Yadav & Glatzel 2016, 2017).In the present study, we intend to investigate instabilities and their consequences in models of MWC 137.Instabilities identified in the linear stability analysis will be further considered for nonlinear numerical simulations to find out the final fate of unstable models.Comparing the obtained numerical results with observed variabilities in MWC 137 will allow us to confirm or rule out pulsation as the cause of reported 1.9 d variability.Photometric variability of MWC 137 using the data of TESS is presented in section 2. In section 3, properties of considered models of MWC 137 are described.Linear stability analysis and associated results are presented in section 4. Nonlinear numerical simulation for selected unstable models is mentioned in section 5. Evolutionary status of this star is highlighted in section 6. Discussion is mentioned in section 7 followed by conclusions in section 8.

PHOTOMETRIC VARIABILITY OF MWC 137 FROM TESS
The normalized light curves for MWC 137 were generated using data from sectors 6 and 33 of TESS.The data for both sector 6 and sector 33 was taken in a 2 yr gap with high cadence (120 s) and are displayed in Fig. 1 and Fig. 2, respectively.Further in order to study the TESS time series in the frequency domain, periodograms / Fourier analysis is required.For TESS observations, we have computed Lomb-Scargle periodogram as they are more appropriate for unevenly spaced data.This Lomb-Scargle analysis is done using the Lightkurve package (Lightkurve Collaboration et al. 2018).For each sector, we have detected a dominant period of 1.9 d, as is evident from Fig. 1, 2. The same variability is also reported in a recent paper by Kraus et al. (2021) using PERIOD04 (Lenz & Breger 2005).

MODELS FOR MWC 137
It is challenging to determine fundamental parameters such as mass, luminosity and surface temperature for B[e] type massive stars.To construct the models of MWC 137, we have adopted the average value of surface temperature T eff = 28200 K and luminosity log L/L ⊙ = 5.84 from Kraus et al. (2021) .There is a large uncertainty in the distance of MWC 137 which ranges from 1 kpc to 13 kpc.Kraus et al. (2021) have derived the value of luminosity using the distance 5.15 kpc.This distance is in agreement with the value (5.2 kpc) determined by Mehner et al. (2016).There is a large uncertainty in the mass of MWC 137 hence models having a sequence of mass in the range of 27 to 70 M ⊙ have been considered.This mass range includes the mass 37 +9 −5 M ⊙ derived by Kraus et al. (2021) for this star.The present study is focused on the models of this star having solar chemical composition (Z = 0.02, Y = 0.28, X = 0.70).Additionally, to study the dependency of instabilities on metals, we have also considered models having Z = 0.01 and Z = 0.03.
For constructing the models of MWC 137, we have integrated stellar structure equations -mass conservation, momentum conservation, energy conservation and energy transport -from the surface upto a point where temperature is 10 7 K.With this integration, we obtain envelope models of the star to perform linear stability analysis followed by nonlinear simulations.Since instabilities are mostly found in the envelopes of massive O and B type stars therefore we have restricted our present study in envelopes to avoid the complications of nuclear reactions in the core.Schwarzschild criterion is applied for convection.Standard mixing length theory is used with the mixing length parameter = 1.5 for convection (Böhm-Vitense 1958).Rotation and magnetic fields have been disregarded.For the opacity, OPAL tables (Rogers & Iglesias 1992;Rogers et al. 1996;Iglesias & Rogers 1996) are used.

STABILITY ANALYSIS FOR MODELS OF MWC 137
Linearised pulsation equations to perform the stability analysis for radial perturbation have been adopted from Gautschy & Glatzel (1990a).This set of equations form a fourth order boundary eigenvalue problem.Out of four boundary conditions, two boundary condition are applied at the surface and remaining two are used at the bottom of the envelope.The eigenvalue problem is solved using the Ricatti method as described by Gautschy & Glatzel (1990a).Obtained eigenfrequencies are complex where real part denotes pulsation period and imaginary part indicates damping ( > 0) or excitation ( < 0).The calculated eigenfrequencies are normalised by global free fall timescale.
The outcome of linear stability analysis is generally depicted in a form of 'modal diagram'.In the modal diagram, real and imaginary part of eigenfrequencies are given as a function of fundamental stellar parameter such as mass or surface temperature.Modal diagram for the case of models having solar chemical composition is given in Fig. 3, where real part of the eigenfrequency is mentioned in Fig. 3(left) and imaginary part is given in 3(right).Thick blue lines in Fig. 3(left) correspond to unstable modes where imaginary part of eigenfrequency is negative.All of the considered sequence of models starting from 27 M ⊙ to 70 M ⊙ are unstable.Number of unstable modes is increasing for models having higher luminosity to mass   ratio.We also note mode interaction phenomena which is frequently present in the models of high luminosity to mass ratio.One of the modes is unstable for all the considered models from 27 M ⊙ to 70 M ⊙ .The strength of instability associated with this mode increases for lower mass models.In the vicinity of this mode, another mode is unstable for models having mass below 40 M ⊙ .For the model of mass 29 M ⊙ , two low order unstable modes are showing interaction in real part of eigenfrequeny while their imaginary parts are well separated.Real part of eigenfrequency associated with one of the damped modes is approaching towards zero for model having mass 48 M ⊙ .Linear stability analysis has also been performed for models having following chemical compositions: (i) X = 0.71, Y = 0.28, Z = 0.01 (ii) X = 0.69, Y = 0.28, Z = 0.03 The modal diagram for models having Z = 0.01 is given in Fig. 4. In this case, models only in the mass range between 27 M ⊙ to 41 M ⊙ and 59 M ⊙ to 70 M ⊙ are unstable.Number of unstable modes and strength of instabilities have reduced compared to models having Z = 0.02.All the considered higher order modes are damped.
For the models having Z = 0.03, modal diagram is given in Fig. 5.In this case, all the considered models in the mass range 27 M ⊙ to 70 M ⊙ are unstable.The modal diagram for Z = 0.03 is qualitatively similar to the Z = 0.02.However, for the case of Z = 0.03, one additional mode is unstable for models having mass above 62 M ⊙ .Strength of instabilities associated with the unstable modes are slightly more than the models with Z = 0.02.
Outcomes of the linear stability analyses mentioned in the modal diagrams (Fig. 3 to Fig. 5) are based on the models of MWC 137 having surface temperature of 28200 K (log T eff ≈ 4.45).However, the previous studies have shown that the value of log T eff for this star is in the range of 4.40 to 4.49 (Kraus et al. 2021).To study the presence and strength of instabilities as a function of surface temperature, we have also carried out linear stability analysis for the models having mass of 45 M ⊙ with solar chemical composition and surface temperature (log T eff ) in the range of 4.40 to 4.49.The result is presented in the Fig. 6 where real and imaginary part of eigenfrequencies are plotted as a function of surface temperature.Similar to the modal diagram of Fig. 3, blue lines in the real part and negative imaginary part indicate unstable modes.From this figure, we infer that all the considered models of MWC 137 in the surface temperature (log Teff) range of 4.40 to 4.49 are unstable and pulsation period associated with unstable mode of different models is almost same.
Similar to the ambiguity in the determination of surface temperature, luminosity of the star is also subject to possible errors mostly due to distance determination.Kraus et al. (2021) have determined the luminosity log L/L ⊙ = 5.84 with an error of ± 0.13 for this star.Earlier studies have shown that the high luminosity to mass ratio can induce several dynamical instabilities such as strange mode instabilities (Gautschy et al. 1990;Kiriakidis et al. 1993).To examine the luminosity dependency in the present case we have selected models with solar chemical composition having mass of 45 M ⊙ , log T eff = 4.45 and luminosity in the range of log L/L ⊙ = 5.68 to log L/L ⊙ = 6.00.The outcomes are given in the modal diagram mentioned in Fig. 7 where real and imaginary parts of the eigenfrequencies are plotted as a function of luminosity.In total, two modes are unstable in the considered luminosity range.One of the low order modes is unstable in all the models while the other mode is only unstable in models having log L/L ⊙ > 5.92.For both unstable modes, strength of instabilities is increasing in models having higher luminosity to mass ratios which is expected and in agreement with previous studies (see e.g.Glatzel 1994;Saio et al. 1998).
For models having solar chemical composition (Z = 0.02), periods associated with different modes as a function of mass are given in Fig. 8.In this figure, thick blue lines indicate unstable modes.A dotted horizontal line represents the observed period of 1.9 d.The vertical pink columns in Fig. 8, represent two different mass range (31 to 34 M ⊙ and 43 to 46 M ⊙ ) where unstable modes have period close to 1.9 d.Periods associated with two unstable modes for models having mass close to 33 M ⊙ and 45 M ⊙ are exactly matching with the observed period of 1.9 d.Model having mass of 33 M ⊙ has two unstable modes with a period of 1.9 d and 3.1 d.Strength of instability associated with the mode having period 3.1 d is more than the mode with the period of 1.9 d (see Fig. 3).Therefore unstable mode with the period of 3.1 d is likely to dominate over the other unstable mode.Periods of all the unstable models are below 4 d in the considered mass range.

INSTABILITIES IN NON-LINEAR REGIME
Linear stability analysis has revealed that several low order modes are unstable in the considered three sets of models for MWC 137.
To determine the final fate of unstable models, we have followed instabilities in nonlinear regime using the similar set of equations as described by Grott et al. (2005) for selected models of MWC 137.The adopted numerical scheme is energy conservative which is essential for the reliability of simulation's outcomes.To handle the shock waves, artificial viscosity is introduced (see e.g.Grott et al. 2005) with the value of parameter 10.Sensitivity of the numerical outcomes on the artificial viscosity parameters are discussed in detail by Grott et al. (2005) and Yadav et al. (2018).During the non-linear simulation of O and B type stars, the used value (10) of artificial viscosity is found to be adequate for smearing out discontinuities and suppressing Gibbs's phenomena (Yadav & Glatzel 2016;Yadav et al. 2018).For nonlinear numerical simulations, four boundary conditions are used.Two of them are imposed at the bottom of the envelope for which zero velocity and a constant flux are considered.The remaining two are applied at the surface where the gradient of compression is set to zero and no heat is stored.This choice of outer boundary conditions will enable the shock waves to pass without being reflected (for a detailed discussion the reader is referred to Grott et al. 2005).Outcome of nonlinear simulation for two of the considered models having mass of 45 M ⊙ and 33 M ⊙ are given in Fig. 9 and 10, respectively.The outcome of these two models are presented here as these models show unstable mode in linear stability analysis with the period of 1.9 d.In Fig. 9, variation in the radius 9(a), velocity 9(b) and temperature 9(c) associated with outermost grid point, change in the bolometric magnitude 9(d), time integrated acoustic energy 9(e) and error in the energy balance 9(f) are given as a function of time.Evolution of instabilities starts from an initial value 10 −4 cm/s for velocity.Code picks up the instability from this value of numerical noise without external perturbation and enter into the phase of exponential growth.After 50 days velocity amplitude saturates with a value close to 157 km/s (Fig. 9b).The final value of velocity amplitude is approximately 22 percent of the escape velocity of this model.After entering into the non-linear regime, the mean radius of the model has increased to 2.8 × 10 12 cm from the initial hydrostatic value of 2.3 × 10 12 cm (Fig. 9a).Due to this inflation of the radius, a drop in surface temperature is expected.However, due to ambiguity in the outer boundary condition, the temperature associated with the outermost grid point may not necessarily represent surface temperature of the star (Fig. 9c).Variation in the bolometric magnitude is clearly showing a period of 1.9 d (Fig. 9d).This 1.9 d period is also present in the variation profile of radius, velocity and temperature.The value of time integrated acoustic energy is of the order of 10 40 erg (Fig. 9e) which is four order of magnitude larger than the error in the energy balance.This time integrated acoustic energy (4 R 2 vp, where p and v stand for pressure and velocity perturbations at the model's outer boundary, and R is the radius of that boundary) is representing the mechanical energy lost from the system in the form of acoustic waves.As described by Yadav & Glatzel (2017), during a pulsation period, phases of incoming and outgoing fluxes are present.While integrating for a complete cycle, the outgoing energy is slightly higher than the incoming energy.Resultantly, the integrated acoustic energy is increasing as a function of time after the onset of the regular pulsation cycle in the simulation.The slope of the integrated acoustic energy can be used to determine the mean mechanical luminosity (see also, Grott et al. 2005).In the performed non-linear numerical simulation, the error in the energy balance is several order of magnitude smaller than time integrated acoustic energy and kinetic energies.It represents the energy conservation which is essential in the modelling of stellar instabilities and pulsation (see e.g.Grott et al. 2005;Yadav et al. 2018).Similar to the model of 45 M ⊙ , evolution of radius, variations associated with the outermost grid point of velocity and temperature for a model having mass of 33 M ⊙ are given in Fig. 10.In the nonlinear regime, we note a substantial inflation as the radius becomes more than seven times the initial hydrostatic radius (Fig. 10a).Velocity variations (Fig. 10b) show that the code picks-up the instability from the numerical noise and it goes through the phase of exponential growth.Due to the inflation, as noted in the radius, temperature associated with the outermost grid is dropping and attaining a value less than 7000 K (Fig. 10c).Due to unavailability of opacity table for further calculations, numerical simulation had to be terminated.Linear stability analysis has shown that this model has two unstable modes with periods of 1.9 d and 3.1 d.During the phase of oscillatory exponential growth, a period of 3.1 d is present which corresponds to the more strongly unstable mode.
As pointed out earlier, outermost grid point used for the numerical simulation may not necessarily represents the photosphere of the star.Therefore for direct comparison with observation, variation of quantities in the vicinity of photosphere is more useful.For a model having mass of 45 M ⊙ , variation in velocity, radius and temperature associated with the grid point in the vicinity of photosphere as a function of normalised time is given in Fig. 11 for three pulsation cycles.In addition to these quantities, variation in the bolometric magnitude of the outermost grid point is also mentioned in Fig. 11.For this purpose, the grid point is selected where T 4 (where is Stefan-Boltzmann constant and T represents temperature) is acquiring an equivalent value of the total flux.The location of the grid point satisfying the adopted criteria is changing during the pulsation cycles.Period of 1.9 d is also present at this grid point in the variation profile of velocity, radius and temperature (see Fig. 11).Temperature and radius variations are in the range of 25300 K to 29800 K and 2.35 × 10 12 cm to 2.95 × 10 12 cm, respectively.Variation in the bolometric magnitude is approximately in the range of -0.15 to +0.10 during these three cycles.The shape of the magnitude variation profile of Fig. 9 and Fig. 11 is the same as these variations are connected with the outermost grid point.The pulsation period of 1.9 d is present in both the cases.

EVOLUTIONARY STATE OF MWC 137
Linear stability analysis and non-linear numerical simulation have shown that 1.9 d variability can be explained by stellar pulsation as found in the models having mass of approximately 45 M ⊙ .Evolutionary stage of many B[e] type stars including MWC 137 is uncertain (Mehner et al. 2016;Kraus 2019;Kraus et al. 2021).To determine the evolutionary state of MWC 137, we have plotted evolutionary tracks with initial masses of 44 M ⊙ , 46 M ⊙ , 48.2 M ⊙ , 50 M ⊙ , 52 M ⊙ and 54 M ⊙ in Fig. 12 using solar chemical composition with the help of MADSTAR 1 which is based on the stellar evolution code Evolve ZAMS (Paxton 2004).The Evolve ZAMS (EZ) can evolve the stellar models having masses in the range of 0.1 to 100 M ⊙ with 1 http://user.astro.wisc.edu/~townsend/static.php?ref=ez-web several possible values of metallicity.It creates one-dimensional, non-explosive, non-relativistic and non-rotating stellar models using OPAL opacity tables.For more details on the adopted numerical scheme and parameters used for the stellar evolution, we refer to Paxton (2004) and references therein.The evolutionary track of 48.2 M ⊙ is passing through the location of MWC 137 after the over all contraction phase.The location of MWC 137 is marked using a red circle with the errors in the luminosity and surface temperature in Fig. 12.The location of this circle is fixed using the average value of luminosity and surface temperature of this star as suggested by Kraus et al. (2021).Due to ambiguity in the distance of this star, error in the luminosity is relatively high and consequently leading to uncertainty in the mass determination using evolutionary tracks.
Considering the location of MWC 137 including the error bars in the Hertzsprung-Russell (HR) diagram and evolutionary tracks for different massive stars, it is evident that MWC 137 is either a post main sequence star or about to enter in the post main sequence evolution.We infer from the evolutionary track of the model having initial mass 48.2 M ⊙ that the star has just crossed the overall contraction phase and currently H-buring is taking place in shell around the He-core.
The derived age for the star from this evolutionary track is 3.9 × 10 6 yr which is in the range as estimated by Kraus et al. (2021) using rotating and non-rotating stellar models.Present analysis is limited to the non-rotating evolutionary tracks hence the inclusion of rotation may influence the derived value of the age and the structure of the star.Number of unstable modes and pulsation periods associated with these modes are generally affected by the parameters such as mass, surface temperature, luminosity and chemical composition.For linear stability analysis, additional models of MWC 137 have been included to consider the errors present in the determination of surface temperature and luminosity of this star (Fig. 6 and Fig. 7).Variation in the period of the unstable mode as a function of surface temperature for models with 45 M ⊙ is represented by a secondary horizontal axis created using light brown shaded region in Fig. 12. Similarly changes in the period associated with the most unstable mode as a function of luminosity is shown by a secondary vertical axis through a shaded light green region in Fig. 12.We infer from Fig. 12 that for a given mass, models with higher luminosity and lower surface temperature are having longer periods (see also Fig. 6 and Fig. 7).

DISCUSSION
In the present study, we have analyzed the TESS light curve of MWC 137 observed during sector 6 and 33 using the package 'Lightkurve' (Lightkurve Collaboration et al. 2018).In agreement with Kraus et al. (2021), we have also found the dominant period of 1.9 d in the data of both sectors for MWC 137.To understand the 1.9 d variability in this star, we have constructed three different set of models in the mass range between 27 to 70 M ⊙ and performed linear stability analysis.Linear stability analysis has revealed that the low order radial modes are unstable.In model-sequence having metallicity Z = 0.02 and Z = 0.03, all the models of MWC 137 are unstable while in case of Z = 0.01, models having mass in the range of 42 M ⊙ to 58 M ⊙ are stable.The number of unstable modes in models of Z = 0.01 is smaller compared to Z = 0.02 and Z = 0.03 which is in accordance with earlier studies where higher metallicity leads to more unstable modes in models of massive stars (Kiriakidis et al. 1993;Shiode et al. 2012).For the case of solar chemical composition, the unstable mode having period of 1.9 d is also present in all the models having mass of 45 M ⊙ and surface temperature in the range of 4.40 < log T eff < 4.49 (Fig. 6).Therefore the instability present in these models is not being influenced by the ambiguity present in the measurement of surface temperature for this star.Models of MWC 137 having higher luminosity tend to have more unstable modes and the strength of instabilities associated with different unstable modes are increasing with luminosity (Fig. 7).Mode interaction phenomena are also evident in all three modal diagrams of the considered models.Generally mode interaction is frequently present in models having high luminosity to mass ratios (more than 10 3 in solar unit) while the mode interaction in the real part of eigenfrequencies is rare for models of less luminous stars such as scuti stars (Yadav 2019).For the case of scuti stars, real part of the eigenfrequency associated with different consecutive modes are approximately equispaced (see e.g.Yadav 2019) while the models of Wolf-Rayet stars, luminous blue variables and other massive stars show frequent mode interaction in their modal diagrams (see e.g.Kiriakidis et al. 1996;Glatzel & Kiriakidis 1993).Strange mode instabilities (Glatzel 1998;Saio et al. 1998) are often present in mode interaction scenarios.Since unstable modes in the considered models of MWC 137 are also exhibiting mode interaction phenomena therefore we suspect that the found instabilities may be associated with strange modes.However non-adiabatic-reversible (NAR) approximation (Gautschy & Glatzel 1990b) will be required to confirm the presence of strange mode instabilities which is beyond the scope of present study.
In the linear stability analysis, models of 45M ⊙ with solar chemical composition has an unstable mode with a period of 1.9 d.This theoretically obtained period is in agreement with the dominant variability as found in TESS data for the star MWC 137 in the present study and by Kraus et al. (2021).To determine the final fate of unstable models, we have performed non-linear numerical simulation for a model of this star having mass 45 M ⊙ and solar chemical composition.In the non-linear regime, the instabilities lead to finite amplitude pulsation with a period of 1.9 d for the considered model (see, Fig. 9).During the evolution of instabilities into the non-linear regime, error in the energy balance should be smaller than the time integrated kinetic and acoustic energies (see e.g.Grott et al. 2005;Glatzel & Chernigovski 2016).Non-linear numerical simulations have been carried out carefully where error in the energy balance is four order of magnitude smaller than the time integrated acoustic energy.Therefore the obtained finite amplitude pulsation with a period of 1.9 d in the model of MWC 137 is resulting from physical instabilities.In addition to the dominant 1.9 d period, MWC 137 also has other variabilities with periods in the range of 3 to 8 d (Kraus et al. 2021) likely to be associated with non-radial modes.The variations of quantities mentioned in Fig. 9 and Fig. 11 are resulting from the radial instabilities therefore the direct comparison of variation profiles with observations needs to be done cautiously.
Evolutionary state of MWC 137 is uncertain.In earlier studies, this star is classified as a pre-main sequence star (see e.g.Hillenbrand et al. 1992;Berrilli et al. 1992) while in recent studies this star is classified as a post main sequence star (see e.g.Mehner et al. 2016;Kraus et al. 2017Kraus et al. , 2021)).With the help of evolutionary tracks, we find that this star is either a post main sequence star or about to start post main sequence evolution.Evolutionary track which is passing through the average value of surface temperature and luminosity in HR diagram is indicating that the star is burning hydrogen in a shell around helium core.Evolutionary tracks are indicating that the mass of MWC 137 is more than 40 M ⊙ which is in agreement with the model's mass of 45 M ⊙ showing the pulsation of 1.9 d in linear stability analysis and non-linear numerical simulation.Suggested mass range of 10 to 15 M ⊙ by Mehner et al. (2016) appears to be quite low while the mass of 45 M ⊙ is within the range as proposed by Kraus et al. (2021).
Using the travel distance of 1.2 pc and jet velocity of 650 km/s, Mehner et al. (2016) has estimated that the emanated jet from the central star is at least 1800 yr old.Similarly the authors have also derived the minimum age of 20 000 yr for the large nebula surrounding the central star MWC 137.The value of log T eff for MWC 137 ranges between 4.40 to 4.49 (Hillenbrand et al. 1992;Esteban & Fernández 1998;Alonso-Albi, T. et al. 2009).Considered stellar evolutionary models having surface temperature in this range are also going through the phase of overall contraction where radius and surface temperature decreases and increases, respectively.Overall contraction phase for the model having mass 48.2 M ⊙ is of the order of 50 000 yrs which is longer than the minimum age of the jet and nebula of the central star.Formation of disk around stars and emergence of jets are complex processes which are poorly understood.Generally jet is found in the astrophysical systems having accretion, magnetic field and rotation (Bogovalov & Tsinganos 1999;Coffey et al. 2008).
The location of overall contraction phase is also partially coinciding with the position of MWC 137 in the HR diagram.By comparing the involved timescales of overall contraction phase, minimum age of nebula and jet, we speculate that the formation of disk and jets may be linked with the overall contraction phase in this star.To establish this link, further theoretical and observational studies are required.

CONCLUSIONS
TESS light curve analysis followed by linear stability analysis and non-linear numerical simulations have been performed for the star MWC 137.The outcome of the present study is summarized in the following four points: (i) Similar to the Kraus et al. (2021), we also find a dominant variability with a period of 1.9 d in TESS data of MWC 137.This period is present in both sector 6 and sector 33 of TESS in which this star is observed.
(ii) Linear stability analyses performed in models of MWC 137 show that the low order radial modes are excited as found in case of other B-type stars (Yadav & Glatzel 2017).Frequent mode interactions in the modal diagram is indicating the presence of strange mode instabilities in models of this star.
(iii) Linear stability analyses followed by nonlinear numerical simulations show that a model with mass of 45 M ⊙ has an excited mode with a period of 1.9 d.Therefore observed dominant variability can be explained using low order radial pulsation mode.
(iv) Evolutionary tracks passing through the location of MWC 137 in the HR diagram is indicating that this star is likely in post main sequence phase or about to end the main sequence evolutionary phase.
Mikulski Archive for Space Telescopes (MAST) located at the Space Telescope Science Institute (STScI).

Figure 1 .
Figure 1.Lightcurve (left) and Periodogram (right) of MWC 137 observed under sector 6 by TESS.Variability in the lightcurve and a dominant period 1.9 d are present.

Figure 8 .
Figure 8. Period associated with different modes are plotted as a function of mass for models of MWC 137 having solar chemical composition (z = 0.02).Unstable modes are represented by blue lines.Horizontal dotted line is indicating the dominant period of observed variability.Vertical brown columns represent the location of models having period close and equal to 1.9 d.Same as Fig. 3, value of log T eff (K) = 4.45 and luminosity log L/L ⊙ = 5.84 have been used for these models.

Figure 9 .Figure 10 .
Figure 9. Evolution of instability in a model of MWC 137 having mass of 45 M ⊙ and solar chemical composition: Variation in radius (a), variation in velocity (b) and temperature (c) associated with outermost grid point are plotted as a function of time.Changes in the bolometric magnitude (d), time integrated acoustic energy (e) and error in the energy balance (f) are given as a function of time.Instabilities lead to finite amplitude pulsation with a period of 1.9 d which is clearly visible in the profile of bolometric magnitude.

Figure 11 .
Figure 11.Variations in velocities (a), radius (b), temperature (c) associated with the grid point close to the photosphere and bolometric magnitude (d) of the outermost gridpoint as a function of normalised time for the model with mass of 45 M ⊙ .

Figure 12 .
Figure 12.Evolutionary tracks of stellar models having solar chemical composition.The position of MWC 137 including errors in the surface temperature and luminosity is indicated by orange dot.Initial masses associated with the evolutionary tracks are given with the corresponding color code.Periods associated with the most unstable mode for the models having mass of 45 M ⊙ are also mentioned for the range of errors in surface temperature and luminosity.