{'fraction_non_alphanumeric': 0.041501052139350014, 'fraction_numerical': 0.027436556896458696, 'mean_word_length': 4.256085456350144, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Upper Mass-Loss Limits and Clumping in the Intermediate and Outer Wind Regions of OB stars

We probe the radial clumping stratification of OB stars in the intermediate and outer wind regions (r>~2 R*) to derive upper limits for mass-loss rates, and compare to current mass-loss implementation. Together with archival multi-wavelength data, our new far-infrared continuum observations for a sample of 25 OB stars (including 13 B Supergiants) uniquely constrain the clumping properties of the intermediate wind region. We derive the minimum radial stratification of the clumping factor through the stellar wind, fclmin(r), and the corresponding maximum mass-loss rate, Mdotmax, normalising clumping factors to the outermost wind region (clfar=1). The clumping degree for r>~2 R* decreases or stays constant with increasing radius for almost the whole sample. There is a dependence on luminosity class and spectral type at the intermediate region relative to the outer ones: O Supergiants (OSGs) present a factor 2 larger clumping factors than B Supergiants (BSGs). The maximum clumping of roughly 1/3 of the OB Supergiants occurs close to the wind base (r<~2 R*) and then decreases monotonically. This contrasts with the more frequent case where the lowermost clumping increases towards a maximum, and needs to be addressed by theoretical models. Additionally, the estimated Mdotmax for BSGs is at least one order of magnitude lower than theoretical values, whereas for OSGs our results and predictions agree within errors. Assuming values of clfar=4-9 from hydrodynamical models would imply a reduction of mass-loss rates included in stellar evolution models by a factor 2-3 for OSGs and by factors 6-200 for BSGs below the first bi-stability jump. This implies large reductions of mass-loss rates applied in evolution-models for BSGs, independently of the actual clumping properties of these winds, and a thorough re-investigation of BSG mass-loss rates and their effects on stellar evolution.


Introduction
Wind outflows from OB Supergiants are the most widely studied examples of the radiatively driven-wind physical phenomenon (e.g., Friend y Abbott 1986, Pauldrach et al. 1986, Puls et al. 1996. These radiation-driven winds were first described theoretically by Lucy y Solomon (1970) and Castor et al. (1975), hereafter CAK, assuming a stationary, homogeneous and spherically symmetric outflow. By comparing observed spectral lines with expanding Non-Local Thermodynamic Equilibrium (NLTE) model atmospheres (e.g., cmfgen, Hillier y Miller 1998;fastwind, Santolaya-Rey et al. 1997 andPuls et al. 2005; powr, Gräfener et al. 2002) 1 , it is possible to derive both physical stellar and wind parameters by means of quantitative spectroscopy (e.g., Najarro 1995, Puls et al. 1996, Herrero et al. 2002, Urbaneja et al. 2003, Berlanas et al. 2018, Sander et al. 2014, Martins et al. 2019. However, despite the initial success of theoretical predictions based on stationary outflows (Vink et al. 2000, Kudritzki 2002, Puls et al. 2003, it is well established nowadays that stellar winds from massive stars are time-dependent and structured in velocity and density, dis-playing small-scale inhomogeneities (see review by Puls et al. 2008, Hamann et al. 2008, Sander 2017. These inhomogeneities change the structure of the atmosphere and wind, affecting the quantitative spectroscopic-based diagnostics used to obtain crucial stellar and wind parameters of massive stars. Hydrodynamical wind simulations have shown that the presence of strong instabilities 2 in the line-driven wind leads to formation of small-scale regions of very high densities (Owocki et al. 1988, Feldmeier 1995. These dense 'wind clumps' can be described using either the fractional volume of dense gas, the volume filling factor ( f v , Abbott et al. 1981), or via a clumping factor ( f cl , Owocki et al. 1988): where ρ 2 and ρ are the mean of the squared density, and the mean density, respectively, and where under the assumption of a negligible amount of gas in between the dense clumps one has f cl = f −1 v . Time-dependent simulations further show that the clumping factor across the wind is not homogenous, but presents radial stratification, f cl = f cl (r). Theoretical predictions agree on O Supergiants (OSGs) in that the clumping structure first increases rapidly in the wind-acceleration region, and then starts to fade again in regions far away from the star (Runacres y Owocki 2002. However, depending on the approximations and physical conditions used in the simulations (e.g., treatment of the source function, photospheric perturbations or dimensionality), the onset of the clumping is predicted further away (Runacres y Owocki 2002, Dessart y Owocki 2003 or closer to the photosphere (Sundqvist et al. 2011, Sundqvist y Owocki 2013, Driessen et al. 2019, which can affect formation of critical wind diagnostics as well as the spot where maximum clumping is achieved. The actual physical conditions close to the base wind (r 2 R * ) represent an open challenge in line-driven wind theory.
It has been suggested than a turbulent photosphere (e.g. Cantiello et al. 2009), may be lead to the formation on clumps at the base of the wind. Of course, this process would affect quantitatively the density structure of the wind, at least close to the stellar photosphere. However, it is not clear how the influence of a turbulent photosphere could persist far out in the wind, and there is not theoretical simulations to probe. Therefore, for the purposes of this work hydrodynamical LDI theory is used, since it has been widely tested successfully.
Since the presence of a density structure affects the matterradiation interaction across the wind, a key consequence of wind-clumping regards its impact upon the effective opacity of the wind (see detailed overview in Sundqvist y Puls 2018). For processes scaling linearly with ρ, the mean opacity of the clumped medium is the same as for the homogeneous wind model, whereas for processes scaling with ρ 2 , the mean opacity is enhanced by a factor f cl . Empirically, this means that clumping differently affects the spectral diagnostic used to derive wind parameters. Assuming that all clumps are optically thin for a given mass-loss rate leaves unaltered diagnostic X-ray lines, electronscattering wings, and scattering resonance-lines (ρ-dependent; e.g., C iv, P v), whereas it causes an opacity-enhancement of recombination lines or free-free continuum (ρ 2 -dependent; e.g., Hα, Br γ , mid-/far-infrared and radio continua). Therefore, if clumping is not appropriately taken into account, inconsistencies in mass-loss rate estimates between different diagnostics 2 The line-deshadowing instability (LDI). will arise (Fullerton et al. 2006, Cohen et al. 2010, Sundqvist et al. 2011. In particular, it is well established now that massloss rates of OB-stars derived from the classical diagnostic Hα line are overestimated by a factor f cl if unclumped models are used in the analysis. Moreover, if clumps become optically thick this leads to additional light-leakage through porous channels in between the clumps. Porosity can occur either spatially or in velocity-space due to Doppler shifts in the rapidly accelerating wind. Studies of velocity-space porosity have shown that clumps indeed very easily become optically thick in UV resonance-lines and that this effect is critical to include when analysing such P-Cygni line formation (Oskinova et al. 2007;Hillier 2008;Sundqvist et al. 2010;Surlan et al. 2013). On the other hand, such velocityspace porosity does not affect the continuum diagnostics studied in this work. And as shown by Sundqvist y Puls (2018), under typical OB-star wind conditions effects of spatial porosity upon the formation of mid-/far-infrared (MIR and FIR) and radio continua should be negligible. Nonetheless, also these diagnostics are subject to uncertainties related to wind-inhomogeneities, in particular the level of clumping in the inner to intermediate (MIR and FIR) and outermost (radio forming region) wind parts. These uncertainties, however, can be reduced by combining suitable diagnostics, effectively mapping the wind at different radial distances, from close to the base (V-band, Hα), over intermediate regions (Brα, MIR/FIR continua), to the outermost region (radio continuum). A consistent analysis would here place tight constraints on the radial clumping-stratification. Moreover, if used in combination with ionisation and velocity law diagnostics (resonance lines; e.g., P v), or X-ray emission (Cohen et al. 2010), then absolute clumping factors and mass-loss rates could be uniquely derived.
First efforts have been performed by Puls et al. (2006) (hereafter, P06) and Najarro et al. (2008Najarro et al. ( , 2011, by means of Hα, Brα and IR/mm/radio-continuum emission models for a sample of Galactic hot O stars. These studies show that whereas O Giants seem to have similar clumping degrees throughout the wind, OSGs show, on average, a clumping factor of about four times larger in the inner wind than in the outermost one. Thus, the mass-loss rate estimates for OSGs derived from Hα and smooth wind models should be scaled down by a minimum factor 2. The empirical clumping stratification derived by P06 and Najarro et al. (2011) suggested that the clumping degree starts to decrease at r ≈ 2 -6 R * . This overall agrees rather well with the predicted clumping stratification at r 20 R * by hydrodynamics LDI simulations of OSGs (Sundqvist et al. 2011;Sundqvist y Owocki 2013;Driessen et al. 2019). However, due to the scarcity of reliable continuum observations of the analysed OB stars at FIR wavelengths, the clumping factor in the intermediate wind region (2 R * r 15 R * ) remains poorly constrained.
To investigate the radial clumping stratification of OB stars at r 2 R * we obtained Herschel-PACS 70 and 100 µm fluxes (and 160 µm flux for the brighter sources) for a carefully selected sample of OB stars. The sample consists of 25 OB Galactic Supergiant, Giant and Dwarf stars, covering spectral types from O3 to B3 -4. Targets were selected to cover a wide OB-star parameter space, allowing us to analyse the behaviour of wind-clumping stratifications as a function of luminosity class and spectral type. Moreover, the large number of OB Supergiants in the sample allows for a deeper, statistically more significant analysis than previous work. In particular, our analysis represents a first attempt to derive the average properties of the clumping structure of B Supergiants (BSGs) by means of a multi-wavelength analysis, carefully comparing to the corresponding OSG results. In addition, the analysis allows us to derive reliable upper limits to mass-loss rates, and to compare these with the existing mass-loss rates recipes (Vink et al. 2000, 2001, hereafter V00 andV01) usually used in evolutionary tracks.
This represents a unique opportunity to test (at least in a relative way) theoretical mass-loss predictions for BSGs across the so-called bi-stability jump. Since the quantitative mass-loss rates across this jump are critical for stellar evolution modelling (e.g., Vink et al. 2010, Keszthelyi et al. 2017), this has a quite important impact also upon our general knowledge about the massivestar life-cycle.
The paper is organised as follows. Section 2 summarises the obtained Herschel-PACS observations and flux estimates of the sample, as well as the collected observations at different wavelengths from the literature. A description of the methodology, assumptions and procedures can be found in Section 3, whereas the performed analyses and corresponding results are presented in Section 4. We discuss our findings in Section 5, and present conclusions and a summary in Section 6.

Observations and data reduction
Although the main goal of this work is to constrain maximum mass-loss rates and degrees of clumping in the intermediate wind region of OB Supergiants, our initial 27 star sample also contains a few OB Giants and Dwarfs, including one LBV, two confirmed binaries, four early B-Hypergiants (eBHGs) and one magnetic star (see Table 1). This sample allowed us not only to analyse clumping stratification by luminosity class and spectral type, but also to investigate clumping properties of some peculiar objects for the first time. Finally, although blue supergiant stars are known to display photometric and spectroscopic variability, suggested to be linked to stellar pulsations, our sample has been carefully selected to encompass Blue Supergiants in long stable stages (Clark et al. 2012), with the exception of HD 198478, an extreme case of spectroscopic variability 3 , with an important scatter in the stellar and wind parameters for this source (see Table A.1). Nonetheless, we included this star in the sample in order to further analyse the uncertainty in our estimated values due to the use of certain sets of stellar and wind parameters (see discussion Appendix A.2), and to account for the effect of possible strong variability not detected so far for the rest of stars in our sample.
In addition, we collected archival data in the literature for V-band, near-infrared (NIR), MIR, mm and radio fluxes, when available. In Table 1 we summarise the sample of stars observed with Herschel-PACS and the data collected for this work. PACS column refers to our FIR flux observations. A description of the data reduction and the photometry processing of Herschel-PACS observations can be found in the following subsection, and the collected data are summarised in the successive sections.

FIR observations
FIR flux observations for our sample of 27 OB stars (Table 1) at 70 and 100 µm were taken with the Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) onboard the Herschel spacecraft (Pilbratt et al. 2010) in photometer mode (PI: Rubio-Díez, ID: OT1_mrubio_1, OT2_mrubio_2). A few of the objects in our sample have additional FIR-IRAS 60 and 100 µm data available, however most of them are upper limits (see references in Table 1).
Observations spans from December 2011 to March 2013, and were done using the mini-scan map observing mode. In this mode, two scan maps are taken along the two array diagonal directions (110 and 70 degrees) at a constant speed of 20 /s in parallel lines. The nominal spatial resolution for this scan velocity is 3.2 /pixel. Since our targets are point sources, a scan length of 2.5 arcmin was enough to cover the sky region of interest. The exposure times were estimated to reach a S/N ≥ 10, using theoretical emission flux estimates. Additionally, PACS operates simultaneously in the two sidebands characteristic of heterodyne receivers. Thus, for each source in our sample a total of eight scans were obtained, two mini-scan maps at 70 and 100 µm (blue and green band), and four scans at extra band 160 µm (red band). The PACS point spread function (PSF) for these photometric bands have a full-width at half-maximum (FWHM) default value of 5.2 , 7.7 and 12 , respectively.
The maps were processed using the map reconstruction task PhotProject with the Multi-resolution median transform (MMT) deglitching method implemented in HIPE (Herschel Interactive Processing Environment, Ott 2010). A second deglitching grade was applied. Subsequently, the processed mini-scan maps were combined into a final map ( Figure 1). We observed that whereas for the faint sources a second deglitching does not makes major changes in the final maps, for our brightest sources this step improved the accuracy of the recovered fluxes. Therefore, for the photometry in our sample stars we decided to use the second deglitching processed maps. Since the observing strategy focused on flux observations at the blue and green bands, the quality of the final extra red maps are generally poor and just a few of sources have been detected at this band, most of them giving upper limits. Table 2 displays the observed sources and the bands where they were detected. Only two sources, of the initial sample, have not been detected or reliably detected, HD 190603 and CygOB2#7, respectively. Thus, our final sample is conformed of a total of 25 OB stars (Tables 3 and 5).

Herschel/PACS Photometry
To perform the point source aperture photometry we used Hyper (Traficante et al. 2015). This routine was initially designed for FIR point source photometry in complex sky regions in the framework of the Herschel infrared survey of the Galactic plane (Hi-Gal). Because of its modularity and versatility, this code can be easily adapted to our Herschel-PACS observations, and it has been successfully used by different authors already (e.g., Benedettini et al. 2015, Svoboda et al. 2016, Li et al. 2018, Paulson y Pandian 2020. Hyper combines 2D multi-Gaussian fitting with aperture photometry to provide reliable photometry in regions with variable background, and in crowded fields. The 2D Gaussian fitting takes into account the beam observations and the source elongation to estimate the region over which to integrate the source flux, i.e., it computes the PSF of the source onto the final map. The background was locally evaluated and removed for each source using polynomial fits of various orders. The code selected the background based on the lowest r.m.s of the individual residual maps. For most of the sources, a polynomial with order not larger than 2 was used. In addition, in the case of blended sources, Hyper would perform a simultaneous multi-Gaussian fit of the main source and its companion(s), subtracting the modelled companion(s) from the science target. This was no the case for any of our targets. Finally, although this code allows simul- Table 1. List of the OB stars observed with Herschel-PACS (IP: OT1_mrubio_1/OT2_mrubio_2) and the photometric data compiled from the literature used in this work (see references below).

Source
The error in the estimated final flux includes the uncertainty in the maps calibration from PACS (10%), the error from the photometric procedure (from Hyper outputs) and from the aperture and colour corrections applied to the extracted flux. These corrections are needed to i) account for the missing flux due to the finite aperture; and ii) convert monochromatic flux densities in the PACS data products, which refer to a constant energy spectrum, to the true object Spectral Energy Density (SED) flux densities at the reference band wavelengths 70, 100 and 160 µm, using the tabulated aperture and colour correction factors to the extracted PACS flux as 1/ f apc f cc 4 .
Finally, attending to the S/N we flagged flux values with quality flags (qF λ ) as: 'A' for detections with S/N ≥ 10, and 'B' for detections with 5 ≤ S/N < 10. In addition, measured fluxes with S/N ≤ 5 or affected by image artifacts, are flagged as upper limits (U). In Table 2

V/NIR and MIR observations
The V, NIR and MIR magnitudes used in this work are listed in Table 1 and were mostly collected from the literature (see references within). All sources are observed in the V and NIR bands, most of them also at MIR wavelengths. Since these observations were made with different instruments (different photometric systems) we performed an absolute flux calibration. We used Vega as a calibrator to convert the magnitudes into absolute fluxes by extrapolating the visual absolute flux calibration of Vega to a specific wavelength by model atmosphere (Kurucz models) and calculating the absolute fluxes of Vega in the different photometric systems. A detailed description of the absolute flux calibration procedure can be found in Puls et al. (2006).

Millimeter and Radio Observations
Millimeter and radio observations were collected from literature (see Table 1 and references therein). Most of our sample stars were observed with the Very Large Array (VLA) at 2, 3.5 and 6 cm or with the Australia Telescope Compact Array (ATCA) at 3.6 and 6 cm in different epochs. In addition, several objects were also observed at 1, 3, 13, 20 and 90 cm (J-VLA) and at 21 cm (e-MERLIN). Only for one object, HD 53138, no radio observations are available, whereas for 7 of the 25 targets only upper limits could be obtained.
With respect to the sub-mm and mm flux continuum measurements, we collected from the literature those obtained at 0.7, 0.85, 1.3 and 1.35 mm, when available. Only 7 stars in our sample have millimeter observations (Submillimetre Common-User Bolometer Array, SCUBA). Table 2. Herschel/PACS 70, 100 and 160 µm fluxes for our sample stars. Sources are listed as in Table 1. The 1σ error of the flux includes the uncertainty in flux calibration from PACS, the error in the photometric procedure and the error arising from the aperture and colour corrections (see Sec. 2). The quality of each of these flux densities is designated by the quality flags (qF λ ) and is either of high quality (qF λ = A; S/N ≥ 10), of moderate quality (qF λ = B; 5 ≤ S/N < 10), or an upper limit (qF λ = U; S/N ≤ 5 and/or confused detection).

Source
F 70µm ± σ (Jy) qF 70µm F 100µm ± σ (Jy) qF 100µm F 160µm ± σ (Jy) qF 160µm This work uses distance determinations based on GAIA DR2 parallaxes (Luri et al. 2018, Gaia Collaboration et al. 2016, 2018. We only used parallaxes where the associated errors are lower or similar to 15% (8 out of 25 objects; marked with an asterisk in Table 3), otherwise the obtained distance becomes unreliable when compared to previous estimates (see below). The differences between GAIA DR2 and older distance estimates ranges from 0.01 kpc to ≈ 1 kpc, with an average of 0.47 kpc. The largest deviations correspond to HD 169454, HD 80077 and CyOB2#12, whose differences with the new values are ∼ 0.6, 1.05 and 0.9 kpc, respectively. In particular, for CygOB2#12 the distance derived from GAIA DR2 parallaxes (0.85 kpc) puts the star much closer than its association, Cygnus OB2 (1.75 kpc). This casts some doubts about its physical properties (Nazé et al. 2019). For this source we provide the results of the analysis for both distances (two entries in Tables 6, 7 and 9), although only the result pertaining to the larger, association distance (1.75 kpc) is considered in the discussion 5 .
The case of HD 80077 is different. Since there were always severe doubts about this star being a membership of the Pismis 11 cluster (located at 3.6 kpc; Marco y Negueruela 2008), and given the lower uncertainty in the parallaxes by GAIA DR2, we keep this distance (2.55 kpc) for the analysis.
For the rest of the sample, initially we used photometric distances from the literature (see Tables 3 and 5), whose uncertainties range between 0.1 and 0.4 kpc 6 , and locate most of the stars in their associated cluster. For HD 152236 (τ 1 Sco) we assumed a distance of 1.64 kpc (Sana et al. 2006), used in previous studies (Clark et al. 2012).
In Table 3 we present the final used distances and the derived stellar radii and extinction parameters for the sample as described in Section 3. For HD 66811 (ζ Pup), in order to be consistent with previous analyses by P06 (see Table 5), we used two different distances: a commonly used shorter distance (0.46 kpc), assuming ζ Pup is part of Vela OB2 association, and a larger one (0.73 kpc, Sahu y Blaauw 1993), which corresponds to the star being runaway. Although we provide the results of the analysis for both distances (two entries in Tables 6, 7 and 9), the discussion in this work refers to the shorter distance, unless otherwise specified. Table 3. Final used distances and associated errors, reddening parameters and derived stellar radii. Sources are sorted as in Table 1. Asterisked values in the distance column correspond to GAIA DR2 distances, and the two values for HD 66811 and CygOB2#12 refer to different alternatives for these sources (see Sec. 2.4), quoted with their associated asymmetric errors. Some parameters have been modified from the initial values found in the literature (see Table 5

Modeling
To investigate the wind clumping stratification and mass-loss rates of the 25 OB stars in our final sample, we modeled their observed spectral energy distribution (SED) from V/NIR to radio wavelengths. We used the interactive procedure developed by P06, which is based on continuum emission models (photospheric plus wind emission), assuming a spherically symmetric wind calibrated against full NLTE stellar atmosphere models, and accounting for optically thin clumping. Although the basic method thus neglects potential effects from porosity (see Sec. 1), these should be relatively small for the diagnostics and spectral ranges considered in this work (see Sundqvist y Puls 2018 for a discussion). Below we present a summary of the approach to model SEDs depending on the stellar and wind parameters of the objects. An in-depth description and verification of the method can be found in P06.

Infrared and radio flux emission
The infrared and radio fluxes are calculated using the approximations described by Lamers y Waters (1984). Following P06, we use: a) Wind velocity law.
The electron temperature was computed using Lucy's temperature law (Lucy 1971) with a lower temperature cut-off at 0.5 T eff . c) Ionisation equilibrium. For all the objects in our sample, hydrogen is considered to be completely ionised, whereas the helium ionisation structure used in the modeling of the SED depends on the temperature of the source and the wavelength domain (for the rationale see P06). Specifically, for the majority of the sample (17 000 K T eff 35 000 K) helium is considered singly ionised in the radio regime whereas for the coolest (T eff 17 000 K) and hotter objects ( T eff 35 000 K), it is assumed to be neutral and fully ionised, respectively. In the NIR, helium is assumed to be fully (T eff 32 500 K) or singly (13 000 K T eff 32 000 K) ionised. A particular treatment is required for P Cyg (HD 193237) whose ionised helium structures depart from the above standard scalings with effective temperature (Najarro y Figer 1998). As such, for P Cyg we assumed that helium is fully recombined in the outer parts of the wind (radio domain) and singly ionised in the inner ones (NIR domain).

d) Photospheric input fluxes.
For λ < 1 µm, where the stellar photosphere dominates, the emitted flux is modeled using Kurucz's fluxes. On the other hand, for λ > 1 µm, where the wind starts to dominate the resulting flux, a black body emission model is used.

e) Wind Clumping treatment.
Clumping is included using the following approach: all material in the wind is redistributed into clumps which are over-dense with respect to the average density, we assume optically thin clumps (see above), and further that the inter-clump medium is effectively void (Abbott et al. 1981, Schmutz 1995. Under these assumptions, the spatial mean density ( ρ = M/4πr 2 υ) and mean squared density can be expressed as a function of the volume filling factor ( f v ) as: where ρ + denotes the density inside the clump. Thus, the clumping factor ( f cl ) as described by Eq.1, becomes the inverse of the volume filling factor, describing the over-density of the clumps as compared to the mean density. Within this approach, the opacity depends only on the matter and the physical processes inside the clumps (recombination, scattering, absorption, etc.). Since the mean opacity of the processes with opacities depending on ρ 2 is enhanced by a factor f cl Article number, page 7 of 35 A&A proofs: manuscript no. AA202040116 (see Sec. 1), the optical depth invariant for thermal emission, Q (Lamers y Waters 1984, P06), is modified as: Therefore, since the opacities of all processes considered in this analysis depend on ρ 2 (IR/mm/radio emission), the effects of clumping can be included in the models by multiplying the opacities evaluated for the mean wind by a clumping factor. Thus, the theoretical-emitted flux (Panagia y Felli 1975, Wright y Barlow 1975, for such processes might be extended according to: Consequently, the emitted flux is sensitive to changes in clumping factors, allowing us to constrain the clumping structure of the wind by fitting emission models to flux observations for a given source. Since clumping is not constant throughout the wind but presents a radial stratification, f cl (r/R * ), it is to be evaluated in different regions of the wind.  For that objective, the wind is divided in different regions with corresponding clumping factors as defined in Table 4.
Typical boundaries are r in = 1.05 R * , r mid = 2 R * , r out = 15 R * and r far > 50 R * , which roughly agree to the formation zones of Hα and NIR (Region 1 -2), MIR/FIR (Region 3), mm (Region 4) and radio (Region 5). Note that i) the clumping at the base of the stellar wind (1 < r/R * < r in ) is set to 1 within the fitting procedure, and ii) that the above-mentioned limits for the wind regions in Table 4 can be adapted within the fitting procedure when necessary (see below). Figure 2 displays a sketch of the defined wind regions as a function of the spectral energy distribution.
The parametrisation adopted in this work is designed to empirically constrain the clumping factor in different radial wind zones. Therefore, the derived clumping factors should be considered as average values describing the global behaviour of clumping throughout the wind. Since all used diagnostics depend on ρ 2 , it is not possible to derive absolute values for clumping factors and mass-loss rates, but only relative ones. What we really obtain is the scaling invariant Q ∝Ṁ f cl for a given radius interval. To derive absolute values, a simultaneous analysis of ρ-dependent diagnostics (e.g. resonance lines, see Sec. 1) would be required to break the degeneracy, which is out of the scope of this work. However, since all derived (minimum) clumping factors and (maximum) mass-loss rates can be scaled via Q , it is still possible to carry out important comparisons with theoretical predictions as well as with previously derived empirical values.
It is worth noting that our fitting-results are independent of the individual values for distance and stellar radius, and that they can be easily scaled to different values as long as Q' and (R * /d) remain conserved (see Eq. 9). Of course, changes in distances (and thus in stellar radius) affect directly the derived absolute values of mass-loss rates, but it does not affect the behaviour of the clumping structure (e.g. HD 66811 and CygOB2#12, see Table 6).
All comparisons with other empirical studies and theoretical predictions were performed for the same values of distance and stellar radius, for a given source. Thus, for comparisons with other empirical studies, their mass-loss rates were scaled to the distance and stellar radius used in our analysis, whenever necessary. Similarly, stellar parameters used to compute theoretical predictions forṀ (Section 5.1.2) are the same as those used to derive the empiricalṀ max estimates. Regardless of variations in absolute values of mass-loss rates with changing distance to the objects, the general behaviour found in this work will thus be unaffected, and the ratio between empirical and theoretical massloss rates is remains conserved. For example, we estimate that a 40% increase in distance 7 leads to a decrease in empirical to theoretical mass-loss rates of around 30%. In other words, such changes in distance would not change the overarching trends and conclusions found in this work.

Fitting Procedure
We compare continuum emission models with the observed SED for the 25 OB stars in the sample. In Table 5 we summarise the stellar and wind parameters of the sample from the literature (references therein), used as input values in our simulations. The best-fit model for each source is obtained from a simple maximum likelihood method (χ 2 ) as follows: First, we de-reddened the observed flux using the extinction law provided by Cardelli et al. (1989). By comparing the ob- Table 5. Sample of analysed OB stars in this work and the initial set of stellar and wind parameters from the literature ('ref' column) used as input values in our simulations. log g indicates gravities without centrifugal correction. For HD 66811 (ζ Pup) we provide two entries, based on different distances (P06).

Source
Spec.Type Secondly, we computed the stellar radius. For a given distance, the initial value of the stellar radius is adapted to match the VJHK de-reddened observed fluxes. Since flux is diluted by (R * /d) 2 (see Eq.9), the ratio between radius and distance represents a scaled flux factor, which has to be conserved for any distance. In Table 3 we present the final used distances, and the corresponding derived stellar radiii and de-reddening parameters.
Thirdly, we estimated mass-loss rates required to reproduce the continuum distribution across the electromagnetic spectrum. Since the optical depth and the diluted flux scale withṀ/R * 3/2 , this ratio has to be conserved for a given (R * /d).
Initially, all clumping factors in our simulations are set to the minimum value, f in cl = f mid cl = f out cl = f far cl = f min cl =1 (unclumped wind) and the mass-loss rate is adapted to reproduce the observed fluxes throughout the wind (1.05 R * < r < 100 R * ). Note that the emitted flux decreases on average with increasing wavelength and, therefore, radio thermal emission represents the lowest absolute flux values for a certain source 8 . As such, the de-rivedṀ represents the maximum possible mass-loss rate (Ṁ max ) consistent with radio observations assuming a minimum clumping condition ( f far cl = f min cl = 1). In the cases where radio fluxes are well-defined but somehow larger than at shorter wavelengths, the estimatedṀ max relies instead on the consistent observed fluxes at those wavelengths (usually the FIR domain). The same is true for those sources whose radio observations are not well constrained, although in this cases the adoptedṀ max is only an upper limit.
Finally, the radial stratification of the clumping is derived by adapting the values of the clumping factors to match the continuum emission. The best-fit solution for each source corresponds to the emission model (Ṁ max , R * , d and f in cl , f mid cl , f out cl , f far cl ) providing the minimum value of the statistics fitting function χ 2 . Unlike P06, we follow two different approaches: i) Fixed-regions approach. Here, we fix the defined radial boundaries (see Table 4) for all sources. This allows us to investigate consistently the global behaviour of the clumping stratification for our sample, comparing the derived average clumping factors across the wind regions. ii) Adapted-regions approach. In these simulations, the limits of the wind regions can be adapted, in addition to the clumping factors -when necessary-to secure the best possible fit compared to fixed-region approach. This allows us to probe also the boundaries defined above for the different wind regions.

Additional considerations
In order to be consistent, besides setting r out = 15 R * in the fixedregions approach in our simulations, we also probe r out = 10 R * , since this limit was used by P06 when analyzing low or intermediate density winds (where, in their case, the profile shape of Hα served as a proxy for wind-strength.) Additionally, due to the relative scarcity of observations in the sub-mm regime, we computed the maximum clumping values for Region 4 still compatible with flux emission at shorter and longer wavelengths ( f out max in Tables 6 and 7). For the sources previously analysed by P06, and with new distance estimates from GAIA (asterisked values in Table 3), we first tested the de-rivedṀ max and clumping stratification by these authors against the obtained FIR flux values and the additional measured fluxes at MIR and radio ranges found in the literature since P06. In a second step, we scaled the stellar radius and the mass-loss rate to the new distance according to Eq.8 and derived the best-fit solution to all data as described above.

Prototypical examples
In this section we briefly describe three prototypical examples of theṀ max and clumping factor estimates through the fitting procedures in the fixed-region approach. They cover different possibilities with respect to the flux values, and how arising peculiarities are reflected in the table entries displaying the results for the whole sample (Tables 6,7,9). From top to bottom Figure 3 displays observed and synthetised fluxes for HD 151804, HD 193237 (P Cyg), and HD 24398 (ζ Per). For a detailed discussion of every star in the sample, see Appendix A.

HD 151804
This is one of the most straightforward cases, where the radio flux is a unique, determined value providing a well con-strainedṀ max . ForṀ max = 6.4 × 10 −6 M yr −1 , the synthesised flux perfectly matches the observed values across the wind without the need for clumping ( f in cl = f out cl = f far cl = 1), except in the intermediate wind region (Region 3). The synthesized fluxes with all clumping factors set to unity are displayed as a magenta-dotted line. To reproduce the observed fluxes in the FIR, we increase to f mid cl = 3.2 (best-fit, solid line). Although the best-fit for this star has f out cl = 1, given the observational gap in the mm-regime a maximum value of the clumping degree in Region 4 f out max = 8 (blue dashed-line) is still consistent with the data at shorter wavelengths. This f out max is typed between parentheses in the corresponding entry for this star in Table 6.

HD 193237 (P Cyg)
This known LBV is well sampled throughout the wind. It is well established that the observed variability on timescales on the order of months at mm and radio wavelengths (Abbott et al. 1981, van den Oord et al. 1985, Bieging et al. 1989, Scuderi et al. 1998, Ofek y Frail 2011, Perrott et al. 2015 is not related to P Cyg being a non-thermal source, but rather to changes in the ionisation stage of the outer wind (recombined outermost wind model, Najarro et al. 1997; see also Pauldrach y Puls 1990). In addition, the observed larger IRAS and SCUBA flux measurements at 60 and 850 µm may be due to spatial resolution effects (such as contamination by the characteristic LBV nebula surrounding P Cyg), whereas the lower value at 6 cm observed by Bieging et al. (1989) has not been confirmed by newer measurements. In view of this, we estimatedṀ max for P Cyg for being consistent with all measured fluxes from 2 to 21 cm, regardless of the upper IRAS 60 µm and SCUBA 850 µm values and the lower flux at 6 cm, obtainingṀ max = 12.8 × 10 −6 M yr −1 .
With this mass-loss rate the synthesised flux underestimates the data at shorter wavelengths (magenta-dotted line). Thus, the clumping degree needed to be increased to reproduce the observed SED, yielding f in cl = 3, f mid cl = 3.5, f out cl = 5 and f far cl = 1 (best-fit, solid line). To match the larger well-defined fluxes observed in the radio-regime the clumping degree in Region 5 had to be increased to f far cl = 1.5 (blue-dashed line). The corresponding clumping factors for these two models for f far cl are separated by a forward slash (/) in Table 6, with the number on the left referring to the best-fit solution. The well-sampled mm-regime of this star further allows us to constrain the clumping degree in Fig. 3. From Top to Bottom, observed and best-fit fluxes vs. wavelength in the fixed-regions approach for HD 151804, HD 193237 (P Cyg), and HD 24398 (ζ Per). Solid lines correspond to the best-fits; magentadotted, blue-dashed and orange-dashed-dotted lines correspond to different models (see text).
Region 4 to f out max = f out cl = 5.
HD 24398 (ζ Per) Only an upper limit of the flux at radio wavelengths is available for this star. An estimate ofṀ max by matching this limit is not consistent with FIR flux values, since any model with suchṀ max would overestimate the FIR flux (magenta-dotted line). Instead, an upper limitṀ max 0.18 M yr −1 derived from the FIR flux values is consistent with the radio upper flux limit. In this case the observed SED is perfectly matched by the minimum clumping degree throughout the whole wind, f in cl = f mid cl = f out cl = f far cl = 1 (best-fit, solid black line). It is possible, however, to still match the upper radio flux limit if increasing up to f far max = 5 (parenthesed value in Table 6; blue-dashed line). Although the best-fit for this star is for f out cl = 1, given the observational gap at λ > 100 µm, f out max = 45 (orange-dashed-dotted line) is still consistent with the data (also between parentheses in Table 6).

Analysis and results
In this section we present the derived maximum mass-loss rates, the corresponding minimum radial clumping stratification, and their associated uncertainties.
The results obtained in the fixed-and the adapted-regions approaches (see previous section) are summarised in Tables 6 and  7. A detailed description of the fits for the individual objects can be found in Appendix A.

General findings
Overall, we find that the minimum χ 2 of the fits in the fixed-regions approach is comparable to that obtained with the adapted-regions approach. As such, unless specified otherwise, the following general findings refer to the best-fit solutions in the fixed-regions approach (Table 6). A definiteṀ max was derived for 17 of the 25 sources. For the remaining 8 objects, we only provide upper limits ofṀ max consistent with the observations, due to lack of detections in the radio regime.  Figure 4 displays log Q = log (Ṁ max /R * 1.5 ) as a function of spectral type and luminosity class. We can see that log Q decreases with luminosity class and spectral type, and that the values departing from this trend correspond to the non-standard sources in our sample, i.e. confirmed binaries, magnetic stars, LBVs and eBHGs (see coloured symbols in Figure 4). In particular, this figure shows that non-standard OB stars seem to be displaced towards higher log Q values and later spectral types when compared to standard OB Supergiants 9 .
Due to the low number of OB Giants and Dwarfs in our sample, it is not possible to clearly observe a similar behaviour as that seen for the OB Supergiants. Overall, the general trend and the estimated values ofṀ max and log Q for the complete sample appear consistent. Indeed, the large values obtained for some of the non-standard objects agree with some previous estimates (Clark et al. 2012, Crowther et al. 2006. For instance, for the eBHG star HD 152236, we derivedṀ max = 6.2 × 10 −6 M yr −1 , a value close to the (unclumped) mass-loss rates obtained by these authors, 6.33 × 10 −6 M yr −1 and 6.0 × 10 −6 M yr −1 , respectively.
The derived minimum clumping structures derived for r 2 R * are displayed in Table 6. Only one source, CygOB2#11, remains unconstrained in the intermediate wind region ( f mid   cl 3) since only upper flux limits were detected at 70, 100 and 160 µm. In Figure 5 we present the ratio of the clumping factors in the intermediate and outer part of the wind ( f mid cl / f out cl ; left panels) and in the intermediate and outermost wind regions ( f mid cl / f far cl ; right panels) as a function of the wind density (log Q ), for the OB Supergiants (top and middle panels, respectively) and for the OB Giants and Dwarfs (bottom panels). From top to bottom, Figure 5 shows that the clumping degree decreases or remains nearly constant with increasing radius throughout the wind regardless of luminosity class or wind density. Overall, the minimum clumping degree in Region 4 ( f out cl ) is similar to that in the outermost Region 5 ( f far cl ). OB Giants and Dwarfs are the sources with generally the lowest, most homogeneous clumping across all wind regions. OB Supergiants, on the other hand, display a more radially structured clumped wind.
There are four exceptions to these general trends, where the clumping structure seems to increase with increasing distance to the photosphere. These are P Cyg (HD 193237), Cy-gOB2#12, ι Ori (HD 37043) and τ Sco (HD 149438). Note that all these are non-standard objects: respectively, a LBV, an eBHG with a recently measured variable radio flux at 21 cm (Morford et al. 2016), a confirmed binary, and the well-known magnetic B Dwarf τ Sco.

Clumping in the intermediate wind of OB Supergiants
The relatively large number of OB Supergiants in the sample allows for a statistically significant analysis of the average radial clumping structure of these stars at r 2 R * . To do this, we use the fixed-regions approach and remove from the analysis those OB Supergiants with a peculiar nature and/or whose derived clumping properties significantly depart from the general trend (HD 193237, Cyg0B2#12, HD 149404, HD 198478; see Table 6). Due to the large uncertainty in PACS flux estimates for HD 194279 at 70 and 100 µm, the scarcity of data at mm wavelengths and the lack of evidence of it being a non-thermal emitter, we consider the first derived solution to be more likely for this target (first entry in Table 6; for discussion see Appendix A.2). Figure 6 displays the individual (top panels) and the average (bottom panels) minimum radial clumping stratification for the OSG and BSG subsamples (left and right, respectively). The OB Supergiants show a similar or larger clumping degree in the intermediate wind region than in the outermost ones (Regions 4 and 5). We estimate that for OSGs the average minimum clumping factor for the intermediate wind region is f mid cl ∼ 4.0 ± 1.5, whereas for BSGs it is f mid cl ∼ 2.1 ± 1.1 (errors represent standard deviations).
Concerning the outermost wind regions, due to the scarcity of sub-mm and mm observations it is possible to derive precise clumping properties in Region 4 for only 7 sources (boldfaced names in Tables 6 and 7). For the rest of the objects, f out cl = 1 generally provides the best-fit solution, but larger values may also be possible (see f out max col. in Table 6). We compute an average minimum clumping factor for OSGs and BSGs of f out cl = 1 (+2) and 1 (+1.3), respectively. Unlike for f mid cl , the number in parentheses here do not represent formal errors, but are instead rough indications based on the average ratio f out max of those 7 OB Supergiants whose clumping properties are fully constrained.

Clumping factor uncertainties
The mass-loss rate and clumping factor derivations have been made using literature-values of β and υ ∞ .
As discussed above, the synthesized and observed fluxes depend mostly on the local quantitiesṀ f cl (for a given R * ). Assuming an (almost) perfect fit, the derived clumping factors scale inversely with the mass-loss rate as δ f cl / f cl ≈ −2 δṀ/Ṁ. Taking into account that uncertainties inṀ max depend directly on the errors in radio fluxes (temporal variability and/or intrinsic errors) with 2 δṀ/Ṁ ≈ 1.5 δF ν /F ν , it follows that the uncertainty in clumping factors also scales inversely with flux as δ f cl / f cl ≈ −1.5 δF ν /F ν . We computed an average uncertainty ∼ 15% forṀ max due to flux uncertainties, which leads to an average uncertainty in derived clumping factors of ≈ 30%. For those sources with uncertainties inṀ max above average, the errors in the clumping factors can reach 50%.
Other factors affecting the derivedṀ max , and thus f cl , are the assumed helium fraction and the ionisation stage. We performed several simulations varying both parameters. We found that a factor 2 difference in the helium fraction leads to a variation iṅ M max of ∼ 16%, quite similar to the uncertainty found (∼ 15%) when the ionisation state in the region probed in the mm and radio domain is different than assumed. From these simulations, we estimate that the average uncertainty of the clumping factors derived in this work, due to uncertainties inṀ max , is in between 30% to 40%.

Clumping dependence on β
Whereas the clumping degree in the outermost wind regions (r 15 R * ) does not depend on the exponent of the velocity field, β, the opposite is true for the inner wind (r 2 R * ). Since the β and υ ∞ used in our analysis come from a range of different studies (see references in Table 5), they are not really coherently derived. Thus, some analysis regarding the impact of β on the derived clumping factors is warranted, and carried out below using two experiments. In the first one, we computed the boundaries of the validity interval of β for which a given best-fit model still Table 6. Maximum mass-loss rates and minimum clumping factors for the sample as derived in the fixed-regions approach (see Sec. 3). Sources are sorted as in Table 1. R * corresponds to the estimated stellar radius for the assumed distance (see Table 3), whereas β and υ ∞ are taken from the literature (references in Table 5). Clumping factors denoted by f in cl , f mid cl , f out cl , f far cl , and r in = 1.05 R * , r mid = 2 R * , r out = 15 R * and r far > 50 R * are corresponding wind-region (cwr) boundaries. The maximum possible clumping factor in Region 4 still consistent with the data denoted by f out max . Boldfaced names refer to stars with the best constrained clumping factors for all wind regions. Values separated by a forward slash (/) indicate alternate well-defined solutions, with the best-fit solution being always the number on the left. Values in parentheses refer to upper limits of the maximum clumping factor in the cwr (see Sec. 3.4  (1) Value interval forṀ max and clumping factors derived from the stellar parameters in the literature (see Table 5 and Appendix A.2 for further discussion.) (a) Upper limit ofṀ max , derived from upper radio-flux limits.  Table 7. Maximum mass-loss rates, minimum clumping factors and alternative boundaries of the defined wind regions (adapted-regions approach, see Section 3) providing the best-fit models for some objects in the sample. Column descriptions and values in parentheses or separated by forward slashes as in Table 6. r in indicates values different from the nominal, adopted r in = 1.05 R * , and r far is always set to 50 R * . reproduced the data with a goodness of fit within 10% of the corresponding minimum χ 2 . The resulting β-intervals are presented in Table 8 as with β 0 being the used β in our simulations. We found that the limits of the intervals are not symmetric with respect to the used β (except in four objects). Moreover, our best-fit solutions are less sensitive to a decrease than to an increase of β, regardless of luminosity class, spectral type or "nature" of the source (Figure 7). Thus, our best-fit solutions are valid for average changes of β of |δβ − /β 0 | ≈ 24% and |δβ + /β 0 | ≈ 16%.
In the second experiment we computed new fitting models. We  ferent values of β, until a best-fit solution was found. This test estimates the errors of the derived clumping factors as: with δβ = β − β 0 . In Table 9 we present the results of these simulations for f cl (β ), with β = (β 1 , β 2 ), where β 1 β − , β 2 β + , and β 0 is the used β in our simulations. Note that for some sources the values of the alternative β 1 and β 2 are larger (or lower) than the typical limits derived from spectral line fitting (though always consistent with the theoretical limit β > 0.5). Nonetheless, our approach allows us to check the validity of the analysis as a function of β. From this experiment we found that, as expected, i) large clumping factors, f in,mid cl , correspond to low values of β and vice versa, i.e. δ f cl ∝ −δβ ; ii) changes in β do not uniformly (or symmetrically) affect all the defined clumping factors, and also depend on luminosity class. Thus, for OB Supergiants average values |δβ | ∼15% -30% lead to average uncertainties of about 15% -60% in f in cl and ∼ 15% in f mid cl . On the other hand, for the sources with weaker winds in our sample, i.e OB Giants and Dwarfs, average values of |δβ | ∼15% -40% translate into average uncertainties of about 10% -15% for f in cl , and up to 8% for f mid cl ; and iii) for a few sources, an increase of β cannot be compensated by a decrease of clumping factors, since the clumping factors for the used β reach unity, so that instead a decrease iṅ M max is required (last column in Table 9). Note that the objects with weaker winds and upper limits iṅ M max also require a clumping increase in the outermost wind region, f far cl ( f cl max ) > 1, in order to be able to consistently reproduce upper radio-flux limits. Our results agree with the average uncertainties estimated for f in cl and f mid cl by P06 (∼ 30% and ∼ 20%, respectively) by means of Hα profiles and IR/radio fluxes analyses in parallel.
Considering these experiments, the average error due to β uncertainties for the estimated values of f mid cl is ∼ 30%, even in the extreme case that the actual values of the wind acceleration parameter β would be significantly different from those adopted in this work.

Discussion
The previous section showed that the clumping structure can be well described relative to the outermost wind in Region 5 (radio domain). Assuming for now minimum values f far cl = f min cl = 1, Article number, page 15 of 35 A&A proofs: manuscript no. AA202040116 A key finding of our analysis is that clumping at r 2 R * presents similar radial stratification regardless of the strength of the wind, always fullfilling the condition f mid cl f out cl f far cl = f min cl (see Tables 6 and 7).
The specific behaviour from the intermediate to outermost wind seems to depend on luminosity class and spectral type, though, on average, the decrease in f cl from the intermediate to outermost regions is steeper for OSGs (a factor 4) than for BSGs (a factor 2). Moreover, the clumping properties of the few OB Dwarfs and Giants in our sample show a smoother, more homogeneous behaviour than OB Supergiants, with similar clumping degrees at the intermediate and outer wind regions for most sources.
This varying clumping degree for different luminosity classes and spectral types is supported not only by larger values of relative clumping factors, but also by the fact that the clumping stratification seems to be confined to a narrower region in OSGs and some BSGs, as revealed by the use of our adaptedregions approach (see Sec. 3). For the majority of the sample the boundaries used in the fixed-regions approach seem to describe the radial clumping stratification rather well. However, for some  Fig. 7. Absolute ratio between the semi-widths (δβ + , δβ − ) of the validity interval of β as a function of spectral type for the best-fit solutions derived in the fixed-regions approach in this work (see Table 6). Symbols and colours as in Fig. 4. The circled symbol indicates that the ratio was divided by 3 for display purposes.
OSGs in our sample the best-fit solution was found after adapting mainly r mid , and in a few cases also r in and/or r out . This is also the case for half of the BSGs (one B0 I and all eBHGs) and for one of the O Giants (see Table 7). Changes in the extension of the wind regions are balanced by changes in the clumping factors, thus the overall clumping stratification is rather similar but more precisely described. Moreover, in all simulations we observed that the quality of the fits barely changed when varying r out = 15 R * to r out = 10 R * , with the exception of OSGs, for which r out = 12 -15 R * always provided the best-fit. Keeping in mind that the average minimum clumping degree at this radius is similar to the clumping degree in the outermost wind region, r > 50R * , this implies that the clumping-degree of these sources drops quickly after reaching its maximum value at r ≈ 1.1 -6 R * .
There are a few exceptions to this general behaviour, involving non-standard sources such as a binary system (HD 37043) and variable thermal emitters (P Cyg, CygOB2#12). The other binary source in our sample, HD 149404, seems to follow the general trend, however the characteristic radio flux variability associated to binarity could not be tested for, since only one flux measurement at 3.6 cm is available. Therefore, this object could eventually behave as the other non-standard objects in the sample. This discrepancy in the clumping behaviour between standard and non-standard objects suggests that those sources showing an increasing clumping degree from the intermediate to outermost regions do not reflect intrinsic wind-clumping differences. Instead, this might rather be explained by different physical conditions in the outer wind when compared to standard sources (colliding winds, binarity, magnetic fields, ionisation changes, etc), which can significantly modify the flux emission in the mm and the radio regimes. If such a correlation was confirmed by further studies, an analysis of radial clumping stratifications in the intermediate to outermost wind regions would help to discriminate between objects of different natures (e.g., standard or peculiar). Table 9. Clumping factors and maximum mass-loss rates (in units of 10 −6 M yr −1 ) derived for alternative values (β 1 and β 2 ) of the used β in our simulations in the fixed-regions approach (Table 6; labeled here as 'used'). Entries in the column labeled 'Ṁ max (β 2 )' indicate sources for which only a value lower than the estimated maximum mass-loss rate for the used β is consistent with all observations (see text); in these cases, the outermost wind must be clumped as indicated by the column labeled ' f far cl (β 2 )'. Usually f out cl did not require further modification except for the asterisk-marked sources (see table foot text). Superscript-letters and values in parentheses or separated by forward slashes as in Table 6.

Comparison with empirical and theoretical studies
Our derived clumping structures at r 2 R * agree well with previous studies. For instance, P06 found a similar dependency of the clumping degree with luminosity class when comparing the intermediate and outermost wind regions. In addition, for those stars in common with their sample, our flux estimates from PACS observations at 70, 100, and 160 µm allowed us to here better constrain clumping in the intermediate wind region, lowering P06's previous upper-limit estimates for this region. Moreover, we find qualitatively similar results for the clumping structure to those reported by Najarro et al. (2011) andClark et al. (2012), using a different clumping parametrisation . Overall, these authors find a sharp decrease in the clumping degree for OB Supergiants beyond r 1.5 -2 R * , while OB Giants and Dwarf stars present a roughly constantpossibly unclumped-value throughout the wind. If our clumping factors in the outermost wind region (radio regime) are normalised to a similar value as those in Najarro et al. (2011) and Clark et al. (2012) the wind structure of the common targets can be compared quantitatively. For HD 66811 (ζ Pup) and HD 36861, the agreement is excellent, despite the different luminosity class and clumping properties derived for these two stars. For HD 152236 (ζ 1 Sco), and CygOB2#12, our absolute values of clumping factors are very similar in the inner and outermost wind regions, and a factor 3 lower at r = 2 R * . Finally, for the other two OB Supergiants in common, HD 30614 and HD 37128 (α Cam and Ori, respectively), the derived clumping values differ considerably by more than a factor 6 at r = 2 R * , but they converge quickly for larger radii. Nevertheless, this apparent discrepancy might be due to the lack of reliable FIR flux observations in those studies, which resulted in poorly constrained degrees of clumping in this wind region. Interestingly, in all sources in common with Najarro et al. (2011) andClark et al. (2012), the best-fit model is found in the adapted-regions approach. This would imply that the physical conditions in the wind of these stars change considerably over small distance increments, suggesting a more confined or structured wind. On the other hand, Sundqvist et al. (2011) analysed the clumping structure in λ Cep by means of radiation-hydrodynamic (RH) simulations and empirical, stochastic models (including effects of optically thick clumping). Here the agreement with our work is again very good once we scale our clumping factors accordingly, although we find here that f cl reaches its maximum at a slightly larger radius (r 2 R * ) when compared to their simulations (r ≈ 1.2 -1.5 R * ).
The above results further allow us to compare to theoretical predictions from LDI simulations. Overall, our results agree quite well with 1D simulations for OSGs by Sundqvist y Owocki (2013) and Driessen et al. (2019) that account for limb-darkening and/or photospheric perturbations (see also dis-Article number, page 17 of 35 A&A proofs: manuscript no. AA202040116 cussion below). The clumping factor in these models typical peaks at ∼ 1.5 -2 R * and decreases beyond, in general agreement with our empirical findings here. Moreover, our empirical study also tentatively agrees with recent simulations by Driessen et al. (2019), which show that the winds of BSGs should be overall less clumped than OSGs. However, all these LDI simulations only reach r = 10 R * and would need to be extended to higher radii in order to confirm this.

The innermost wind region
There are also a few interesting aspects of the clumping structure in the inner-to-intermediate wind region transition. We find that for several OB Supergiants in our sample (13/20), the clumping degree seems to increase from the inner to the intermediate wind region ( f in cl f mid cl ). This agrees with P06, who found the same trend for all their sample but one star, HD 66811 (ζ Pup). Such a trend implies that the maximum f cl would occur at r ∼ 2 -6 R * , depending on the source. However, this is not really a general trend, since for a non-negligible fraction of our sample (7/20) we find the opposite behaviour ( f in cl > f mid cl ). In other words, the exceptional behaviour of ζ Pup found by P06 is shared among several stars in our sample (and also those studied by Najarro et al. 2008, Clark et al. 2012and Sundqvist et al. 2011. This would imply that in these stars the clumping onset occurs very close to the base of the wind (r 1.2 R * ), and that the maximum clumping is achieved slightly below the intermediate wind region (r ∼ 1.5 -2 R * ).
This discrepancy could, however, be influenced by the βf in cl degeneracy problem (see Sec. 4.4). Lowering/increasing the values of β requires a corresponding increase/decrease in clumping factors, which is larger in the inner region than in the intermediate one (see Table 9). Therefore, for stars with f in cl < f mid cl to become reversed (i.e., f in cl > f mid cl ), a very large decrease of β would be needed, and vice versa. It is reasonable to assume some uncertainties in the used β values, however it seems unlikely that the estimated values of β are off by such a large amount in so many stars. That is, the existence of two different trends in our sample may point to intrinsic differences in the efficiency of the different mechanisms governing the onset of clumping.
As also discussed above, these observed trends may be compared to the different theoretical LDI simulations by Sundqvist y Owocki (2013) and Driessen et al. (2019). Namely, depending on the different initial conditions of the simulations, and whether the modeled star is an OSG or BSG, the onset, the peak, and the overall predictions for the radial stratification of clumping may differ. For example, for high-luminosity objects it is possible that the radiative acceleration will exceed gravity already in deep sub-surface layers (e.g., at the so-called 'iron opacity-bump'), which might then trigger a turbulent atmosphere that may also affect quantitative predictions for clumping factors (Cantiello et al. 2009;Jiang et al. 2015). Effects on clumping factors from a perturbed photosphere were investigated by Sundqvist y Owocki (2013) (see also Feldmeier et al. 1997), who indeed found that this can affect predictions both for the inner-wind clumping and for the clumping stratification as a function of radius. As mentioned above, the empirical results for the outer wind regions found here seem overall consistent with these simulations. We note, however, that since these models had an outermost wind radius ∼ 10 R * it remains to be investigated how far out in the wind effects from a potentially turbulent photosphere could persist.
As such, the observational constraints obtained here will provide a very sound base for future theoretical models targeting more detailed investigations of this radial behaviour. Observationally, to understand and to properly describe the physical mechanisms that determine the clumping degree in the innermost wind region, further investigations using multi-wavelength analyses, continuum and line fits, including both ρ-dependant (resonance lines) and ρ 2 -dependant (Hα and NIR lines + V to radio continuum emission) clumping diagnostics are required.

Upper-limit mass-loss rates and comparison to theory across the bi-stability jump
In the following, we compare our empirical estimates of maximum mass-loss rates (Ṁ max ) to the theoretical predictions by V00 and V01 (Ṁ Vink th ), since they are the mass-loss rate recipes most commonly used in key applications such as stellar evolution 10 These authors provide simple recipes to estimate massloss rates for various ranges of effective temperatures, depending on the so-called first and second 'bi-stability jumps'. These jumps occur in the models when iron recombines first from Fe iv to Fe iii (first jump) and then to Fe ii (second jump); since the lower iron ionisation states have more efficient driving lines, crossing these jumps from the hot sides result in significant increases in predicted mass-loss rates.
The recipe suggested in V00 (and V01), switches from 'Fe iv' to 'Fe iii' in the range T jump1 eff = 27.5 -22.5 kK and then further to the 'Fe ii' branch around T jump2 eff ∼ 12.5 -18.5 kK, depending on the mean wind density (see eqn. 6 in V00, and discussion in their section 5.3). This means that the central jump temperature that defined the first and the second bi-stability jump regions are T jump1 eff ∼ 25 kK and T jump2 eff ∼ 15 kK, respectively. However, studies by Lamers et al. (1995), Crowther et al. (2006), Markova y Puls (2008), Petrov et al. (2014) and Vink (2018) indicate that the first jump, if any, is around 22 -20 kK, whereas the second jump, according to Petrov et al. (2016), might be below 9 kK. Since our OB Supergiant sample covers effective temperatures from 39 000 K to 13 700 K (O4 -B3/B4) we can empirically investigate the mass-loss behaviour across the predicted bi-stability jumps using our data.
Let us note already here that, very recently, Björklund et al. (2020) presented new theoreticalṀ predictions for O-stars. These are based on steady-state hydrodynamics and NLTE line radiative transfer in the co-moving frame (instead of a prescribed β velocity law and NLTE Sobolev-Monte-Carlo line-transfer computations, as in V00 and V01), and consistently find loweṙ M values than V00 and V01. However, since these models have not yet been extended to BSGs, we opt here to compare our findings only to those by V00 and V01 (Ṁ Vink th ), rather than using some proxy such as the wind-momentum luminosity relation.
We compare our empirical resultsṀ max with theoretical pre-dictionsṀ Vink th as they are used in the often-cited grids of evolutionary models in codes such as Geneva (e.g. Ekström et al. 2012, Yusof et al. 2013, Bonn (e.g. Brott et al. 2011, Köhler et al. 2015 and MESA (Modules for Experiments in Stellar Astrophysics; Paxton et al. 2011).
First, we specifically use the definitions of T jump1 eff and T jump2 eff from equations 15 and 6 from V01 and V00, respectively, and the mass loss recipe corresponding to the location of the source in the temperature/bi-stability jump space is applied (see V00 and V01). This implementation corresponds to the bi-stability jump regions and the mass-loss rates recipe as outlined above, which is indeed used in the Geneva code (hereafter Geneva approach). For a second test, we adopt a similar approach to that used in the MESA and Bonn codes (hereafter MESA approach): T jump1 eff is obtained as a function of the density of the wind via metallicity (eqn. 14 and 15 by V01), T jump2 eff is set to 10 kK, and the mass loss prescriptions for the 'hot' and 'cold' side of the first jump are used (eqns. 24 and 25 by V01). For T eff between 27.5 and 22.5 kK, we use the same interpolation method than in MESA, and developed by Brott et al. (2011). Finally, to investi-gateṀ Vink th as a function of the temperatures of the first and second jumps reported by several authors (see above), we carry out a third test, following the MESA approach but fixing also T jump1 eff at 22 kK (hereafter Fixed-jumps approach). It is worth to mention that, whereas the Bonn and MESA codes define T jump2 eff below 12.5 -10 kK -switching to different mass loss predictions 11 -and T jump1 eff depending on the wind density of the source via metallicity, the Geneva code defines them via Γ e (eqn. 5 by V00). Since our sample are galactic OB stars, this means thatṀ Vink th estimates will be barely affected by which definition of T jump1 eff we use in this analysis. 12 Figure 8 shows the empirical to theoretical mass-loss ratio (Ṁ max /Ṁ Vink th ) as a function of T eff for the OB Supergiants subsample, for the three different implementations of V00 & V01 described above. The non-thermal and the variable thermal emitters (HD 149404, P Cyg, CygOB2#12) were excluded of this analysis.
There is a clear trend in the sample, showing higher levels of discrepancy betweenṀ max andṀ Vink th for lower effective temperatures, regardless the used definition for the temperature of the bi-stability jumps. For the hotter OSGs,Ṁ max agree within errors (Section 4.3) with theṀ Vink th recipes. But as T eff is lowered reaching the BSG regime, the difference between empirical and theoretical mass-loss rates increases significantly. Indeed, for the coolest BSGsṀ max is up to almost 2 orders of magnitude lower than the prediction when using the Geneva approach, and close to 1.5 orders of magnitude when using the MESA approach.
At first glance, these very big discrepancies for luminous BSGs may seem surprising. However, theṀ max values derived in this work are (on average) actually not that different than many others obtained by various studies present in the literature. For example, Haucke et al. (2018) derivedṀ max for a sample of BSGs by means of Hα fitting using (unclumped) fastwind models. For the 4 stars that are present in both samples, having T eff = 17 500 -25 000 K, we obtain an av-erageṀ Haucke max /Ṁ thiswork max ≈ 0.8. Moreover, inspection of their table A.1, where they also list rates derived by a few other studies, seem to not reveal any major systematic discrepancies. Similarly, Crowther et al. (2006) also used Hα to derivė M max (using unclumped cmfgen models); here, for the 9 overlapping stars, all within T eff = 15 500 -29 000 K, we find an averagė M Crowther max /Ṁ thiswork max ≈ 1.0. As such, although the scatter is large it seems unlikely that the big discrepancies between observations 11 The switch to a different wind density scenario and mass loss prescription in these codes following different criteria as T eff reaching the defined T jump2 eff . 12 For our sample, T jump1 eff ≈ 25 kK (26.1 kK-24.1 kK) when eqn. 5 and eqn. 4 (V00) are used, and T jump1 eff ≈ 25.9 kK when eqn. 14 and eqn. 15 (V01) are applied instead.    Fig. 8. From Top to Bottom, empirical to theoretical mass-loss rates ratio, in logarithmic scale, as a function of effective temperature for the OB Supergiants subsample. Empirical mass-loss rates correspond to thė M max derived in this work. Theoretical mass-loss rates,Ṁ Vink th , correspond to the mass-loss rates computed via recipes from V00 & V01 for different definitions of the temperatures of the jumps (see Sec. 5.2); top: Geneva approach, middle: MESA approach, and bottom: Fixed-jumps approach (see Sec. 5.2). Different colours indicate at which side of the bi-stability jumps the sources are located. Arrows and symbols as in Figure 4. and theory found for BSG mass-loss rates should be caused by any systematic effect specific for our study here (even though considerable uncertainties exist for various individual objects, see Table 6 and discussions in Appendix A).
As just one specific example, for HD 198478 with T eff = 17 500 K we derive an interval (see Table 7, and discussion in Appendix A)Ṁ max = 0.14 -0.38 × 10 −6 M yr −1 . This  Fig. 9. Top: Empirical-maximum (Ṁ max ) and theoretical (Ṁ Vink th ) wind performance numbers, η =Ṁ υ ∞ /(L * /c), as a function of T eff , for our OB Supergiants subsample, in the Geneva (magenta) and MESA (orange) approaches. Dashed and dotted-dashed lines correspond to theoretical predictions for a source with log L/L = 5.75 and M * = 45 M for solar metallicity respectively based on Geneva and MESA implementations of V00 and V01. The shadowed regions represent the first (27.5 -22.5 kK) and second (18.5 -12.5 kK) bi-stability jump zones as defined by V00. Bottom: Same as top panel, but showing theoretical wind performance numbers in the Geneva (magenta) and Fixed-jumps (light green) approaches, and a dotted-dashed line marking the theoretical model based on the Fixed-jumps implementation of V00 and V01. The shadowed region represents the first bi-stability jump zone (24.5 -19.5 kK) for a central temperature of the jump of 22 kK. In both panels the early B-type Hypergiants in our sample are indicated with squares, whereas arrows indicate upper limits ofṀ max . Finally, arrows, dotted lines, and symbols as in Figure 8. agrees well with the intervalṀ max = 0.12 -0.41 × 10 −6 M yr −1 obtained by Markova y Puls (2008) from their unclumped Hα analysis. Direct application of Geneva approach discussed above to HD 198478 places this star at the cool side of the second bistability jump (see dotted-lined source in Figures 8 and 9); as such, even when using the empirical upper-limit υ ∞ = 470 km s −1 from Table 6, the recipe predicts a very largeṀ Vink th = 6.38 × 10 −6 M yr −1 , causing the large discrepancy visible in Figure 8 (see top panel). Even if the star is placed in the 'Fe iii' region, i.e. when MESA and Fixed-jumps approaches are used, the recipe would still return a rate,Ṁ Vink th = 5.025 × 10 −6 M yr −1 , that is more than an order of magnitude higher than the empirical upper limit (middle and bottom panel in Figure 8). This shows that the reason we find large discrepancies betweenṀ max andṀ Vink th for BSGs, but not for OSGs, is because of the large predictedṀ Vink th increase when crossing the bi-stability jump regions.
In fact, when inspecting the so-called wind-performance number η =Ṁ υ ∞ /(L * /c) for our sample (Figure 9) we find a gradually decreasing trend -albeit again with significant scatterwith decreasing T eff , and do not see any evidence for sudden increases, or a secondary local maximum (e.g. Benaglia et al. 2007), at any of the predicted bi-stability limits. This is quite similar to the findings by Crowther et al. (2006), Markova y Puls (2008) and Haucke et al. (2018), who also did not find any empirical evidence for a sudden mass-loss increase with decreasing T eff . Note that, although the difference between theoretical predictions and empirical values is reduced when setting the central temperatures for the first and second bi-stability jumps at ∼22 kK and ∼10 kK respectively, the observed discrepancies for B Supergiants are still considerable.
These large discrepancies may thus point toward a problem with current theoretical mass-loss rate predictions around the first bi-stability jump, regardless of its exact location, rather than toward issues related to empirical mass-loss determinations for BSGs derived by different authors and methods. This suggests that the basic mechanism responsible for increasing the massloss rate at the bi-stability jump (essentially recombination of iron) needs to be revisited in future work, preferably by locally consistent hydrodynamic calculations that do not assume a predescribed β velocity law. This is important also in view of the quite big role this bi-stabilityṀ increase plays in general stellar evolution modelling (see e.g., Vink et al. 2010;Keszthelyi et al. 2017).

Impact of clumping in scaling the derived upper-limit mass-loss rates
As mentioned above, the empirical values of the maximum OSG mass-loss rates derived from radio fluxes agree (within uncertainties) with the theoretical predictions by V00 and V01 (see also P06). However, 'real' agreement would only be true if the outermost wind of OSGs were unclumped. We note again that ourṀ max estimates were obtained normalising the clumping stratification to the outermost wind region by setting f far cl = f min cl = 1. Therefore, if instead this region is affected by clumping, downward corrections become necessary. Hydrodynamical wind models extending all the way to the radio-emitting regions suggest typical factors of about f far cl ≈ 4 -9 (Runacres y Owocki 2002. These models, however, were calculated assuming an unperturbed lower boundary at the stellar surface, and so it remains to be seen whether effects from, e.g., a turbulent photosphere could persist all the way out to the radio region (see also discussion in Sect. 5.1.3). Assuming for now the outer-wind values from Runacres y Owocki (2002, thė M max estimates of this paper would need to be scaled down by a factor ≈ 2 -3. We re-emphasise that these typical f far cl ≈ 4 -9 only refers to the outermost wind. Clumping factors predicted for the inner wind regions are typically significantly higher; for example, in the O-star simulations by Driessen et al. (2019) f cl peaks around r/R * ∼ 1.5 − 2 and the average value for the inner wind (r/R * = 1 − 2) is f cl ≈ 15 − 20. Thus, O-star mass-loss rates derived from unclumped models of inner-wind diagnostics (e.g. using optical recombination lines like Hα) should be scaled down more than corresponding radio rates. This is consistent with the empirical finding here that clumping factors on average are higher in the inner wind regions than in the outer ones. Typical reductions of ≈ 2 − 3 would agree well with the new theoretical mass-loss rate predictions for O-stars by Björklund et al. (2020) (see also Sander et al. 2017 andSundqvist et al. 2019), as well as with various multi-wavelength empirical studies such as Najarro et al. (2011) (see figs. 12 and 13 in Björklund et al. 2020). As demonstrated by Figure 8, however, the situation is very different for BSGs, whereṀ Vink th can be overestimated by as much as 2 orders of magnitude. Of course, if the outermost wind of BSGs should be significantly clumped, this discrepancy would become even larger. For example, assuming the same clumping degree in the outermost wind as for OSGs (but see Driessen et al. 2019), current theoreticalṀ predictions by V00 & V01 (the standard rates used in most evolutionary models), would be overestimated by factors 6 -200 in case of the Geneva implementation.

Summary and conclusions
We have constrained clumping properties of the intermediate stellar winds for a sample of 25 OB Supergiants, Giants and Dwarfs stars, by using Herschel-PACS flux measurements at 70, 100 and 160 µm. For the analysis, we followed the approach developed, and further tested, by P06, assuming optically thin clumping. Additional available continuum observations in the literature from optical to radio wavelengths further allowed us to perform a consistent analysis and to derive robust estimates for the (minimum) wind-clumping structures at r 2 R * , as well as maximum mass-loss rate estimates ( f far cl = f min cl = 1). Accurate radio emission measurements allowed us to derive preciseṀ max values and to infer robust clumping structures for 17 of the 25 stars in our sample. In addition, mm flux observations of 7 of these stars led to better constraints of the outer wind region (15 r/R * 50), providing a precise determination of the radial stratification of the clumping at r 2 R * . Our major findings can be summarised as follows: -Our findings agree well with the empirical clumping properties at r 2 R * derived by Najarro et al. (2011) and Clark et al. (2012) following a different parametrisation. In addition, our results overall support the hydrodynamical OSG models by Sundqvist y Owocki (2013), where clumping starts to decrease at r ≈ 2 -6 R * , and, tentatively, the recent 1D LDI-simulations of OSG and BSG winds by Driessen et al. (2019), which predict lower amounts of clumping in BSGs.
-We found that for 8 OB Supergiants in our sample f in cl > f mid cl . This significantly extends the findings of P06 in just one star of their sample (ζ Pup), and is in agreement with the empirical clumping properties by Najarro et al. (2011) and the theoretical predictions for OSGs by Sundqvist et al. (2011) andSundqvist y Owocki (2013). This suggests that such a behaviour, rather than being an exception, could imply the existence of two trends characterised by different physical conditions at the base of the wind.
-For OSGs the derived upper-limit mass-loss rates,Ṁ max , agree with the theoretical predictions by V00 and V01 for unclumped winds, whereas BSGs show a discrepancy which severely increases with decreasing effective temperature: the estimatedṀ max start to differ from theoretical recipes in the predicted first bi-stability transition zone, and up to 1.5-2 per SCUBA value, which does not help to properly constrain the clumping degree for Region 4. However, both 160 µm and SCUBA fluxes have large errors 13 , preventing us from drawing any further conclusions. In the following, the described clumping structure refers to the fixed-regions approach. The best-fit solution is achieved for the minimum possible clumping degree for the inner and outermost wind regions, and a larger clumping for the intermediate one, f in cl = f out cl = f far cl = 1 and f mid cl = 7 (solid line in Figures A.1 and A.4). In Region 4, although a minimum clumping factor provides the best-fit solution, the clumping degree could range from 1 to 8, with a maximum value also ranging from 5 to 12. Thus, with f out cl = 1 or f out cl = 8, the lower (solid line in Figures A.1 and A.4) or the upper (magenta-dotted line in Figures A.1 and A.4) SCUBA value are respectively matched. f out max = 5 -12 is still consistent with the 160 µm and with both SCUBA fluxes, respectively. Only extra flux observations at sub-mm and mm spectral ranges can help derive definite estimates for wind region 4. For the new distance estimate from GAIA, d = 0.62 kpc, the corresponding scaled maximum mass-loss rate isṀ max = 1.3 × 10 −6 M yr −1 , and the same clumping properties described above (Table 6) provide the best-fit solutions ( Figure A.4).

HD 152408
The flux observations at 2 and 6 cm are well determined, thusṀ max is definite. In the corresponding panel in Figure A.4, we present our best-solution fitting all the data (magenta-dotted line) and the best possible one in the fixed-regions approach (solid line), respectively. The maximum mass-loss rate is 9.5 × 10 −6 M yr −1 one of the largest values in the sample. For the fixed-regions approach the best solution is obtained 13 The spatial resolution of the PACS photometer at 160 µm is lower than at 70 and 100 µm, and the PSF-FWHM is large enough to enclose extended emission from its surrounding. Indeed, a close inspection of PACS images reveals a strong background emission at 160 µm. by setting f in cl = f out cl = f far cl = 1, and f mid cl = 3.25 (Table 6). However, the best-fit solution for HD 152408 is achieved for the adopted-regions approach by slightly increasing to f mid cl = 3.8, provided that the inner wind region is extended to r mid = 3.5 R * ( Table 7). The SCUBA measurement (1.3 mm) further allows us to constrain the maximum possible value of the clumping degree in Region 4 to f out max = 3.
HD 151804 (V973 Sco) Discussed in detail in the main text.
HD 30614 (α Cam) This source was analysed by P06, and we therefore show P06's best solution together with our best-fits models. It can be observed in the fit-diagram that P06's solution underestimates fluxes at MIR and FIR wavelengths, and that although theiṙ M max is still consistent with the additional measured fluxes at 6 cm 14 , the additional radio flux at 21 cm points to a slightly largerṀ max . We found that a value ofṀ max = 1.75× 10 −6 M yr −1 with f in cl = 6, f mid cl = 3.3 and f out cl = f far cl = 1, provide the best possible solution in the fixed-regions approach (solid line; Table 6). However, the best-fit solution (magenta-dotted line) for the estimatedṀ max is achieved by extending to r mid = 4 R * , so that f in cl = 6 and f mid cl = f out cl = f far cl = 1 (Table 7). In all simulations the maximum value of f out cl is constrained by f out max =2.

Appendix A.2: B Supergiants
HD 37128 ( Ori) There are well determined available radio fluxes for this source, thereforeṀ max is definite. Foṙ M max = 1.25 × 10 −6 M yr −1 , we present the two fit-solutions consistent with all the data, corresponding to the fixed-(solid line) and the adapted-regions (magenta-dotted line) approach, respectively ( Figure A.5). For the fixed-regions approach, the best possible solution corresponds to a 1 ( Table 6). The best-fit solution, however, is obtained in the adapted-regions approach by increasing a factor 2 the clumping degree in the intermediate wind region, f in cl = 5, f mid cl = 8 and f out cl = f far cl = 1, provided that r out = 8 R * (Table 7).

HD 38771 (κ Ori)
Only one radio observation, an upper limit at 6 cm, is available for this source. Thus, only an upper limit forṀ max is provided. The panel corresponding to this star in Figure A.5 shows the fit-solutions for both the fixed-(solid line) and the adapted-regions (magenta-dotted line) approach. We es-timatedṀ max 0.7 × 10 −6 M yr −1 . For this upper limit, the best possible solution in the fixed-regions approach corresponds to a significantly clumped inner wind, f in cl = 8, and f mid cl = f out cl = f far cl = 1 ( Table 6). The best-fit solution is achieved in the adapted-regions approach, with a considerable large clumping degree in the inner wind region, f in cl = 20, and again f mid cl = f out cl = f far cl = 1, provided that r mid = 1.55 R * ( Table 7). Despite of the lack of flux observations from 100 µm to radio wavelengths, f out max = 8 is still consistent with all data. Only mm and well-determined radio flux measurements can provide definiteṀ max and the clumping structure of the outer wind for this source.

HD 154090 (κ Sco)
Only one radio measurement is available at 3.5 cm, but it is well determined and thusṀ max is definite. Figure A.5 displays fit-solutions for both the fixed-(solid line) and the adapted-regions approach (magenta-dotted line). Foṙ M max = 1. × 10 −6 M yr −1 , a moderate clumping degree in the inner and the intermediate wind regions, f in cl = 2 and f mid cl = 4.5, and f out cl = f far cl = 1, provide the best possible fit solution in the fixed-regions approach ( Table 6). The best-fit solution for this source is achieved in the adapted-regions approach, with a moderate, and almost constant, clumped wind restricted to r 15 R * , f in cl = 4, f mid cl = 3.5 and f out cl = f far cl = 1, provided that r in = 1.4 R * instead of r in = 1.05 R * (Table 7). In the absence of sub-and mm observations, f out max = 8 is still consistent with the 70 and 100 µm and 3.5 cm fluxes. Only additional flux observations at mm can better constrain wind region 4.

HD 193237 (P Cyg)
Discussed in detail in the main text.
HD 24398 (ζ Per) Discussed in detail in the main text.

HD 169454
This source is one of the four eBHGs in our sample. The available radio fluxes at 2 and 6 cm are well determined, and thusṀ max is definite. In Figure A.5 we present fit-solutions for the fixed-(solid line) and the adapted-regions (magenta-dotted line) approach. ForṀ max = 10.4 × 10 −6 M yr −1 , the best possible solution in the fixed-regions approach is provided for a weak intermediate clumped wind, f mid cl = 2.5, and an unclumped inner and outer wind, f in cl = f out cl = f far cl = 1 (Table 6). However, a slightly better solution is obtained in the adapted-regions approach, by increasing the clumping degree only in the intermediate wind region, f mid cl = 4.5, provided it begins at r mid = 4 R * (Table 7). In all simulations, the single mm flux estimate (at 1.3 mm) constrains the maximum clumping degree in Region 4, still consistent with all data, f out max = 2.5 HD 152236 (ζ 1 Sco) This star is another of the four eBHGs in our sample. For the well determined radio measurements, we obtaiṅ M max = 6.2 × 10 −6 M yr −1 . Figure A.6 displays the derived solutions for the fixed-(solid line) and the adapted-regions approach (magenta-dotted line). For the fixed-regions approach, we obtain a weak clumping degree in the intermediate wind region, f mid cl = 2.8, and f in cl = f out cl = f far cl = 1 (Table 6). However, a slightly better solution is obtained in the adapted-regions approach with a marginally larger clumping factor in the intermediate wind region, provided that it is restricted by r out = 10 R * (Table 7).
HD 41117 (χ 2 Ori) All available radio measurements are well determined, and thereforeṀ max is definite. ForṀ max = 1.8 × 10 −6 M yr −1 , the best-fit solution is achieved in the fixed-regions approach for f in cl = 7.3, f mid cl = 1.35 and f out cl = f far cl = 1 (solid line in the corresponding panel in the Figure A.6). The lack of observations in the sub-and mm regime lead to f out max = 6 consistent with FIR and radio observations at 2 cm.

HD 194279
Our measured fluxes at 70 and 100 µm are lower than theoretical predictions for an homogeneous, unclumped wind, whereas the available, well-determined radio fluxes at 3.5 and 6 cm are consistent with those. We derived two solutions intending to fit consistently the FIR and/or radio fluxes (corresponding fit-diagram in Figure A.6). The described clumping structures below correspond to the fixed-regions approach, which provides the best-fit solution for this source.
A first solution (solid line), matching radio observations, is provided forṀ max = 2.12 × 10 −6 M yr −1 and for a constant f in cl = f mid cl = 1 (first entry Table 6). However, it can be observed that although this solution reproduces the radio fluxes, it only hardly matches the flux measurement at 70 µm, and overestimates the one at 100 µm.
A second, and better solution (magenta-dotted line), is achieved for a much lowerṀ max = 0.505 × 10 −6 M yr −1 . For this value, the minimum clumping degree at r 15 R * and a strong one at r 15 R * , perfectly match all observations, i.e. f in cl = f mid cl = 1 and f out cl = f far cl = 20 (second entry in Table 6). The second solution fits all flux observations throughout the wind and provides the best-fit for this source. However, if this is the case, HD 194279 is the only single star and radio thermal emitter in our sample with such a high clumping degree in the outer wind regions (r 15 R * ). On the other hand, if HD 194279 turned out to be a non-thermal emitter, then this clumping structure would follow the same trend found for the other non-thermal emitters in our sample as HD 37043 (SB) or CygOB2#12. Due to the scarcity of flux observations of HD 194279 at large wavelengths, the high uncertainties derived in the estimated FIR fluxes and the fact that the two radio measurements are consistent with thermal emission, we consider the constant clumping degree solution to be more likely (first entry in Table 6). Only re-observing at FIR and radio wavelengths can help discriminating between these scenarios.

HD 198478
Due to its stellar pulsation activity this star has variable spectral profiles, so that the stellar and wind parameters obtained via optical spectroscopic analysis show substantial scatter (Crowther et al. 2006, Markova y Puls 2008, Kraus et al. 2015. This is not due to the different approaches to derive stellar parameters, but to the phase of activity of the source when spectra where taken. In order to study how different input values affect the esti-matedṀ max and clumping structure, we use the range of input parameters as provided by Markova y . In Table A.1 we present our fit solutions in the fixed-regions approach for each combination of υ ∞ , β and Y He by these authors. The two diagrams in Figure A.2 display the corresponding solutions present in Table 6. For clarity, in Table 6 we indicate the full range for derivedṀ max and clumping factors, and Figure A.6 presents fitsolutions corresponding to the lowest (magenta-dotted line) and the largest (solid line)Ṁ max .
Only two radio flux measurements at 3.5 and 6 cm are available (Scuderi et al. 1998), and both are upper limits. Therefore, the resultingṀ max is an upper limit. The minimum χ 2 obtained for the different solutions in Table A.1 are only marginally different (solid and magenta-dashed lines in Figures A.2 and A.6). Note that although the absolute values of clumping factors change, the overall clumping properties throughout the entire wind are conserved, and the largest difference in f in cl and f mid cl is a factor of 1.33. RegardingṀ max , the largest difference, a factor 2.7, arises from simultaneously changing υ ∞ and Y He (first and last entry in Table A.1). However, for a given velocity field (υ ∞ and β) and a factor of 2 different helium-fraction,Ṁ max and clumping factors vary only by ∼ 15% and less than 12%, respectively. Moreover,   Figure A.2). T eff is in kK; log g is the gravity without centrifugal correction; YHe denotes the helium fraction by number used in the simulations; and υ ∞ , R * andṀ max are in units of km s −1 , R and 10 −6 M yr −1 , respectively. f out max gives the maximum possible clumping factor in Region 4 still consistent with the data. log Q = log (Ṁ max /R * 1.5 ) is in units of M yr −1 and R , and, for clarity, the last column indicates the label for the solution displayed in Figures A.2 and A.6.

Fig. A.2.
Observed and best-fit fluxes vs. wavelength for HD 198478 depending on the input parameters used in the simulations (see Table A.1 and text in Appendix A.2). Colours and symbols as in Figure Figure A.6). Therefore, we are confident that the presentedṀ max ranges and clumping properties are reliable.
In summary, qualitatively the observed SED for HD 198478 can be perfectly reproduced with a strong clumping degree in the innermost wind region, moderate clumping in the intermediate one, and the minimum clumping degree in the outermost ones, with a range of mass-loss rates (maybe a reflex of mass-loss episodes due to stellar pulsations).

HD 80077
The available radio measurements at 3.5 and 6 cm for this eBHG are well determined, and thusṀ max is definite. The best-fit model is obtained in the fixed-regions approach (solid line in corresponding panel in Figure A.6). ForṀ max = 3.45 × 10 −6 M yr −1 a weak clumped inner and intermediate wind, f in cl = 2.5 and f mid cl = 1.8, and an unclumped outer one, f out cl = f far cl = 1 (Table 6) perfectly match the SED of this source. Due to the lack of observations at sub-and mmwavelengths we estimated that f out max = 6 is still consistent with the observed 70 and 100 µm, and 3.5 cm fluxes (magenta-dotted line). Only flux observations at mm regime will better constrain this wind region.

HD 53138 (o 2 CMa)
There are no mm or radio observations available for this source. Therefore, our analysis only provides an upper limit for the maximum mass-loss rate consistent with the observed fluxes at the largest possible wavelengths (at 70 and 100 µm). In the corresponding panel in Figure A.6, we present the solutions for both the fixed-and adapted-regions approaches (solid and magenta-dotted line, respectively).
ForṀ max = 1.8 × 10 −6 M yr −1 , the fixed-regions approach yields f in cl = 2, and f mid cl = f out cl = f far cl = 1 ( Table 6). The best-fit solution, however, is achieved in the adapted-regions approach by slightly increasing the clumping degree in the inner wind region, f in cl = 3 and f mid cl = f out cl = f far cl = 1, provided that r mid = 1.3 R * (Table 7). Due to the lack of observations of this source at λ > 100 µm, f out max = 18 at r 15 R * , and f far max = 50. Only mm and radio observations will definitively constrain thė M max and the clumping structure of the outer wind of this source.

CygOB2#12
This is an exotic eBHG 15 which, in addition, is suspected to be a non-thermal or variable thermal source; several authors have reported short term variations in the measured radio fluxes (Bieging et al. 1989, Scuderi et al. 1998, and more recently Morford et al. (2016) detected flux variations of ∼ 14 days at 21 cm. However, the nature of this variability, whether it is due to changes in the state of ionisation of the wind like in P Cyg, or by the object being a non-thermal source (maybe a binary), still remains to be unveiled.
In the corresponding panel in Figure A.6, we present the best possible solutions in the fixed-regions (solid and magentadotted lines) and the adapted-regions approaches (blue-dashed and orange-dashed-dotted lines). CygOB2#12 is one of the sources in our sample for whichṀ max is constrained by infrared fluxes instead of by radio ones. We noted that the mass-loss rate required to match the estimated low flux at 21 cm,Ṁ max ∼ 2.4 × 10 −6 M yr −1 , overestimates by large the observed fluxes at MIR. From MIR and FIR fluxes, we deriveḋ M max = 1.02 × 10 −6 M yr −1 .
For thisṀ max , the fixed-regions approach yields f in cl = 1, f mid cl = f out cl = 10 and f far cl ranges from 5 to 15 (Table 6) to consistently reproduce the variable radio flux values (solid and magenta-dotted lines, respectively). The best-fit solution, however, is achieved in the adapted-regions approach for a larger clumping degree in the intermediate wind region, f in cl = 1, f mid cl = 12, f out cl = 10 and f far cl = 5 -15 (blue-dashed and orangedashed-dotted lines, respectively), provided that r mid = 2.5 R * and r out = 8 R * (Table 7). In all simulations the flux at 1.3 mm constrains f out max = 12.

HD 24912 (ξ Per)
The clumping stratification of this source was also analysed by P06. These authors derived two possible solutions from fitting the variable upper limits of the radio flux estimates at 0.7, 2, 3.5 and 6 cm. In Figure A.3 we show these best-fit solutions from PO6. The first solution (solid line) is consistent with the larger radio upper flux limits, and the second one (magentadashed line) with all radio upper flux limits. It can be observed that both solutions by P06 overestimate our new flux estimates at 70 and 100 µm and, moreover, that the additional and well determined fluxes at 3.6 cm and 21 cm used in our analysis seem to confirm variable radio emission from this source. Note that although these new radio observations are slightly larger (which can be due to either calibration effects or real radio variability of the source), they are also consistent with the upper flux limits previously detected at the radio regime. Moreover, since the radio fluxes estimated for other sources in our sample (λ Cep and α Cam) by Schnerr et al. (2007) are consistent with the existing radio observations by different authors at different wavelengths, we are also confident in the estimated values provided by the authors for this source.
We obtainṀ max = 1.4 × 10 −6 M yr −1 . In order to fit the aforementioned variable radio fluxes we derived several solutions for the clumping structure. The corresponding fitdiagram in Figure A.7 displays our fit-solutions in the fixedregions (solid and magenta-dotted lines) and adapted-regions approaches (blue-dashed and orange dashed-dotted lines). In the fixed-regions approach an almost constant and moderate clumping degree for the entire wind, f in cl = f mid cl = 3.5 and f out cl = f far cl = 3 (second value in Table 6), fits perfectly the upper flux limits and the well determined measurement at 21 cm (magenta-dotted line in Figure A.7), whereas a moderate clumped wind restricted to r 15 R * and an unclumped outer one, f in cl = f mid cl = 3.5 and f out cl = f far cl = 1 (first value in Table 6), are consistent with all radio fluxes but the flux value at 21 cm (solid line in Figure A.7).
The best-fit solution for this source, however, is achieved in the adapted-regions approach for a much more structured clumped wind. For the estimatedṀ max , a significantly clumped inner wind, f in cl = 5, and a constant f mid cl = f out cl = f far cl = 1 (first value in Table 7) provide the best-fit solution consistent with all radio data but the flux value at 21 cm (blue-dashed line in Figure A.7), with r in = 1.1 R * and r mid = 4 R * . To fit the larger upper radio flux limits and the flux at 21 cm (orange dashed-dotted lines in Figure A.7), the clumping degree for the outer wind regions have to be increased to f out cl = f far cl = 3 (second value in Table 7). The lack of observations in the mm regime gives f out max = 5. Finally, the new measured FIR fluxes for this source enforce a certain clumping degree within 1.05 -1.1 r/R * 2 -4. The observed variability at the MIR (from 4.6 to 11.4 µm) and radio regimes may be explained by means of either clumping or some other physical process such as co-rotational interaction zones or a wind compressed equatorial region, as suggested by P06. In particular, a potential magnetic field interacting with the stellar wind leads to wind confinement and channelling (ud-Doula y Owocki 2002, Owocki et al. 2008, and in presence of relativistic electrons, non-thermal synchrotron emission arises and may dominate the radio spectrum under certain conditions (Trigilio et al. 2004). Another possibility is that ξ Per is a binary system, or both magnetic field and colliding winds are present as some numerical simulations suggest (van Loo et al. 2005). However, positive magnetic field detections for ξ Per have not been reported so far, and there is no observational evidence of a close companion, although several sources in its surroundings have been detected. Therefore, the nature of this variability remains unanswered, and only mm and additional radio observations can help ascertain it.

HD 36816 (λ Ori A)
This source was previously analysed by P06, and there is a new distance estimate available from GAIA. All radio flux observations are upper limits, and consequently the estimatedṀ max is an upper limit. In the right panel in Figure A.7 we present P06's best-fit solution (magenta-dotted line) and the best-fit model de-rived in this work (solid line), corresponding to the fixed-regions approach.
P06 derivedṀ max 0.4 × 10 −6 M yr −1 (for d = 0.5 kpc) and f in cl = 2, f mid cl = f out cl = f far cl = 1 (with r in = 1.05 R * , r mid = 2 R * ,r out = 10 R * and r far > 50 R * ). Moreover f mid max = 20 and f out max = 2 (blue-dashed line), respectively. The clumping stratification derived by P06 for the upper limitṀ max is consistent with the additional MIR observations available, but slightly overestimates the obtained fluxes at 70 and 100 µm. The best-fit solution, however, is found by decreasing the clumping degree in the inner wind region to f in cl = f mid cl = f out cl = f far cl = 1. There are no well-defined flux observations available at wavelengths larger than 100 µm. Despite that, we can lower the maximum clumping degree for the intermediate and outer wind regions derived by P06 to f mid max = 5 and f out max = 2. Only additional, and well determined observations at mm and radio regimes can provide definiteṀ max and clumping factors in Regions 4 and 5, f out cl and f far cl , respectively. For the GAIA distance used here, d= 0.27 kpc, the scaled upper limit maximum mass-loss rate becomesṀ max 0.16 × 10 −6 M yr −1 , and the derived clumping structure described above (Table 6) provide the best-fit solution.

Appendix A.4: OB Dwarfs
HD 149757 (ζ Oph) All available mm and radio observations for this rapidly rotating source are upper limits, thereforeṀ max is also an upper limit. We noted that if we deriveṀ max to fit the radio upper limits at 6 cm, the flux-model emission is consistent with the observations from V-to NIR-band but overestimates fluxes at MIR and FIR wavelengths. Therefore, the upper limitṀ max is constrained by the fluxes at 70 and 100 µm. The left panel in Figure A.8 displays our best-fit solution, which corresponds to the fixed-regions approach. This best-fit solution (solid line) is obtained byṀ max 0.07 × 10 −6 M yr −1 and a constant f in cl = f mid cl = f out cl = f far cl = 1 ( Table 6). This model reproduces all well-defined flux observations and is consistent with all upper limits in the radio regime.We also obtain f out max = 20 and f far max = 6 (magenta-dotted line). Additional well defined mm and radio fluxes are needed to derive definiteṀ max and clumping properties in the outermost wind region.

HD 149438 (τ Sco)
τ Sco is a confirmed magnetic B-star (Petit et al. 2013), and a very slow rotator (∼ 5 km s −1 ; Donati et al. 2006). Mass-loss rates found in the literature range fromṀ = 0.0013 × 10 −6 M yr −1 (from UV resonance lines; Hamann 1981), toṀ = 0.0614 × 10 −6 M yr −1 (from optical spectroscopy; Mokiem et al. 2005), with a value in-betweeṅ M = 0.02× 10 −6 M yr −1 (from NIR spectroscopy; Repolust et al. 2005). We note that all these derivations used non-magnetic wind models to derive the mass-loss rate, and thus did not take into account the considerable effect the magnetic field has upon the wind structure.
The right panel in Figure A.8 displays our best-fit models for this source, corresponding to the fixed-regions approach. Whereas the upper limits at 1 and 3 cm (Kurapati et al. 2017) are consistent with thermal emission, those at 6 cm (Bieging et al. 1989) and 13 cm (Kurapati et al. 2017) point to non-thermal emission. In view of the SED trend up to 160 µm, we suspect that the likely flux value at 13 cm is much lower than the upper limit provided by Kurapati et al. (2017), and therefore we do not attempt to fit it. Despite that, the upper flux limits from 1 to 6 cm still suggest τ Sco is a non-thermal emitter. Since it is a magnetic star, the observed radio variability might be explained by non-thermal synchrotron emission. However, non-thermal radio emission seems to favour centrifugal magnetospheres instead of dynamical ones (Kurapati et al. 2017), and τ Sco is classified as having a dynamical magnetosphere (Petit et al. 2013). Therefore, the detected variability and its nature still needs to be confirmed by deeper radio observations. Nevertheless, we derived two possible upper limits foṙ M max , attending to variable radio upper flux limits, and being consistent with the well determined fluxes at 70 and 100 µm. A first solution (magenta-dotted line) is obtained for aṀ max 0.315 × 10 −6 M yr −1 and a constant f in cl = f mid cl = f out cl = f far cl = 1 (first entry in Table 6). Although this model perfectly matches the observations at λ 100 µm and the upper limit at 6 cm, it overestimates by large the radio upper flux limits at shorter wavelengths. A second solution (solid line) is achieved for a much lowerṀ max 0.018 × 10 −6 M yr −1 , an extremely large clumping degree in the inner wind region, f in cl = 300, and f mid cl = f out cl = f far cl = 1 (second entry in Table 6). In this case, the model consistently reproduces all observations and the radio upper limits, but the flux at 70 µm only marginally, and slightly underestimates the one at 100 µm. The FIR fluxes and the radio upper limit at 6 cm could be perfectly matched with this value ofṀ max by increasing the clumping degree of the wind at r 2 R * as f mid cl = f out cl = f far cl = 300. However, the radio upper limits at 1 and 3 cm still remains overestimated by large. We caution, however, that none of the above takes into account any effects of the magnetic field upon the wind structure. Additional and well determine flux observations at mm and radio regime, along with proper magnetic modelling, may help reaching conclusive results regarding this source.

HD 149404
This source is a known massive spectral binary (Rauw et al. 2001, Raucq et al. 2016. From the only reliable radio observation at 3.6 cm we obtainṀ max = 8.27 × 10 −6 M yr −1 . The corresponding panel in Figure A.9 displays the best possible solution in the fixed-regions (solid line) adapted-regions approaches (magenta-dotted line). The best possible solution in the fixed-regions approach is obtained for f in cl = f out cl = f far cl = 1 and f mid cl = 5.4 (Table 6). The best-fit solution, however, is achieved in the adapted-regions approach (magenta-dotted line) by increasing the clumping degree in the intermediate wind region, f in cl = f out cl = f far cl = 1 and f mid cl = 7, provided that it starts at r mid = 2.5 R * (Table 7). For both approaches, f out max = 6 is still consistent with the data (blue-dashed line).

HD 37043 (ι Ori)
ι Ori is a known binary (Pittard et al. 2000, Pablo et al. 2017, with consistent observed variable flux at radio wavelengths. This source was previously analysed by P06 as a thermal emitter. In the corresponding panel in Figure A.9, we present the two solutions derived by P06, together with our best-fit model, which corresponds to the fixed-regions approach. P06 derived a first solution (blue-dashed line) consistent with both the upper limit at 2 cm and the well determined radio flux at 6 cm:Ṁ max = 0.8 × 10 −6 M yr −1 , f in cl = f mid cl = f out cl = f far cl = 1 with r in = 1.05 R * , r mid = 2 R * , r out = 10 R * and r far > 50 R * ; and a second one (orange dashed-dotted line) fitting only the lower flux value at 3.5 cm:Ṁ max = 0.25 × 10 −6 M yr −1 , f in cl = 12, f mid cl = f out cl = f far cl = 1 with r in = 1.05 R * , r mid = 1.3 R * , r out = 10 R * , r far > 50 R * .
As shown in the fit-diagram, both models provide a marginal fit of the flux value at 70 µm. However, if HD 37043 is a non-thermal source, the lowerṀ max solution by P06 represents the upper and lower limits of the mass-loss rate and clumping structure for this source, respectively, even in the case that our estimated flux value for this source at 70 µm is wrong. Note that this solution is also consistent with the upper measured radio fluxes by increasing the clumping in the outer wind regions. We found that the best-fit for all the data (solid line) is provided by the lowest value ofṀ max derived by P06 (Ṁ max = 0.25 1 (first values in Table 6). To reproduce the upper radio fluxes (magenta-dotted line), large clumping factors at the outer wind regions are required, f out cl = 5 and f far cl = 15 (second values in Table 6). In all simulations f out max = 10. In this case, our best-fit solution is only slightly better than that obtained for the loweṙ M max and the strong inner clumped wind solution derived by P06.