© The Authors 2026

1. Introduction

The missing mass problem in galaxies has been studied since the 1930s (Babcock 1939), but gained prominence in the second half of the 20th century (van de Hulst et al. 1957), especially since the 1970s thanks to observations of flat rotation curves (RCs) in disk galaxies (e.g. Rubin & Ford 1970; Bosma 1978) that deviated from the expected Keplerian decline. This mass discrepancy (M_tot/M_bar ∼ v_obs²/v_bar² ≫ 1, where v_obs and v_bar are the circular velocities, respectively observed and generated by the baryonic content of the galaxy) led to the hypothesis of dark matter (DM) halos surrounding galaxies (e.g. White & Rees 1978; see also reviews such as Bullock & Boylan-Kolchin 2017; Bosma 2023).

This problem was later also quantified in terms of the local radial acceleration in galaxies (Sanders 1990; McGaugh 2004), which led to the identification of the apparently fundamental radial acceleration relation (RAR; McGaugh et al. 2016; Lelli et al. 2017; Stiskalek & Desmond 2023) that connects the observed acceleration (a_tot = v_obs²/r = |∂Φ_tot(r)/∂r|) to the gravitational acceleration generated by the baryons (a_bar = v_bar²/r = |∂Φ_bar(r)/∂r|) across all radii. The RAR has been empirically calibrated using the SPARC sample (Lelli et al. 2016), and is usually parametrised by the following relation (McGaugh 2008; Famaey & McGaugh 2012):

$$\begin{array}{c}\hfill a_{\mathrm{tot}}=\frac{a_{\mathrm{bar}}}{1-exp(-\sqrt{a_{\mathrm{bar}}/a_0})}·\end{array}$$(1)

This relation exhibits a scatter of ∼0.11 dex (McGaugh et al. 2016), with a characteristic acceleration scale of a₀ = 1.2 ± 0.26 × 10⁻¹⁰ m/s², below which a_bar begins to deviate systematically from a_tot. Desmond (2023) and Desmond et al. (2024) refined this analysis on the SPARC sample, and found an even lower intrinsic scatter:  σ_intrinsic = 0.034 ± 0.001_stat ± 0.001_sys dex.

A few numerical studies have shown that a relation resembling the RAR can occur in the Lambda cold dark matter (ΛCDM) framework using hydrodynamical simulations (Keller & Wadsley 2017; Garaldi et al. 2018; Tenneti et al. 2018; Dutton et al. 2019; Mayer et al. 2023) and semi-empirical analytic models (van den Bosch & Dalcanton 2000; Di Cintio & Lelli 2016; Navarro et al. 1997; Grudić et al. 2020; Paranjape & Sheth 2012; Li et al. 2022). Nevertheless, the detailed results of these studies differ slightly from each other, so that it remains unclear whether some of them do actually reproduce in detail the properties of the observed RAR, including its functional form, very low scatter, and the full diversity of RC shapes to which it applies (Ghari et al. 2019; Desmond 2017).

Alternatively, the RAR may reflect a modification of gravity (or, more generally, of dynamics) rather than unseen matter, which would naturally lead to a very small scatter of the relation. Modified Newtonian dynamics (MOND; Milgrom 1983a,b,c; see also Famaey & Durakovic 2026 for a recent review) introduced a characteristic acceleration scale a₀ ∼ 10⁻¹⁰ m/s⁻², below which Newtonian dynamics break down. Remarkably, MOND predicted the existence and shape of the RAR decades before it was empirically established.

In the standard MOND framework, a₀ is typically considered a universal constant. However, since MOND is a non-relativistic theory, extending it to a cosmological context requires additional assumptions (e.g. Bekenstein 2004; Famaey & McGaugh 2012; Milgrom 2020; Blanchet & Skordis 2024). As a result, there is no clear consensus on whether a₀ should remain constant or evolve with cosmic time in the MOND framework. For example, some studies consider a₀ as a constant parameter, while others propose that it might scale with the Hubble parameter, H(z), or follow alternative evolutionary paths (e.g. Hossenfelder & Mistele 2018; Maeder 2023; Del Popolo & Chan 2024; Gueorguiev 2024).

In contrast, within ΛCDM, it remains unclear whether the RAR can be fully reproduced as it would have to arise from the complex interplay between baryonic and DM distributions, and any z-evolution of the relation would then reflect the evolving properties of these components. In any case, the degree of z-evolution of a₀ in hydrodynamical simulations remains uncertain (e.g. Garaldi et al. 2018; Tenneti et al. 2018; Mayer et al. 2023), likely due to differences in the feedback models used.

While these theoretical frameworks offer different perspectives on the redshift evolution of the RAR, observational studies beyond z = 0 remain scarce, except for Vărăşteanu et al. (2025), who investigated the RAR up to z ∼ 0.08 and found a larger a₀ value than at z = 0. For this letter, we studied the RAR at intermediate redshifts up to z ∼ 1.44, offering new insights into the evolution of a₀ over cosmic time. This analysis is based on the results presented in Ciocan et al. (2026, hereafter Paper I), where we conducted a disk–halo decomposition using GALPAK^3D (Bouché et al. 2015) on a large sample of star-forming galaxies (SFGs), using deep data from the MUSE Hubble Ultra Deep Field survey (MHUDF, Bacon et al. 2023) offering high signal-to-noise (S/N) data ranging from S/N ∼ 10 to > 100.

2. Sample selection and methodology

To investigate the RAR at intermediate z, we selected from the parent sample in Paper I, 79 SFGs with M_★ > 10^8.8 M_⊙. This sample has the same minimum stellar mass (i.e. is mass-complete) over 0.33 < z < 1.44 and is part of the MHUDF (Bacon et al. 2023); further details are provided in Appendix A.

In Paper I we performed a detailed 3D disk–halo decomposition on the parent sample of SFGs presented in Fig. A.1, and tested six different DM halo profiles. For the present study, we used the generalised profile of Di Cintio et al. (2014, hereafter DC14) since they provided the best fits compared to other halo models, including the Navarro-Frenk-White (NFW; Navarro et al. 1997) profile. We note that Di Cintio & Lelli (2016) showed that the RAR is more readily reproduced in models with mass-dependent DM density profiles (DC14) than with universal NFW profiles, which struggle to match the observed shape of the relation, particularly in low-M_★ galaxies.

We computed a_bar and a_tot from the best-fitting DC14 disk–halo models, propagating uncertainties from the Markov chain Monte Carlo (MCMC)-derived velocity components. To minimise resolution-driven biases, measurements within r < 2 kpc (corresponding to approximately one spatial resolution element, i.e. ∼0.2″) were excluded (see Appendix B).

While our analysis relies on parametric modelling, we performed independent consistency checks of the disk–halo decomposition in Paper I, and found consistent results. Additionally, we assessed potential systematic uncertainties in the stellar mass estimates by deriving stellar mass-to-light ratios for our sample. As discussed in Appendix C, these are lower than those assumed at z = 0 in SPARC, but fully consistent with the expected decrease with z from independent studies (Drory et al. 2004).

3. Results

3.1. RAR at 0.33 < z < 1.44

Figure 1 presents the RAR at 0.33 < z < 1.44, using data points from 79 individual model RCs. As illustrated, the RAR at intermediate z is reminiscent of the RAR observed at z  =  0 (McGaugh et al. 2016). Our sample, unlike SPARC, predominantly probes the low-acceleration regime as it is largely composed of low-mass, DM-dominated SFGs lacking prominent bulges.

Fig. 1. RAR for the MHUDF sample. The purple curve shows the best-fit RAR at 0.33 < z < 1.44 using Eq. (1), while the red and green curves show the McGaugh et al. (2016) relation at z = 0 and the Vărăşteanu et al. (2025) relation at z < 0.08, respectively. The black dotted line indicates the 1:1 relation. The histogram in the inset displays the residuals with respect to our fit. The cross in the upper left corner indicates the mean 1σ uncertainties.

We fit the data with Eq. (1), leaving a₀ as a free parameter. For the fits, we used the Python package Roxy (Bartlett & Desmond 2023), which implements the marginalised normal regression (MNR) method. This approach simultaneously accounts for uncertainties in both x and y, as well as for the intrinsic scatter in the relation. We applied the MNR to the accelerations in base 10 logarithmic space, and we adopted a uniform broad prior on log a₀, ranging from −15 to 5, and on the intrinsic scatter between 0 and 3 dex. We show the best-fit curve in purple in Fig. 1, and as the purple data point in Fig. 2, which corresponds to a characteristic acceleration scale of

$$\begin{array}{c}\hfill a_0{|}_{z\sim 1}=2.{38}_{-0.10}^{+0.12}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2,\end{array}$$(2)

Fig. 2. Best-fit a0 values obtained by fitting Eq. (1) on the full sample, using the resolved RAR tracks derived under different modelling assumptions. The purple data point shows a0|z ∼ 1 obtained using the DC14 halo profile applied uniformly to the full sample (Sect. 3.1). The cyan point shows a0|z ∼ 1 derived from a galaxy-by-galaxy best-fitting ΛCDM model, where each galaxy is assigned the DM profile that maximises its Bayesian evidence (Appendix D). The blue point shows a0|z ∼ 1 recovered from the MOND-based framework (Appendix E). The red star indicates the fiducial a0|z ∼ 0 from McGaugh et al. (2016), while the green star shows the result from Vărăşteanu et al. (2025).

with the errors denoting the 95% confidence intervals (CI) from our MCMC fits. This value is significantly higher, by ∼19σ, than the canonical value of a₀|_z = 0 = 1.2 ± 0.26 × 10⁻¹⁰ m/s² inferred for the SPARC sample (McGaugh et al. 2016, shown as the red curve in Fig. 1, and as the red star in Fig. 2), and ∼5σ higher than the value of a₀|_z < 0.08 = 1.69 ± 0.13 × 10⁻¹⁰ m/s² inferred by Vărăşteanu et al. (2025) for the MIGHTEE-HI sample (shown as the green curve in Fig. 1, and as the green star in Fig. 2). At face value, this would imply that disks are more submaximal at higher z than at z = 0 in agreement with the large DM fractions found in Paper I.

The histogram in the inset of Fig. 1 shows the residuals with respect to the best-fit RAR at 0.33 < z < 1.44, which are tightly peaked, with a mean offset of μ = 0.016 dex and a standard deviation of σ = 0.19 dex, similar to the intrinsic scatter (σ_intr. ∼ 0.17 dex). The scatter that we recover for the intermediate-z sample is ∼ × 1.5 higher than that inferred by McGaugh et al. (2016) for the SPARC sample, which is of the order of ∼0.11 dex. The larger scatter that we measure is most likely related to the broad z range that we probed, as discussed in the next section, and to the poorer data quality at higher z.

As mentioned above, our analysis relies on the parametric models from Paper I, where we found that the DC14 density profile provides the best representation of the data. To assess the robustness of our results, we explored alternative DM profiles (see Appendix D) and also refitted the data within a self-consistent MOND framework (see Appendix E). The resulting a₀|_z ∼ 1 values from these various assumptions are always larger than the z = 0 value and agree within the errors with the results presented above, as shown in Fig. 2 by the cyan and blue data points.

3.2. Evolution of the RAR with z

To carefully trace the potential evolution of the RAR with cosmic time, we applied a quantile-based binning of the galaxy redshift distribution. Specifically, we divided the sample into four quantile intervals containing approximately the same number of data points and, for each interval, we fit the RAR with Eq. (1) to extract the best-fit characteristic acceleration scale. The fits were performed following the same procedure as in Sect. 3.1. We adopted this equal-population binning to stabilise the Bayesian fits by ensuring comparable statistical weight in each bin; we note that the bins correspond to adjacent, non-overlapping z intervals.

We show the resulting a₀(z) values in Fig. 3 as the black data points. As illustrated, the characteristic acceleration systematically increases with redshift, rising from ∼1.99 × 10⁻¹⁰ m/s² in the lowest-z bin to 2.71 × 10⁻¹⁰ m/s² in the highest.

Fig. 3. Redshift evolution of the characteristic acceleration scale. In each quantile-based z bin (blue bars), we fit the RAR with Eq. (1) to derive a0, with uncertainties shown as grey error bars. The star marks the SPARC z = 0 value. The purple line shows the predicted evolution, using the best-fitting parameters from the global z-dependent fit to the RAR for the full sample (i.e. using Eq. (1) where a0 was substituted with Eq. (3); see Eq. (4)). The red line shows the same, but including SPARC in the fit. The coloured shaded region shows the 1σ errors for the fits.

We evaluated the intrinsic scatter in each redshift bin and found an increase with z, from σ_intr. = 0.13 dex in the lowest-redshift bin–below the full-sample value–to σ_intr. = 0.19 dex in the highest. Nevertheless, the scatter in the individual z bins remains larger than that observed in the SPARC sample, likely reflecting the degradation in data quality and spatial resolution with increasing z, and the lower statistics.

To complement this binned analysis, we also performed a global parametric fit to quantify the z-dependence of the RAR by introducing a z-dependent acceleration scale of the form

$$\begin{array}{c}\hfill a_0(z)=a_0(0)+a_1·z,\end{array}$$(3)

where the coefficient a₁ captures the evolutionary trend with z, following the approach of Vărăşteanu et al. (2025). We then refit the RAR for the entire sample using Eq. (1), where we substitute a₀ with the z-dependent formulation (Eq. (3)), treating both a₀(0) and a₁ as free parameters in the fits. We obtain

$$\begin{array}{c}\hfill a_0(0)=1.0_{-0.04}^{+0.04}\times {10}^{-10}\mathrm{m}{s}^{-2},a_1=1.{59}_{-0.10}^{+0.11}\times {10}^{-10}\mathrm{m}{s}^{-2},\end{array}$$(4)

with the errors denoting the 95% CI. We note that including the SPARC sample (Lelli et al. 2017) in the fits yields similar results.

We show this linear relation, using the best-fit values of a₀(0) and a₁ from the fit to the MHUDF RAR (Fig. 1) as the purple line in Fig. 3. The red line shows the same, but with the SPARC sample included in the fit. In Appendix F we present the resulting posterior distributions of the fitted parameters, and discuss the small offsets between the binned estimates and the global fit.

We emphasise that this linear parametrisation of a₀(z) (Eq. (3)) is intended as a simple phenomenological description to quantify the presence and sign of a possible z-dependence rather than to provide a physically motivated description. If the real z-dependence of a₀ is non-linear or more complex, the resulting evolutionary trend could be affected by systematic biases.

In conclusion, the results presented in this section suggest that while the RAR persists out to z ∼ 1.44, the value of a₀ increases with z, potentially reflecting changes in the baryon–DM coupling, in the feedback efficiency, or in modified gravity over cosmic time.

4. Discussion

The results presented above indicate that the RAR evolves with z. We quantified this evolution with a simple, linear model (Eq. (3)), obtaining $a_1=(1.{59}_{-0.10}^{+0.11})\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$. Our result can be compared to the z ≤ 0.08 study of Vărăşteanu et al. (2025), who found a₁ = (4.47 ± 1.88) × 10⁻¹⁰ m s⁻² from a smaller sample. While their analysis provided only mild evidence (∼ 2.4σ) for a z-evolution, our broader z-baseline and larger sample yield a more statistically significant detection (at a ∼30σ level). The two measurements are statistically consistent within ∼ 1.5σ.

In the context of ΛCDM, our results show a roughly comparable z-trend as reported by Mayer et al. (2023) and Keller & Wadsley (2017), although the robustness of some of these results to reproduce the observed low-z RAR has been questioned (Milgrom 2016). For example, Mayer et al. (2023) find that a₀ grows by a factor of ∼3 from z = 0 to 2, slightly less than the factor of ∼4 increase inferred here over the same z range. In the context of modified gravity, our measured a₀(z) is faster than that of H(z) (Milgrom 1983a), and is in tension with the evolution predicted by the scale invariant vacuum paradigm (Gueorguiev 2024).

An evolution of the RAR with z should be reflected in the baryonic Tully–Fisher relation (bTFR). Current bTFR measurements based on H I detections at z ≲ 0.4 (Ponomareva et al. 2021; Jarvis et al. 2025) show little to no evidence for evolution, but are limited by small samples, narrow z baselines, and large uncertainties. Our z-dependent RAR fit (Eq. (4)) predicts a bTFR zero-point shift of ΔZP ≈ −0.2 dex between z ∼ 0 and z ∼ 1. At higher z, Übler et al. (2017) find a decrease in the bTFR zero-point between z ∼ 0 and z ∼ 0.9 (ΔZP = −0.44), though their M_bar estimates rely on scaling relations and neglect M_H I, while Jeanneau et al. (2026) found no evolution of the bTFR at z ∼ 1, accounting for neutral gas from scaling relations. Looking ahead, deep H I surveys with the Square Kilometre Array will provide critical tests of the potential evolution of both the bTFR and the RAR with z.

5. Conclusions

In this study we analysed the RAR at intermediate redshifts (0.33 < z < 1.44), for a mass-complete sample of 79 SFGs with 8.8 < log(M_★/M_⊙) < 11, using 3D forward modelling. To do so, we used the results obtained in Paper I from the disk–halo decomposition. Our analysis has yielded several key findings:

• The RAR for our z ∼ 1 sample exhibits a higher acceleration scale $a_0{|}_{z\sim 1}=2.{38}_{-0.10}^{+0.12}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$ compared to the canonical value inferred by McGaugh et al. (2016) for the z = 0 SPARC sample (Fig. 1).

• A larger value of a₀|_z ∼ 1 is found regardless of the model used to derive the total and the baryonic accelerations (Fig. 2), including various ΛCDM and MOND models for self-consistency (Appendices D and E).

• By fitting Eq. (1) independently in four redshift bins, we recover a systematic increase of the best-fit a₀ with z (Fig. 3).

• Fitting Eq. (1) to the full sample by substituting a₀ with the redshift-dependent formulation, a₀(z) = a₀(0)+a₁ ⋅ z, yields $a_1=(1.{59}_{-0.10}^{+0.11})\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$, providing evidence for an increase in a₀ with redshift (Figs. 3 and F.1).

Some hydrodynamical simulations and modified-gravity frameworks predict a possible redshift evolution of a₀, although the magnitude and functional form vary across models. At present, theoretical expectations remain too diverse to make definitive comparisons. Our results, therefore, motivate future efforts to refine predictive models and explore the physical mechanisms that could shape any evolution of the RAR with cosmic time.

Acknowledgments

We thank the anonymous referee for constructive comments that improved the manuscript. JFr acknowledges stimulating discussions on the topic with S. Trujillo-Gomez. This work is based on observations collected under ESO programmes 094.A-0289(B), 095.A-0010(A), 096.A-0045(A), 096.A-0045(B) and 1101.A-0127. This research made use of the following open source software: Astropy (Robitaille et al. 2013), numpy (van der Walt et al. 2011) and matplotlib (Hunter 2007). BC, NB and JFe acknowledge support from the ANR DARK grant (ANR-22-CE31-0006). HD is supported by the Royal Society University Research Fellowship grant 211046.

References

Appendix A: Sample of intermediate-z SFGs

For this study, we used the sample of SFGs with RC decomposition from Paper I, which are drawn from the MHUDF survey (Bacon et al. 2023). A comprehensive description of the initial sample selection is provided in Paper I; here we briefly summarise the key criteria. The parent sample consists of intermediate-z SFGs with high S/N (∼10 to > 100) in either [O II] λ3727, Hβ, [O III] λ5007, Hα, which are well resolved, have inclinations i > 30°, are rotationally supported (v_max/σ > 1), and show no evidence for ongoing mergers. For the present analysis, we further restrict the sample to galaxies classified as ‘regular’ in Paper I. This classification was based on a visual inspection of rest-frame optical/IR morphologies and kinematic maps, selecting systems showing no obvious signs of gravitational interactions (e.g. tidal features, asymmetries, or disturbed velocity fields). In addition, we only selected galaxies that exhibit small residuals in the kinematic modelling, i.e. which have good fits, quantified by log(𝒵) < 15000, where log(𝒵) is the model evidence.¹

To robustly investigate the redshift evolution of the RAR, it is crucial to work with a sample defined by a uniform stellar-mass distribution across the redshift range considered, i.e. imposing the same minimum M_★ at all redshifts. Figure A.1 indicates that the MHUDF survey is complete down to a stellar-mass of log(M_★/M_⊙) = 8.6 over 0.2 < z < 1.44, whereas the sample from Paper I is complete to log(M_★/M_⊙) = 8.8 due to the additional S/N criteria that were imposed. We therefore adopt a stellar-mass cut of log(M_★/M_⊙) > 8.8, yielding a final sample of 79 SFGs (blue-red data points in Fig. A.1) used in this study. These galaxies span stellar masses and redshifts in the ranges 8.8 < log(M_★/M_⊙) < 11 and 0.33 < z < 1.44, respectively.

Fig. A.1. Stellar masses vs redshift for all MHUDF SFGs. The dashed black line indicates the mass completeness limit of the survey, which is reasonably complete for 0.2 < z < 1.5 down to log(M★/M⊙) = 8.6. The dashed red line shows the completeness of the parent sample from Paper I. The grey points show all SFGs observed in the MHUDF survey, the red-outlined points correspond to the galaxies analysed in Paper I, and the blue points indicate the galaxies studied in this work. The histograms on the top and right display the redshift and stellar-mass distributions, respectively, colour-coded to match the main panel.

Appendix B: Methodology

Based on the disk-halo decomposition from Paper I,

$$\begin{array}{c}\hfill v_{\mathrm{c}}^2(r)=v_{\mathrm{DM}}^2(r)+v_{\mathrm{disk}}^2(r)+v_{\mathrm{HI}}(r)|v_{\mathrm{HI}}(r)|(+v_{\mathrm{bulge}}^2(r)),\end{array}$$(B.1)

we calculated the baryonic acceleration due to the baryonic mass distribution at every radius as

$$\begin{array}{c}\hfill a_{\mathit{bar}}(r)=v_{\mathrm{disk}}{(r)}^2/r+v_{\mathrm{HI}}{(r)}^2/r(+v_{\mathrm{bulge}}{(r)}^2/r)\end{array}$$(B.2)

and the total acceleration as

$$\begin{array}{c}\hfill a_{\mathit{tot}}(r)=v_{\mathrm{c}}{(r)}^2/r.\end{array}$$(B.3)

Here v_disk(r), v_HI(r), and v_bulge(r) denote the contributions of the stellar disk, neutral atomic gas, and bulge² to the RC, respectively, while v_DM(r) represents the contribution of the DM halo. The circular velocity, v_c(r), is corrected for pressure support (asymmetric drift) following Dalcanton & Stilp (2010), such that v_c²(r) = v_⊥²(r)+v_AD²(r), where v_⊥(r) is the observed rotation velocity and v_AD(r) is the pressure support correction. We note that all velocity terms are intrinsic, i.e. corrected for any instrumental effects, including seeing.

We briefly summarise here the main ingredients of the 3D disk–halo decomposition (see Paper I for full details). The DM halo was parametrised using the best-fitting DM density profile from Paper I, namely DC14 (Di Cintio et al. 2014). For the neutral gas component, we used a parametric model which assumes a constant H I surface mass density, leading to $v(r)\propto \sqrt{{\mathrm{\Sigma }}_{\mathrm{HI}}r}$. The bulge, if present, was modelled as a Hernquist (1990) spheroidal component. The disk contribution was parametrised from the observed ionised gas surface brightness profiles using a Multi-Gaussian Expansion modelling (Emsellem et al. 1994), whereas the normalisation of v_disk(r) was given by M_★. Importantly, M_★ itself was not fixed by photometric priors. Instead, it was dynamically inferred in our 3D forward modelling. Within the DC14 framework, the disk-halo degeneracy is mitigated by fitting for log(M_★/M_halo) together with the virial velocity. This allows the stellar mass to be constrained directly by the kinematics.

For illustration, we show in Fig. B.1 the 3D disk–halo decomposition for two galaxies from our sample. We note that all catalogues and data products from our disk–halo decomposition, including the RCs, can be found on the DARK website.³ For MXDF 1266 (left, M_★,SED = 10^9.05 M_⊙), the RC is dominated by the DM component, indicating a strongly submaximal disk compared to local higher-mass spirals (e.g. Lelli et al. 2016). On average, we find that most of our sample has DM fractions larger than 50% within R_e (i.e. they are submaximal), as shown in Paper I. In contrast, the higher-mass galaxy, MOSAIC 905 (right, M_★,SED = 10^10.3 M_⊙), is dominated by the stellar component within 2R_e.

Fig. B.1. Examples of the 3D disk-halo decomposition for the two galaxies from the sample. The solid black line represents the circular velocity vc(r) corrected for pressure support. The solid green curve represents the RC of the DM component. The solid orange and blue curves represent the disk and cold gas components, respectively. The light shaded regions show the 95% CI. The dotted brown line represents the stellar component measured using the multi-Gaussian expansion modelling of resolved stellar mass maps obtained from pixel-by-pixel SED fitting, while the dash-dotted green curve shows the DM component obtained independently from a 1D disk–halo decomposition. All velocities are intrinsic, i.e. corrected for inclination and instrumental effects, including seeing (PSF). The black dotted line shows the physical extent of a MUSE spaxel.

The 1σ errors in the velocity components were derived from the 68% CI of the MCMC chains. Specifically, for each model parameter, we computed the 16th and 84th percentiles of the posterior distributions, providing the lower and upper bounds of the CI. To estimate the uncertainties in the baryonic and total accelerations, we propagated the errors from the velocity components using standard error propagation techniques.

To minimise the impact of observational uncertainties, particularly in regions where the rotation and acceleration profiles are poorly resolved, we excluded measurements in the innermost regions (r < 2 kpc) from the final fits to the RAR. Due to the finite spatial resolution of the MUSE instrument (0.2″spaxels, black dotted lines in Fig. B.1), measurements in these innermost regions are affected by resolution limitations that might bias our final results. A similar approach was taken in Vărăşteanu et al. (2025). We note, however, that excluding measurements within the central 2 kpc for both the MHUDF and SPARC (Lelli et al. 2017) samples changes the inferred a₀ by less than 10%.

Appendix C: Robustness tests and systematic uncertainties on disk masses

In this appendix, we assess the robustness of our results with respect to the main modelling assumptions. We emphasise that our measurement of the RAR at intermediate z is not strictly identical to the traditional local-Universe approach, as it relies on 3D forward modelling. While this technique represents the state of the art for intermediate- to high-z galaxy dynamics, it necessarily introduces a degree of model dependence, stemming from assumptions about the halo profile, disk geometry and gas distribution, among other factors. In this framework, both the total acceleration, a_tot(r), and the baryonic acceleration, a_bar(r), are derived from the same underlying model rather than independently, as is typically done at z ∼ 0.

Nonetheless, we note that in Paper I, we presented a series of consistency checks supporting the robustness of our disk–halo decomposition. Specifically: (i) The total circular velocity, v_c(r), was independently measured using the traditional 2D line-fitting algorithm CAMEL (Epinat et al. 2012), yielding results in good agreement with those obtained from the 3D parametric modelling (rendering v_c(r) effectively model-independent). (ii) The stellar disk contribution, v_disk(r), was independently derived via a Multi-Gaussian Expansion (Emsellem et al. 1994) modelling applied to resolved stellar-mass maps obtained from pixel-by-pixel SED fitting across ∼ 19 Hubble Space Telescope and James Webb Space Telescope bands, again showing good agreement (as illustrated by the brown dotted curves in Fig. B.1). (iii) For the unknown neutral gas contribution, v_HI(r), we explored a range of plausible H I surface-density profiles–given their lack of direct constraints–and found that our results are consistent across various assumed profile shapes. (iv) Finally, we performed a fully independent 1D disk–halo decomposition by computing the baryonic contribution to the RC. This was achieved by solving the Poisson equation for the surface-density distribution of the stellar and cold gas components (i.e. by solving Eq. 4 of Casertano 1983). This 1D approach yielded results consistent with those from the 3D modelling (as illustrated by the green dash-dotted curves in Fig. B.1).

Given the central role of M_★ in our results, it is important to note the dynamically inferred stellar masses are in good agreement with independent SED-based estimates (Fig. 11 of Paper I), showing that the potential systematic uncertainties in M_★ are negligible. In addition, we computed the rest-frame K-band mass-to-light ratios (M/L_K) of our sample, using photometric data from the Hubble Space Telescope and James Webb Space Telescope (see Paper I), and compared them to independent measurements at similar redshifts from Drory et al. (2004). Figure C.1(left) shows that the M/L_K of our sample are fully consistent with typical values reported at higher redshifts accounting for the mass dependence, and are systematically lower than those adopted for SPARC at z = 0, which have M/L_K ∼ 0.6 in the K-band (corresponding to M/L ∼ 0.5 at 3.6 μm, McGaugh et al. 2016).

Fig. C.1. Left: Average mass-to-light ratios (M/L) in the K-band as a function of redshift for the MHUDF sample split into two M★ bins, the low-M★ bin (yellow) and high-M★ bin (orange), along with the evolution inferred for a sample of more than 5000 K-selected galaxies in Drory et al. (2004). The red star shows the M/L-K ratio adopted for the SPARC sample. Right: Impact of systematic disk mass offsets on the inferred a0(z). The black points correspond to the fiducial fits shown in Fig 3. The blue points show the inferred a0(z) assuming a global shift of +0.2 dex in the mass budget of the sample, while the green points show the same, but when assuming a global shift of +0.45 dex . The red star and dotted lines show the canonical a0(z = 0) inferred for the SPARC sample.

To further assess the impact of systematic uncertainties in the disk mass estimates on the RAR, we also explore how global systematic shifts in the M_★ would affect the inferred a₀(z). Specifically, we recompute the RAR after applying uniform offsets to M_★ across the sample. We find that reconciling our measurements with the canonical value a₀(0) = 1.2 × 10⁻¹⁰ m s⁻² at all redshifts would require systematically larger stellar masses, with offsets ranging from ∼ + 0.2 dex in the lowest-redshift bins to ∼ + 0.45 dex in the highest (Fig. C.1, right). As argued above and discussed in Paper I, such large shifts are not supported by independent consistency checks.

We note that in order to improve our modelling, direct constraints on the gas components will be required. However, as discussed in Paper I, given typical molecular gas fractions of ∼30-50% at z ∼ 1 (Freundlich et al. 2019), this would introduce a systematic uncertainty of ∼0.2 dex in the total disk mass.

Finally, in Appendices D and E, we explored the RAR from RC decompositions using different parametric models in the ΛCDM and MOND frameworks, respectively, finding consistent results with those presented in the main text. Together, these tests demonstrate that our results are not driven by any single modelling assumption, but instead, most likely, reflect a genuine physical trend.

Appendix D: Robustness tests: RAR from different DM halo models in ΛCDM

To ensure that our conclusions on the redshift evolution of the RAR are not driven by methodological choices, we investigated how variations in our modelling assumptions influence the recovered accelerations a_bar(r) and a_tot(r). In Paper I we tested six DM density profiles in our disk–halo decomposition: (1) DC14 (Di Cintio et al. 2014), which links the halo profile shape to the stellar-to-halo mass ratio; (2) NFW (Navarro et al. 1997); (3) Burkert (Burkert 1995; (4) Dekel–Zhao (Freundlich et al. 2020); (5) Einasto (Navarro et al. 2004); and (6) coreNFW (Read et al. 2016). Our analysis showed that DC14 provides the best population-level fit, and we refer to the Bayesian model comparison in Paper I for details.

In this section, we adopted a galaxy-by-galaxy approach: for each system, we selected the profile that yields the best Bayesian evidence in the pairwise comparisons. Each galaxy was therefore assigned its own preferred model (e.g. DC14 for some, NFW for others). Using these best-fitting models for each galaxy, we recomputed a_bar(r) and a_tot(r) (as detailed in Appendix B). The resulting RAR tracks are shown as coloured points in Fig. D.1. We then refited the RAR for the full sample (Eq. 1), treating the characteristic acceleration scale as a free parameter (as in Sect. 3.1); the best-fit relation is shown as the purple curve. Additionally, in Fig. 2, we also show the best-fit value of a₀|_z ∼ 1 obtained in this framework as the cyan data point. We find

$$\begin{array}{c}\hfill a_0{|}_{z\sim 1}=2.{61}_{-0.09}^{+0.13}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2,\end{array}$$(D.1)

Fig. D.1. RAR of our intermediate-z sample, with abar and atot derived from the model RCs corresponding to each galaxy’s best-fitting ΛCDM halo profile. The coloured points indicate the halo profile used in the disk–halo decomposition. The cross in the lower right corner indicates the mean 1σ uncertainties on the data points. The dotted black line shows the 1:1 relation. The purple curve represents the best-fit RAR, with $a_0{|}_{z\sim 1}=2.{61}_{-0.09}^{+0.13}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$ (i.e. as in Fig. 1).

a result that is consistent, within the errors, with the value derived in the main text (Eq. 2), where the DC14 profile is applied uniformly to the full sample (purple data point in Fig. 2). Figure 2 illustrates that regardless of the model used to derive the total and the baryonic accelerations, the a₀|_z ∼ 1 value inferred from the RAR is offset from the z = 0 canonical value.

Following the same procedure presented in Sect. 3.2, we fit the RAR in this framework using the z-dependent acceleration scale (Eq. 1 substituting a₀ with Eq. 3), and obtain

$$\begin{array}{c}\hfill a_0(0)=1.{05}_{-0.05}^{+0.05}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2;a_1=1.{63}_{-0.12}^{+0.13}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2,\end{array}$$(D.2)

which agrees within 1σ with the results presented in the main text when using DC14 uniformly for the whole sample.

Overall, these tests demonstrate that our results are fairly robust against the choice of DM profile: whether using DC14 uniformly or adopting the best-fitting model for each galaxy individually, the recovered redshift evolution of the RAR remains consistent. In the next section, we discuss the RAR derived from the MOND framework.

Appendix E: Robustness tests: RAR from MOND models

To investigate the RAR in the MOND framework, we first need to refit our MHUDF data directly with MOND parametric models, and then re-analyse the RAR from these new fits as in Sect. 3. To refit our data, we extended our 3D modelling framework, GALPAK^3D (Bouché et al. 2015), to MOND (Milgrom 1983a) ⁴ using the interpolation functions defined in Sect. 3.3 of Famaey & Durakovic (2026). Specifically, we used two interpolation function families (Milgrom 1983c, McGaugh et al. 2016),

$$\begin{array}{cc}& \mathrm{n}-\mathrm{family}:{\nu }_n(x)={\left[\frac{1+\sqrt{1+4{(x/a_0)}^{-n}}}{2}\right]}^{1/n}and\hfill \end{array}$$(E.1)

$$\begin{array}{cc}& \delta -\mathrm{family}:{\nu }_{\delta }(x)={(1-exp(-{(x/a_0)}^{\delta /2}))}^{-1/\delta },\hfill \end{array}$$(E.2)

which relate the total and baryonic accelerations as

$$\begin{array}{c}\hfill a_{\mathrm{tot}}=\nu (a_{\mathrm{bar}})a_{\mathrm{bar}}.\end{array}$$(E.3)

Briefly, in our implementation within GALPAK^3D, the Newtonian baryonic acceleration a_bar(r) is computed directly from the disk, gas and, where applicable, bulge parametric models described in Paper I (and briefly detailed in Appendix B, and thus, does not depend on any MOND-specific parameters). The code then evaluates the total acceleration a_tot(r) using Eq. E.3. Again, corrections for pressure support, following Dalcanton & Stilp (2010), were applied to the RCs.

The free parameters of the model are the same as in Paper I, with the addition of the MOND-specific parameters, namely n or δ, depending on the interpolation function, and the acceleration scale, a₀. We adopt the same priors as in Paper I for all shared parameters. For the n-family model (Eq. E.1), we use a flat prior 0.7 < n < 1.3 , while for the δ-family model (Eq. E.2), we use a similar prior: 0.7 < δ < 1.3. These values are motivated by Desmond et al. (2024), who showed that the best fits to galaxy RCs were systematically close to n ∼ 1 or δ ∼ 1. We note that δ = 1 immediately gives back the usual form of Eq. 1.⁵ We find that the two MOND models (E.1 and E.2) perform similarly well, with a Bayes factor −2 < Δ log(𝒵) < 2, meaning that they are indistinguishable, and we therefore quote results for the δ-family (Eq. E.2, with δ ∼ 1, i.e. effectively Eq. 1) in the remainder of this appendix.

Before discussing the RAR in the MOND framework, we attempted to quantify the redshift evolution of a₀(z) from the individual fits to our MUSE data cubes with: (i) a₀ free in each galaxy, with a broad flat prior 10⁻¹¹ < a₀ < 10⁻⁹; (ii) a₀ fixed to the z = 0 canonical value; and (iii) with the redshift evolution a₀(z) (Eq. 3), where a₁ is a free parameter for each galaxy (with a broad flat prior −1 × 10⁻⁹ < a₁ < 1 × 10⁻⁹) and a₀(0) is kept fixed to the z = 0 value.

Regarding (i), a₀ was treated as a free parameter in the kinematic modelling of the individual MUSE cubes, and we do not recover a single, well-defined characteristic acceleration for the sample. Instead, the inferred a₀|_z ∼ 1 values display broad galaxy-to-galaxy variations, ranging from 1.3 × 10⁻¹¹ to 9.6 × 10⁻¹⁰ m/s², with a median value of a₀|_z ∼ 1 = 2.27 × 10⁻¹⁰ m/s² (close to the value inferred in the main text in the ΛCDM framework - Eq. 2). It is worth noting that similar results were obtained by Rodrigues et al. (2018), who used Bayesian inference on the SPARC and THINGS samples to show that the probability of a single, universal acceleration scale common to all galaxies is effectively zero (but see also McGaugh et al. (2018) on the role of inclination and M/L uncertainties and priors on this statistical result).

We tested whether the acceleration scales inferred from our fits to the individual data cubes are consistent with being at or below the canonical SPARC value, by performing a one-sided Kolmogorov–Smirnov (KS) test. The fiducial value was modelled as a Gaussian distribution centred on a₀ = 1.2 × 10⁻¹⁰ m/s² with a width equal to the measurement error 0.26 × 10⁻¹⁰ reported in McGaugh et al. (2016). Our measurement uncertainties were incorporated by generating 10000 Monte Carlo realisations of the sample. In every realisation, the null hypothesis is rejected at a confidence level of 99%, with a median statistic of D = 0.55 and a median p − value = 5.5 · 10⁻²⁰, indicating that the inferred values are systematically higher than the canonical a₀(0) value.

Regarding (ii), when a₀ is fixed to a₀(0) = 1.2 × 10⁻¹⁰ m/s² in the kinematic modelling, we found that the fits are worse, with the Bayes factor indicating positive evidence (Δlog(𝒵) < −2) for the models with a₀ as a free parameter. We find that the model with free a₀ (i.e. model i) performs as well or better than that for a fixed a₀ (i.e. model ii) for ≳ 70% of the sample. Similar results were obtained when comparing the model with z-varying a₀ (i.e. model iii) to the one with fixed a₀. We also note that the dynamically inferred M_★ using the model with fixed a₀ are overestimated by ∼0.28 dex with respect to the ones obtained in the main text using DC14, and are likewise elevated compared to independent SED-based estimates (although it does bring the M/L closer to those assumed in the SPARC database at z = 0).

Regarding (iii), we show the resulting a₁ values from the individual fits in Fig. E.1 (solid histogram). This figure reveals that the inferred a₁ values are predominantly positive, suggesting a positive correlation between the characteristic acceleration scale and redshift, and that the peak of the a₁ distribution is near the value (vertical line) obtained with the z-dependent fit to the RAR discussed in the main text (see Sect. 3.2: Eq. 4, Fig. 3). While individual a₁ values carry measurement uncertainties, we quantify the significance of the overall positive trend using a one-sample KS test (described below), which incorporates these uncertainties.

Fig. E.1. Distribution of the MOND acceleration parameter a1. The purple histogram shows the recovered a1 values from individual galaxy fits. The vertical red line marks the reference value obtained from the RAR in the main text, $a_1=(1.{59}_{-0.10}^{+0.11})\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$, using the DC14 framework (Eq. 4).

To quantify the statistical significance of this apparent offset from zero, we tested whether the a₁ values are consistent with being drawn from a zero-mean Gaussian distribution. This was achieved by performing a one-sample KS test against a standard normal null cumulative distribution function. Our measurement uncertainties were incorporated by generating 10000 Monte Carlo realisations of the sample. In every realisation, the null hypothesis is rejected at the 99% confidence level, with a median KS statistic of D = 0.30 and a median p-value of 1.2 × 10⁻⁵. This indicates that the distribution of a₁ is inconsistent with a zero-mean Gaussian.

Taken together, these results indicate that the individual fits of our intermediate-z galaxies do favour a positive a₁ value when fitted in a MOND framework, which agrees with the evolution inferred from the RAR when the data is fitted with DM halos (Sect. 3.2), hinting at a positive correlation between the characteristic acceleration scale and redshift.

Finally, for consistency, we re-analyse the RAR from these new fits, using the results from the MOND model (i), i.e. with a₀ free in each galaxy and δ ∼ 1, following the methodology presented in Sect. 3. First, we derive a_bar(r) and a_tot(r) as described in Appendix B from the individual model RCs. The resulting RAR is shown in Fig. E.2. We then fit the resolved RAR tracks of our sample with Eq. 1, leaving a₀ as a free parameter (using the same methodology as in Sect. 3.1), and we show the result as the purple curve in Fig. E.2 and as the blue point in Fig. 2:

$$\begin{array}{c}\hfill a_0{|}_{z\sim 1}=2.{19}_{-0.10}^{+0.12}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2.\end{array}$$(E.4)

Fig. E.2. RAR of our intermediate-z sample, where abar and atot are derived from the RCs modelled in the MOND framework. The cross in the lower right corner indicates the mean 1σ uncertainties on the data points. The dotted black line shows the 1:1 relation. The purple curve represents the best-fit RAR, with $a_0{|}_{z\sim 1}=2.{19}_{-0.10}^{+0.12}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$ (i.e. as in Fig. 1). The cyan curve shows the relation, as inferred using the median a0 value from the MOND fits to the individual data cubes, namely a0|z ∼ 1 = 2.27 × 10−10 m/s2.

The recovered a₀|_z ∼ 1 value is in fairly good agreement with the values of a₀|_z ∼ 1 obtained previously when the data were fitted with DM halos (Eq. 2, purple data point in Fig. 2), and is larger than the canonical z = 0 value (red star in Fig. 2). In Fig. E.2, we also plot as the cyan curve Eq.1 using for the characteristic acceleration scale the median value inferred for the sample from the individual MOND fits on the data cubes, i.e. a₀|_z ∼ 1 = 2.27 × 10⁻¹⁰ m/s², which close to the a₀|_z ∼ 1 inferred from the fit to the RAR.

Following the procedure presented in Sect. 3.2, we also fit the RAR extracted from the MOND framework using the z-dependent acceleration scale (Eq. 1 substituting a₀ with Eq. 3), and obtain

$$\begin{array}{c}\hfill a_0(0)=1.{03}_{-0.05}^{+0.05}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2;a_1=1.{20}_{-0.10}^{+0.10}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2.\end{array}$$(E.5)

To test whether this result is consistent with the values inferred from individual MOND fits, we also performed a direct linear regression of the best-fit a₀ values from each galaxy as a function of z, which yields a slope $a_1=1.{42}_{-0.89}^{+0.94}\times {10}^{-10}\mathrm{m}{\mathrm{s}}^{-2},$ and an intercept $a_0(0)=1.{11}_{-0.51}^{+0.39}\times {10}^{-10}\mathrm{m}{\mathrm{s}}^{-2}$, in fairly good agreement with the z-dependent fit to the RAR presented above (Eq. E.5), and with the values obtained in the main text (Eq. 4). Together, these results reinforce the presence of a systematic increase of the characteristic acceleration scale with redshift as we found previously when the data were fitted with DM halos (Sect. 3.2 and Appendix D).

We end this appendix by noting that these MOND fits underperform compared to the DM-based models (Appendix D) for over 60% of the sample using the Bayesian model comparison test as in Paper I.

Appendix F: Quantifying the z-evolution of the RAR

To complement the z-dependent RAR fit presented in Sect. 3.2, we show the full posterior distributions of all fitted parameters in Fig. F.1. The fit was performed on the whole sample using the MNR method implemented in Roxy (Bartlett & Desmond 2023), by replacing a₀ in Eq. 1 with the redshift-dependent form defined in Eq. 3. In addition to the primary parameters, a₀(0) and a₁, the corner plot displays the intrinsic scatter of the relation, σ_intr, and the Gaussian hyperparameters μ_gauss and w_gauss. Here, μ_gauss and w_gauss represent the mean and standard deviation of the Gaussian prior on the true log(a_bar) values. The hyperparameters are also marginalised over when inferring a₀(0) and a₁. The 2D contours reveal no significant correlations, confirming that a₀(0) and a₁ are well constrained.

Fig. F.1. Corner plot showing the posterior distributions of the parameters from the redshift-dependent fit to the RAR for the full sample. The acceleration scale is replaced by a0(z) = a0(0)+a1 ⋅ z in the MOND-inspired interpolation function (Eq. 1). The Gaussian hyperparameters μgauss and wgauss define the mean and width of the Gaussian prior on the true log(abar) values, as implemented in Roxy, while σintr represents the intrinsic scatter of the RAR. The units of a0 and a1 are 10−10 m/s2, while for μgauss, wgauss, and σintr the units are in dex. The contours indicate the 68% and 95% confidence levels.

This linear relation, using the best-fit values of a₀(0) and a₁ is shown as the purple line in Fig. 3. We note that this global fit is obtained from a single likelihood analysis of all individual data points and is not a regression through the a₀ values obtained from the binned analysis, i.e. the black data points shown in Fig. 3. Because data points enter the fit with their full measurement uncertainties, higher-z measurements typically carry less statistical weight, and small offsets between the global relation and the binned estimates are therefore expected.

All Figures

Fig. 1. RAR for the MHUDF sample. The purple curve shows the best-fit RAR at 0.33 < z < 1.44 using Eq. (1), while the red and green curves show the McGaugh et al. (2016) relation at z = 0 and the Vărăşteanu et al. (2025) relation at z < 0.08, respectively. The black dotted line indicates the 1:1 relation. The histogram in the inset displays the residuals with respect to our fit. The cross in the upper left corner indicates the mean 1σ uncertainties.

In the text

Fig. 2. Best-fit a0 values obtained by fitting Eq. (1) on the full sample, using the resolved RAR tracks derived under different modelling assumptions. The purple data point shows a0|z ∼ 1 obtained using the DC14 halo profile applied uniformly to the full sample (Sect. 3.1). The cyan point shows a0|z ∼ 1 derived from a galaxy-by-galaxy best-fitting ΛCDM model, where each galaxy is assigned the DM profile that maximises its Bayesian evidence (Appendix D). The blue point shows a0|z ∼ 1 recovered from the MOND-based framework (Appendix E). The red star indicates the fiducial a0|z ∼ 0 from McGaugh et al. (2016), while the green star shows the result from Vărăşteanu et al. (2025).

In the text

Fig. 3. Redshift evolution of the characteristic acceleration scale. In each quantile-based z bin (blue bars), we fit the RAR with Eq. (1) to derive a0, with uncertainties shown as grey error bars. The star marks the SPARC z = 0 value. The purple line shows the predicted evolution, using the best-fitting parameters from the global z-dependent fit to the RAR for the full sample (i.e. using Eq. (1) where a0 was substituted with Eq. (3); see Eq. (4)). The red line shows the same, but including SPARC in the fit. The coloured shaded region shows the 1σ errors for the fits.

In the text

Fig. A.1. Stellar masses vs redshift for all MHUDF SFGs. The dashed black line indicates the mass completeness limit of the survey, which is reasonably complete for 0.2 < z < 1.5 down to log(M★/M⊙) = 8.6. The dashed red line shows the completeness of the parent sample from Paper I. The grey points show all SFGs observed in the MHUDF survey, the red-outlined points correspond to the galaxies analysed in Paper I, and the blue points indicate the galaxies studied in this work. The histograms on the top and right display the redshift and stellar-mass distributions, respectively, colour-coded to match the main panel.

In the text

Fig. B.1. Examples of the 3D disk-halo decomposition for the two galaxies from the sample. The solid black line represents the circular velocity vc(r) corrected for pressure support. The solid green curve represents the RC of the DM component. The solid orange and blue curves represent the disk and cold gas components, respectively. The light shaded regions show the 95% CI. The dotted brown line represents the stellar component measured using the multi-Gaussian expansion modelling of resolved stellar mass maps obtained from pixel-by-pixel SED fitting, while the dash-dotted green curve shows the DM component obtained independently from a 1D disk–halo decomposition. All velocities are intrinsic, i.e. corrected for inclination and instrumental effects, including seeing (PSF). The black dotted line shows the physical extent of a MUSE spaxel.

In the text

Fig. C.1. Left: Average mass-to-light ratios (M/L) in the K-band as a function of redshift for the MHUDF sample split into two M★ bins, the low-M★ bin (yellow) and high-M★ bin (orange), along with the evolution inferred for a sample of more than 5000 K-selected galaxies in Drory et al. (2004). The red star shows the M/L-K ratio adopted for the SPARC sample. Right: Impact of systematic disk mass offsets on the inferred a0(z). The black points correspond to the fiducial fits shown in Fig 3. The blue points show the inferred a0(z) assuming a global shift of +0.2 dex in the mass budget of the sample, while the green points show the same, but when assuming a global shift of +0.45 dex . The red star and dotted lines show the canonical a0(z = 0) inferred for the SPARC sample.

In the text

Fig. D.1. RAR of our intermediate-z sample, with abar and atot derived from the model RCs corresponding to each galaxy’s best-fitting ΛCDM halo profile. The coloured points indicate the halo profile used in the disk–halo decomposition. The cross in the lower right corner indicates the mean 1σ uncertainties on the data points. The dotted black line shows the 1:1 relation. The purple curve represents the best-fit RAR, with $a_0{|}_{z\sim 1}=2.{61}_{-0.09}^{+0.13}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$ (i.e. as in Fig. 1).

In the text

	Fig. E.1. Distribution of the MOND acceleration parameter a1. The purple histogram shows the recovered a1 values from individual galaxy fits. The vertical red line marks the reference value obtained from the RAR in the main text, $a_1=(1.{59}_{-0.10}^{+0.11})\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$, using the DC14 framework (Eq. 4).
In the text

Fig. E.2. RAR of our intermediate-z sample, where abar and atot are derived from the RCs modelled in the MOND framework. The cross in the lower right corner indicates the mean 1σ uncertainties on the data points. The dotted black line shows the 1:1 relation. The purple curve represents the best-fit RAR, with $a_0{|}_{z\sim 1}=2.{19}_{-0.10}^{+0.12}\times {10}^{-10}\mathrm{m}/{\mathrm{s}}^2$ (i.e. as in Fig. 1). The cyan curve shows the relation, as inferred using the median a0 value from the MOND fits to the individual data cubes, namely a0|z ∼ 1 = 2.27 × 10−10 m/s2.

In the text

Fig. F.1. Corner plot showing the posterior distributions of the parameters from the redshift-dependent fit to the RAR for the full sample. The acceleration scale is replaced by a0(z) = a0(0)+a1 ⋅ z in the MOND-inspired interpolation function (Eq. 1). The Gaussian hyperparameters μgauss and wgauss define the mean and width of the Gaussian prior on the true log(abar) values, as implemented in Roxy, while σintr represents the intrinsic scatter of the RAR. The units of a0 and a1 are 10−10 m/s2, while for μgauss, wgauss, and σintr the units are in dex. The contours indicate the 68% and 95% confidence levels.

In the text