On attempting to automate the identification of mixed dipole modes for subgiant stars

The existence of mixed modes in stars is a marker of stellar evolution. Their detection serves for a better determination of stellar age. The goal of this paper is to identify the dipole modes in an automatic manner without human intervention. I use the power spectra obtained by the Kepler mission for the application of the method. I compute asymptotic dipole mode frequencies as a function of coupling factor and dipole period spacing, and other parameters. For each star, I collapse the power in an echelle diagramme aligned onto the monopole and dipole mixed modes. The power at the null frequency is used as a figure of merit. Using a genetic algorithm, I then optimise the figure of merit by adjusting the location of the dipole frequencies in the power spectrum}. Using published frequencies, I compare the asymptotic dipole mode frequencies with published frequencies. I also used published frequencies for deriving coupling factor and dipole period spacing using a non-linear least squares fit. I use Monte-Carlo simulations of the non-linear least square fit for deriving error bars for each parameters. From the 44 subgiants studied, the automatic identification allows to retrieve within 3 $\mu$Hz at least 80\% of the modes for 32 stars, and within 6 $\mu$Hz at least 90% of the modes for 37 stars. The optimised and fitted gravity-mode period spacing and coupling factor agree with previous measurements. Random errors for the mixed-mode parameters deduced from Monte-Carlo simulation are about 30-50 times smaller than previously determined errors, which are in fact systematic errors. The period spacing and coupling factors of mixed modes in subgiants are confirmed. The current automated procedure will need to be improved using a more accurate asymptotic model and/or proper statistical tests.


Introduction
The internal structure of stars has been derived with great details with the advent of the space missions CoRoT and Kepler (Michel et al. 2008;Chaplin et al. 2010). The detection of pressure modes and mixed modes in red giants provided a wealth of information about the internal structure of these stars (Mosser et al. 2014;Bedding et al. 2011) and their internal dynamics (Mosser et al. 2012a). Subgiant stars bridges the gap between the evolved red giants and the Main Sequence (MS) stars. Subgiant stars being closer to MS stars can help to infer informations about the evolution of their interior and dynamics. For subgiant stars, Deheuvels et al. (2014) also measured their core rotation using these mixed modes. The properties of these mixed modes are such that they are easily detectable because they propagate like pressure modes until the surface of the star and are excited by turbulent convection just like the pressure modes (Goldreich & Keeley 1977;Houdek et al. 1999;Dupret et al. 2009)t; but these modes also propagate like gravity modes in the stellar core Unno et al. (1989). Given these properties, these modes have also been searched in the Sun with no positive results .
The detection of these mixed modes in subgiants stars provides also an important marker of the age of these stars. The subgiants evolve away from the MS by burning hydrogen in shells after exhaustion on the hydrogen in the core. The time it take for these stars from the core exhaustion to the Helium burning phase is called the First Giant Branch (Vassiliadis & Wood 1993). The time spent on this branch last about 15% to 30% of the time spent on the MS (Vassiliadis & Wood 1993). In addition, due to the presence of an Helium core, there is large increase in the Brunt-Vaïsala frequency which makes possible the existence of a propagation zone for gravity modes (Dziembowski et al. 2001;Christensen-Dalsgaard 2004). It leads to the existence of the so-called mixed modes that have a pressure mode character in the outer stellar regions, and a gravity mode character in the stellar core (Dziembowski et al. 2001;Christensen-Dalsgaard 2004). In addition, it is clear that since the change in the Brunt-Vaïsala frequency is so quick, the frequencies of the mixed modes change even faster because of their sensitivities to that frequency . Therefore the detection of mixed modes in subgiants is a powerful tool for getting a precise stellar age (Lebreton et al. 2014a), even for getting a precise estimate of the mass (Deheuvels & Michel 2011).
A major challenge of the ESA PLATO mission is to derive stellar ages to better than 10% (Rauer et al. 2014) for the sample 1 of stars brighter than m V = 11 magnitude and with an instrumental noise lower than 50 ppm.hr −1/2 provided with the 24 cameras (Marchiori et al. 2019;Moreno et al. 2019). The derivation of the stellar age is obtained by comparing theoretical mode frequencies to observed mode frequencies (Lebreton et al. 2014b). The number of stars in the P1 sample for which an asteroseismic age with a precision can be derived to better than 10% is larger than 8000. The measurement of the observed mode frequencies can be achieved using the following steps: detecting the stellar oscillation mode envelope, inference of the large separation of the modes (related to the diameter of the star), tagging of the degree of the modes and mode frequency fitting using either Maximum Likelihood Estimators or Markov Chain Monte Carlo estimators. All these steps can be automated for fitting a large number of stars (See Mosser & Appourchaux 2009;Hekker et al. 2010, and references therein). These automated method can be easily applied to MS stars having no mixed modes (Fletcher et al. 2011).
It was already identified by Appourchaux et al. (2004) that the identification of mixed modes for a large amount of in subgiants could be a challenge. For red giants having mixed modes, Vrard et al. (2016) devised a technique for automatically deriving the gravity mode period spacing allowing to fit the mixed modes. Their technique is applied to stars for which the large separation is smaller than 18 µHz but they use guess for the gravity mode period spacing derived from Mosser et al. (2014); their technique is then not entirely automated. For subgiant stars, Bedding (2012) devised a replicated echelle diagrams allowing to visually identify dipole mixed modes. For subgiants stars with a large separation typically higher than 30 µHz (Mosser et al. 2014), an automated detection procedure has been recently devised by Corsaro et al. (2020) for power spectra having the l = 1 mode ridge slightly distorted. For more complicated case, the mixed mode frequencies of the subgiant stars of Appourchaux et al. (2012a) were all identified by hand. In the case of PLATO, the manual identification becomes impractical for the mere 3000 subgiants of the P1 sample of PLATO, therefore requiring an automated procedure for that aim.
The goal of this paper is to present a way to automatically provide guess for the frequencies of dipole mixed modes in subgiants. The first section recalls how the frequencies of these modes can be derived from an asymptotic expression based upon the work of Shibahashi (1979) and of Mosser et al. (2015). The second section explains how one can build a diagnostic of the location of the mixed mode frequencies in a power spectrum derived from the Fourier transform of a stellar light curve (or from stellar radial velocity). The third section details how the diagnostic can be optimised for finding the main characteristics of the parameters describing the location of the frequencies of the dipole mixed modes. The next three sections provides the results for about 44 subgiant and a few early red giants using the optimisation procedure and the fitting of available frequencies. I finally discuss the scientific implications and conclude.

Asymptotic mixed dipole mode frequency
The approach to the derivation of the asymptotic mixed mode frequency is based upon the work of (Mosser et al. 2015). The methodology is repeated here for completeness. The frequencies of the mixed dipole modes is asymptotic because the frequencies of the pressure modes and the gravity modes are assumed to follow the asymptotic description of Tassoul (1980). A mixed mode has by definition a dual character of being both a pressure mode and a gravity mode. The main equation describing the mixed mode character is related to the following continuity equation provided by Shibahashi (1979): 6. Compute difference of phase tangent as: δ(ν) = tan θ p − q tan θ g (Eq. 8) 7. Compute sign of δ. Since the derivative of δ(ν) wrt to ν is always positive, δ(ν) is an increasing function of ν, then the only transition through zero is when δ changes sign from negative to positive. When it occurs, it means that it crossed the zero line, i.e. we have a solution for a mixed-mode frequency. 8. The final detection of the location of mixed-mode frequency is done by computing the first difference of the sign of δ.
When there is a crossing from negative to positive the first difference is then 2 (+1-(-1)).
This pseudo code can easily be implemented in any language.

Collapsing of mixed-mode peak power
The identification of dipole mixed modes in subgiant stars rely on a visual assessment of the location of the mixed modes.
As mentioned earlier, the use of replicated echelle diagram (Bedding 2012) or just plain experience of a regular or irregular echelle diagram is sufficient for most application, i.e for tens of stars. For the PLATO mission, I had to devise an automated procedure that would make the identification more immune to human intervention. The novelty of the method developed in this section lies in the use of all the information available in the power spectrum for the monopole and dipole modes coupled with the asymptotic description introduced in the previous section. The information used allows to derive a figure of merit that is related to the location of the dipole mode frequencies in the power spectrum. The idea behind the derivation of the guess frequencies for the dipole mixed modes is based upon the visualisation of thé echelle diagram (Grec 1981) applied to a power spectrum of stellar time series. If one co-adds the power along each segment line of theéchelle, one obtains the so-called collapsed power. An example of collapsed power is shown on Fig. 1 where one can see two main doublets: one doublet for the l = 0 − 2 modes located at 0 µHz; one doublet for the l = 1 − 3 modes located at +63 µHz. The other patches of power located at -35 µHz and +18 µHz seen in Fig. 1 are due to the l = 4 modes and the l = 5 modes, respectively. Fig. 1 is computed using data from the Luminosity Oscillations Imager (LOI) seeing the Sun as a star (Appourchaux et al. 1997). The additional sensitivity to the higher degree modes are due to the LOI pixels resolving the surface of the Sun (ibidem). Fig. 2 shows the collapsed for a KIC11137075, a star with l = 1 mixed modes already analysed by Tian et al. (2015). It is obvious that the power of the l = 1 mixed modes is not concentrated at all compared to the l = 0−2 mode pairs. The idea is then to compute l = 1 mixed mode frequencies resulting from solving Eq. (8) and to collapse the power around these frequencies.
The power is then extracted around an l = 1 mixed frequency using a window of the size of the large separation (±∆ν/2), then the power is summed up over all l = 1 mixed modes found from solving Eq. (8). Here the total power is not averaged over the number of mixed modes because the division would prevent to have a large total power when more mixed modes correspond to real peak of power. In addition, since I am interested in the location of the dipole mode frequencies, I compensate for the intrinsic fall out of mode power away from ν max by dividing the original power spectrum by the inverse of a Gaussian mode envelope whose full width at half maximum is 2/3 of the full mode frequency range defined as [ν low ,ν high ]. This compensation per- Fig. 1. Superposed power for the LOI data: power averaged over 15 modes centered on the l=0 modes: (top, green) for a fixed large separation of 136 µHz; (bottom, black) for an optimised large separation of 134.7 µHz linearly dependent upon frequency. The original power spectrum, smoothed to 10 bins, is taken from one year of data of the Luminosity Oscillations Imager seeing the Sun as a star (Appourchaux et al. 1997). mits to give a lower weight to the modes close to ν max , and a higher weight to the modes on the edge of the p-mode range.

Optimisation of collapsed power
The computation of collapsed power is performed as explained in the previous section. Then I need a way to optimise the location of the l = 0 mode frequencies and of the l = 1 mode frequencies. The next idea is to compute the mean collapsed power within 1 µHz of the center located at 0 µHz. I use this mean power as a figure of merit which is then maximised. The procedure is done in the following steps: 1. Get a stellar power spectrum, smoothed it over 10 bins 2. Select a frequency range for the p-mode power as [ν low ,ν high ] 3. Divide the power spectrum by a proxy of the Gaussian envelope. 4. Optimise for the l = 0 mode frequencies: (a) Set a starting value for ∆ν, ǫ p and α for the l = 0 mode frequencies (Eq. 4) (b) Compute collapsed power and get the figure of merit for l = 0. (c) Optimise (∆ν, ǫ p and α) for getting the highest figure of merit for l = 0. 5. Optimise for the l = 1 mode frequencies: (a) Use the previously determined (∆ν, ǫ p and α) as fixed value. (b) Set d 01 for getting the l = 1 p-mode frequencies (Eq. 3) (c) Set starting values for P 1 and ǫ g for getting the l = 1 g-mode frequencies (Eq. 7) (d) Set starting value for q (e) Get l = 1 mixed mode frequencies from solving Eq. (8) (f) Compute collapsed power and get the figure of merit for l = 1. (g) Optimise (d 01 , P 1 , ǫ g , q) for getting the highest figure of merit for l = 1.
The maximisation of the l = 0 figure of merit (Step 2c) is performed using a regular steepest-descent method since it is clear Table 1. Table of optimised parameters. First column: Kepler Identification Catalog (KIC) number; second column: duration of observation; third column: percentage of mode within ±3 microHz of the asymptotic dipole frequency; fourth column: source of available fitted frequencies; fifth column: frequency of maximum poser; sixth column: large separation; seventh column: parabolic coefficient of the l = 0 mode frequencies; eighth column: small separation between the l = 0 and the l = 1 mode frequencies; ninth column: g-mode phase; tenth column: p-mode phase; eleventh column: coupling factor; twelfth column: dipole period spacing. that there is only one minima that is very local and very close to the guess value of ∆ν, ǫ p and α.
The maximisation of the l = 1 figure of merit (Step 5g) is a lot more complicated. Regular algorithms based on steepest descent fails because of the many local extrema. An algorithm that is able to search for the extremum of local extrema must be used for finding a proper solution. Different algorithms were tried such as the following ones: -Random search -Gaussian process prediction -Genetic algorithm The random search is a brute force method which consist in doing 1 million evaluation of the l = 1 figure of merit taken from the four parameters taken at random in a hypercube for the 4 parameters to be searched. This algorithm was used for a similar problem involving solar g-mode detection (Appourchaux & Corbard 2019). The algorithm do find the global extremum but is very slow due to the large number of function evaluation.
Another approach tested was to use to Gaussian Processes (GP) for evaluating / forecasting the figure of merit in the hypercube (See Frazier 2018, for a tutorial). An instructive application of GP can be found in Jones et al. (1998). The application of GP is most useful when the figure of merit to be computed takes a large amount of time in the hyperspace. Instead of computing the figure of merit on a very fine grid of the hyperspace for finding the extremum, the idea is to forecast with a small number of The EI function can be much faster to compute on a fine grid of the hyperspace that the original figure of merit. Although this is a promising approach, the transfer of the optimisation problem from the figure of merit to EI works on a small number of parameters (say less than 3), but fails for 4 parameters because the EI has local maxima that are Dirac like which are then easily missed unless the computing grid is extremely fine. The original figure of merit has already large variation and does not slowly vary; in hyperspace the EI is then similar to an empty space with small hypersphere sparsely distributed.
The most useful algorithm for the current problem is based upon a genetic approach. The optimisation is based upon the Shuffled Complex Evolution (SCE) algorithm, which mixes evolutionary and simplex algorithm Duan et al. (1993). This algorithm allows to search globally for the best maximum, and then search locally for the highest maximum. This algorithm was used by Thi et al. (2010) for fitting absorption profile observed in the infrared with European Southern Observatory's Very Large Telescope. The fit employs an Interactive Data Language (IDL) software, that can be found at Thi's home page (www.astrowing.eu), which I use for the optimisation of the figure of merit. Owing to the random character of the genetic algorithm, I also choose to run the SCE in several steps as follows: = (P max 1 + P min 1 )/2, ∆ P1 = (P max 1 − P min 1 )/2 3. Optimisation 10 times using SCE with the following parameters: same range for d 01 , ǫ g , q, P 1 = [P mean This is more efficient than to run a million evaluation of the figure of merit since about it results in 10 times less evaluations in total. Figures 3, 4 and 5 show results of the optimisation for three typical cases. The (P 1 , q) maps show the difficulties for finding the maximum of the figure of merit. The typical cases are as follows: -For period spacing larger than 200-300 s: a maximum clearly located for q below 0.2. Lower local maxima are present but can clearly excluded (Fig 3) -For q larger than 0.2: Several local maxima located along parallel ridges along q. -For P 1 less than 200 s: several isolated local maxima.
The boundaries for these categories are by no means really strict but just rough indication. Some stars may fall in one category or the other while still being out of the boundaries defined. Still the three maps are typical of the various maps observed. From these maps, it is obvious that the optimisation is not a simple problem as many local minima may prevent to find the global minimum. Table 1 gives the parameters resulting from the optimisation for each star.

Results from optimisation: period spacing and coupling factor
I used two criteria that are useful to separate the red giants from the subgiants: the evolution criteria and the density of mixed modes. Mosser et al. (2014) provided a boundary for the evolution criteria that subgiants fulfill (∆ν/36.5µHz) (P 1 /126 s) > 1. The density of mixed modes is to the first order equivalent to the ratio of gravity-mode order to the pressure modeorder (N ∼ n g /n p ); Mosser et al. (2012b) gave this density as N = ∆ν/(P 1 ν 2 max ). Figure 12 shows the comparison with the measurement made using the asymptotic relation as in Mosser et al. (2014), and using avoided crossing as in Li et al. (2020) for the period spacing. I note that the period spacing as derived from the asymptotic relation of Shibahashi (1979) cannot be directly compared with the formulation used by Deheuvels & Michel (2011) based on avoided crossing for harmonic oscillators. For example, Benomar et al. (2013) and Li et al. (2020) provide several examples of echelle diagrams providing the location of the frequency of the gravity modes (γ modes) using the formalism of Deheuvels & Michel (2011). It is clear that the resulting frequencies of the mixed modes are greatly affected close to the frequency location of the gravity γ modes, i.e. the mixed modes are further away than the regular l = 1 ridge; on the other hand the mixed modes are closer to the original ridge in between the gravity γ mode frequencies. In the formalism of Shibahashi  (1979), the situation is opposite: close to a gravity mode, the mixed mode frequency is close to the original dipole mode frequency since θ g ∼ 0, then θ p ∼ 0; but when θ g ∼ π/2, then θ p ∼ π/2 which means that the mixed mode frequency can be anywhere within a ∆ν. In the formalism of Shibahashi (1979), the gravity modes have a g-mode phase offset of 1/2 compared to the gravity γ mode of Benomar et al. (2013). Figure 7 shows the comparison with the measurement made by Mosser et al. (2017) for the coupling factor (supplemented by data shown in Mosser et al. (2017) but not published; Mosser, 2020, private communication). The comparison shows very good agreement with the work of Mosser et al. (2017) and Mosser et al. (2014). The period spacings from the optimisation also agree with the results of Deheuvels et al. (2014). The increase of the coupling factor for subgiant stars is due to a very thin zone where mixed modes are evanescent (Pinçon et al. 2020). This is the strong coupling assumption leading high trans-mission between the gravity mode cavity and that of the pressure mode cavity (Takata 2016).
It is also quite clear that there is a simple relation between the evolution criteria and the density of mixed modes. Since ν max depends on the large separation, it can be shown that the density of mixed modes is inversely proportional to power of the evolution criteria; depending the state of evolution the exponent factor varies from -1 (for subgiants) to -1.5 (for red giants). A better relation can be derived but the end result is that the evolution criteria as devised by Mosser et al. (2014) is sufficient for understanding the complexity of the mixed-mode pattern.
6. Results from optimisation: d 01 and ǫ p Figure 8 shows the relative small separation obtained from the optimisation compared to previous measurements obtained by Benomar et al. (2013). The results are comparable with the    Fig. 9. ǫ p as a function of the effective temperature. (Black disks) From the optimisation procedure; (Blue disks) From the theoretical model of zero age MS stars given by Li et al. (2020). relative small separation provided by Stello (2012) for a 1 M ⊙ star. The small separation d 01 is much less powerful in terms of asteroseismic diagnostic compared to the small separation d 02 that can provide estimate or proxy of stellar ages (Lebreton & Montalbán 2009;Appourchaux et al. 2015). Figure 9 shows ǫ p as a function of effective temperature compared to a theoretical model from Li et al. (2020). This (ǫ p , T eff ) diagram is used for as a diagnostic for identifying the l = 0 − 2 modes vs the l = 1 modes for stars with high effective temperature for which the mode linewidth is larger than the small separation δ 02 (White et al. 2011). In the case of subgiant stars, the identification is eased because of the presence of mixed modes. The dependence of ǫ p with high effective temperature for the subgiant stars is similar to that of White et al. (2011).

Results: Asymptotic frequencies vs fitted frequencies
For 35 stars out of the 44 stars, we have fitted frequencies available from Campante et al. (2011) , Appourchaux et al. (2012b), Deheuvels et al. (2014), Tian et al. (2015) and Li et al. (2020). We use these frequencies for comparing with the location of the asymptotic dipole frequencies that could have been used as guess. For the remaining 9 stars, I visually checked whether a mode was present close to the asymptotic dipole frequency. Table 1 provides the percentage of frequencies (fitted or visual) that are within 3 µHz of the asymptotic frequencies. The automatic identification allows to retrieve within 3 µHz at least 80% of the modes for 32 stars, and.within 6 µHz at least 90% of the modes for 37 stars. Figure 10 shows the best and worst matching for stars with fitted frequencies. Figure 11 shows the best and worst matching for stars with visual frequencies. Figures 15 to  20 give the comparison for all the other stars with fitted frequencies.
It is known that the asymptotic approximation is not applicable to mixed modes in subgiants (Deheuvels et al. 2014). This is mainly due to the dipole mode frequency being of the same order of magnitude that of the Brunt-Vaïsala frequency. In addition, Ong & Basu (2020) mentioned that the asymptotic approximation fails for low n g and high n p which is the case when the density of mixed mode N is less than 1. For most of the stars, the agreement is quite remarkable (Left of Figure 10). Nevertheless a success rate of 80% is not satisfactory for an automatic guess determination. It means for instance that out of the 3000 subgiants from P1 sample of PLATO, 900 stars will have between 50% and 80% of their dipole mode frequencies properly fitted. 7. Results from least square fit: period spacing, coupling factor and d 01 The optimisation procedure does not provide a direct access to the error bars for each parameter. I could have derived the error bars from a Monte Carlo simulation of the optimisation process but this would have taken a vary long time; in addition it would have been not really pertinent because the original goal of the optimisation was to obtain guess frequencies for proper extraction of the dipole mode frequencies.
Therefore I choose to fit directly the asymptotic frequencies to the fitted frequencies. In that case, it is feasible to make Monte Carlo simulation of the fitting by using the error bars provided for the fitted frequencies by several authors. As outline previously, the asymptotic model is far from perfect. For several stars, the deviation of the asymptotic frequencies to the fitted frequencies can be very large amounting to more than 1000 σ ! In other words, the fit is plagued with systematic errors due to the model itself. Therefore, for the fitting the model using the fitted frequencies, I choose to minimise unweighted least square fit. For computing error bars on the 4 parameters defining the location of the dipole mixed modes, I ran 100 Monte-Carlo simulations of the least square fitting using the fitted frequencies and their error bars as input randomised input. For the least square fit, I also use the SCE algorithms with the following ranges: d 01 = [−0.2, 0.2], ǫ g = [0., 2.0],q=[0.,1.0], P 1 = [0.6P optim , 1.4P optim ]; where P optim is the period spacing from the optimisation procedure. In contrast, Buysschaert et al. (2016) used also non-linear least square fit but using a grid search on three parameters (ǫ g , q and P 1 ) Table 2 provides results of the fit with their error bars. Note that only the dipole mode frequencies were fitted in the process, the parameters defining the l = 0 mode frequencies were assumed to be the same as for the optimisation. Figure 12 provides the result of the fit compared to the optimisation for the coupling factor and the period spacing, which shows a good agreement between the two determinations. Figure 13 provides the errors from the Monte-Carlo simulation for the coupling factor and the period spacing. In this Fig. 10. Echelle diagram of the amplitude spectra with: dipole frequencies fitted on the power spectra (Red circles), l = 0 frequencies from the optimisation (Orange diamonds); asymptotic dipole frequencies from the optimisation (White diamonds) and dipole frequencies from fitting the asymptotic model on the fitted frequencies (Red diamonds). (Left) for a case with 100% of asymptotic frequencies matching fitted frequencies. (Right) for a case with 55% of asymptotic frequencies matching fitted frequencies. Fig. 11. Echelle diagram of the amplitude spectra with visual frequencies for the dipole modes (Orange circles) and the asymptotic dipole frequencies from the optimisation (White / Black crosses). (Left) for a case with 100% of asymptotic frequencies matching visual frequencies. (Right) for a case with 80% of asymptotic frequencies matching visual frequencies.. figure, I added the error bars derived by Mosser et al. (2014) which are rather large compared to the error bars derived from the Monte-Carlo simulation. A simple explanation to these large error bars is that the many solutions providing the period spacing can be mistreated as an error while it is a plain systematic error. Multiple solutions for the period spacing arise when for given g-mode frequency at n g , P 1 can be expressed such that: (n g + ǫ g )P 1 = (n g + ǫ g ± 1)(P 1 + ∆P 1 ) leading to ∆P 1 = ± P 1 n g + 1 + ǫ g I can replace (n g + 1 + ǫ g ) by P max /P 1 , where P max is the period corresponding to ν max . Then I have: This is a rather crude approximation since the next g-mode frequency (n g +1) should have also been taken into account. It just gives a simple explanation for multiple solutions. The same equation as Eq. (11) was given by Vrard et al. (2016) but with a different meaning. As a matter of fact, this equation does not provide an estimate of an error, as stated by Vrard et al. (2016), but the separation between multiple solutions, or pseudoperiodic solutions. Figures 3, 4 and 5 give good examples of the pseudo-periodic solutions that are roughly separated according to Eq. (11). As shown by Figure 13, the error derived by Mosser et al. (2014) are very close to the systematic error provided by Eq. (11). This explains why my error bars for the period spacing are about 30-50 times smaller than these systematics errors. Figure 14 shows the relative small separation obtained from the frequency fit compared to the optimisation and to previous measurements obtained by Benomar et al. (2013). The results with the fit are not very different from the optimisation procedure.

Discussion
Several alleys for improvement of the method are envisaged that use a combination of an improved asymptotic description together with detection scheme similar to that of Corsaro et al. (2020). The improvement should lead to an automatic identify-    Benomar et al. (2013) . ing guess with a success rate better than 99%. The application of the optimisation procedure for the l = 2 and l = 3 mixed modes can be envisaged in a very similar fashion.
The current computation of the figure of merit does not take into account the variation in mode height due to the mode inertia, as was observed by Benomar et al. (2014) for a couple of subgiant stars. The extension of the method to red giants would require to take into account the variation of the mode height with the dipole mode frequency. Dupret et al. (2009) showed that many dipole mixed modes are strongly attenuated, and would not contribute to any power in the figure of merit. Benomar et al. (2014) and Mosser et al. (2015) provided an asymptotic expression for the mixed mode inertia that can be used for both subgiant and red giants. Last but not least, the rotational splitting is not included in the current model but would be necessary as the rotation in red giants can be similar to the mixed mode spacing Mosser et al. (2015).

Conclusion
I presented a new way of deriving the coupling factors and period spacing of dipole mixed modes. The results obtained agree with previous measurements made by Mosser et al. (2014Mosser et al. ( , 2017. I also obtained error bars on the parameters defining the dipole mixed modes using Monte-Carlo simulations. The automated procedure is far from being perfect as it can provide for 32 stars out of 44 stars, 80% of the dipole mode frequencies to within ± 3 µHz of the dipole peak power. The improvements of the method for subgiant stars, and the extension to red giant stars are discussed.  of fitted parameters. First column: KIC number; second column: relative small separation between the l = 0 and the l = 1 mode frequencies, and its error bars; third column: g-mode phase; and its error bar; fourth column: coupling factor and its error bar; fifth column: dipole period spacing and is error bar; sixth column: rms distance between the fit and the observations.