Revised Architecture and Two New Super-Earths in the HD 134606 Planetary System

Multiplanet systems exhibit a diversity of architectures that diverge from the solar system and contribute to the topic of exoplanet demographics. Radial velocity (RV) surveys form a crucial component of exoplanet surveys, as their long observational baselines allow for searches for more distant planetary orbits. This work provides a significantly revised architecture for the multiplanet system HD 134606 using both HARPS and UCLES RVs. We confirm the presence of previously reported planets b, c, and d with periods of 12.0897−0.0018+0.0019 , 58.947−0.054+0.056 , and 958.7−5.9+6.3 days and masses of 9.14−0.63+0.65 , 11.0 ± 1, and 44.5 ± 2.9 Earth masses, respectively, with the planet d orbit significantly revised to over double that originally reported. We report two newly detected super-Earths, e and f, with periods of 4.31943−0.00068+0.00075 and 26.9−0.017+0.019 days and masses of 2.31−0.35+0.36 and 5.52−0.73+0.74 Earth masses, respectively. In addition, we identify a linear trend in the RV time series, and the cause of this acceleration is deemed to be a newly detected massive companion with a very long orbital period. HD 134606 now displays four low-mass planets in a compact region near the star, one gas giant further out in the habitable zone, an additional companion in the outer regime, and a low-mass M dwarf stellar companion at large separation, making it an intriguing target for system formation/evolution studies. The location of planet d in the habitable zone proves to be an exciting candidate for future space-based direct imaging missions, whereas continued RV observations of this system are recommended for understanding the nature of the massive, long-period companion.


INTRODUCTION
Exoplanet discoveries have revealed a vast range of planetary architectures, many of which differ significantly from that seen in the solar system (Ford 2014;Li et al. Winn & Fabrycky 2015;Horner et al. 2020;Kane et al. 2021;Mishra et al. 2023a,b).The expansion of exoplanet discovery space into the outer regions of planetary systems has been strongly enabled by the improving precision and increasing baseline of radial velocity (RV) surveys (Fischer et al. 2016).Such surveys allowed increased RV sensitivities towards the larger semimajor axis regime and provided insights into the prevalence of Jupiter analogs (Wittenmyer et al. 2011(Wittenmyer et al. , 2020a;;Fulton et al. 2021;Rosenthal et al. 2021), indicating that, for orbits of 1-10 AU, ∼6% of solar-type stars harbor giant planets (Wittenmyer et al. 2016).Further out, although companion orbits cannot be fully resolved due to incomplete RV sampling, their presence can usually be hinted in the form of a linear or quadratic acceleration on top of already recovered RV signals.
At the other end of the spectrum, another benefit brought by some of these RV surveys is their capabilities of probing the low-mass short-period planetary regime if such surveys were conducted on a high cadence schedule with high precision spectrographs.Smaller planets induce lower RV semiamplitudes from the stellar reflex motion close to the noise floor.As such, these planets often require a lot more sampling over many orbital periods to gradually build up their statistical significance.With long running time and high precision, some of the RV surveys, or multiple ones combined, opened the door for discoveries of multiple planets with drastically different orbital properties within the same system.These multi-planet systems provide an essential contextual platform from which to evaluate formation scenarios within a statistical framework, and determine common versus rare outcomes of planetary architectures (Chambers et al. 1996;Lissauer et al. 2011;Tremaine & Dong 2012).Furthermore, planetary systems with bright, nearby host stars are essential targets for direct imaging surveys that enable spectroscopic studies of planetary atmospheres to be undertaken (Kane 2013;Clanton & Gaudi 2016;Kopparapu et al. 2018;Nielsen et al. 2019;Stark et al. 2020;Kane 2022;Quanz et al. 2022).Thus, the discovery of nearby multi-planet systems, particularly those that harbor planets of different characteristics, are of enormous value in determining their contribution to exoplanet demographical studies and the potential for significant follow-up opportunities.
The HD 134606 system presents itself as one such exciting target.The system was reported by Mayor et al. (2011, M11 hereafter) to host three exoplanet candidates using RV data acquired from the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph (Pepe et al. 2000).The planetary system was presented as part of a larger catalog of exoplanet candidates detected from HARPS observations.At the time, the three planets were reported to have orbital periods of 12, 59, and 459 days.Over 12 years since the original report, we dive back into the system equipped with a much extended HARPS observation baseline complemented by RVs from the University College London Echelle Spectrograph (UCLES) (Diego et al. 1990) to fully explore the system architecture and to check the validity of the originally reported candidates.Aliases and other false positives present in the analysis of RV data can often result in incorrect periodic signature extraction, and significantly impact the inferred system architecture (Dawson & Fabrycky 2010).In addition, it is likely that further planets may be present in the system that remain undetected due to the precision of the RV data (Kane et al. 2007;Laliotis et al. 2023;Newman et al. 2023).In these cases, a planet injection/recovery methodology may be applied to determine the limits to which currently known planets represent the true system inventory (Howard & Fulton 2016).
In this paper, we present a significantly revised architecture for the HD 134606 system, including a total of five planets and indications of a long-period additional companion.Section 2 provides descriptions of the various observations conducted and data sources used in our work, including RV, imaging, astrometry, and photometry.In Section 3, we present a full analysis of the system, including updated stellar parameters, a new RV model, and constraints on the additional companion within the system.Section 4 describes our study of stellar activity signals within the data, using both spectral activity indicators and photometry to explore the possibility of false positive signals within the RV data.Section 5 discusses the various periodic signals in the data, aliases, dynamical stability, and the prospects for direct imaging.Our results are summarized in Section 6, along with suggestions for future work.The original discovery of the system by M11 utilized 113 RV observations by the HARPS spectrograph spanning a total of 2,548 days from July 2004 to July 2011.Since then, HARPS continued its monitoring of HD 134606 until May 2017.The additional observations significantly increased the number of RVs for the star to 219 and extended the baseline coverage to 4,677 days, or roughly 13 years.The extended coverage of the HARPS dataset includes a fibre upgrade for the instrument that resulted in a discontinuous jump in the RV time series around June 2015 (Lo Curto et al. 2015).This full HARPS RV dataset was later re-analyzed and re-reduced as part of the new reduction for all the public HARPS spectra by Trifonov et al. (2020) using the SpEctrum Radial Velocity AnaLyser (SERVAL) pipeline (Zechmeister et al. 2018).The new reduction with the SERVAL pipeline corrected several systematics including nightly zero-point RVs and average intra-night drifts that slightly improves the precision of the HARPS RVs compared to those previously derived.The high observation cadence of HARPS and the improved precision thanks to the new reduction allowed for detection and characterization of possible short-period and lowsemiamplitude exoplanet signals.The newly reduced HARPS RVs were published in the HARPS RV database (Trifonov et al. 2020) and we make use of those data in our work.

UCLES
In addition to the HARPS data, we obtain RVs for HD 134606 as part of the Anglo-Australian Planet Search (AAPS) (Tinney et al. 2001) taken by the UCLES spectrograph mounted at the Anglo-Australian Telescope.AAPS is one of the longest running "legacy" RV surveys that provided a long temporal baseline for many targets in the program that allowed for characterization of very long orbital period giant planets.Due to the nature of the AAPS program that is optimized for sampling long period orbits, the cadence of UCLES RVs is low and is thus not useful for detecting short period planets in this case.Nevertheless, HD 134606 was observed by UCLES 66 times from April 1998 to July 2015 for a total of 6,315 days of coverage.The data reduction procedure can be found in Butler et al. (1996).Combined with the HARPS observations, an RV baseline of almost 7,000 days, or a little over 19 years in total was established that could potentially reveal the presence of previously undetected long period companions.We present a portion of the UCLES RVs in Table 1.

HARPS
A common problem many RV searches face is the confounding effect of stellar activity on the RV time series.Commonly, stellar rotation signals due to active regions or dark spots on the stellar surfaces may introduce RV semiamplitudes similar to those of small mass planets that would last for multiple stellar rotation periods (Robertson et al. 2015;Giles et al. 2017;Robertson et al. 2020) whereas stellar magnetic cycles induced by the stellar dynamo effect could exhibit long-term periodic signals resembling induced RV by long-period gas giant planets that may last from several years to decades (Meunier et al. 2010;Dumusque et al. 2011;Costes et al. 2021).There have been numerous cases where previously identified exoplanet candidates were later refuted because the RV signals were actually of stellar activity origins or at their harmonics (Robertson et al. 2014;Robertson & Mahadevan 2014;Robertson et al. 2015;Kane et al. 2016a;Lubin et al. 2021;Simpson et al. 2022).It is therefore of the utmost importance that any RV observations should have at least one accompanying stellar activity indicator to disentangle stellar activity induced RV variations from those due to genuine exoplanet candidates.
The re-reduction of the HARPS RVs from the Trifonov et al. (2020) HARPS RV database not only improves the RV precision, but also provides many activity indicators for each star that could be useful for singling out activity signals.In this work, we make use of all the available activity indicators from the database, namely, Hα index, chromatic index (CRX), differential line width (dLW), Na I D index (NaD1), Na II D index (NaD2), cross correlation function (CCF) bisector inverse slope span (BIS), CCF full width at half maximum (FWHM), and CCF line contrast.The latter three activity indicators were derived by the HARPS Data Reduction Software (DRS) whereas the rest of the indicators were obtained from the re-reduction using the SERVAL pipeline.Briefly, Hα, NaD1, and NaD2 indexes measure emissions in these spectral lines that are known to be sensitive to activity.CRX provides information of the wavelength dependence of the RVs extracted from each echelle order and dLW measures the variations in line widths of spectral lines.Details regarding the indicators from the SERVAL spectral analysis can be found in Zechmeister et al. (2018).
Unfortunately the HARPS RVs provided by Trifonov et al. (2020) do not include measurements of S HK index values for the spectra, which are known to correlate well with stellar chromospheric activity.In order to monitor stellar magnetic variability of HD 134606 via chromospheric emission in the Ca II H & K lines, we compute the Mt.Wilson S HK index (Vaughan et al. 1978;Duncan et al. 1991) for the HARPS spectra.S HK is essentially a ratio of flux in the H & K line cores to that in reference continuum bands on either side of the lines.We measure the calcium index using the 1D HARPS spectra as produced by the ESO data reduction pipeline 1 .S HK values are computed using the ACTIN 2 software package (Gomes da Silva et al. 2018Silva et al. , 2021)), and use the pyrhk extension to calibrate the index to the original Mt.Wilson scale.The full S HK time series for the HARPS spectra contains a number of extreme outliers, all of which come from spectra with low signal-to-noise (SNR) ratios.We have thus limited our analysis to spectra with SNR ≥ 25 in the spectral order containing the H & K lines.

UCLES
The limited spectral coverage of the UCLES spectrograph makes it impossible to provide simultaneous coverage of the Iodine region for precise extraction of RVs from the stellar spectra and the Ca II H & K region for derivation of stellar activity S-index measurement.In order to provide the RVs with an accompanying stellar activity indicator for the UCLES data, we measure the equivalent widths (EW) of Hα absorption lines in the UCLES spectra to detect variations in the EW as a proxy to stellar activity.We follow the similar approach of the Hα EW analysis described in Robertson et al. (2014) and Wittenmyer et al. (2017), using an automated algorithm for normalizing the continuum and identifying telluric contamination close to the Hα line.Although employing the Hα line profile as a stellar activity indicator for a G-type star is non-ideal, it has been shown by Gomes da Silva et al. (2022) that the use of the Hα line for activity studies can be optimized with a proper choice of bandwidth for extracting line profile changes.Gomes da Silva et al. (2022) pointed out that the usual choice of a broader 1.6 Å bandwith can sometimes lead to degradation of activity signals due to the inclusion of flux in the line wings, whereas a narrower 0.6 Å bandwidth for Hα can circumvent the issue and maximize the correlation between Hα and Ca II H&K.We therefore adopt the recommendation of using a 0.6 Å bandwidth for calculating the Hα EW to identify stellar activity signals in the UCLES data for HD 134606.

Speckle Imaging
In an effort to search for and/or rule out nearby stellar companions HD 134606, we carried out a speckle imaging observation using the Zorro instrument at Gemini South on July 24, 2021.Zorro is a dual-channel, dualplate-scale imager that is able to observe simultaneously in two bands obtaining diffraction limited images with a field-of-view (FOV) of 6.7 arcseconds (Scott et al. 2021).The simultaneous imaging is achieved by the use of a dichroic beamsplitter which splits collimated beam at the wavelength of 675 nm, with red channel light passed through and blue channel light reflected by it.The light in each channel is then recorded onto an electronmultiplying charged coupled device where two speckle patterns are simultaneously recorded using a sequence of 1000 short 60 ms exposures.These images were combined using Fourier analysis techniques, used to produce reconstructed speckle images, and examined for stellar companions (Horch et al. 2011;Howell et al. 2011;Scott et al. 2018).Detailed instrument description, calibration, and data reduction process can be found in Scott et al. (2021).At the time of our observation, the blue channel was not used for speckle imaging due to instrument alignment issue and we only obtained data from the red channel with filter of central wavelength and bandwidth of 832 nm and 40 nm, respectively.

Astrometry
We take our absolute astrometry for HD 134606 from the Hipparcos-Gaia Catalog of Accelerations, or HGCA (Brandt 2018(Brandt , 2021)).The HGCA cross-calibrates Hipparcos (esa 1997;van Leeuwen 2007) and Gaia (Gaia Collaboration et al. 2021) astrometry onto a common reference frame with recalibrated uncertainties suitable for orbit fitting.The HGCA provides a proper motion of HD 134606 from Hipparcos near Jyr 1991, a proper motion from Gaia near Jyr 2016, and a long-term proper motion given by the positional difference between Hipparcos and Gaia divided by their time baseline.Discrepancies between these three proper motions signify acceleration in an inertial reference frame.
The HGCA astrometry of HD 134606 suggests nonlinear motion, albeit at low significance, with a χ 2 value of 8.4 for a model of constant proper motion.This is primarily due to a slight disagreement between the Gaia and the long-term proper motions, by far the most precise of the three measurements.Section 3.5 derives the implications of the absolute astrometry for the companions to HD 134606 and their orbits. 2.5.Photometry

TESS
We utilize time-series photometry obtained by the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) in order to search for the presence of transiting exoplanets and photometric stellar activity.HD 134606 was observed by TESS during Sectors 12, 38, and 39.We obtain 2-min cadence TESS light curve photometry from the Mikulski Archive for Space Telescopes (MAST).The light curve products include both the simple aperture photometry (SAP) and pre-search data conditioning simple aperture photometry (PDC-SAP), which were processed by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016).We additionally remove any data that are flagged as being poor in quality or that are > 5σ outliers relative to the root mean square (RMS) of the light curve.

ASAS-3
Long term photometric monitoring allows identification of long period transiting planets as well as long period activity cycles.To check for any long term variation from photometry, we obtain the All Sky Automated Survey -3 (ASAS-3) V band data using the online catalogue portal.We pick data taken by aperture 3 given they yield the lowest scatter.Any data that are not of the best quality (quality grade A) are filtered out along with any outliers that are rejected after applying a 4σ clip.The final ASAS-3 dataset results in 479 data points spanning from 2,452,441 to 2,455,111 BJD for a duration of 2,670 days, or about 7.3 years.

Stellar Parameters
HD 134606 is a G6IV spectral type star and a previous publication on its age and radius indicate the possibility of it being near the end of its main sequence stage (Takeda et al. 2007).However, the metalrich nature of the star extends its main sequence lifetime to ∼12.5 Gyrs, as verified via isochrones from the MESA Isochrones & Stellar Tracks (MIST) (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015;;Choi et al. 2016;Dotter 2016).Here we derive new and updated stellar parameters for the HD 134606 system.We start by providing one of the high signal-to-noise HARPS spectra as an input to SpecMatch-Emp in order to constrain fundamental stellar parameters such as stellar effective temperature (T eff ) and metallicity ([Fe/H]).Details of how SpecMatch-Emp works can be found in Yee et al. (2017).In short, the provided spectrum is calibrated and matched against every other star available in the library and SpecMatch-Emp picks the five closest matches according to chi-squared statistics.A new best matching spectrum is then constructed using the five selected library spectra through linear combination.The set of coefficients of linear combination that minimized chisquared when compared against the input spectrum is then chosen to output a weighted average of stellar parameters of the five selected stars along with their associated uncertainties.
Derived T eff and [Fe/H] values from SpecMatch-Emp are used as Gaussian priors for EXOFASTv2, which is used to derive a precise and consistent set of stellar parameters through modeling the spectral energy distribution of the star with archival broadband photometry from Gaia (Gaia Collaboration et al. 2018), Two Micron All Sky Survey (Skrutskie et al. 2006), and Wide-field Infrared Survey Explorer (Wright et al. 2010) in combination with the aforementioned MIST stellar evolutionary models, all of which are included within EXOFASTv2.Additional constraints provided to EXOFASTv2 include a Gaussian parallax prior from Gaia Data Release 3 (DR3) (Lindegren et al. 2021a;Gaia Collaboration et al. 2023) with a bias correction provided by Lindegren et al. (2021b), along with an upper limit on the V-band extinction using the galactic reddening maps from Schlafly & Finkbeiner (2011).We provide derived stellar parameters from the converged EXOFASTv2 fit in Table 3.

Radial Velocity Solutions
We utilize all of the available HARPS and UCLES data for our RV analysis.The joint dataset from the two sources provides an ideal opportunity to study both potential short and very long period signals at the same time thanks to the two long duration surveys contributing to a very long observation baseline as well as the high cadence observations conducted by HARPS.Special treatment is needed for both HARPS and UCLES data before the analysis.In June 2015, a new set of optical fibers was installed on HARPS, resulting a discontinuous jump in the RV time series (Lo Curto et al. 2015).To account for the velocity offset caused by the fiber upgrade, we treat the pre-upgrade (HARPS1) and postupgrade HARPS (HARPS2) data as two separate instruments.For UCLES, there appears to be an upward velocity offset in the RVs around 2,455,500 BJD.After analyzing the RV time series of other stars observed by UCLES in the AAPS program, the velocity shift at the same timestamp appears to be a common feature, suggesting the discontinuity in the RVs is probably not of astrophysical origin.However, there were no known upgrades or maintenance made to the spectrograph that would change the instrumental profile.The exact cause of the upward velocity shift in the RVs unfortunately is still unknown and is likely either instrumental or atmospheric.To err on the side of caution, we decide to treat UCLES data before and after 2,455,500 BJD as two separate instruments, UCLES1 and UCLES2, respectively.
We analysed four separate datasets from two instruments using RVSearch (Rosenthal et al. 2021), an iterative planet searching tool to search for periodic signals within the combined RV dataset.RVSearch fits sinusoids to different fixed periods within the defined period search grid and utilizes the change in Bayesian Information Criterion (∆BIC) to measure of the goodness-of-fit of models.If there is a signal returned by RVSearch having a significance above certain empirical false alarm probability (FAP) threshold, which in our case we chose to be 0.1%, then RVSearch carries out a maximum a posteriori (MAP) fit to refine the Keplerian model of that particular signal.If there are multiple signals in the datasets, ∆BIC is once again used to determine if an n+1-planet model should be preferred over a previous n-planet model.This search algorithm has the advantage over the traditional Lomb-Scargle periodogram in that we are able to fit for instrumental velocity offsets, stellar jitters, linear or quadratic trend as part of the fitting process.In addition, already searched signals or known planet parameters are free to vary during the search and fitting process of additional signals, allowing the overall model to reach a higher maximum likelihood and therefore more precise orbital parameter results returned by the search (Rosenthal et al. 2021).At each search iteration, the calculation of the FAP threshold follows a similar methodology as in Howard & Fulton (2016) where a linear fit is modeled to the periodogram histogram in log scale against its ∆BIC power values.The 0.1% empirical FAP is then obtained by extrapolating the linear fit to a ∆BIC value where only 0.1% of the periodogram values are beyond this threshold.We set the algorithm to search between a minimum period of 0.8 days and a maximum period of five times the observing baseline of our combined data, which is around 35,000 days, while allowing linear and quadratic trend to be included in the search.The result of the planet search is shown in Figure 1.The sub-panels are shown in the order of when each planet was detected by RVSearch rather than in order of orbital period.The search process successfully recovers two of the signals originally reported in M11: 12-day and 59-day, with FAP of each signal being 4.31×10 −16 and 5.49×10 −21 , respectively.The third signal in M11 appears to manifest as a new signal with over double the originally reported period of ∼960 days (FAP = 4.57×10 −27 ), suggesting the possibility that either M11 or our search algorithm here has found an alias rather than the true signal.Besides the three previously reported signals, the search detects two new low amplitude variations both with semiamplitude around 1 m s −1 , yet with great significance: one at ∼27 days with FAP = 2.63×10 −9 and the other at ∼4 days with FAP = 8.63×10 −8 .Some peaks appear to be significant around the 1-day periodicity but disappear after we fit for the strongest signal during each iteration, indicating those are 1-day aliases of the signals we recovered.On top of all recovered signals, a linear trend appears to be significant with a large amplitude, and is indicative of an additional long period signal associated with either a physical companion or a magnetic cycle.The significance of all returned signals show a monotonic upward trend according to the running periodogram in panel (m), especially during the HARPS observing period, implying that these signals are unlikely to arise from pure noise and are not only temporary and limited to specific seasons.
The MAP fit result for each returned signal from RVSearch is then passed to the RV modeling toolkit RadVel (Fulton et al. 2018) as an initial guess to sample Keplerian model posteriors and estimate parameter uncertainties using Markov Chain Monte Carlo (MCMC).We use P , T c , √ ecosω, √ esinω, and K as the fitting basis, where P , T c , e, ω, and K are orbital period, time of inferior conjunction, eccentricity, argument of periastron, and RV semiamplitude, respectively.All parameters including linear trend, quadratic trend, instrumental jitters, and instrumental offsets are allowed to vary freely during the fitting process except for eccentricity and semiamplitude where the e prior is forced to be less than 1 and K is forced to be positive.Once again, we treat the UCLES and HARPS datasets before and after the velocity offsets as different instruments.The MCMC chain is carried out with 8 independent ensembles in parallel where each has 80 walkers that can take up to 15,000 steps.Chain convergence is evaluated under the following criteria: Gelman-Rubin statistic (< 1.01), minimum autocorrelation time factor (> 40), maximum relative change in autocorrelation time (< 0.03), and minimum independet draws (> 1000).The chain successfully converged before the walkers could reach the maximum allowed steps and we present the 16%, 50%, and 84% solution along with the maximum likelihood result for both the fitted orbital parameters as well as derived physical parameters of the planetary candidates in Table 4.
UCLES data are extremely helpful in providing a longer observation baseline.However, it is worth noting that due to the limited precision of UCLES, the five signals recovered in Figure 1 are completely driven by the HARPS data.When running RVSearch on the UCLES data alone, no signal can be returned other than the linear trend.This is expected considering UCLES data exhibit lower precision and greater scatter: mean error of 1.60 and 1.74 m s −1 for UCLES1 and UCLES2, respectively, compared to 0.72 and 0.70 m s −1 for HARPS1 and HARPS2.In addition, UCLES suffers significantly higher instrumental jitters than the HARPS data, with jitters for UCLES1 around 2.98 m s −1 and UCLES2 around 5.20 m s −1 , whereas those for HARPS1 and HARPS2 are about 1.80 and 0.70 m s −1 , respectively.For that reason, we carry out another fitting process using only the HARPS1 data.HARPS2 data are not combined with HARPS1 data in this case to simplify the fitting since we are able to fit with two fewer free parameters at a cost of only five fewer data points.The MCMC chain converges quickly and we provide solutions of the 5-planet candidate model using the HARPS1 data alone in addition to those from the combined dataset including all the UCLES and HARPS data in Table 4.The RV solutions and uncertainties from models utilizing all the data and only HARPS1 data are consistent with each other.Taking into consideration of similar model result but with much longer observation baseline, we opt for the solution with all the data as the model we employ for the rest of this work.Model comparison based on the Akaike Information Criterion (AIC) shows all five planetary candidate signals are favored in the final solution along with a linear trend, which bears a significance of about 17σ.However, there is no great statistical evidence for a quadratic trend in either of the solutions.Models with fewer free parameters containing planets less than five all have ∆AIC > 47, indicating these models are not supported by our data and are completed ruled out.

RV Injection-Recovery Test
The sensitivity of RV data gives information as to the possible presence of additional companions that may lurk in the system below the current detection threshold.Such an analysis is an important component of a complete system characterization (Wittenmyer et al. 2020b;Dalba et al. 2021;Li et al. 2021).Here, we characterize sensitivity of our RV data by carrying out injection-recovery tests, where we inject synthetic signals into the RV data and attempt to recover these injected signals using RVSearch.We inject 3,000 synthetic planets with orbital period and minimum mass drawn from loguniform distributions and orbital eccentricity from the beta distribution with parameters from Kipping (2013).The result of the injection-recovery tests are then used to compute the search completeness contour in the minimum mass versus semimajor axis space.As before, due to the higher precision HARPS data have, we carry out a separate test with only the HARPS1 data in addition to the one with all the HARPS and UCLES RVs and we show the results as two panels in Figure 2.
The five planetary candidates are shown in solid black dots and are all above the 50% completeness curve, with the candidate planet e at 4.32 days sitting closest to the curve given it has the lowest semiamplitude among all signals and being closest to the noise floor of our data.Both panels show a similar level of sensitivity to low amplitude signals thanks to the presence of HARPS data.The combined RV (right panel) show higher sensitivity to higher mass companions at longer periods compared to the HARPS only data (left panel) due to the longer baseline of UCLES data and their sampling optimized for long period companion search.In both cases, potential planets with similar or larger masses than the five in Table 4 appear to be recoverable from our data within 2 au but were not detected, indicating no more such planetary presence in the inner region of the system.According to the right panel, presence of a gas giant planet is very unlikely within 10 au orbital radius and the ability to recover very long period companions is fairly limited beyond the baseline of the combined dataset (∼19 years, or a∼7 au), where the detection probability drops below 50% for a Jupiter-mass planet at a little over 10 au.Therefore, the possibility of one very long orbital period companion hinted by our RV linear trend cannot be ruled out.

Speckle Imaging Constraint
Reconstructed images derived from the speckle imaging reduction process mentioned in section 2.3 are used to examine for presence of stellar companions around HD 134606.Similar to the process presented in Horch et al. (2011) and Howell et al. (2011), we estimate the location of a secondary peak under the assumption of the presence of a companion and modeled the fringe pattern using the position of the secondary.The signal of the tentative companion is then inspected by various metrics in the procedure and also manually to check whether the secondary peak is due to a genuine star or is because of a noise spike or cosmic ray.Only a true companion detec-

Li et al.
tion would be retained.For HD 134606, no companion is detected in our speckle imaging data using the red channel of the instrument, suggesting that either there is no companion within the FOV of the instrument, or maybe a companion is simply not bright enough to be detectable.In this case, it is useful to derive a relative detection limit to answer the question of how bright a potential companion could be while remaining below the detection sensitivity of the instrument.
To get a sense of the limiting magnitudes of our imaging data as a function of separation from the host star, we examine the distribution of all local maxima and minima in the background of the image by drawing five concentric annuli centered on the primary star, each with width of 0.2 ′′ centered at radii of 0.2, 0.4, 0.6, 0.8, and 1.0 arcsec.We treat each maximum as a potential stellar signal and compute the mean peak values within each annulus and the corresponding magnitude difference.Since the distributions of maxima and minima are essentially mirror images of each other with similar absolute mean and standard deviation values, we average the values from both maxima and minima when calculating the magnitude difference.The detection limit as a function of separation is then estimated using a conservative 5σ limit, which we define as being 5 times the standard deviation of the mean peak values plus the mean peak values in each annulus.The detection limit derived for each annulus in terms of instrumental magnitude difference (∆m) is used as a conservative upper limit above which stars should be detected, thus providing a constraint of the possible undetected low mass companions nearby.We show the detection limit curve of our speckle imaging observation from the red channel of the instrument in Figure 3.
The derived 5σ ∆m upper limit is then used to place upper limits on masses of potential companions as a function of separation similar to what has been done in previous works (Kane et al. 2014(Kane et al. , 2019;;Dalba et al. 2021).We first iterate through all the available stellar spectral types in the Pickles spectral library (Pickles 1998) and calculate respective luminosity ratios between the primary and the companion, and therefore, estimate instrumental magnitude differences from the given spectral type of the companion (∆m ′ ) with the Zorro 832 nm red band filter.The calculated ∆m ′ were then compared to the ∆m limits we derived from our speckle imaging observation at each separation and we pick the spectral type that yields the closest matching magnitude difference as the star that would yield the magnitude upper limit at that separation.For the selected spectral type of the companion at each location, we once again compute luminosity ratios and magnitude differences, this time in the apparent V-band magnitude (∆m V ).Using the known V-band magnitude of the primary star and the distance to the system, calculated absolute magnitude values M V are then used along with the mass-luminosity relations (Henry & McCarthy 1993) to yield the final upper companion mass limit as a function of angular separation using the speckle imaging data.The derived mass upper limit curve is shown in Figure 4.

Joint RV+Astrometry Constraints
HD 134606 has a Renormalised Unit Weight Error value around 1.022 from Gaia DR3, suggesting the single-star model is a good fit to the astrometric observations.However, such a binarity flag is only sensitive to binaries that contribute a significant fraction of the light and/or are on a period shorter than a few times the Gaia mission life.Therefore, additional large mass companions are possible at larger orbital separations even if Gaia indicates no binarity for the star.HD 134606 has a detectable RV acceleration and a detectable astrometric acceleration between Hipparcos and Gaia, though the astrometric acceleration is less than 3σ significant.The RVs themselves show no deviation from a linear trend over 19 years of monitoring (see section 3.2).For the following analysis, we therefore assume that the star's acceleration is constant over the entire time baseline monitored by UCLES, HARPS, Hipparcos, and Gaia.
The Hipparcos-Gaia Catalog of Accelerations includes HD 134606 as a low-significance astrometric accelerator; its χ 2 value with respect to constant proper motion is 8.4 with two degrees of freedom (2.4σ Gaussian equivalent The purple sensitivity curve for speckle imaging is an upper limit, while the blue curve is a measurement for joint RV and astrometry constraint.RV observational baseline is indicated by the vertical green dashed line.Any potential companion is likely to be either a high mass giant planet, brown dwarf, or a very low mass star.significance).We compute the difference in proper motion between Gaia and the Hipparcos-Gaia mean proper motion separately in R.A. and Decl.We divide this difference by one-half the time baseline between Hipparcos and Gaia, i.e., by 12.43 years for R.A. and by 12.45 years for Decl.We add the proper motions in R.A. and Decl. in quadrature and use standard propagation of uncertainties to obtain an astrometric acceleration of 5.4 ± 1.9 µas yr −2 .Finally, we use the measured Gaia parallax of 19.14 ± 0.08 mas to convert this proper motion acceleration to a physical value of 0.68 ± 0.24 m s −1 yr −1 .
The accelerations in astrometry and in RV together constrain the ratio M sec /ρ 2 , where ρ is the projected separation.We define r to be the physical separation of the star and its companion and ϕ to be the angle between this separation vector and the plane of the sky.We then combine a RV = GM sec r 2 sin ϕ (1) where a RV and a ast are the accelerations along and perpendicular to the line of sight, respectively, to obtain We do not use standard propagation of uncertainties because the signal-to-noise ratio on the astrometric acceleration is modest and the true distribution for GM sec /ρ 2 will be significantly non-Gaussian.Instead, we use Monte Carlo with Gaussian uncertainties on a RV and a ast to obtain (at 16/50/84% quantiles), or a 95% confidence interval of 12-68 M Jup /arcsec 2 .These constraints are shown as shaded intervals in Figure 4.

Photometric Variability
We search for photometric variability in the 2min TESS light curves following the variability analysis of Fetherolf et al. (2023), which we describe briefly here.Variability is searched on short-timescales (< 13 days) in the TESS light curves using a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) and auto-correlation function.
The Lomb-Scargle periodogram searched for periodic variability on timescales of 0.01-1.5 days and 1-13 days in each individual TESS Sector (12, 38, and 39) in both the SAP and PDCSAP photometry.Photometric variations in the SAP photometry are consistent with known timings of spacecraft thruster firing events, and these variations have been removed in the PDCSAP photometry.The PDCSAP photometry is consistent with being flat (normalized LS power < 0.01), with an RMS of 235, 180, and 190 ppm in Sectors 12, 38, and 39, respectively.Therefore, we do not detect photometric variability the TESS light curves.We also search for the presence of transiting planets by investigating the phase folding the concatenated TESS light-curve on the measured radial-velocity orbital periods for each planet in the HD 134606 system.However, we do not detect any evidence for transit events in the TESS light curves at a limit of 50 ppm.
We carry out a similar search for the ASAS-3 photometry.The scatter and precision of the ASAS-3 data are unfortunately too poor for transit detection of our planetary candidates here and are therefore not useful for this purpose.However, running the data through a Generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kürster 2009) does reveal signatures with ∼1200-day and ∼4600-day periods within the ASAS-3 photometry that may indicate a long term stellar activity variation, which we will discuss in section 4.

STELLAR ACTIVITY
An essential step of all RV analyses as mentioned in section 2.2 is to check whether any of the signals recovered from the studied RV time series are susceptible to contamination by stellar activity, as signals such as stellar rotation, magnetic cycles and their aliases are often misidentified as exoplanet candidates.Here, we use activity indicators based on activity-sensitive spectral lines from both the HARPS and UCLES spectra as well as photometry from TESS and ASAS to identify signals of stellar activity origin.

HARPS
We take the activity indicators from the HARPS database that includes all re-reduced RVs and their associated activity indicators, from both the SERVAL pipeline and ESO's default DRS pipeline (Trifonov et al. 2020).Provided HARPS indicators include Hα index, CRX, dLW, NaD1 index, NaD2 index, BIS, FWHM, and CCF line contrast.Details of these indicators can be found in section 2.2 of Zechmeister et al. (2018).In addition, we derived our own S HK index from the HARPS re-reduced spectra and calibrated to the Mt.Wilson scale (see section 2.2).The Ca II H&K measurement later became available as part of the second version of HARPS RVBank (Perdelwitz et al. 2023) and we use it to double check our result from the S HK .We once again only utilize HARPS1 activity data since we do not want to introduce an unwanted vertical offset between data coming from before and after the fiber upgrade, especially considering there are only 5 data points for HARPS2.We detrend activity indicators that show a significant trend (see section 5.1) before passing them through a GLS periodogram to uncover any periodic features.Search range is limited between 2 and 4,500 days and we show the result in Figure 5.For each indicator, we mark the 0.1% FAP threshold as indicated by the horizontal dashed line.Signals with powers extending above these horizontal lines have high significance and less than 0.1% FAP.Additionally, the five planetary candidate signals we recovered in section 3.2 are overlaid on each periodogram for reference.At first glance, there appears to be a peak occurring at 1-year periodicity in many of the indicators such as the Hα index, CRX, dLW, NaD2, and CCF contrast, of which the Hα, dLW, and NaD2 data provide great signal power.This is likely due to the consistent sampling during each observing season when the target was visible.At longer periods, S HK index, CRX, and dLW all have strong peaks below or close to 0.1% FAP, probably due to a stellar magnetic cycle.These peaks occur at periods comparable or longer than the HARPS1 baseline, which is just less than 4,000 days.Results of Ca II H&K from Perdelwitz et al. ( 2023) and our own S HK are consistent with each other, except the long-period weak exhibit a much weaker power, possibly due to different methods for flux extraction.One more activity signal shows up at a little over a 100-day period and is visible in Hα, dLW, and NaD1.The spectral window function (Dawson & Fabrycky 2010), which is useful in identifying artificial periodicities incurred by the sampling frequency, displays no strong peaks for HARPS other than the 1-year signal, confirming the origin of it being the consistent yearly observations.The window function power rises towards longer periods are due to the testing frequency range extending beyond the HARPS observation baseline, which is a common behavior.No other significant HARPS activity peaks are found and none of the strong signals in the periodogram coincide with any of the five planetary signals.

UCLES
For the UCLES data, S HK index was unobtainable due to the limited coverage of the UCLES spectrograph.We therefore utilize the Hα EW we derived with a Å bandwidth to best use Hα as a proxy to S-index for this G-type star.As mentioned in section 3.2, there appears to be a systematic velocity offset near 2,455,500 BJD in the RV time series.We therefore split the EW dataset into before and after the offset, the same way we treated the RVs, to avoid the offset leading to a spurious activity signal.We follow the same procedure as for the HARPS activity indicators and we show the pre-offset UCLES EW periodogram and spectral window function in Figure 5 with a baseline around 4,500 days along with other HARPS indicators.Two peak stand out with less than 0.1% FAP in the UCLES EW periodogram.The peak at ∼3,000-day periodicity appear to be indicative of the long-term magnetic cycle that some of the HARPS activity indicators detected, confirming its identity.The other peak detected by Hα EW at around 29 days, which is also evident in the UCLES window function, is likely due to the observing sampling frequency affected by the moon's orbital period.Additionally, the 365-day periodicity that is strong in many of the HARPS indicators is also apparent.Activity periodicity search for the UCLES post-offset Hα EW dataset reveals a negative result and all peaks detected in the pre-offset activity data do no overlap with the five candidate signals.
It is worth mentioning the choice of using a 0.6 Å bandwidth recommended by Gomes da Silva et al. ( 2022) for deriving our Hα EW time series values indeed yield better result.The signals at both 29-day and 3,000-day display stronger power in the periodogram with 0.6 Å bandwidth data compared to 1 Å data.With the use of a traditional 1.6 Å bandwidth, the two signals are not even recovered.In this case, the use of a wider bandwidth for extracting the Hα data provides a poor result, likely due to additional flux from the line wings degrading the signal as suggested by Gomes da Silva et al. (2022).The use of a 0.6 Å bandwidth works well here, but it remains to be seen whether all studies should employ this strategy or its usage should be determined case by case.

Photometry
Photometry is often useful as another diagnosis of stellar activity by looking at the photometric light curves and identifying variability (Meunier & Lagrange 2019;Fetherolf et al. 2023;Simpson et al. 2022).We study time series photometry from TESS and ASAS-3.TESS analysis in section 3.6 indicates no obvious variability present in the light curve and therefore no stellar rotation or other activity signals were detected in TESS photometry.The ASAS-3 data, however, appears to exhibit a long term cycle.Two peaks are detected, one around 1,200-day periodicity and the other at the upper end of the observing baseline.The higher peak has a longer period extending beyond the baseline of the ASAS-3 photometry and therefore the period cannot be precisely located.However, by extending the search range to longer periods allows us to broadly estimate the period of this signal of ∼4,600 days (Upper panel of Figure 6).Subtracting this signal removes the other peak at 1,200-day as well, suggesting the 1,200-day signal might be an alias of the longer period signal.Although the location of the higher peak is poorly defined, the origin of it could very likely be the same magnetic cycle that is detected in other spectroscopic activity indicators.Although the window function power rises again towards longer periods, this is again the same window function behavior as we saw earlier for the HARPS data where the testing range is going beyond the observing baseline.

Rotation Period
Neither spectral activity indicators nor photometry show any periodicities at short periods, indicating no short term activity signals persisting over the duration of the long observing baseline we have.Rotationallymodulated signals from starspots, however, are not coherent, and degrade over a few stellar rotation periods (Robertson et al. 2015(Robertson et al. , 2020;;Giles et al. 2017;Lubin et al. 2021).Therefore, the power in rotation signals are typically lost over a long period of time if spots do not consistently occur on the stellar surface.To double check for the presence of a stellar rotation signal, we divide the HARPS1 dataset into seasons and attempt to retrieve any short-term periodic signatures.UCLES data are not utilized here since their sparse sampling were more optimized for the long-period planet search.We follow a similar procedure as before and use GLS to search for periodicities within the range of each season, ∼200-250 days, for each one of the HARPS activity indicators as well as for ASAS photometry.None of the periodograms return significant peaks and therefore we conclude the stellar rotation signal was not detectable in our spectral and photometric data, possibly due to the reduced level of activity at a relatively senior stage the star is at in its lifetime right now.
Given the non-detection of stellar rotation in our data, we attempt to estimate the rotation period of HD 134606 using the pyrhk function as a final check for stellar activity and possible aliases.The estimation calculates the chromospheric rotation period using activity-rotation relations from Mamajek & Hillenbrand (2008), where the Rossby number is fitted against log R ′ HK for a sample of 169 solar-type stars.Rotation period is then derived from the Rossby number using convective turnover time relation as a function of B-V color (Noyes et al. 1984).For HD 134606, the time-averaged log R ′ HK is -5.1, a fairly quiet star, and the stellar rotation period is estimated to be 42.0 ± 3.9 days.This period and its aliases do not overlap with any of the five planetary candidates in Table 4.However, given log R ′ HK of the inactive stars (log R ′ HK < -5.0) in the sample from Mamajek & Hillenbrand ( 2008) is only very weakly correlated with Rossby number, the derived rotation period should not be taken with trust and is therefore only used for a precautionary check for RV false positives.

Nature of the Recovered Signals
Although the stellar activity checks we performed in section 4 do not yield any significant periodicities that coincide with any one of the five recovered signals from Table 4, there remains some questions regarding the validity of the five planets: 1) there are two seasons within HARPS CCF contrast or FWHM data that show a weak ∼30-day periodicity, indicating the need to further scrutinize our planet f due to its proximity to the lunar cycle period despite the activity signals having FAP greater than 0.1% in the periodograms; 2) M11 reported a period for planet d that is nearly half of what we have recovered, which suggests that either M11 or this work is seeing an alias of the actual signal; 3) planets b, c, e, and f appear to be near integer ratios of each other, which could potentially be a sign of aliasing; 4) there is a significant linear trend present in the HARPS dLW, CCF contrast, and FWHM data, therefore making the physical companion nature of the RV linear trend questionable.Given these concerns, we perform final checks within the dataset to further verify the significance or nature of our recovered planet signals.Since all five planetary signals are driven by the HARPS data, we utilize only the HARPS1 data here to simplify the checking procedure.
5.1.1.Planet f 's 26.9-day Signal Planet f's 26.9-day period is close to the 29-day lunar cycle.To make it more worrisome, the period almost overlaps with the 1-year alias of the lunar cycle, making its planetary nature highly suspicious.Given that only two of the spectral activity indicators show a weak signal that may be attributed to the lunar cycle in two of the seasons (from BJD 2,456,311.9to 2,456,516.5, and from BJD 2,454,879.89to BJD 2,455,025.57),we create a synthetic dataset consisting of white noise with timestamps and error bars for the HAPRS1 dataset.The velocity scatter is drawn from a Gaussian distribution with standard deviation being the quadrature sum of original data point errors and RMS of the five-planet fit from section 3.2.Then, we inject a 29-day signal into the seasons in question with similar error and scatter.A periodogram analysis is then run with RVSearch on the entire synthetic data set and we find the injected signal is not recoverable with < 0.1% FAP unless the variation amplitude is increased to ∼5 m s −1 .In all scenarios, there is no indication of an 1-year alias of the lunar cycle being produced without the 29-day signal being recovered first, if that 29-day cycle is recoverable at all.This suggests that the time sampling and the lunar cycle within the season in question is extremely unlikely to have produced the 26.9-day period.As an additional check to see if such a periodicity is localized, we remove all the RVs within these two seasons and run another periodogram search.In this case, all five planets are once again recovered with FAP being 1.19×10 −12 , 8.26×10 −9 , 1.10×10 −14 , 1.24×10 −6 , and 4.07×10 −5 for planets b, c, d, e, and f, respectively.These signal significance are still well below 0.1% FAP albeit with less power considering two seasons of data are subtracted.The failure to falsify the planetary nature of the 26.9day period combined with the fact we do not see the lunar cycle in any other HARPS activity indicators or seasons suggest such a signal is likely genuine.However, future high precision high cadence RV observations are recommended to support or refute this claim.

Planet d's True Periodicity
M11 reported planet d to have a period of 459.3 days, which is nearly half of the orbital period we have recovered in our analysis.However, we could only see a weak power of such a 459-day signal in our periodogram (see panel (f) in Figure 1) and the signal is completely gone after we fit for the 959-day signal.We attempt to recover the 459-day periodicity by running the periodogram search only on the epochs used by M11.However, similar result is obtained where only the peak corresponding to the 959-day signal is recovered with FAP = 9.52×10 −5 with a peak containing weak power at around 459 days.We further create a synthetic dataset with an injected 459-day signal with parameters from M11 using the HARPS1 timestamps and errors to check whether such a signal could be retrieved.The periodogram is able to successfully detect the injected signal, with a weaker power around 1000-day periodicity.Our analysis suggests that the 459-day and 959-day periodicities could indeed be aliases of each other, depending on which one is the true signal.Considering the 459day signal is recoverable if it is really present but none of our periodogram searches returned the 459-day signal, we therefore conclude the 959-day periodicity is the true period of this planet and what M11 recovered was an alias of the 959-day signal.Indeed, a model run with epoch and orbital parameters from M11 cannot reach MCMC chain convergence.Therefore, results derived from M11 are likely due to incomplete posterior sampling of orbital parameters.

Aliasing or Independent?
Orbital parameters we derived for the HD 134606 system show orbital periods of planets b, c, e, and f are close to integer ratios of each other.Such orbital configuration often times are targets of dynamical studies for their possibility of being in a resonance chain.Before such analysis can be conducted, one must rule out the case where they could be aliases of each other.We proceed by creating five new datasets, each dataset having one of the planetary signals subtracted from the original HARPS1 dataset using parameters from Table 4.These datasets are then fed into a RVSearch to see if signals that are not subtracted can still be recovered.Indeed, all planets other than the one subtracted for the particular dataset are recovered with FAP well below 0.1% FAP and we can safely assume the five signals in Table 4 are not aliases of each other.

Cause of the RV Linear Trend
As mentioned in section 3.2, a highly significant (17σ) linear trend is observed in the RV timeseries, which could be indicative of a very long period high mass companion.However, such a linear trend should be scrutinized when such a similar trend is also observed in activity indicators, which here is the case for HARPS dLW, CCF contrast, and FWHM.Both dLW and FWHM exhibit a positive trend while contrast shows a negative trend (see upper panels of Figure 7).The continued increase in line width and decrease in contrast point to an apparent long term drift in the line shape.One possible explanation is that a long-term magnetic cycle is driving such a change in the line shape.However, at the same time we do not see any change in the line skewness in the BIS indicator nor do we see a corresponding linear trend in the S-index, which usually is known to correlate well with magnetic cycles in G-type stars.Interestingly, numerous previous publications have observed similar long term drift in the HARPS CCF FWHM and contrast (Lo Curto et al. 2015;Benatti et al. 2017;Dumusque 2018;Barbato et al. 2019;Costes et al. 2021;Li et al. 2022) and such a drift has been attributed to the HARPS instrumental long-term de-focusing issue which was improved during the 2015 fibre upgrade.To remove the effect of the systematic drift, we apply two different polynomials below in Equations 6 and 7 to FWHM and contrast, respectively, based on relations derived in Costes et al. (2021) for G-type stars.The result can be seen in the bottom panels of Figure 7 where the linear trends in both FWHM and contrast are completely removed after applying the drift correction.No more trend features persist in the FWHM and contrast data and we conclude that the linear trend we observe in the RV timeseries is indeed due to a physical long orbital period companion rather than a magnetic cycle.As suggested by our joint astrometry and RV analysis in section 3.5, the companion is likely to be a high mass giant planet, brown dwarf, or a very low mass star.Our speckle imaging result in section 3.4 yields non-detection in the vicinity of the host star, ruling out presence of bright stellar companions within ∼1.2 ′′ of the star, or a little over 30 au of projected separation.The Washington Double Star Catalog however lists HD 134606 having a stellar companion of spectral type M3V at a projected angular separation of ∼57 ′′ that shares similar proper motion as the host star we are interested in.Assuming 57 ′′ being the minimum angular separation between HD 134606 and the M dwarf companion at a distance of 26.8 pc from Table 3, we estimate the minimum orbital distance between the pair to be ∼1500 au.At this distance, the induced RV line-of-sight acceleration due to the M dwarf with a mass of 0.37 M ⊙ for a M3V star (Pecaut & Mamajek 2013) using equation 1 can be estimated to be ∼0.03m s −1 yr −1 , or ∼8×10 −5 m s −1 d −1 , which is far smaller than the RV linear trend acceleration of 3.34×10 −3 m s −1 d −1 that we derived from Table 4. Consequently, the distant M dwarf companion itself is not capable of producing the observed RV linear trend, and it is extremely likely that an additional substellar companion is present in the outer regime of the HD 134606 system.Continued RV observations of this target are therefore needed over long term to resolve the orbit of this newly identified long period companion.

Dynamical Analysis
The inner four planets b, c, e, and f orbit closer to each other in the inner regime of the system and have orbital periods near integer ratios.To explore the dynamics of the system, we conduct an N-body simulation using the REBOUND package (Rein & Liu 2012).The system is initialized using our derived stellar mass value from Table 3 and the maximum likelihood orbital parameters for the planets from Table 4 with planetary orbits assumed to be coplanar.The system is integrated using the symplectic integrator WHFast (Rein & Tamayo 2015) with a time step of 1/20 of the inner most planet's orbital period, or 4.8 hours, consistent with the recommended step size from Duncan et al. (1998) to ensure proper time resolution orbit sampling.The integration duration is set to 10 million years and we record the eccentricity values for all the planets every 100 years.The system is dynamically stable over the duration of the simulation as shown by the eccentricity variation of planets in Figure 8.However, the inner four planets appear to be experiencing high levels of mutual interaction as evident by the mild eccentricity variation of planet c from near circular to up to 0.1 eccentricity.Planet b, e, and f exhibit higher amounts of eccentricity variation, with planet e showcasing variation up to 0.4.Given the low mass of planet e, its proximity to the G-type star, and the age of the system, the observed eccentricity variation is unlikely to have inherited from the formation of the planetary system since tidal circularization would have already dampened the orbital eccentricity of planet e.Therefore, the observed moderate eccentricity of planet e opens up a possible scenario where the planet's orbit may have been perturbed by a dynamical event either internal or external to the system.If planet e is indeed present, the 0.4 eccentricity variation could be the key to the intriguing dynamical history of this system.Based on our dynamical result, planet e exhibits orbital eccentricity of currently observed value of 0.2 or smaller for 58% of the sampled time steps, suggesting observation of the current eccentricity or smaller is not a low probability event.However, more high cadence follow-up observations are needed to further validate the presence of this planet and refine its orbital eccentricity (presently with 0.14 uncertainty) such that the current system configuration can be used for future dynamical studies with less ambiguity.
The orbital periods of the inner four planets are nearinteger ratios of each other.Specifically, planet e (∼4 days) and b (∼12 days)'s periods are near ratios of 1:3, and planet b (∼12 days), f (∼27 days), and c (∼59 days)'s periods are close to a 1:2:4 ratio.To study the dynamical interaction between these pairs, we conduct a separate N-body simulation for a duration of 1,000 years and record evolution of resonant angles for these pairs every 0.1 year.Unsurprisingly, given the amount of deviation from the perfect period integer ratios, none of the resonant angles exhibit evidence of libration and the four inner planets therefore are not in orbital resonance with each other, consistent with previous results that planets do not show much preference to be near mean-motion resonances (Fabrycky et al. 2014).

Direct Imaging Prospects
Future space-based direct imaging missions such as the Habitable Worlds Observatory will be targeting ter-  restrial exoplanets within the habitable zone (HZ) of nearby stars (Kasting et al. 1993;Kane & Gelino 2012;Kopparapu et al. 2013Kopparapu et al. , 2014;;Kane et al. 2016b;Hill et al. 2018Hill et al. , 2023)), opening doors to direct atmospheric retrieval of low mass planets as well as direct detection of such planets that are hard or impossible to be discovered using other detection methods.Here, we estimate the feasibility of directly imaging planets in the HD 134606 system using a first-order flux ratio estimation with telescope configuration of the Habitable Exoplanet Observatory (HabEx) as a proxy to the next generation Habitable Worlds Observatory as recommended by the Decadal Survey of Astronomy and Astrophysics 2020 (National Academies of Sciences, Engineering, and Medicine 2021).We follow the methodology presented in Li et al. (2021) that used the required 1×10 −10 contrast as the detection limit for the HabEx mission with starshade (Gaudi et al. 2020) and taking into consideration the transmittance profile of starshade for the flux ratio estimation.Detectability of planets in an imaging mission also depends on their angular separations.The inner working angle (IWA) of HabEx, which is the angular size of the starshade as seen from the telescope, is taken to be 70 mas.In addition, we employ a rela-tively optimistic angular separation limit, IWA 0.5 , which is defined as the angular radius from the starshade center where the transmittance is 50%.For HabEx, IWA 0.5 is estimated to be at 56.4 mas when averaged across all wavelengths and can serve as the minimum angular separation where exoplanet detections can be made (Gaudi et al. 2020;Li et al. 2021).Our flux ratio estimation assumes Lambertian reflectance models and selects relatively conservative geometric albedo values of 0.2 in the visible band for all the planets.Planetary radius values are estimated using mass-radius relationship from Chen & Kipping (2017).Using stellar and orbital parameters from Tables 3 and 4 along with an assumed inclination value of 60 • , we present the estimation in Figure 9.
The planetary masses and radii are adjusted according to the inclination value and the system is confirmed to be stable at this orbital inclination.Clearly, the inner four planets are too close to the star and are completely blocked by the starshade as indicated by the two shaded disks representing IWA and IWA 0.5 in the right panel.However, the outermost planet d appears to be bright enough and with its orbit wide enough to be detected by imaging, especially near maximum angular separation.Maximum brightness of this planet is achieved around ∼4.7×10 −9 planet-to-star flux ratio at locations when the planet is near the edge of the IWA (east of the star), making it a promising candidate for future direct imaging missions.The conclusion of direct imaging feasibility of planet d still stands if we assign even lower geometric albedo values and more edge-on inclination angles.Interestingly, as indicated by the left panel in Figure 9, over half of the planet's orbit is in the optimistic HZ (OHZ) and a significant portion of it in the conservative HZ (CHZ), the boundaries of which are defined in Kopparapu et al. (2013Kopparapu et al. ( , 2014)).Although planet d is unlikely to be habitable considering its relatively large mass and size, it could still serve as an interesting target for future studies of habitability of exomoons around giant planets in the HZ.As mentioned previously in section 5.1, the origin of the observed RV linear trend is due to a physical companion in a very long period orbit.According to Figure 4 from the combined RV and astrometry analysis, such a companion could either be a very large mass gas giant planet, a brown dwarf, or a very low mass stellar companion, depending on the orbital separation of such an object.Based on the angular separation and derived mass information from the RV and astrometry analysis, we estimate the imaging feasibility of such a companion as a function of orbital semimajor axis by creating a pool of synthetic samples.We provide such an estimate between 7.5 au and 30 au with a step of 0.5 au.At each semimajor axis step, we randomly draw 10,000 exoplanet samples with eccentricity values drawn from a beta distribution (Kipping 2013), orbital inclination from a uniform distribution of cosi, and with argument of periastron, longitude of periastron, and true anomaly all drawn from uniform distributions.For each drawn sample, on-sky angular projection is calculated and this value is passed into an interpolation function based on the median companion mass values as a function of angular separation from the RV and astrometry analysis in Figure 4 to yield an estimated companion mass value at this semimajor axis.The derived mass is then passed into mass-radius relationship from Chen & Kipping (2017) and the flux ratio calculation then followed the same procedure as above using the drawn orbital parameters with two different geometric albedo values A g : 0.1, and 0.5, to take into consideration that the companion could either be a gas giant or a brown dwarf.These values are consistent with albedo models for solar system objects (Roberge et al. 2017) and brown dwarfs (Marley et al. 1999).Since we are using a reflectance only model, we place an upper limit of 80 Jupiter mass on the mass of the samples drawn.The above process is repeated for all 10,000 samples at each semimajor axis location.If for certain samples the derived mass exceeded the upper limit, we discard such samples and redraw from the distributions such that all 10,000 samples are below the stellar mass limit.Finally, we obtain a distribution of 10,000 flux ratio values for each of the orbital locations from 7.5 to 30 au.At each semimajor axis location, we take 95% confidence interval along with the median value from the estimated flux ratio distribution as the brightness range of the companion as a function of different orbit sizes.The result is shown in Figure 10.The horizontal dashed and dotdashed lines represent the required 1×10 −10 and optimistic 4×10 −11 contrast limit of HabEx.At an orbital distance of ∼30 au, ∼50% of the samples could potentially be directly imaged assuming the required contrast limit and A g of 0.5.For A g of 0.1, the sample portion drops to < 20%.As expected, planet-to-star reflected light flux ratio of this companion greatly diminishes the farther away the possible location of the companion is from the star and the detectability of such a companion through direct imaging in the visible band is low unless it orbits closer to the star.However, as the expected companion mass increases as a function of separation (see Fig. 4), the companion would be massive enough at large separations (beyond 20-30 au) at which thermal emission would dominate and therefore be visible in the infrared wavelength.At even larger distances (beyond 60-70 au), the companion nature is likely to be a low mass star and would be self-luminous.

Planet Naming Convention
According to the policy adopted by NASA Exoplanet Archive regarding the inclusion of exoplanets in the archive, results of a study must be peer-reviewed and accepted for publication in literature.The HD 134606 system was originally reported to host three planets by M11.Although the M11 paper was not refereed and only a table of orbital parameters without further analysis was given, we respect the b, c, and d letter assignment to the 12-day, 59-day, and 959-day (originally 459 days by M11) planets.Because of this, despite the 4day and 27-day planets having shorter orbital periods, we assign them with letters e, and f, respectively, instead of reassigning designations for all planets according to their orbital distances from the star.This designation is consistent to the ordering presented in Table 4.

CONCLUSIONS
In this work, we revisited one of the multi-planet systems observed by HARPS.Assisted by the continued HARPS observations since the M11 paper, the addition of UCLES data, and the re-reduction of the HARPS data by Trifonov et al. (2020), we are able to significantly revise the architecture of this system with updated orbital parameters and newly discovered planets.In total, 285 RV observations spanning a total of over 19 years were utilized in the RV modeling process, resulting in a total of five planets detected in the RV time series with great statistical significance.Out of the five detected planets, we confirmed planets b and c from M11 and revised planet d's periodicity to 959 days and deemed the result from M11 was likely an aliasing issue Li et al.
due to insufficient posterior sampling of orbital parameters.Two new super-Earths e and f were discovered with periods of 4 and 27 days, respectively.Additionally, we identified a significant linear trend in the RV residual and confirmed its origin most likely to be a new sub-stellar companion.Major steps were undertaken in an attempt to falsify the planetary nature of the recovered signals through analysis of various stellar activity indicators and RV sampling.No stellar rotation signatures can be recovered from spectroscopy or photometry and the star was deemed quiet.None of the suspicious periodicities in the activity time series would rule out the likelihood of the five RV signals being planets.
Multiple detection methods were employed in this study.Although both photometry and imaging yielded negative results, photometry helped confirm the quiet nature of the star whereas imaging placed upper limit constraints on the mass of the newly identified companion in the outer regime of the planetary system.Using combined RV and astrometry analysis, we were able to derive mass estimates as a function of angular separation, which later was used to estimate the feasibility of directly imaging this distant companion causing the RV linear trend and astrometric acceleration.The HD 134606 multi-planet system is now characterized to host 5 planets, with the inner four being low mass super-Earths and mini-Neptunes orbiting in a relatively compact inner region, with another gas giant in the HZ.Further out, a sub-stellar companion resides in a very wide orbit unresolved by RVs and astrometry and an additional M-dwarf star at a sufficiently large distance that we see no evidence of it in the RVs and astrometry.The different characteristics of orbiting bodies within the system prompts interesting questions regarding the possible formation scenarios and the pathways it took to reach the current observed architecture.The eccentricity variation of the inner four planets may provide some clues on the evolution and recent dynamical history of the system.Dedicated long-term high-precision RV monitoring of HD 134606 is required to further verify the inner planets and resolve the orbital characteristics of the sub-stellar companion to hopefully get a glimpse of the true nature of this object.Direct imaging of this companion is certainly possible with either reflected light or thermal emission from this companion.Planet d in addition offers itself as an exciting target for future imaging missions thanks to its relatively large size, favorable planet-to-star flux ratio, and orbital location being in the habitable zone.HD 134606 provides plenty of excitement and opportunities for future studies and follow up observations.ACKNOWLEDGEMENTS The authors would like to thank the anonymous referee for the valuable feedback.Z. L. wishes to thank Mirek G. Brandt for the many ideas regarding the astrometric analysis of this target, Howard Isaacson for the discussion on RVs and stellar activities, and Alex Venner for the conversation on HARPS RVs and systematics.P. D. acknowledges support by a 51 Pegasi b Postdoctoral Fellowship from the Heising-Simons Foundation.This work is based in part on data acquired at the Anglo-Australian Telescope.We acknowledge the traditional custodians of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past and present.Dynamical simulations in this paper made use of the REBOUND code which is freely available at http://github.com/hannorein/rebound.This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.This research has also made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory.

Figure 1 .
Figure 1.RVSearch result of the HD 134606 system.Upper panels show the best fit along with residuals to the model.Individual fit to each of the five signals are shown in the lower panels along with their periodograms.Panel (m) shows running periodogram for the five signals.No more significant signal is found in the residual as shown in panel (n).
Note-Model based on the entire UCLES and HARPS data measures a linear trend of −0.00340 ± 0.00020 m s −1 d −1 with the maximum likelihood value of -0.0034 m s −1 d −1 .Linear trend measurement from the model based on the HARPS1 data alone provides a consistent value of −0.00333 ± 0.00018 m s −1 d −1 with the maximum likelihood value of -0.00335 m s −1 d −1 .ω values in each row of the table are those of the star, not of the planets.a The best-fit model we employ for the rest of this work.

Figure 2 .
Figure 2. RV completeness contour of HD 134606 from the injection-recovery tests.Tests conducted with only HARPS1 are shown in the left panel and those carried out with all the RV data are in the right panel.Blue dots are injected signals that were successfully recovered whereas red dots were not.The black line represents the 50% detection probability, with redder color showing lower probability.

Figure 3 .
Figure3.Speckle imaging detection limit in terms of instrumental magnitude difference as a function of separation from the primary star.No bright companion was detected.Only the red channel was used for observation as the blue channel suffered an instrument alignment issue.

Figure 4 .
Figure4.Detection sensitivity curves from speckle imaging and joint RV+astrometry in terms of Jupiter masses as a function of separation from the central star.The purple sensitivity curve for speckle imaging is an upper limit, while the blue curve is a measurement for joint RV and astrometry constraint.RV observational baseline is indicated by the vertical green dashed line.Any potential companion is likely to be either a high mass giant planet, brown dwarf, or a very low mass star.

Figure 5 .
Figure 5. GLS periodogram of activity indicators for HARPS and UCLES data.Window functions for each are in black.Horizontal dashed lines are 0.1% FAP and vertical solid lines represent period locations for the five planet candidates.No peaks in this periodogram coincide with any of the candidates' periods.

Figure 6 .
Figure 6.GLS periodogram of ASAS-3 photometry.Horizontal dashed lines are 0.1% FAP thresholds.Upper and middle panel shows the periodogram before and after the subtraction of the 4,600-day peak, respectively.Lower panel shows the window function of the photometry.

Figure 7 .
Figure 7. HARPS CCF FWHM and contrast data.Upper panels show the original time series of the two activity indicators and the lower panels show data that have been corrected for the systematic drift due to the instrument de-focusing issue.The linear trends seen in the original FWHM and contrast time series are completely instrumental.

Figure 8 .
Figure 8. Dynamical simulation result of the HD 134606 system with five planets.Eccentricities for each planet were recorded for a simulation duration of 10 million years.The inner planets experience significant eccentricity variations.

Figure 9 .
Figure 9. Direct imaging visibility of the five planets in the HD 134606 system.Left panel: top down view of the system with CHZ (bright green) and OHZ (dark green) overlaid.Right panel: sky projected view of the system with 60 • orbital inclination.IWA and IWA0.5 are represented by the bright and dark shaded disks, respectively.Orbits are color coded if they are above the 1×10 −10 flux ratio with 60 • inclination.Inner four planetary orbits are too compact to be shown with detail on this scale.

Figure 10 .
Figure10.Flux ratio estimate as a function of semi-major axis of the long period companion based on random sampling of orbital parameters.The horizontal dashed line and the dot-dashed line represent the required and optimistic contrast ratio limit of future space-based direct imaging missions.

Table 4 .
System parameters of the five planetary candidates in HD 134606.