© The Authors 2026

1 Introduction

µ Herculis (µ Her) is a nearby quadruple star system located at a distance of 8.3 pc. Its hierarchical architecture is defined by a wide A–BC configuration with a semi-major axis of ∼1000 au. The primary subsystem (A) has a semi-major axis of 20.4 au and contains the bright (V = 3.42) G5IV subgiant µ Her Aa (HD 161797; HIP 86974) and its M-dwarf companion, µ Her Ab. The secondary subsystem (BC) is a tight M-dwarf pair with a semi-major axis of 11.7 au, composed of primary component µ Her C (GJ 695 C) and secondary component µ Her B (GJ 695 B). Figure 1 shows the geometric architecture of the quadruple as derived in our analysis.

The sky position of µ Her Aa has been tracked since antiquity: its earliest digitised entry appears in Ptolemy’s Almagest(137 AD) (Lequeux 2014), and the first measurement usable in our analysis dates to 1750 (Lequeux 2014). This nearly three-century record spans a substantial fraction of the wide A–BC orbit and several complete BC revolutions, while the Aa–Ab pair has been monitored for a few decades with modern techniques. Together, these data capture all geometric aspects of the system’s motion. To recover the full architecture, we assemble three complementary tracers: relative astrometry maps sky-projected orbital ellipses, radial velocities resolve the line-of-sight component, and absolute astrometry anchors the trajectory on the sky. The combination of them allows us to reconstruct all relevant orbits in a single coherent framework that naturally accounts for proper motion, orbital motion, and parallax–perspective effects, yielding precise, model-independent dynamical masses under minimal assumptions.

Earlier efforts already hinted at this possibility. Heintz (1994) published orbital parameters for the Aa–Ab pair based on decades of astrometric monitoring of µ Her Aa, and although only his final solution is available, it broadly agrees with the modern picture. Despite statements in Heintz (1991) that the original measurements had been archived, extensive searches assisted by the Swarthmore Friends Historical Library¹ did not recover them. More recently, tentative Aa–Ab orbits were derived by Roberts et al. (2016), using radial velocities and relative astrometry, and by Feng et al. (2021), combining HIPPARCOS and Gaia data. Crucially, the system has since passed through periastron, providing the curvature needed to anchor the orbital geometry and refine the component masses to high precision. The BC subsystem, discovered in 1856 (Clark & Dawes 1857) and observed continuously thereafter, was last characterised by Prieur et al. (2014), who obtained a complete orbital solution but not individual masses.

µ Her Aa is of particular interest because it exhibits clear solar-like oscillations (Grundahl et al. 2017). With more than 120,000 radial-velocity (RV) measurements made possible by the robotic nature of the Hertzsprung Stellar Observations Network Group Telescope (SONG; Grundahl et al. 2017; Andersen et al. 2019), a detailed asteroseismic analysis that rivals those of the best-studied stars is underway. Combined with interferometric measurements of its radius (Baines et al. 2018) and the dynamical mass derived here, µ Her Aa will become one of the most precisely characterised stars known. The independently determined dynamical mass makes µ Her Aa an exceptionally clean benchmark for testing and calibrating asteroseismic scaling relations and stellar models.

Fig. 1 To-scale overview of the orbital configuration and naming scheme of the µ Her system. Components are colour-coded as Aa (green), Ab (red), B (yellow), and C (blue). The black curve traces the total sky motion of the system’s CM over a 40-year interval, combining a linear term (dashed grey line) with the annual parallax. Star symbols mark the positions on January 1 2026. Solid lines show the projected forward motion of each component over the next 40 years, and of the wide A–BC orbit over the next 400 years. To visualise the orbital inclination, segments closer to the observer are rendered in white and cast a shadow onto the sky plane; the transition between the white and grey segments marks the line of nodes. The inner A and BC subsystems are presented in greater detail in the right-hand panels.

2 Observations and data

For the Aa–Ab subsystem, we assembled five complementary datasets: (i) relative astrometry, (ii) primary-star radial velocities, (iii) absolute astrometry from HIPPARCOS, (iv) absolute astrometry from Gaia, and (v) ground-based absolute astrometry. For the B–C subsystem and the wide A–BC pair, we relied on (i) Gaia astrometry, (ii) published relative astrometry measurements, and (iii) a single spectroscopic measurement. To facilitate reproducibility, we compiled all observational data used in this work into a single, machine-readable archive².

2.1 Doppler data

To trace the reflex motion of µ Her Aa, we used high-cadence spectra of µ Her obtained over 12 years with the spectrograph at the Hertzsprung SONG Telescope operating at Observatorio del Teide. Originally acquired for asteroseismology, these data provide exceptionally dense coverage for orbital modelling. For the present analysis, we therefore compressed each night’s time series to a single RV measurement by taking the median of all exposures. This procedure averages over solar-like p-mode oscillations and mitigates intra-night instrumental drift (Chaplin et al. 2019). Uncertainties were estimated from the scaled median absolute deviation within each night. A detector upgrade in January 2024 increased SONG’s resolving power but introduced a zero-point shift. We therefore treated pre- and post-upgrade measurements as separate instruments with independent offsets and jitter terms in the model (see Sect. 3).

We supplemented the SONG data with archival RVs from the Hamilton and APF spectrographs at Lick Observatory (Rosenthal et al. 2021), from HIRES at Keck Observatory (Teklu et al. 2025), and from CORAVEL at Haute-Provence Observatory (Duquennoy et al. 1991). We further incorporated historical measurements from Küstner (1908), Lunt (1918), Wilson (1953), Abt (1973), Campbell (1913), Adams & Joy (1923), Spencer Jones (1928), Campbell (1928), Shajn & Albitzky (1932), Harper (1933), and Young (1939). Table 1 summarises the number of RV data points used in our analysis per facility, and their uncertainties. An RV offset was applied to each unique instrument.

Doppler coverage of the BC subsystem is significantly sparser. The sole literature value (Fouqué et al. 2018) is a blended measurement lacking an epoch designation, rendering it unsuitable for orbital modelling. To secure a resolved snapshot of the binary, we observed the subsystem on December 3 2025 (BJD 2461012.31) using the FIES spectrograph at the Nordic Optical Telescope (Djupvik & Andersen 2010). This observation successfully disentangled the components, yielding v_C = −11.33 ± 0.23 km s⁻¹ and v_B = −6.23 ± 0.23 km s⁻¹.

2.2 Relative astrometry

Historical measurements of the relative orbits of the A–BC and BC systems were drawn from the Washington Double Star Catalogue (Mason et al. 2013). For measurements for which formal uncertainties were not originally reported, we performed a preliminary fit to characterise the residuals for each dataset, and then adopted the root mean square of these residuals as the measurement uncertainty in the final analysis.

For the Aa–Ab pair, we adopted the archival measurements compiled by Roberts et al. (2016) and added two additional observations of our own. The first of these observations was obtained June 8, 2023 with Keck-II and its NIRC2 AO camera, using the 2.2 µm K_cont filter to obtain 56 short exposures (t_tot = 0.48 s) in a subarray containing the binary. The images were reduced and two-source PSF fits were measured following Kraus et al. (2016), where the best empirical PSF templates were adopted from the 1000 archival single-star K_cont images that were closest in time. The astrometry was corrected for optical distortion following Service et al. (2016). We obtained the second observation on September 17, 2025 with Robo-AO-2 Baranec et al. (2024), using its natural guide star adaptive optics mode. The target was observed in the H band for a total of 5 minutes using 3 ms exposures. Each exposure was corrected for detector bias and sky-background before sorting by image quality. Using the best 10% of frames, aligned to the image position of Aa, we created a final co-added image. We determined the stellar locations using the Aperture Photometry Tool.

2.3 Hipparcos

The HIPPARCOS one-dimensional abscissae epoch data are referred to as intermediate astrometric data (IAD) and exist in multiple versions and reductions, as is described in Brandt et al. (2021). We used the latest version, extracted from the 2014 Java tool³. Only the bright component µ Her Aa has HIPPARCOS data, and its HIPPARCOS ID is HIP 86974.

2.4 Gaia DR3

In Gaia DR3 (Gaia Collaboration 2023), three of the four components of µ Herculis are resolved and have astrometric solutions in the gaia_source table⁴: µ Her Aa (Gaia DR3 4594497769766809216), µ Her B (Gaia DR3 4594497834189213184), and µ Her C (Gaia DR3 4594497838483177984). Our use of these catalogue parameters, and of the Gaia Observation Forecast Tool (GOST)⁵ for scan angles and epochs, is described in Sect. 3.

2.5 Ground-based absolute positions

We also incorporated ground-based absolute astrometry for µ Her Aa. Modern milliarcsecond-level positions were taken from the U.S. Naval Observatory Bright-star Astrometric Database (UBAD; Munn et al. 2022). To extend the temporal baseline, we extracted additional positions from the Harvard DASCH project (Grindlay et al. 2012), computing median coordinates in 10 year bins to suppress plate-to-plate systematics. Finally, we included historical catalogued positions from Boss (1937), Dyson (1913), Argelander (1903), Lequeux (2014), Fricke et al. (1963), and Fricke et al. (1988) that help to trace the long-term motion of the system across the sky.

3 Methods

3.1 Overview

We modelled all spectroscopic and astrometric data simultaneously in a forward-modelling framework comprising four elements: the centre-of-mass (CM) motion of the system across the sky, and the Keplerian orbits of the Aa–Ab subsystem, the B–C subsystem, and the outer A–BC system. The CM motion is described by the five standard astrometric parameters (position, proper motion, and parallax), while the orbits add the reflex displacements of each star about their respective centres of mass. All of these dynamical components and their relative strengths are shown in Fig. 1.

As is outlined in Sect. 1, the radial velocities, relative astrometry, and absolute astrometry each constrain complementary views of the same dynamics. Our model combines all of them in a single joint fit. At each measurement epoch, our model generates the corresponding predicted value for the relevant data type. We quantified the agreement between predictions and data using a χ² statistic that incorporates both the values and uncertainties of the measurements. From this, we built the likelihood function, which was then explored with an MCMC sampler to map the parameter space and obtain the posterior probability distribution. The full list of model parameters and their definitions is given in Table 2. Our adopted parametrisation was chosen for its close mapping to the datasets.

3.2 Conventions and nomenclature

All astrometric quantities are expressed in the International Celestial Reference System (ICRS), and our model is anchored to the HIPPARCOS reference frame, meaning that we adopt t₀ = 1991.25 jyear, and define the position offsets, $\text{Δ}{\alpha }_0^{*}$ and ∆δ₀, relative to the published HIPPARCOS reference position (α_H, δ_H); see Table 3.

When an orbit is viewed through astrometry alone, it is projected onto the plane of the sky. This flattening erases the distinction between the near and far sides of the orbit, producing a (Ω, ω)↔(Ω + π, ω + π) degeneracy. Radial velocity data provide the missing line-of-sight information, breaking this degeneracy and allowing the full three-dimensional geometry of the orbit to be determined. For purely astrometric or purely spectroscopic studies, the precise definitions of Ω and ω often matter little, since the degeneracy prevents a unique solution in any case. In our joint analysis, however, this degeneracy is finally lifted, making it crucial to state our conventions unambiguously. Unfortunately, the literature is inconsistent in how Ω and ω are defined (Householder & Weiss 2022), and if left unspecified, the very advantage gained by combining RVs and astrometry could be undone.

To avoid confusion, we explicitly state our conventions, matching RV convention 2 of Householder & Weiss (2022):

• Co-ordinate system: a right-handed sky-plane frame ($\widehat{X}$, Ŷ, $\widehat{Z}$) with $\widehat{X}$ pointing east (along +α^∗), Ŷ pointing north (along +δ), and $\widehat{Z}$ directed away from the observer.

• Radial velocity model: the RV model uses the argument of periastron of the star being modelled, rather than that of its companion.

• Ascending node: defined as the point where the orbit crosses Z = 0 while moving away from the observer.

• Longitude of the ascending node: Ω is measured as a position angle (from north through east) to the line of nodes.

• Scan angle: we adopt the Gaia convention for the scan angle, ψ, measured from north to east (Holl et al. 2023). HIPPARCOS uses a different definition, which we converted via ψ = π − ψ_hip (Brandt et al. 2021).

• Orbits: the relative orbit refers to the ellipse traced by the secondary around the primary, with the primary centred at (0, 0). We refer to the barycentre of the Aa–Ab orbit as A, and the barycentre of the B–C orbit as BC.

• Subscripts: unsubscripted quantities (a, ω, K, M) refer to the relative orbit or total system values, as appropriate. Subscripts 1 and 2 refer to the primary and secondary, respectively. Using these definitions, the relative orbit is in phase with the barycentric orbit of the secondary, such that ω = ω₂ = ω₁ + π.

• Astrometric notation: p denotes the vector of the five astro-metric parameters. Superscripts indicate the reference frame (p^H for HIPPARCOS, p^G for Gaia, and p^H←G for Gaia rotated into the HIPPARCOS frame). Subscripts denote the entity (e.g. p_CM for the CM). A hat, $\widehat{p}$, denotes a model-implied ‘Gaia-like’ solution.

3.3 Stellar mass derivation

To express the individual component masses in terms of our fitted parameters, we start from Kepler’s third law in astronomical units and years: $$M=\frac{a^3}{p^2},$$(1)

where M = M₁ + M₂ is the total mass of the binary. The barycentric semi-major axes a₁ and a₂ satisfy a = a₁ + a₂, and the CM condition M₁a₁ = M₂a₂. These relations imply $$\begin{array}{cc}{M}_1=M\frac{a_2}{a},& M_2=M\frac{a_1}{a}.\end{array}$$(2)

Substituting M from Eq. (1) and eliminating a₂ via a₂ = a − a₁, we obtain $$\begin{array}{cc}{M}_2=\frac{a_1{a}^2}{p^2},& M_1=\frac{(a-a_1)a^2}{p^2},\end{array}$$(3)

so that both component masses are written solely in terms of the fitted quantities a, P, and a₁.

3.4 Forward model

The following subsections detail the construction of the likelihood for each data type. To keep this overview concise, we defer second-order effects such as perspective acceleration, frame rotation, and the distinction between the photocentre and primary motion to Appendix A.

We begin with the motion of the system’s CM, expressed in the tangent plane as $$\left[\begin{array}{c}\Delta {\alpha }_{\text{CM}}^{\text{*}}(t)\\ \Delta {\delta }_{\text{CM}}(t)\end{array}\right]=\left[\begin{array}{c}\Delta {\alpha }_{0,\text{CM}}^{*}\\ \Delta {\delta }_{0,\text{CM}}\end{array}\right]+\left[\begin{array}{c}{\mu }_{{\alpha }^{*},CM}\\ {\mu }_{\delta ,CM}\end{array}\right](t-t_0)+\left[\begin{array}{c}{\Pi }_{{\alpha }^{*}}(t)\\ {\Pi }_{\delta }(t)\end{array}\right]{\varpi }_{\text{CM}},$$(4)

where ${\text{Π}}_{{\alpha }^{*}}(t)$ and Π_δ(t) are dimensionless parallax factors describing the geometric shift due to the observatory’s changing vantage point. They trace the shape of the annual parallax ellipse in the local (α^∗, δ) basis and depend on the observatory ephemerides and the target’s sky position. The overall amplitude of this apparent motion is set by the parallax ϖ_CM.

We model each star’s sky-projected position as the sum of this CM path and two Keplerian reflex motions: one due to the orbit between the A and BC barycentres (the wide A–BC orbit), and one due to the corresponding inner orbit. To describe the orbital geometry, we use the standard Campbell elements {P, τ, e, ω, i, Ω, a}. A simple approach is to first model the relative orbit and later scale it to obtain the barycentric orbits of the primary and secondary. The instantaneous linear displacement due to orbital motion in the two tangent-plane directions and along the line of sight is (Leclerc et al. 2023) $$\left[\begin{array}{c}\Delta {\alpha }_{\text{orb}}^{*}(t)\\ \Delta {\delta }_{\text{orb}}(t)\\ \Delta Z_{\text{orb}}(t)\end{array}\right]=D\left[\begin{array}{c}\cos (\nu +\omega )\sin \Omega +\sin (\nu +\omega )\cos \Omega \cos i\\ \cos (\nu +\omega )\cos \Omega -\sin (\nu +\omega )\sin \Omega \cos i\\ \sin (\nu +\omega )\sin i\end{array}\right]\text{\hspace{0.17em},}$$(5)

with $$D=\frac{a(1-e^2)}{1+e\cos \nu (t)}.$$(6)

The true anomaly, ν(t), is obtained by solving Kepler’s equation, $$\frac{2\pi }{P}(t-\tau )=E(t)-e\sin E(t).$$(7)

for the eccentric anomaly, E(t), and then converting via $$\nu (t)=2\arctan \left[\sqrt{\frac{1+e}{1-e}}\tan \frac{E(t)}{2}\right]\text{.}$$

To obtain the secondary’s position, we simply scale the relative co-ordinates of Eq. (5) by a₂/a. The primary’s position is given by the same co-ordinates scaled by −a₁/a, where the minus sign reflects the 180° phase offset between the primary and relative orbits (ω₁ = ω + π). We are now ready to assemble the full CM-plus-orbital sky path of the Aa component: $$\left[\begin{array}{c}\Delta {\alpha }_{Aa}^{*}(t)\\ \Delta {\delta }_{Aa}(t)\end{array}\right]=\left[\begin{array}{c}\Delta {\alpha }_{CM}^{*}(t)\\ \Delta {\delta }_{CM}(t)\end{array}\right]-\frac{1}{{\varpi }_{CM}}\left(\frac{a_{1,A}}{a_A}\left[\begin{array}{c}\Delta {\alpha }_{orb,A}^{*}(t)\\ \Delta {\delta }_{orb,A}(t)\end{array}\right]+\frac{a_{1,A-BC}}{a_{A-BC}}\left[\begin{array}{c}\Delta {\alpha }_{orb,A-BC}^{*}(t)\\ \Delta {\delta }_{orb,A-BC}(t)\end{array}\right]\right).$$(8)

Here, the Aa component is the primary of the Aa–Ab subsystem, and A is the primary of the A–BC subsystem; hence the orbital co-ordinates are scaled by −a_1,A/a_A and −a_1,A–BC/a_A–BC, respectively. The factor ${\varpi }_{\text{CM}}^{-1}$ converts the orbits from linear to angular size. By explicitly separating the CM motion from the orbital motion, we obtained an independent parallax for the CM, rather than relying entirely on Gaia and HIPPARCOS catalogue solutions in which unmodelled orbital curvature can leak into the fit and bias the parallax.

To model HIPPARCOS and Gaia data, we project the sky-plane displacement in Eq. (8) onto the along-scan direction, ψ: $$w(t)=\sin \psi \Delta {\alpha }_{\text{Aa}}^{*}(t)+\cos \psi \Delta {\delta }_{\text{Aa}}(t).$$(9)

The corresponding projected parallax factor is ${\text{Π}}_{\psi }={\text{Π}}_{{\alpha }^{*}}\sin \psi +{\text{Π}}_{\delta }\cos \psi$. Finally, to model spectroscopic data, we require the radial velocity, i.e. the time derivative of the orbital Z-displacement, Eq. (5): $$v_x(t)=K_x[\cos (\nu (t)+{\omega }_x)+e\cos {\omega }_x],$$(10)

where K_x denotes the RV semi-amplitude (with x = 1 for the primary and x = 2 for the secondary), $$K_x=\frac{2\pi a_x\sin i}{P\sqrt{1-e^2}}.$$(11)

With this forward model in place, we can evaluate the predicted observables for RV, relative astrometry, HIPPARCOS, Gaia, and ground-based datasets and construct the joint log-likelihood: $$\ln {\mathcal{L}}_{\text{tot}}={\displaystyle \sum _k\ln }{\mathcal{L}}_{\text{RV},k}+\ln {\mathcal{L}}_{\text{RelAst}}+\ln {\mathcal{L}}_{\text{Hip}}+\ln {\mathcal{L}}_{\text{Gaia}}+\ln {\mathcal{L}}_{\text{GB}}.$$(12)

Each term is outlined below.

3.5 Radial velocity

Spectroscopic Doppler measurements constrain the orbital parameters through RV time series. For µ Her Aa, we have measurements, v_k,i, from multiple instruments, k, at epochs, t_i. Using Eq. (10), the model-predicted RV of Aa at time t is $$v_{\text{Aa,model}}(t)={\gamma }_k+v_{1,\text{A}}(t)+v_{1,\text{A-BC}}(t),$$(13)

where γ_k is an instrument-specific velocity zero-point, and v_1,A(t) and v_1,A–BC(t) are the primary’s reflex velocities induced by the inner Aa–Ab and wide A–BC orbits, respectively. For each instrument, k, the likelihood contribution from its RV time series is $$\ln {\mathcal{L}}_{\text{RV},k}=-\frac{1}{2}{\displaystyle \sum _i\left[\frac{{(v_{k,i}-v_{\text{Aa,model}}(t_i))}^2}{{\sigma }_{k,i}^2}+\ln (2\pi {\sigma }_{k,i}^2)\right]},$$(14)

where the sum runs over i measurements, and the total variance is $${\sigma }_{k,i}^2={\sigma }_{k,i,\text{obs}}^2+{\sigma }_{\text{jit},k}^2.$$

Here, σ_k,i,obs is the formal uncertainty of measurement i from instrument k, and σ_jit,k is a per-instrument jitter term that captures additional noise, including stellar variability and instrumental systematics.

3.6 Relative astrometry

Relative astrometry measures the secondary’s sky-projected position relative to the primary. The two observables are the angular separation, ρ, and the position angle of the secondary, θ. From the orbital displacement in Eq. (5), the model predictions for these observables are $${\rho }_{\text{model}}(t)=\sqrt{{[{\varpi }_{\text{CM}}\Delta {\alpha }_{\text{orb}}^{*}(t)]}^2+{[{\varpi }_{\text{CM}}\Delta {\delta }_{\text{orb}}(t)]}^2},$$(15) $${\theta }_{\text{model}}(t)=\mathrm{atan}2(\Delta {\alpha }_{\text{orb}}^{*}(t),\Delta {\delta }_{\text{orb}}(t)).$$(16)

where ‘atan2’ denotes the two-argument arctangent function. Unlike the ordinary arctangent, it uses the signs of both coordinate offsets to return an angle in the correct quadrant, ensuring a continuous position angle over the full 0–2π range. For measurements ρ_i and θ_i with uncertainties σ_ρ,i and σ_θ,i, the log-likelihood is $$\ln {\mathcal{L}}_{\text{RelAst}}=-\frac{1}{2}{\displaystyle \sum _i\left[\frac{{({\rho }_i-{\rho }_{\text{model}}(t_i))}^2}{{\sigma }_{\rho ,i}^2}+\frac{{({\theta }_i-{\theta }_{\text{model}}(t_i))}^2}{{\sigma }_{\theta ,i}^2}\right]}.$$(17)

Because position angles are circular quantities, the residuals, θ_i − θ_model(t_i), are wrapped into the interval (−180°, 180°] before squaring, ensuring continuity across the 0°/360° boundary.

3.7 Absolute astrometry

Absolute astrometry measures sky positions in absolute coordinates defined relative to background stars. In this work, we used three independent sources of absolute astrometry: HIPPARCOS, Gaia, and ground-based measurements. Each is treated with a dedicated likelihood function, as described below.

For an unresolved binary, such as µ Her Aa–Ab, the position inferred from absolute astrometry is not the position of either star individually, but corresponds to an effective image centroid (the photocentre) that lies somewhere between the two components. When the secondary is much fainter than the primary, this effective position remains close to the primary but is slightly displaced towards the companion. As the system orbits, this displacement reduces the observed amplitude of the primary’s apparent barycentric motion compared to its true orbit. For clarity, when describing the modelled sky path of µ Her Aa below, we first express it in terms of the primary’s barycentric position alone, neglecting this effect. In Appendix A, we then describe how we model small corrections to this position that arise from the mixing µ Her Aa and µ Her Ab’s light.

3.7.1 Hipparcos

For HIPPARCOS, we used IAD, which consist of the residuals of the one-dimensional along-scan abscissa measurements with respect to the published five-parameter (5p) solution. Each entry provides the epoch, t_i, scan angle, ψ_i (after conversion to the Gaia convention), parallax factor, Π_ψ,i, measured residual, w_H,res,i, and formal uncertainty, σ_H,i.

Our first step towards creating the HIPPARCOS likelihood function is to use the published residuals, w_H,res,i to reconstruct the abscissa measurements themselves, w_H,i: $$w_{\text{H},i}=w_{\text{H,ref}}(t_i)+w_{\text{H,res},i},$$(18)

where the reference abscissa is $$w_{\text{H,ref}}(t)=\sin \psi \left[\Delta {\alpha }_{0,\text{H}}^{*}+{\mu }_{{\alpha }_{\text{H}}^{*}}(t-t_0)\right]+\cos \psi \left[\Delta {\delta }_{\text{H}}+{\mu }_{\delta ,\text{H}}(t-t_0)\right]+{\Pi }_{\psi }{\varpi }_{\text{H}}.$$(19)

Because our model is expressed in the HIPPARCOS reference frame and at the HIPPARCOS reference epoch, the offsets $\text{Δ}{\alpha }_{0,\text{H}}^{*}$ and ∆δ_0,H vanish by definition. The model-predicted abscissa, w_H,model(t), then follows from Eq. (9). With the model abscissae and reconstructed measurements, we are ready to create the HIPPARCOS contribution to the likelihood, $$\ln {\mathcal{L}}_{\text{Hip}}=-\frac{1}{2}{\displaystyle \sum \frac{{(w_{\text{H},i}-w_{\text{H,model}}(t_i))}^2}{{\sigma }^2{}_{H,i}}},$$(20)

3.7.2 Gaia DR3

Gaia DR3 does not provide the underlying epoch astrometry for each source, only a linear, 5-p astrometric solution containing position, parallax, and proper motion. We therefore cannot fit the Gaia measurements directly in the same way as the HIPPARCOS IAD. Nevertheless, the published DR3 solution still encodes strong geometric information: the sky path predicted by our model must pass close to the reported (α, δ) at the DR3 reference epoch, and the local tangent to that path must align with the quoted proper-motion vector.

To exploit these constraints, we ask what Gaia would have reported had our dynamical model been exactly correct. To answer that, we mimicked Gaia’s measurement and fitting procedure on the model-predicted positions to obtain an astrometric solution that is directly comparable to the published one. To reconstruct the Gaia measurements, we used the DR3 observation epochs and scan angles from the Gaia Observation Forecast Tool. Using these epochs, we evaluated the model sky position and projected it onto Gaia’s along-scan direction using the scan angles (Eq. (9)). This yields a set of synthetic along-scan abscissa measurements, w_G,model(t_i), which we fitted with a linear least-squares 5p model. In other words, we fit a straight-line motion to data that in reality contain slight orbital curvature, exactly as Gaia does in its single-star solution. The fit returns the model-predicted Gaia parameters $${\widehat{p}}_G^{H\leftarrow G}=(\widehat{\Delta {\alpha }_{0,G}^{*}},\widehat{\Delta {\delta }_{0,G}},\widehat{{\mu }_{{\alpha }_G^{*}}},\widehat{{\mu }_{{\delta }_G}},\widehat{{\varpi }_G}).$$(21)

By construction, ${\widehat{p}}_G^{H\leftarrow G}$ shares the same biases as the published solution, $p_G^{H\leftarrow G}$, whenever our model correctly reproduces Gaia’s observations. We built the Gaia likelihood term from the difference between the published and fitted 5p solutions, $$\Delta p=p_G^{H\leftarrow G}-{\widehat{p}}_G^{H\leftarrow G},$$(22)

with the covariance matrix, C, derived from the published uncertainties and correlation coefficients: $$\ln {\mathcal{L}}_{\text{Gaia}}=-\frac{1}{2}\Delta p^{\top }{C}^{-1}\Delta p,$$(23)

3.7.3 Ground-based astrometry

For ground-based catalogues, the published observables are direct (α_i, δ_i) positions at epochs t_i. We transformed them into HIPPARCOS-frame offsets, where our model operates: $$\begin{array}{l}\Delta {\alpha }_{\text{GB},i}^{*}=({\alpha }_i-{\alpha }_{\text{H}})\cos {\delta }_{\text{H}},\\ \Delta {\delta }_{\text{GB},i}=({\delta }_i-{\delta }_{\text{H}}).\end{array}$$(24)

Unlike HIPPARCOS or Gaia, ground-based catalogues generally do not provide precomputed parallax factors. We therefore computed the tangent-plane factors, Π_α(t) and Π_δ(t), from Earth’s position with respect to the Solar System barycentre, r_⊕(t) = (X_⊕, Y_⊕, Z_⊕), evaluated from standard Solar System ephemerides (accessed in practice via the Astropy package). To evaluate these factors, we specified the target’s reference (i.e., HIPPARCOS) direction (α_H, δ_H). The parallax factors are therefore $${\Pi }_{{\alpha }^{*}}(t)=X_{\oplus }(t)\sin {\alpha }_{\text{H}}-Y_{\oplus }(t)\cos {\alpha }_{\text{H}},$$(25) $${\Pi }_{\delta }(t)=X_{\oplus }(t)\cos {\alpha }_{\text{H}}\sin {\delta }_{\text{H}}+Y_{\oplus }(t)\sin {\alpha }_{\text{H}}\sin {\delta }_{\text{H}}-Z_{\oplus }(t)\cos {\delta }_{\text{H}},$$(26)

With these definitions, we used Eq. (8) to derive model positions at each epoch, from which we constructed the likelihood contribution from ground-based astrometry: $$\begin{array}{l}\ln {\mathcal{L}}_{\text{GB}}=-\frac{1}{2}{\displaystyle \sum _i\frac{{(\Delta {\alpha }_{\text{GB},i}^{*}-\Delta {\alpha }_{\text{GB,model}}^{*}(t_i))}^2}{{\sigma }_{{\alpha }^{*},i}^2}}\\ -\frac{1}{2}{\displaystyle \sum _i\frac{{(\Delta {\delta }_{\text{GB},i}-\Delta {\delta }_{\text{GB,model}}(t_i))}^2}{{\sigma }_{\delta ,i}^2}}.\end{array}$$(27)

3.8 Sampling and priors

We sampled the posterior with an affine-invariant MCMC implemented in the emcee package (Foreman-Mackey et al. 2013), evaluating the full joint likelihood at every step. Derived quantities were recomputed at each iteration so that their uncertainties were propagated naturally through the chain. To impose uninformative priors on the orbital parameters, we sampled in the transformed variables cos(i), log(P), log(σ_jit), and $(\sqrt{e}\cos \omega ,\sqrt{e}\sin \omega )$, which ensures uniformity over the physically relevant domains and avoids boundary artefacts.

In addition to these largely uninformative choices, we imposed a dynamical-stability prior on the wide A–BC orbit. Without such a prior, the relatively weak constraints on the outer orbit admit highly eccentric configurations in which the A and BC subsystems pass so close at periastron that the inner binaries would almost certainly be disrupted.

A definitive assessment of long-term stability would require direct N-body integrations of the full quadruple. Given the strong hierarchy of timescales, however, carrying out such integrations for every MCMC sample is computationally prohibitive. An alternative would be to map the stability landscape across the relevant parameters and interpolate during sampling, but resolving small stable islands in the high-dimensional parameter space at a sufficient resolution would be similarly prohibitive.

Instead, we adopted the analytical, empirically calibrated stability criterion of Mardling & Aarseth (2001), formulated for hierarchical triples. We treated the tighter BC pair as a single outer perturber with point mass M_BC on the wide orbit around the A barycentre, and evaluated the criterion for each MCMC draw. The criterion requires that the outer periastron distance, r_p = a_A,BC(1 − e_A,BC), exceed the critical value $r_p^{\text{crit}}$, $$r_p^{\text{crit}}\equiv 2.8a_{\text{A}}\frac{{(1+q)}^{2/5}{(1+e_{\text{A,BC}})}^{2/5}}{{(1-e_{\text{A,BC}})}^{1/5}}\left(1-\frac{0.3}{\pi }{\Phi }_{\text{A/A,BC}}\right),$$(28)

where q = M_BC/(M_Aa + M_Ab). Because the Mardling & Aarseth (2001) criterion was calibrated for hierarchical triples and we applied it here by approximating the BC pair as a point-mass perturber, it should be regarded as a conservative guide rather than an exact stability boundary for the full quadruple. In particular, configurations near the nominal threshold could still be stable once the full four-body dynamics are taken into account. We therefore encoded the criterion as a soft stability prior for which comfortably hierarchical solutions are essentially unaffected, while draws with very small periastron separations are progressively downweighted. Specifically, we added to our log-prior the term $$\ln p\propto -10\ln \left[1+\exp \left(\frac{1-r_p/r_p^{\text{crit}}}{0.05}\right)\right],$$(29)

We set the transition width to 5% in $r_p/r_p^{\text{crit}}$ to keep sampling smooth, and chose the overall strength such that clearly unstable configurations are sufficiently suppressed.

4 Results

Our joint analysis establishes a precise and internally coherent dynamical solution for the µ Her system. A summary of the fitted parameters is provided in Table 4. From these, we derived the following component masses: $$\begin{array}{cc}{M}_{\text{Aa}}=1.134\pm 0.007M_{\odot },& M_{\text{Ab}}=0.2286\pm 0.0006M_{\odot },\\ M_{\text{B}}=0.417\pm 0.005M_{\odot },& M_{\text{C}}=0.445\pm 0.005M_{\odot }.\end{array}$$(30)

To verify the robustness of our results, we assessed the agreement between the different observational techniques by comparing our joint solution to models constrained by individual data subsets. We began by examining the role of relative astrometry.

4.1 Relative astrometry

Figures 3, 4, and 5 illustrate the constraints from relative astrometry for the A, BC, and A–BC systems, respectively. The solutions derived solely from relative astrometry are consistent with our joint model, but the latter provides significantly tighter constraints on the orbital parameters, especially for A. In the case of the wide A–BC system, the measurements trace the separation between the photocentres of the A and BC binaries. Because both photocentres shift due to their respective reflex motions, the observed trajectory deviates from a simple Keplerian ellipse.

4.2 Radial velocity

Figure 6 summarises the µ Her Aa RV data and models. Considered in isolation, the RVs leave the orbital period, P_A, periastron epoch, τ_A, and semi-amplitude, K_1,A, only weakly constrained, although the shape parameters, e_A and ω_A, are relatively well determined from the velocity curve. The constraints strengthen once the angular scale of the orbit is set by relative astrometry. In the joint fit, the RVs fix the primary’s barycentric semi-major axis, a_1,A, in addition to refining e_A and ω_A. For the RV-only solution, we allowed a linear RV drift as a function of time to absorb the contribution from the A–BC system.

In the highest-precision APF and SONG data, the RV residuals show a coherent, low-amplitude structure that we ascribe to rotation-modulated stellar activity; a detailed analysis will be presented in Kjeldsen (in prep.). We modelled this with per-instrument jitter terms at the level of a few meters per second, which absorbed the activity signal without materially biasing the orbital elements. Where the APF and SONG time series overlap, the residual patterns are similar, supporting an astrophysical rather than an instrumental origin (see Fig. 6).

For the BC subsystem, our single spectroscopic observation of µ Her BC exhibits asymmetric absorption lines. These asymmetries indicate that, at the epoch of observation, the brighter C component was approaching along the line of sight, while the B component was receding. To disentangle the contributions from B and C and recover individual radial velocities, we computed broadening functions (Rucinski 1999) based on model atmospheres from Kurucz (1979), and fitted two Gaussian profiles to the resulting distribution. We adopted uncertainties from the diagonal elements of the covariance matrix from the least-squares minimisation. The broadening function and extracted velocities are displayed in Fig. B.1.

4.3 Absolute astrometry

The absolute astrometry data trace the sky-plane trajectory of A’s photocentre. This composite path encodes the combined effects of proper motion, parallax, and the inner and outer orbits, and therefore constrains the full geometry of the model. In practice, the most direct contribution from absolute astrometry is the system parallax. Apart from this role, the absolute astrometry data is broadly analogous to that of the Aa radial velocities: it constrains the orientation of the Aa–Ab orbit on the sky and sets the scale of the primary’s barycentric semi-major axis, a_1,A. Much of this constraining power comes from the Gaia DR3 proper motion, which effectively measures the local tangent of µ Her Aa’s trajectory.

4.3.1 Parallax

Our joint analysis yields a system parallax of ϖ_CM = 120.069 ± 0.089 mas. This result reconciles the independent measurements of Aa from HIPPARCOS with the Gaia parallaxes for components Aa, B, and C. Moreover, the published Gaia DR3 5p solution (transformed into the HIPPARCOS frame, $p_G^{H\leftarrow G}$) agrees with the solution implied by our dynamical model ${\widehat{p}}_G^{H\leftarrow G}$, demonstrating their mutual consistency (see Table 3).

Fig. 2 Sky-plane trajectory and residuals for the absolute-astrometry solution. Top: photocentre trajectory from our joint fit. The black curve comprises proper motion, parallax, and orbital reflex motion. The dashed grey line shows the underlying linear CM motion. Right: residuals. From top to bottom these panels show (I) HIPPARCOS along-scan abscissae, (II) Gaia along-scan abscissae, (III) ground-based positions in α∗, and (IV) ground-based positions in δ. Bottom: orbital reflex motion of µ Her Aa relative to the barycentre A. We display random draws from the posterior distribution of the absolute-astrometry-only solution in green and the joint solution in black. The pink cone represents the proper motion from the Gaia solution (3 σ uncertainty); the black cones represent the orbital tangent vectors corresponding to the joint solutions.

4.3.2 Orbital parameters

The top left panel of Fig. 2 displays the total projected sky path, while the bottom panel subtracts the parallax, proper motion, and outer orbit to isolate the reflex motion of Aa. To test the internal consistency of the absolute astrometry, we performed two validation fits and compared them to our joint solution. First, we solved the Aa–Ab orbit while omitting absolute astrometry data entirely, yielding the result M_Aa = 1.134±0.009 M_⊙. Second, we omitted radial velocities and let absolute astrometry alone determine a_1,A. The resulting solution, also shown in Fig. 2, yields M_Aa = 1.149 ± 0.011 M_⊙.

4.3.3 BC component

For components B and C, the available absolute astrometry is limited to their Gaia DR3 catalogue solutions. After subtracting both the CM proper motion and the contribution from the wide A–BC orbit, the residual proper motions represent the instantaneous sky-plane orbital motion of B and C about their shared barycentre. This information, in turn, allows us to determine the individual masses of the BC subsystem.

4.4 Final parameter estimates

To understand what limits the precision on M_Aa, we decomposed its variance into contributions from each fitted parameter. Because many parameters are correlated in the joint posterior, this decomposition must respect their covariance structure, and we therefore worked directly with the MCMC samples. From these, we computed the covariance matrix, C_xx, of the model parameters and the vector of cross-covariances, $C_{x{M}_{\text{Au}}}$, between each parameter and the derived mass. The linear relation $$\beta =C_{xx}^{-1}{C}_{xM{L}_{Aa}}$$(31)

then gives an effective gradient: it describes how small changes in each parameter translate into changes in M_Aa near the posterior mean. The contribution of parameter i to the total variance of the mass is $$S_i={\beta }_i{\left(C_{x{M}_{Aa}}\right)}_i.$$(32)

Table 4 lists these S_i as normalised percentages, providing an error budget for the dynamical mass of Aa. This exercise reveals that P_A is the main driver of uncertainty to M_Aa, followed by a_A. This result is to be expected, since our data do not cover a full orbit of the A subsystem, in either RV or relative astrometry. Corner plots showing the posterior correlations among the orbital and astrometric parameters are provided in Figs. B.2, B.3, B.4, and B.5.

Fig. 3 Relative astrometry of the µ Her A subsystem. The plot is centred on the primary µ Her Aa. Data points are colour-coded by epoch, progressing from purple to yellow. To illustrate the constraints provided by relative astrometry alone, we display 600 random draws from the corresponding posterior distribution (red traces). The solid black lines trace the orbit derived from our full joint fit, with the median solution highlighted in white. The bottom panels show residuals versus time in separation, ρ, and position angle, θ.

Fig. 4 Relative astrometry of the BC subsystem. The plot is centred on the primary of that system, µ Her C. The format is the same as Fig. 3, but with sample orbits coloured orange.

4.5 Mutual inclinations and spin-orbit angles

By combining radial velocities with absolute and relative astrometry, we lifted the degeneracy between the ascending and descending nodes for all orbits in the system. With the full three-dimensional geometry accessible, we could determine the mutual inclinations, Φ, between the orbital planes and thereby map the spatial architecture of the quadruple. The mutual inclination between any two orbits follows from $$\cos \Phi =\cos i_1\cos i_2+\sin i_2\cos \left({\Omega }_1-{\Omega }_2\right),$$

from which we inferred Φ_A/A–BC = 44.94 ± 3.22° for the inner A subsystem relative to the outer A–BC orbit, Φ_BC/A–BC = 140.00 ± 2.24° for the BC subsystem relative to A–BC, and Φ_A/BC = 122.27 ± 0.12° between the two inner binaries; see also Table 4.

Additionally, the stellar inclination of µ Her Aa was measured from asteroseismology by Grundahl et al. (2017) as $i_{\star }={63}_{-10}^{+9}\circ$. This value is fully consistent with our orbital inclination for the Aa–Ab pair, i_A = 62.53 ± 0.07°, and therefore provides no evidence of a misalignment between the stellar spin axis and the orbital axis, although an azimuthal misalignment remains possible. Empirical constraints on spin–orbit alignment for binaries in the 3–100 au range (Justesen & Albrecht 2020) remain sparse, but by analogy with trends observed at smaller separations (Marcussen et al. 2024), spin-orbit alignment in µ Her A would not be unexpected.

Fig. 5 Relative astrometry of the outer A–BC system. The plot is centred on µ Her Aa and the format is the same as in Fig. 3. Randomly drawn orbits are coloured cyan.

5 Discussion

5.1 Formation and dynamical history

The present-day orbital architecture of µ Her is markedly non-coplanar. The inner Aa–Ab orbit is inclined by ∼ 45° with respect to the wide A–BC orbit, while the BC subsystem is retrograde at ∼140°. With an outer semi-major axis of ∼1100 au, the system falls between the proposed formation pathways of core fragmentation and disc fragmentation, as is discussed by Offner et al. (2023). A more eccentric and misaligned architecture is expected from core fragmentation, whereas a more circular and aligned system is characteristic of disc fragmentation. The observed geometry naively suggests the former pathway. However, the close periastron passages in the wide orbit are expected to drive secular effects, such as chaotic eccentric Kozai–Lidov cycles (Naoz 2016), which can substantially reshape an initially simple configuration and partially erase its formation signature.

An additional clue about the system’s history is provided by the inferred stellar inclination of Aa, which is consistent with spin–orbit alignment in the A subsystem. This suggests that the inner Aa–Ab pair has formed in a shared circumbinary disc, and that the system has likely not experienced strong inclination– eccentricity cycling since formation. Its present orientation may thus still be close to primordial. By contrast, the BC pair is both less massive and more compact, and therefore carries a smaller share of the system’s total angular momentum. In the presence of the highly eccentric wide orbit, long-term secular torques can more readily reorient such a low-angular-momentum inner binary, potentially driving it to the high inclination we observe today.

5.2 µ Herculis as a benchmark system

The precise dynamical mass derived here establishes µ Her Aa as a premier benchmark for stellar physics. By setting the fundamental scale of the system independent of stellar evolution models, we provide a rigorous testbed for validating asteroseismic scaling relations and examining the internal physics of subgiants.

Prior to this work, estimates of the primary mass relied on spectroscopic calibrations or model-dependent asteroseismic scaling, leading to a significant spread in values. Figure 7 compares our model-independent dynamical mass against a representative selection of these literature values.

Our dynamical mass for µ Her Aa, 1.134 M_⊙, falls slightly below the early non-seismic spectroscopic estimate by Fuhrmann (1998) (1.14 M_⊙), while remaining at the upper end of the masses inferred from asteroseismic and evolutionary modelling. In the first modelling efforts based on the initial oscillation detections of Bonanno et al. (2008), Yang & Meng (2010) reported a best-fit mass of 1.00 M_⊙, noting that solutions of up to 1.10 M_⊙ could be obtained under alternative assumptions about the input physics.

With the advent of higher-precision seismic constraints, the inferred masses moved upwards but remained below the early spectroscopic value. The grid modelling by Li et al. (2019) yielded ${1.10}_{-0.02}^{+0.06}{M}_{\odot }$, and the solutions presented by Grundahl et al. (2017) (ASTFIT, BASTA) cluster close to 1.12 M_⊙. Most recently, analyses incorporating TESS data and updated glitch modelling by Lund et al. (2025) and Gupta et al. (2025) have remained in this same regime, with Gupta et al. (2025) deriving ${1.105}_{-0.024}^{+0.058}{M}_{\odot }$.

As is illustrated in Fig. 7, our dynamical mass therefore provides an independent, geometry-driven constraint that is broadly compatible with the seismic consensus, while selecting its high-mass end. Any successful asteroseismic or evolutionary model of µ Her Aa must reproduce both the detailed oscillation spectrum and this dynamical mass.

Beyond its role in calibrating asteroseismology, the µ Her system also offers an opportunity to test models of low-mass stars. Existing spectroscopic and asteroseismic analyses already indicate that µ Her Aa is an old, metal-rich subgiant (e.g. Fuhrmann 1998; Li et al. 2019; Gupta et al. 2025), and the dedicated modelling of the SONG asteroseismic data will deliver precise constraints on its age and heavy-element abundance. Assuming coeval formation, these values can be adopted for Ab, B, and C; combined with our dynamical masses, the three M dwarfs then constitute a small sequence of benchmark stars at fixed age and supersolar metallicity. A dedicated comparison of their radii, luminosities, and spectra with low-mass evolutionary and atmospheric grids would provide a stringent test of M-dwarf models, but we leave such an analysis to future work.

Fig. 6 Radial-velocity data and model fits for µ Her Aa. Top: radial velocities versus time. The points are coloured by instrument (see legend). Note that SONG RVs are nightly medians. Green curves represent random draws from the posterior of the RV-only solution, while black curves depict the joint solution; the median joint model is overlaid in white. The dashed lines trace the long-term trends: the green line indicates the best-fitting linear trend for the RV-only solution, while the black line follows the reflex motion induced by the A–BC orbit in the joint fit. Middle: residuals relative to the median joint model. The error bars include jitter. Bottom: residuals in three-year segments. To visualise the dense time series, residuals are displayed in three-year segments, starting from the year indicated on the left side of the plot and vertically offset by 40 ms−1.

Fig. 7 Comparison of the primary mass, MAa, derived in this work (vertical grey band, showing the 1σ credibility interval) against literature values. The markers represent estimates from non-seismic spectroscopy (Fuhrmann 1998), stellar modelling (Yang & Meng 2010; Li et al. 2019), and asteroseismic analyses (Grundahl et al. 2017; Gupta et al. 2025).

5.3 Future prospects

While the current solution is robust, several avenues remain to further tighten the constraints on this system. The release of Gaia DR4 will provide time-series epoch astrometry, enabling a direct fit to the Gaia measurements that should further improve the precision of orbital elements and the system parallax. At present, the error budget of M_Aa is dominated by the relatively small number of astrometric measurements of the Aa–Ab separation and their modest temporal coverage. Additional high-precision measurements taken over a few years would therefore substantially improve the accuracy of the masses in the A subsystem.

Acknowledgements

This publication includes observations made with the SONG network of telescopes operated by Aarhus University, Instituto de Astrofísica de Canarias, the National Astronomical Observatories of China, the University of Southern Queensland and New Mexico State University. MNL acknowledges support from the ESA PRODEX programme (PEA 4000142995). This research has made use of the SIMBAD database (Wenger et al. 2000), operated at CDS, Strasbourg, France. This research has made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. This research made use of the SONG database SODA (https://soda.phys.au.dk/), operated and maintained at Aarhus University, DK. The Robo-AO-2 system is supported by the National Science Foundation under Grant Nos. AST-1712014 and AST-2509941, the State of Hawaii Capital Improvement Projects, the Mt. Cuba Astronomical Foundation, and by a gift from the Lumb Family. Support for the infrared camera for Robo-AO and Robo-AO-2 was provided by the Mt. Cuba Astronomical Foundation and through the National Science Foundation under Grant No. AST-1106391 Based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku, and the University of Oslo, representing Denmark, Finland and Norway, the University of Iceland and Stock-holm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. We acknowledge the use of the following Python-based software modules: Astropy (Astropy Collaboration 2013), scikit-image (van der Walt et al. 2014).

References

Appendix A Additional effects

In the foregoing description, we intentionally omitted several small effects. In nearby, high–proper-motion systems such as µ Herculis, however, these effects are no longer negligible at the precision we seek. Here, we introduce the additional terms needed to relax those assumptions, explain how each is incorporated into the model, and quantify its contribution. Full derivations are given in the cited references.

Appendix A.1 Photocentre versus primary

When a binary system is not fully resolved, the measured position on the sky corresponds not to either star individually but to a weighted blend of their light, the photocentre. This distinction matters because the secondary’s flux pulls the centroid of the system away from the true position of the primary. In effect, the apparent orbital motion of the primary is diluted by the contribution from the companion.

We capture this dilution through the factor β, which links the photocentre displacement $(\text{Δ}{\alpha }_{\text{orb},0}^{*},\text{Δ}{\delta }_{\text{orb},0})$ to the relative orbital co-ordinates $(\text{Δ}{\alpha }_{\text{orb}}^{*},\text{Δ}{\delta }_{\text{orb}})$ as $$\left[\begin{array}{c}\Delta {\alpha }_{\text{orb},0}^{*}(t)\\ \Delta {\delta }_{\text{orb},0}(t)\end{array}\right]=\left(\beta -\frac{a_1}{a}\right)\left[\begin{array}{c}\Delta {\alpha }_{\text{orb}}^{*}(t)\\ \Delta {\delta }_{\text{orb}}(t)\end{array}\right].$$(A.1)

In the limit of small separations, where the secondary cannot be distinguished in the point-spread function, β takes its maximum value: $$\beta =\frac{f_2/f_1}{1+f_2/f_1},$$(A.2)

with f₁ and f₂ denoting the fluxes of the primary and secondary. At very large separations, β tends to zero. For Gaia, the transition between these regimes occurs at angular separations of roughly 9 mas (fully blended) and 270 mas (fully separated; Lindegren 2022). At intermediate separations, β can be modelled using Eq. 12 of Lindegren (2022).

During the Gaia DR3 observations, µ Her Aa–Ab had a mean separation of ∼1700 mas, well beyond Gaia’s nominal blending limit. Nonetheless, the one-dimensional nature of Gaia’s along-scan measurements means that the effective separation scales by how much the line between the primary and secondary align with scan angle: ∆η = ρ cos(ψ − θ). Thus, even in such a wide system, a small subset of the Gaia measurements still registers a small photocentre offset and are modelled accordingly. We estimate the flux ratio of µ Her Ab to µ Her Aa in the Gaia band as f₂/ f₁ ≈ 1/1000, based on the measured magnitude differences from Roberts et al. (2016). Using this flux ratio, we compute ∆η for each epoch from the Gaia GOST epochs and scan angles. Only four of the 64 modelled Gaia observations are affected by light blending, and even in those cases the induced shift is only of order 0.03 mas.

Appendix A.2 Tangent-plane projection

When comparing absolute positions from different catalogues, we are really comparing two directions on the celestial sphere. These positions are reported in equatorial co-ordinates (α, δ), but our model operates in a local Cartesian system—the tangent plane—centred on the HIPPARCOS reference position. Thus, taking a simple difference such as α_G − α_H implicitly assumes that curved-sky co-ordinates can be treated as flat. For small separations, this approximation is excellent. However, because µ Her has a large proper motion and the HIPPARCOS–Gaia baseline spans 24.75 years, the curvature of the sky becomes non-negligible. We therefore transform all non-HIPPARCOS catalogue positions into the tangent plane using the standard gnomonic projection (Murray 1983). This effectively removes the curvature of the celestial sphere, allowing the corrected positions to be treated linearly in our models. For µ Her, this projection introduces a correction to positions of approximately 0.38 mas over the HIPPARCOS–Gaia baseline.

Appendix A.3 Perspective effects due to radial velocity

µ Her is approaching us with a CM radial velocity of around −16.3 km s⁻¹. As the system moves closer, its parallax increases. Additionally, its apparent proper motion grows too, since unchanged space motion projects differently as perspective changes. These effects are known as perspective acceleration (Lindegren & Dravins 2021). Across the 24.75-year HIPPAR-COS–Gaia baseline, the radial velocity changes the system parallax by ∼6 µas and the total proper motion by ∼0.082 mas yr⁻¹. Accumulated over the baseline, the difference in proper motion displaces the system by ∼1.5 mas. In comparison with Gaia DR3 uncertainties, the effect is negligible for parallax, but comparable to the errors in proper motion and clearly significant for positions (see Table 3). Analogous to the deprojection procedure, we correct the catalogue positions and proper motions for perspective effects once, using a reference radial velocity of −16.3 km s⁻¹. This adjustment removes the geometric distortion caused by the system’s changing line-of-sight distance, without introducing additional parameters into the model itself.

Appendix A.4 Gaia–Hipparcos frame rotation

Although both HIPPARCOS and Gaia are tied to the ICRS, their reference frames are not perfectly aligned. Small residual rotations and spins exist between them. Following the prescription of Fabricius et al. (2021), we therefore rotate Gaia astrometry into the HIPPARCOS frame before combining the datasets. In practice, this amounts to a global rigid transformation of the Gaia 5p solution, ensuring that both catalogues describe the same inertial system. For DR3, the frame-rotation correction corresponds to about ∼0.20 mas yr⁻¹ in total proper motion and ∼5.1 mas in position. The Munn et al. (2022) positions are anchored to the Gaia DR3 reference frame, and we therefore apply the same frame-rotation correction to these measurements before incorporating them into our model.

Appendix B Broadening function and corner plots

Broadening function of the BC component is shown here along with corner plots from our MCMC sampler.

Fig. B.1 Broadening function of µ Her BC. Top: Broadening function of µ Her BC from our single December 2025 epoch. The best-fitting two-Gaussian model is overplotted (blue and yellow components; combined profile in green). Middle: Residuals between the broadening function and the model. Bottom: Radial-velocity curve of the BC subsystem. The two velocities extracted from the broadening function are indicated.

Fig. B.2 Corner plot showing the two-dimensional posterior distributions for a subset of parameters from our full joint fit. The panels display the astrometric parameters together with the derived masses of the four components.

Fig. B.3 Corner plot of subsystem A’s orbital parameters. The figure format is as the same as Fig. B.2.

Fig. B.4 Corner plot of subsystem BC’s orbital parameters. The figure format is as the same as Fig. B.2.

Fig. B.5 Corner plot of subsystem A–BC’s orbital parameters. The figure format is as the same as Fig. B.2.

All Tables

All Figures

Fig. 1 To-scale overview of the orbital configuration and naming scheme of the µ Her system. Components are colour-coded as Aa (green), Ab (red), B (yellow), and C (blue). The black curve traces the total sky motion of the system’s CM over a 40-year interval, combining a linear term (dashed grey line) with the annual parallax. Star symbols mark the positions on January 1 2026. Solid lines show the projected forward motion of each component over the next 40 years, and of the wide A–BC orbit over the next 400 years. To visualise the orbital inclination, segments closer to the observer are rendered in white and cast a shadow onto the sky plane; the transition between the white and grey segments marks the line of nodes. The inner A and BC subsystems are presented in greater detail in the right-hand panels.

In the text

Fig. 2 Sky-plane trajectory and residuals for the absolute-astrometry solution. Top: photocentre trajectory from our joint fit. The black curve comprises proper motion, parallax, and orbital reflex motion. The dashed grey line shows the underlying linear CM motion. Right: residuals. From top to bottom these panels show (I) HIPPARCOS along-scan abscissae, (II) Gaia along-scan abscissae, (III) ground-based positions in α∗, and (IV) ground-based positions in δ. Bottom: orbital reflex motion of µ Her Aa relative to the barycentre A. We display random draws from the posterior distribution of the absolute-astrometry-only solution in green and the joint solution in black. The pink cone represents the proper motion from the Gaia solution (3 σ uncertainty); the black cones represent the orbital tangent vectors corresponding to the joint solutions.

In the text

Fig. 3 Relative astrometry of the µ Her A subsystem. The plot is centred on the primary µ Her Aa. Data points are colour-coded by epoch, progressing from purple to yellow. To illustrate the constraints provided by relative astrometry alone, we display 600 random draws from the corresponding posterior distribution (red traces). The solid black lines trace the orbit derived from our full joint fit, with the median solution highlighted in white. The bottom panels show residuals versus time in separation, ρ, and position angle, θ.

In the text

	Fig. 4 Relative astrometry of the BC subsystem. The plot is centred on the primary of that system, µ Her C. The format is the same as Fig. 3, but with sample orbits coloured orange.
In the text

	Fig. 5 Relative astrometry of the outer A–BC system. The plot is centred on µ Her Aa and the format is the same as in Fig. 3. Randomly drawn orbits are coloured cyan.
In the text

Fig. 6 Radial-velocity data and model fits for µ Her Aa. Top: radial velocities versus time. The points are coloured by instrument (see legend). Note that SONG RVs are nightly medians. Green curves represent random draws from the posterior of the RV-only solution, while black curves depict the joint solution; the median joint model is overlaid in white. The dashed lines trace the long-term trends: the green line indicates the best-fitting linear trend for the RV-only solution, while the black line follows the reflex motion induced by the A–BC orbit in the joint fit. Middle: residuals relative to the median joint model. The error bars include jitter. Bottom: residuals in three-year segments. To visualise the dense time series, residuals are displayed in three-year segments, starting from the year indicated on the left side of the plot and vertically offset by 40 ms−1.

In the text

	Fig. 7 Comparison of the primary mass, MAa, derived in this work (vertical grey band, showing the 1σ credibility interval) against literature values. The markers represent estimates from non-seismic spectroscopy (Fuhrmann 1998), stellar modelling (Yang & Meng 2010; Li et al. 2019), and asteroseismic analyses (Grundahl et al. 2017; Gupta et al. 2025).
In the text

Fig. B.1 Broadening function of µ Her BC. Top: Broadening function of µ Her BC from our single December 2025 epoch. The best-fitting two-Gaussian model is overplotted (blue and yellow components; combined profile in green). Middle: Residuals between the broadening function and the model. Bottom: Radial-velocity curve of the BC subsystem. The two velocities extracted from the broadening function are indicated.

In the text

	Fig. B.2 Corner plot showing the two-dimensional posterior distributions for a subset of parameters from our full joint fit. The panels display the astrometric parameters together with the derived masses of the four components.
In the text

	Fig. B.3 Corner plot of subsystem A’s orbital parameters. The figure format is as the same as Fig. B.2.
In the text

	Fig. B.4 Corner plot of subsystem BC’s orbital parameters. The figure format is as the same as Fig. B.2.
In the text

	Fig. B.5 Corner plot of subsystem A–BC’s orbital parameters. The figure format is as the same as Fig. B.2.
In the text