Jacopo Chevallard

Stéphane Charlot

Peter Senchyna

Daniel P Stark

Alba Vidal-García

Anna Feltre

Julia Gutkin

Tucker Jones

Ramesh Mainali

Aida Wofford

Physical properties and H-ionizing-photon production rates of extreme nearby star-forming regions

Keywords: galaxies: evolution -galaxies: ISM -galaxies: dwarf -HII regions -dark ages, reionization, first stars -methods: data analysis

Measurements of the galaxy UV luminosity function at z 6 suggest that young stars hosted in low-mass star-forming galaxies produced the bulk of hydrogen-ionizing photons necessary to reionize the intergalactic medium (IGM) by redshift z ∼ 6. Whether star-forming galaxies dominated cosmic reionization, however, also depends on their stellar populations and interstellar medium properties, which set, among other things, the production rate of H-ionizing photons, ξ ion , and the fraction of these escaping into the IGM. Given the difficulty of constraining with existing observatories the physical properties of z 6 galaxies, in this work we focus on a sample of ten nearby objects showing UV spectral features comparable to those observed at z 6. We use the new-generation tool to model the UV-to-optical photometry and UV/optical emission lines of these Local 'analogues' of high-redshift galaxies, finding that our relatively simple, yet fully self-consistent, physical model can successfully reproduce the different observables considered. Our galaxies span a broad range of metallicities and are characterised by high ionization parameters, low dust attenuation, and very young stellar populations. Through our analysis, we derive a novel diagnostic of the production rate of H-ionizing photons per unit UV luminosity, ξ ion , based on the equivalent width of the bright [O ] λλ4959,5007 line doublet, which does not require measurements of H-recombination lines. This new diagnostic can be used to estimate ξ ion from future direct measurements of the [O ] λλ4959,5007 line using JWST/NIRSpec (out to z ∼ 9.5), and by exploiting the contamination by Hβ + [O ] λλ4959,5007 of photometric observations of distant galaxies, for instance from existing Spitzer/IRAC data and from future ones with JWST/NIRCam.

INTRODUCTION

The appearance of the first stars and galaxies in the Universe at the end of the Dark Ages marked the beginning of the epoch of reionization (EoR, e.g. Bromm & Yoshida 2011). The EoR appears galaxies were more abundant in the early Universe than at more recent times, potentially providing enough H-ionizing photons to reionize the IGM (e.g. Robertson et al. 2015;Bouwens et al. 2015b). However, this interpretation is subject to major uncertainties: current data do not allow any accurate determination of the abundances of low-mass galaxies during the EoR, because of degeneracies among the luminosity function parameters (e.g. Finkelstein et al. 2015) and the difficulty of probing magnitudes fainter than M UV ∼ -17 (except in lensed fields, e.g. Atek et al. 2015a,b;Ishigaki et al. 2018). Another major source of uncertainty is the relation between the (non-ionizing) UV-continuum emission of EoR galaxies, their production rate of H-ionizing photons, N ion , and the number of these photons escaping from the galaxies, f esc . Direct measurements of the rate of ionizing photons escaping from galaxies (the product f esc × N ion ) are only possible at lower redshifts (e.g. Izotov et al. 2016;Vanzella et al. 2016), since any H-ionizing photon emitted by EoR galaxies would be absorbed by neutral H along the line of sight. Indirect diagnostics on f esc have been proposed, for instance based on the relation between H-Balmer lines and UV-continuum emission (Zackrisson et al. 2013) and on the Lyα-line profile (Verhamme et al. 2015(Verhamme et al. , 2017), but ), but have not yet been tested on EoR galaxies.

In this paper, we focus on the production rate of H-ionizing photons per unit intrinsic monochromatic UV (1500 Å) luminosity, ξ ion = N ion /L UV . This quantity depends only on the properties of the stellar populations emitting the ionizing photons, i.e. mainly age and metallicity. Analytical reionization models often assume a fixed value for ξ ion , typically in the range log(ξ ion /erg -1 Hz) ∼ 25.2 -25.3 (e.g. Kuhlen & Faucher-Giguère 2012;Robertson et al. 2013), although the dependence of ξ ion on stellar population properties would introduce galaxy-to-galaxy variations (both stochastic and systematic). Here, we use the newgeneration spectral analysis tool (Chevallard & Charlot 2016) to measure the production rate of H-ionizing photons of ten extreme nearby star-forming regions found by Senchyna et al. (2017, hereafter S17) to have UV spectra similar to those of z ∼ 6-7 star-forming galaxies. Our simultaneous modelling of several UV and optical emission lines also provide valuable constraints on the properties of the ionizing stars, photoionized gas and attenuation by dust. Beyond the demonstration that the stellar population and photoionization models incorporated into offer accurate fits to the data, this analysis allows us to uncover and calibrate a relation between ξ ion and equivalent width of (Hβ+) [O ] λλ4959,5007. This relation extends from canonical ξ ion values at low line equivalent widths up to nearly an order or magnitude larger for the most extreme galaxies.

DATA AND MODELLING APPROACH

Spectro-photometric observations of Local 'analogues' of high-redshift galaxies

We appeal to a sample of ten nearby analogs of primeval galaxies, for which our team has acquired high-quality UV data from HST/COS (G160M and G185M gratings) and optical spectra from MMT (see S17 for details). These galaxies were selected from a sample of He λ4686 emitters dominated by stellar photoionization, identified by Shirazi & Brinchmann (2012) from the SDSS DR7 (Abazajian et al. 2009). These objects show extremely large equivalent widths of optical emission lines (e.g. [O ] λλ4959,5007 1000 Å), which, although rare in the nearby Universe (only ∼ 0.01 per cent of SDSS galaxies), they become far more common at higher redshifts (e.g. Smit et al. 2014;Rasappu et al. 2016). The ten galaxies in our sample have hard ionizing spectra able to doubly ionize He, i.e. with substantial flux at λ < 228 Å (54.4 eV), although only half show Wolf-Rayet (WR) spectral features. Our follow-up HST observations of these objects succeeded in revealing highionization UV lines (C λλ1548,1551, [C ] λ1907+C ] λ1909, O ] λλ1660,1666) similar to those typically observed in galaxies at z 5 (Stark et al. 2015a(Stark et al. ,b, 2017;;Mainali et al. 2017;Schmidt et al. 2017;Smit et al. 2017), providing a unique reference sample for high-redshift studies.

The location of our galaxies in the '[N ] ' Baldwin et al. (1981, BPT) diagram shown in Fig. 1a is typical of sub-solar metallicity gas (Z/Z ∼ 0.1 -0.7) with a high ionization parameter (see, e.g., fig. 2 of Gutkin et al. 2016). The large O32 values at fixed R23, shown in Fig. 1b, also indicate a high ionization parameter, although these values can also be caused by density-bounded H regions (e.g. Jaskot & Oey 2013;Stasińska et al. 2015). As we show below, ionizationbounded regions are favoured by our model, which well matches the observed [O ] λ6363 fluxes, which should appear weaker for density-bounded regions. Fig. 1a also shows that our galaxies have lower [N ] λ6584/Hα, at fixed [O ] λ4959/Hβ, than those in the KBSS sample of Strom et al. (2017). These ratios are similar to those of the z ∼ 0.3 LyC leakers of Izotov et al. (2016), although they do not reach the extreme [N ] λ6584/Hα ∼ -2.5 of the compact star-forming galaxies of Izotov et al. (2017b). Interestingly, Fig. 1b shows no significant offset between our galaxies, the KBSS sample and the LyC leakers in the O32 vs R23 diagram, while Izotov et al. (2017b) specifically targeted extreme O32 emitters. We also note the presence of a broad (a few × 10 2 km/s) component in several emission lines ([O ] λλ4959,5007, Hβ, Hα), at a level of a few per cent the intensity of the narrow component. Such lowintensity broad components have been observed in other low-mass, low-metallicity, highly star-forming galaxies and linked to energy injection from stellar winds and supernovae (e.g. Izotov et al. 2007;James et al. 2009;Amorín et al. 2012), since Active Galactic Nuclei (AGNs) typically produce broader components ( 10 3 km/s) with intensities comparable to the narrow ones (e.g. Izotov & Thuan 2008;Izotov et al. 2010).

Emission line fluxes measurements

We measure the fluxes of the emission lines listed in Table 1-2 and their associated errors with a custom fitting software described in S17. In summary, for each line, or group of neighbouring lines, we consider a linear function to describe the local continuum and fit the line profile with a Gaussian function. For lines with multiple kinematic components, we consider a second Gaussian function. We do not attempt to subtract the underlying stellar absorptions from the measured line fluxes, since stellar absorptions are self-consistently included in our stellar population + photoionization model (see Section 2.3 below). We rescale the formal errors on the line fluxes measured from SDSS spectra by a factor 1.8, the average error scaling suggested by the MPA/JHU analysis. 1 As in S17, we correct the line fluxes for Galactic extinction using the dust map of Schlafly & Finkbeiner (2011) and assuming the R V = 3.1 extinction curve of Fitzpatrick (1999). Since the predictions of our model include only H-Balmer transitions of principal quantum number n ≤ 6, where n = 6 corresponds to Hδ, we subtract from the measured He λ3889 flux the contribution from the Hζ λ3889 line using the measured Hδ flux, and assuming a case B recombination ratio of Hζ/Hδ = 0.41 valid for a T = 104 K gas (see Table 4.2 of [START_REF] Osterbrock | Astrophysics Of Gas Nebulae and Active Galactic Nuclei[END_REF]). 2 Given the similar (circular) apertures of SDSS fibers and HST/COS (diameters of 3 and 2.5 , respectively), we do not apply any aperture correction to the UV and optical emission lines measurements.

± 0.1 -1.35 -1.46 ± 0.08 -1.39 [O ] λ6363 -2.75 ± 0.3 -2.52 -2.61 ± 0.1 -2.45 -2.41 ± 0.08 -2.39 -2.90 ± 0.2 -2.68 -2.45 ± 0.1 -2.24 [N ] λ6548 -2.53 ± 0.2 -2.49 -1.95 ± 0.09 -1.73 -2.22 ± 0.07 -2.22 -2.23 ± 0.1 -2.16 -1.85 ± 0.09 -1.93 [N ] λ6583 -2.01 ± 0.06 -2.01 -1.41 ± 0.07 -1.24 -1.68 ± 0.05 -1.75 -1.73 ± 0.1 -1.69 -1.32 ± 0.08 -1.45 [S ] λ6716 -1.51 ± 0.03 -1.53 -1.29 ± 0.06 -1.26 -1.37 ± 0.04 -1.41 -1.62 ± 0.1 -1.62 -1.21 ± 0.07 -1.19 [S ] λ6730 -1.68 ± 0.04 -1.64 -1.41 ± 0.06 -1.37 -1.48 ± 0.04 -1.52 -1.76 ± 0.1 -1.73 -1.35 ± 0.

Stellar-population and photoionization modelling

We model the UV-to-optical photometry and emission lines integrated fluxes of our galaxies to constrain their physical properties and production rate of H-ionizing photons. To achieve this, we use , a new-generation spectral analyses tool which incorporates a self-consistent description of the emission from stars (based on the latest version of the Bruzual & Charlot 2003 stellar population synthesis code 3 ) and photoionized gas (using version 13.3, last described in Ferland et al. 2013). In practice, we consider a Chabrier (2003) our galaxies is hence parametrised in terms of the galaxy stellar mass M, age of the oldest stars t (or age of the galaxy) and stellar metallicity Z. Following S17, we account for the effect of dust attenuation by considering the Small Magellanic Cloud extinction curve [START_REF] Pei | [END_REF], and vary the V-band attenuation optical depth τV . The choice of using an extinction curve is justified by the fact that geometry and scattering effects should have only a minor influence on line emission from the small physical regions probed by our data.

We follow Gutkin et al. (2016) to account for (line+continuum) nebular emission. We fix the gas density at n H = 10 2 cm -3 , the typical value observed in local compact star-forming galaxies (e.g. Izotov et al. 2006;Jaskot & Oey 2013), set the interstellar metallicity equal to the stellar metallicity (Z ISM = Z), and vary the ionization parameter log U S , which sets the ratio of H-ionizing-photon to gas density at the Strömgren radius, 4 the dust-to-metal mass ratio ξ d , which sets the depletion of heavy elements onto dust grains, and the carbon-to-oxygen abundance ratio C/O, where C/O = 1 implies a scaled-Solar abundance.

We perform two separate fittings with : one to a combination of UV and optical photometry and Balmer emission lines, to derive the intrinsic UV luminosity of our galaxies; the other to a broad set of UV and optical emission lines, to derive the physical properties of the ionized gas and ionization sources. In the combined photometry+Balmer lines fit, we consider the H-Balmer lines Hα,

exp -τV ∈ [0, 2] V -band attenuation optical depth C/O U ∈ [0.1, 1] Carbon-to-oxygen ratio ξ d U ∈ [0.1, 0.5]
Dust-to-metal mass ratio Table 3. Free parameters of the model adopted in to fit the COS/SDSS photometry + Balmer emission lines [parameters log(M/M ), log (t/yr), log(Z/Z ), log U S and τV ], and to fit the emission line fluxes (all parameters in the table) of the ten nearby star-forming galaxies analysed in this work. The symbol U indicates a uniform distribution in the specified range, N(µ; σ) a Gaussian (Normal) distribution with mean µ and standard deviation σ, truncated in the specified range.

Hβ, Hγ and Hδ, the ugri optical (fiber) photometry from SDSS,5 and a single UV photometric band, which we compute from an integration of the HST/COS spectra. 6 Because of the different filter configurations used when acquiring the COS data, we consider a flat filter covering the region λ ∼ 1620 -1660 Å for the objects SB 36 and SB 182, and a filter covering the range λ ∼ 1480 -1515 Å for the other 8 objects. Both regions sample the rest-frame UV continuum of the galaxies, and are not significantly affected by nebular emission lines or ISM absorption lines. Since Balmer lines alone cannot constrain the several nebular emission parameters, we simplify the model and fix the carbon-to-oxygen abundance ratio to the Solar value [C/O = (C/O) ] and the depletion factor to ξ d = 0.3, to obtain a model described by the 5 free parameters

Θ phot+Balmer = [log(M/M ), log (t/yr), log(Z/Z ), log U S , τV ].
In a separate fitting, we model the broad range of UV and optical emission lines listed in Table 1-2. Besides the H-Balmer lines Hα, Hβ, Hγ and Hδ, we include several O, O + and O 2+ lines, the nitrogen and sulphur doublets [N ] λλ6548,6584 and [S ] λλ6716,6731, the [C ] λ1907+C ] λ1909 line, helium lines, and lines from Ar 2+ , Ar 3+ and Ne 2+ . We fit with the integrated line fluxes, computed from the observed spectra as detailed in Section 2.2, with a model similar to the one described above, but letting the carbon-to-oxygen abundance ratio and depletion factor free to vary. This model is therefore described by the 7 free parameters

Θ lines = [log(M/M ), log (t/yr), log(Z/Z ), log U S , τV , C/O, ξ d ].
The Bayesian approach adopted in allows us to derive the posterior probability distribution of the model parameters Θ phot+Balmer and Θ lines ,

P(Θ | D, H) ∝ π(Θ) L(Θ) , (2.1)
where D indicates the data, H the adopted model, π(Θ) the prior distribution, and L(Θ) the likelihood function. We adopt independent priors for each parameter, uniform for log(M/M ), log(Z/Z ), log U S , ξ d , C/O, Gaussian for log (t/yr), and exponential for τV , and report the prior definitions and ranges in Table 3. We adopt a Gaussian likelihood function with independent errors on each measurement,

-2 ln L(Θ) = i y i -ŷi (Θ) σ i 2 , (2.2)
where the summation index i runs over all observables considered in the fitting, i.e. the Balmer lines and COS+SDSS photometry in one case, the emission lines listed in Table 1-2 in the other case.

In equation 2.2, y i indicates the data, σ i the observational error and ŷi ( Θ ) the model prediction. We adopt the Nested Sampling algorithm (Skilling et al. 2006) in the implementation (Feroz et al. 2009) to sample the posterior probability distribution of the free parameters of our model. We note that this implies that each model spectrum is computed in 'on-the-fly' for each combination of parameters selected by during the sampling of the posterior probability distribution of the model parameters.

RESULTS

The model we have adopted allows us to well reproduce all the observables here considered , i.e. the combination of COS/SDSS photometry and H-Balmer lines, and the emission lines of Table 1-2. In Table 1-2, we also report the emission line fluxes predicted by our model for each object fitted, while we show in Fig. 2 a visual example of the fitting of the object SB 80, an embedded H region at z = 0.01085. We show in Figs 1a,b, the line ratios predicted by for all fitted objects in our sample. Table 1-2 indicates a remarkable agreement between the predictions of our relatively simple model and the data.

Physical properties of extreme star-forming regions

Table 4 lists the physical properties of our sample derived from the -based emission lines fitting. As already shown in S17, our ten galaxies cover a broad range of gas-phase metallicities 12 + log(O/H) ∼ 7.9 -8.5 [log(Z/Z ) ∼ -0.8 to -0.2] and show low dust attenuation values. We derive high ionization parameters in all galaxies [log U S -3], even in the highest metallicity ones, namely SB 179 and SB 191 [12 + log(O/H) ∼ 8.5], which show log U S ∼ -2.60 and ∼ -2.14, respectively. Galaxies showing intermediate metallicities and high ionization parameters are rare in the Local Universe (e.g. among SDSS galaxies, see fig 2 of Carton et al. 2017), while they become more common at higher redshifts, as suggested by the detection of high-ionization metal lines at z 5 (e.g. Stark et al. 2015a;Mainali et al. 2017). The emission lines that we consider in this analysis originate from both refractory (O, C) and non-refractory (N, S, Ar, Ne) elements, and this enables us to place tight constraints on the fraction of metals depleted onto dust grains ξ d . The ξ d values reported in Table 4 are significantly lower than the Solar value ξ d = 0.36, indicating that most metals in these galaxies are in the gas phase and not locked into dust grains.

Fig. 3 compares the gas-phase metallicity 12 + log(O/H) derived with our analysis with the values obtained by S17 using a standard method, the 'direct-T e ' one (see their section As Fig. 3 shows, the agreement between the S17, direct-T e based estimates of 12 + log(O/H) and our values (listed in Table 4) is generally good, although the -based metallicities are slightly larger than the S17 ones, by ∼ 0.05 dex at 12 + log(O/H) 8.3, and by ∼ 0.2 dex for the two galaxies with larger 12 + log(O/H) values, SB 179 and SB 191. These two objects show large ionization parameters, and S17 ignore the contribution from O 3+ when measuring the oxygen abundance, but they evaluate the contribution of O 3+ to be minor by considering a relation between O 3+ and He 2+ /He + (Izotov et al. 2006). The only object that does not follow the trend above is SB 198, for which our metallicity estimate is ∼ 0.3 dex lower than the S17 one. The metallicity estimated for SB 198, however, is the most uncertain, and the statistical significance of the difference among the two methods is ∼ 3 σ, hence more data would be necessary to gain better insight into the origin of this discrepancy.

δ H γ H β H α O i i i ] λ 1 6 6 1 O i i i ] λ 1 6 6 6 C i i i ] λ 1 9 0 8 [ O i i ] λ 3 7 2 7 [ N e i i i ] λ 3 8 6 8 H e i λ 3 8 8 9 [ O i i i ] λ 4 3 6 3 [ A r i v ] λ 4 7 4 0 [ O i i i ] λ 4 9 5 9 [ O i i i ] λ 5 0 0 7 H e i λ 5 8 7 5 [ O i ] λ 6 3 6 3 [ N i i ] λ 6 5 4 8 [ N i i ] λ 6 5 8 3 [ S i i ] λ 6 7 1 6 [ S i i ] λ 6 7 3 0 [ A r i i i ] λ 7 1 3 6 -16 -14 log(F/erg s -1 cm -2 ) Figure 2.
We note that differences between metallicities estimated with photoionization models and with the 'direct-T e ' method, such as those discussed above, have been reported in several recent works (e.g. Blanc et al. 2015;Vale Asari et al. 2016;Pérez-Montero & Amorín 2017). Interestingly, both Blanc et al. (2015) and Vale Asari et al. (2016) find that metallicities estimated with photoionization models are systematically larger (by ∼ 0.2 to 0.4 dex) than those estimated with the direct-T e method, qualitatively consistent with our findings, and that photoionization models can provide metal-licity estimates in better agreement with those derived from metal recombination lines. In Vale Asari et al. (2016) (see their sec. 6.3) it is suggested that the [O ] λ4363 line is at the origin of the measured differences among photoionization-and direct-T e -based metallicities. The [O ] λ4363 line is used to determine the electron temperature in the O 2+ zone, but this line can be affected by temperature fluctuations and non-Maxwellian velocity distributions of the free electrons (Nicholls et al. 2012). These effects can boost the [O ] λ4363 line flux, thus modifying the derived electron temperature and, consequently, the metallicity. In this work, the differences among the photoionization-and direct-T e -based metallicities are, for most objects, significantly smaller than found in previous works. The reason is likely related to the fact that S17 directly measure the T e of both O + and O 2+ , instead of only measuring the temperature of O 2+ and using calibrated relations to derive that of O + . The S17 approach can thus alleviate some of the limitations of the direct-T e method, providing estimates in better agreement with those derived with the photoionization models incorporated into .

Production rate of H-ionizing photons

In this paper, we define ξ ion = N ion /L UV as the ratio between the production rate of H-ionizing photons ( N ion ) and the intrinsic stellar monochromatic UV luminosity (L UV ). The quantity L UV is computed from the pure stellar spectrum considering a flat filter centred at 1500 Å and of 100 Å width (e.g. see Robertson et al. 2013). L UV is the emission from stars that would be observed in the absence of gas and dust in galaxies, i.e. it does not include the effect of absorption and emission from neutral and photoionized gas, nor of dust attenuation. Previous studies have adopted different definitions of the monochromatic UV luminosity, for example, by considering the stellar+nebular luminosity L H UV , corrected for dust attenuation outside H regions but not for dust attenuation inside H regions nor for recombination-continuum emission, or the observed luminosity L UV . For the sake of clarity, we introduce here the symbols ξ H ion and ξ ion to differentiate between our definition of ξ ion and the ratios ξ H ion = N ion /L H UV and ξ ion = N ion /L UV . To compute ξ ion , we combine the monochromatic UV luminosity (L UV ) inferred from the fitting to COS/SDSS photometry + Balmer lines with the rate of H-ionizing photons ( N ion ) obtained by fitting the different UV and optical emission lines. We report in Table 5 the ξ ion values for the ten galaxies in our sample. Our galaxies exhibit log(ξ ion /erg -1 Hz) values in the range ∼ 25.2 to ∼ 25.8, i.e. generally larger then the values adopted in standard reionization models (25.2 to 25.4). We also report in Table 5 the values we derive for ξ H ion and ξ ion , to highlight the impact that different definitions of this quantity can have on the quoted values. In particular, ξ H ion can be either smaller or larger than ξ ion , depending on the relative strengths of dust attenuation inside H regions (which lowers L H UV and raises ξ H ion ; see e.g. object SB 191) and recombinationcontinuum emission (which boosts L H UV and hence lowers ξ H ion ; see e.g. SB 2). Considering the observed (i.e. attenuated) UV luminosity L UV generally implies larger ξ ion values, but the difference between ξ ion and ξ ion depends on the amount of dust attenuation affecting stellar and nebular emission outside H regions.

In Fig. 4, we show the relation between ξ ion and equivalent width of the [O ] λλ4959,5007 doublet. Objects exhibiting larger EW([O ]) also show larger ξ ion . A relation between the two quantities is expected, since large [O ] λλ4959,5007 equivalent widths are typically produced by stellar populations with young ages and sub-solar metallicities, which would also produce copious amounts of H-ionizing photons. Moreover, part of the small scatter visible in Fig. 4 can be ascribed to the relatively narrow range of [O ] λλ4959,5007/Hβ probed by our objects (Fig. 1a). Remarkably, the tight relation in Fig. 4 also implies a limited variation of ionized-gas properties in our objects, since the relation between ξ ion and EW([O ]) depends on the shape of the ionizing spectrum (e.g. stellar age and metallicity), ionization parameter (geometry of the photoionized regions), metal abundances and gas density. The colour coding of the circles in Fig. 4 We fit the relation between ξ ion and EW([O ]) using 'Orthogonal Distance Regression', a linear regression method which allows us to account for errors along both the y-and x-axis. 7 The black line in Fig. 4 show the best-fitting relation

log(ξ ion /erg -1 Hz) = 22.12±0.48+1.15±0.16×log EW [O ]/Å , (3.1)
while the grey bands indicate the 68 per cent credible region around the best-fitting relation. We also fit a relation considering the combined equivalent widths of Hβ + [O ] λλ4959,5007, obtaining a relation with intercept (slope) = 21.95 (1.18) instead of 22.12 (1.15). Such a relation is particularly useful to convert in a physicallymotived way broad-band colour excesses (e.g. IRAC excess, Shim et al. 2011;Stark et al. 2013;Smit et al. 2014) into ξ ion values.

DISCUSSION

The results presented in the previous section allow us to draw several conclusions regarding our modelling approach. The -based fitting of the UV and optical emission lines, and in particular our ability to match simultaneously the emission from different ionization states of the same element (e.g. O + and O 2+ ; Ar 2+ and Ar 3+ ), ions with widely different ionizing energies (e.g. 27.6 eV for Ar 2+ and 41.0 eV for Ne 3+ ) and electron temperature-sensitive ratios (e.g.

[O ] λ4363/[O ] λλ4959,5007), suggests that the ionizing stellar-population prescription and nebular parametrization of our model are suited to interpret observations of extreme star-forming regions. The matching of the different H-Balmer lines validates our adoption of an SMC extinction curve to model dust attenuation, and that of the [S ] λλ6716,6731 line doublet the gas density we have adopted. Also, the ability of our model to match emission lines of different elements (C, N, O, Ne, S and Ar) suggests that these galaxies show similar metal abundance ratios, which we can reproduce by varying only the interstellar metallicity, dust depletion factor and carbon-to-oxygen ratio. We note, however, that in our analysis we did not attempt to fit the He recombination lines (He λ1640 and He λ4686). As we have discussed in S17, the stellar population + photoionization models incorporated in are unable to reproduce the very large equivalent widths of the He λ1640 nebular line (up to ∼ 1.7 Å) observed in the galaxies of our sample that do not show WR spectral features, i.e. galaxies with metallicities 12 + log(O/H) 8. As discussed in Shirazi & Brinchmann 2012, this problem is common to other widely used stellar population models, which fail in providing enough He + ionizing photons for stellar populations with Z/Z 0.2. Several potential sources of additional He + ionizing photons have been proposed in the literature, such as Xray binaries, fast radiative shocks, and low-luminosity AGNs (e.g. Stasińska et al. 2015). As discussed in sec. 5.3 of S17, data on our galaxies do not provide evidence for the presence of X-ray binaries or fast radiative shocks, although deeper X-ray observations and larger grids of shock models (extending to lower metallicities) would be required to rule out such presence. Thanks to the recent incorporation into of the AGN models for Narrow Line Region emission of Feltre et al. (2016), in a future work we will test the possibility of low-luminosity AGNs to explain the observed He λ1640 nebular emission of our non-WR galaxies.

The constraints obtained with our -based modelling of emission lines + photometry on ξ ion have potentially important implications for our understanding of the sources of Hionizing photons in pristine star-forming galaxies. We find that ξ ion can reach values significantly larger than those considered in many reionization models, especially for galaxies exhibiting EW([O ] λλ4959,5007) 1000. Such galaxies are relatively rare at low redshifts, while they become much more common at higher ones. Very large equivalent widths of optical emission lines are routinely inferred from contamination of broad-band Spitzer/IRAC photometry, corresponding to EW(Hα + [N ] λλ6548,6584 + [S ] λλ6716,6731) up to ∼ 400 for galaxies at z ∼ 4.5 (Stark et al. 2013;Smit et al. 2016;Mármol-Queraltó et al. 2016), reaching 1000 for galaxies at z ∼ 5.2 (Rasappu et al. 2016). Smit et al. (2014) derive EW(Hβ + [O ] λλ4959,5007) 1000 for extremely blue IRAC galaxies at z ∼ 7. Observations of high-ionization UV lines provide a consistent picture: while high equivalent widths of UV lines, such as C λλ1548,1551, [C ] λ1907+C ] λ1909 and O ] λλ1660,1666, are rare at low redshift, they have been routinely detected at higher redshifts (e.g. Stark et al. 2015aStark et al. ,b, 2017;;Mainali et al. 2017). The emerging picture is that at z 6 a combination of high star formation rates per unit stellar mass, young stellar populations, low metallicities and high ionization parameters favour the presence of EW 1000 (10) for the brightest optical (UV) lines, which our analysis suggest correspond to values of ξ ion significantly larger than canonical values used in reionization models.

From a theoretical perspective, simulations indicate that the escape fraction of H-ionizing photons from galaxies is timedependent, being tightly coupled to the energy input from stellar feedback (winds, SNe), which opens 'clear' channels through which ionizing radiation can escape (e.g. Wise & Cen 2009;Trebitsch et al. 2017). The presence of broad components in the emission lines of our objects suggests that such feedback mechanisms are at work in star-forming regions with a high production efficiency of H-ionizing photons. If regions with high ξ ion also exhibit large f esc , then short, powerful bursts of star formation in low-metallicity, lowmass galaxies can provide enough H-ionizing photons to achieve the time-averaged value log( f esc ξ ion /erg -1 Hz) = 24.5 (Bouwens et al. 2016) necessary to keep the IGM ionized.8

COMPARISON WITH PREVIOUS WORK

The importance of constraining the production rate of H-ionizing photons to understand the role of star-forming galaxies in cosmic reionization motivated several recent studies in which this quantity has been estimated using different observables. Bouwens et al. (2016) consider star-forming galaxies at z ∼ 3.8 -5.4 in the GOODS North and South fields, and estimate the production rate of H-ionizing photons N ion from the contamination of Hα + [N ] λλ6548,6584 + [S ] λλ6716,6731 to broad-band photometry. In practice, they derive the Hα luminosity by adopting fixed ratios of [N ]/Hα and [S ]/Hα, and estimate the amount of dust attenuation from the UV-continuum slope. They then derive a dust-corrected UV luminosity L H UV . The ξ H ion they obtain span a broad range of values log(ξ H ion /erg -1 Hz) ∼ 24.2 -26, compatible with the values spanned by our galaxies, with the bulk of their galaxies showing values log(ξ H ion /erg -1 Hz) ∼ 25 -25.6, and the bluest and less luminous galaxies exhibiting the highest values. Schaerer et al. (2016) analyse five LyC-leaking galaxies at z ∼ 0.3 with large EW([O ]) 1000, deriving N ion from the (dust-corrected) Hβ luminosity, and L H UV from fitting the UV (GALEX) and optical (SDSS) photometry of the galaxies. As in Bouwens et al. (2016), Schaerer et al. (2016) derive a dust-corrected UV luminosity, but without correcting for the contribution from recombination continuum and dust within H regions. The values they obtain are in the range log(ξ H ion /erg -1 Hz) ∼ 25.1 -25.5, generally lower than the values we derive for galaxies with comparable EW([O ]). Matthee et al. (2017) study a sample of ∼ 600 Hα emitters at z = 2.2 from the HiZELS survey (Geach et al. 2008). They measure the Hα luminosity from a narrow photometric band centred on Hα + [N ] λλ6548,6584, correcting Hα for the [N ] λλ6548,6584 contribution. They derive dust corrections with different methods, and compute N ion from the dust-corrected Hα luminosity. The ξ H ion they derive cover a broad range of values log(ξ H ion /erg -1 Hz) ∼ 24 -26, with the exact values being sensitive to the adopted dust correction method. They find a strong correlation among ξ H ion and EW(Hα), with the highest equivalent width galaxies exhibiting the large ξ H ion values. Izotov et al. (2017a) study a sample of ∼ 14000 compact star-forming galaxies at z ≤ 1 from SDSS DR12. They derive dust corrections from Balmer line ratios, compute N ion from the dust-corrected Hβ luminosity, and obtain the monochromatic UV luminosity from SED fitting to GALEX and SDSS photometry. They derive ξ H ion values in the range log(ξ H ion /erg -1 Hz) ∼ 24.4 -25.6, with the bluest galaxies exhibiting the largest ξ H ion values. Shivaei et al. (2018) consider a large sample of ∼ 700 galaxies at redshift z ∼ 1.4 -2.6 from the MOSDEF survey (Kriek et al. 2015). They derive N ion from the dust-corrected Hα luminosity, and obtain L H UV by fitting the 3D-HST multi-band photometry with the Bruzual & Charlot (2003) population synthesis code. They apply average corrections to N ion to account for a non-zero escape fraction of H-ionizing photons, and find log(ξ H ion /erg -1 Hz) values in the range ∼ 25.00 to ∼ 25.6 (assuming an SMC extinction curve). As found by Bouwens et al. (2016) and Izotov et al. (2017a), they also find that the bluest galaxies exhibit the largest ξ H ion values. Measurements of the production rate of H-ionizing photons in EoR galaxies are difficult because of the challenges in performing spectroscopic observations of such galaxies with existing facilities. Stark et al. (2015b) and Stark et al. (2017) present detections of Lyα, C λλ1548,1551 and [C ] λ1907+C ] λ1909 in 3 galaxies at z 7, and estimate ξ ion from the stellar populations + photoionization modelling of their UV emission lines and photometry. They derive log(ξ ion /erg -1 Hz) ∼ 25.6 -25.8, comparable to the largest values obtained in the present work.

CONCLUSIONS

We have analysed a sample of ten extreme nearby star-forming regions with high-quality UV-to-optical spectra from HST/COS, MMT and SDSS. Our objects are characterised by low metallicities and high ionization parameters, as suggested by the locus they occupe in the [N ]-BPT diagram and in the O32 vs R23 one. Using a combination of stellar population + photoionization models incorporated into the spectral analysis tool, we have fitted a combination of UV-to-optical photometry + Balmer lines, and a broad set of UV and optical emission lines sensitive to the ionizing spectrum, ionization parameter, electrons temperature and density, metal abundances and dust attenuation. This modelling has allowed us to constrain the intrinsic UV luminosity of the galaxies L UV , their production rate of H-ionizing photons, and their physical properties, such as gas-phase and stellar metallicities, ionization parameters, dust depletion factors, stellar ages and dust attenuation.

The ten galaxies of this work are characterised by large ionization parameters (log U S -3), low dust depletion factors (ξ d ∼ 0.1 -0.3), low dust attenuation and very young stellar ages (t 10 Myr in all but one case). As already shown in S17, our galaxies span a broad range of gas-phase metallicities 12 + log(O/H) ∼ 7.9 -8.5. We have compared the 12 + log(O/H) values estimated with with those estimated in S17 with the direct-T e method, finding an overall good agreement among the estimates, although the -based metallicities are systematically larger (by ∼ 0.05 dex) than the S17 ones. Larger differences (∼ 0.2 dex) among the two methods appear at 12 + log(O/H) 8.3. This confirms previous findings that metallicities derived with photoionization models are systematically larger than those derived with the direct-T e method, although the differences found in this work are significantly smaller than those reported in the literature.

The production rate of H-ionizing photons per unit (intrinsic) UV luminosity of our galaxies spans a wide range [log(ξ ion /erg -1 Hz) ∼ 25.2 -25.8] compatible with the values measured in previous works, but with a strong dependence of ξ ion on the equivalent width of the [O ] λλ4959,5007 line, such that high-EW([O ]) galaxies exhibit larger ξ ion values. Such a relation is not completely surprising, since high-EW([O ]) are typically associated with young, metal-poor stellar populations, which are expected to emit large amounts of H-ionizing photons. Energy injection from young stellar populations (i.e. stellar 'feedback') is also expected to create clear channels in the ISM through which ionizing radiation can escape, and it is likely the cause of the low-intensity broad components emission lines observed in our galaxies. The combination of large ξ ion and f esc in short bursts of star formation in low-mass galaxies can hence supply large quantities of H-ionizing photons to the IGM. To enable more accurate determinations of ξ ion in EoR galaxies, we have calibrated a new relation between ξ ion and EW [O ] λλ4959,5007 (EW [O ] λλ4959,5007 + EW Hβ). This relation can be used to derive ξ ion from broad-band observations of galaxies at z 6, e.g. from existing Spitzer/IRAC observations and future JWST/NIRCam ones, and from spectroscopic observations of [O ] λλ4959,5007 with JWST/NIRSpec out to redshift ∼ 9.5.

Figure 1 .

 1 Figure 1. (a) [N ] BPT diagram showing star-forming galaxies at z ∼ 2 -2.6 from the KBSS survey (blue circles, Strom et al. 2017), the 4 LyC continuum 'leaker' at z ∼ 0.3 from Izotov et al. (2016) (red circles), the 5 compact star-forming galaxies at z < 0.07, candidates LyC continuum leakers, from Izotov et al. (2017b), and the ten galaxies presented in S17 analysed in this work (grey diamonds), along with the model predictions obtained with (red diamonds). As a reference, we also show with a grey-scale (linear) density plot galaxies from the MPA/JHU emission line catalogue based on SDSS DR7 data. MPA/JHU fluxes are corrected for Galactic extinction only, while the other fluxes for internal extinction as well. (b) same as (a), but showing the O32 = log ([O ] λλ4959,5007/[O ] λλ3726,3729) vs R23 = log ([O ] λλ4959,5007 + [O ] λλ3726,3729)/Hβ) diagram.

  stellar initial mass function truncated at 0.1 and 100 M , and adopt a constant star formation rate and a constant metallicity. The star formation and chemical enrichment history of 2 The proximity in wavelength of Hδ at λ = 4101.7 Å to Hζ λ3889 allows us to ignore the potential impact of dust attenuation on the estimation of the Hζ flux, since Hδ and Hζ would suffer a comparable amount of attenuation. 3 See appendix A of Vidal-García et al. 2017 for a comparison of the population synthesis code used in this work with the original version of Bruzual & Charlot (2003).

  3.2). S17 measure the electron temperature (T e ) of O + and O 2+ using ratios of auroral ([O ] λλ7320,7330 and [O ] λ4363) and strong ([O ] λλ3726,3729 and [O ] λλ4959,5007) lines, and then compute the abundances of the O + and O 2+ ions using the strong lines.

H

Figure 3 .

 3 Figure 3. Comparison between the gas-phase metallicity 12 + log(O/H) obtained in this work with the tool (y-axis) and by S17 with the 'direct-T e ' method. Colour-coded is the ionization parameter log U S obtained with the analysis. Error bars indicate the 68 per cent central credible interval.

  further indicates the equivalent width of [C ] λ1907+C ] λ1909 measured from the COS spectra. As noted by S17 (their fig. 10), the equivalent widths of [C ] λ1907+C ] λ1909 and [O ] λλ4959,5007 are related to each other, the most extreme C ] emitters showing largest [O ] emission. Additional high-quality rest-frame UV spectra are required to calibrate a relation between ξ ion and C ], but this relation has the potential to allow indirect constraints on ξ ion even for galaxies deep into the EoR.

Figure 4 .

 4 Figure 4. Relation between production rate of H-ionizing photons per unit intrinsic monochromatic UV luminosity (y-axis) and equivalent width of the [O ] λλ4959,5007 doublet (x-axis). Colour-coded is the equivalent width of [C ] λ1907+C ] λ1909, while grey circles indicate objects with no [C ] λ1907+C ] λ1909 detections. Error-bars on the circles indicate the 68 per cent central credible region. The black line indicates the best-fit linear relation between ξ ion and EW([O ]) reported in equation (3.1), while the grey bands indicate the 68 per cent credible region. The red hatched region indicates typical ξ ion values adopted in reionization models.

Table 1 .

 1 UV and optical emission lines fluxes computed from the HST/COS, MMT, and SDSS spectra, along with the fluxes predicted by our model. The model fluxes correspond to the posterior median of the marginal posterior distribution of each line flux. Although in our analysis we fitted the individual flux measurements, in this table, to facilitate the intepretation, we express the fluxes as log [F/F(Hα)], where F is the flux of any emission line and F(Hα) is the observed Hα flux. The complete table is available online.

	log [F/F(Hα)]

Table 2 .

 2 Continues from Table 1.

	07	-1.30

Table 4 .

 4 Example of a fitting of the galaxy SB 80, an embedded H region. Observed fluxes are indicated by grey diamonds, model fluxes by red ones. Note that in most cases the error bars on the observed fluxes are contained within the markers. Physical parameters of the ten objects analysed in this work, obtained by fitting with the emission lines listed in Table1-2. We report the posterior median and 68 per cent central credible interval of stellar age t, gas-phase metallicity 12 + log(O/H), ionization parameter log U S , dust-to-metal mass ratio ξ d and V -band attenuation optical depth τV .

	ID	t/Myr	12 + log(O/H)	log U S	ξ d	τV
	2	1.10 ± 0.1	7.87 ± 0.02	-2.56 ± 0.02	0.15 ± 0.03	0.37 ± 0.01
	36	1.11 ± 0.1	7.97 ± 0.03	-2.77 ± 0.02	0.19 ± 0.04	0.1 ± 0.02
	80	3.10 ± 0.8	8.31 ± 0.03	-2.65 ± 0.03	0.18 ± 0.03	0.17 ± 0.03
	82	2.03 ± 0.1	7.95 ± 0.01	-2.36 ± 0.02	0.28 ± 0.01	0.18 ± 0.008
	110	11.06 ± 4	8.15 ± 0.03	-2.66 ± 0.02	0.35 ± 0.03	0.04 ± 0.02
	111	1.27 ± 0.3	7.90 ± 0.03	-2.68 ± 0.02	0.18 ± 0.04	0.04 ± 0.02
	179	177.99 ± 5	8.49 ± 0.02	-2.60 ± 0.03	0.18 ± 0.03	0.01 ± 0.01
	182	1.07 ± 0.08	8.10 ± 0.03	-2.68 ± 0.02	0.16 ± 0.03	0.08 ± 0.02
	191	3.15 ± 0.2	8.47 ± 0.01	-2.14 ± 0.03	0.12 ± 0.02	0.00 ± 0.005
	198	7.05 ± 2	8.20 ± 0.08	-2.95 ± 0.06	0.12 ± 0.02	0.02 ± 0.02

Table 5 .

 5 Production rate of H-ionizing photons per unit monochromatic UV luminosity [expressed in units of log(erg -1 Hz)]. We report different definitions of this quantity, using the intrinsic stellar luminosity (ξ ion ), the stellar+nebular luminosity (ξ H ion ), and the observed (i.e. attenuated) luminosity (ξ ion , see Section 3.2). The reported values correspond to the posterior median, while the errors indicate to the 68 per cent central credible interval.

	ID	log ξ ion	log ξ H ion	log ξ ion
	2	25.84 ± 0.05	25.81 ± 0.11	26.42 ± 0.10
	36	25.37 ± 0.02	25.36 ± 0.03	25.38 ± 0.04
	80	25.73 ± 0.01	25.69 ± 0.01	26.26 ± 0.03
	82	25.80 ± 0.01	25.84 ± 0.02	25.85 ± 0.02
	110	25.24 ± 0.02	25.26 ± 0.03	25.27 ± 0.04
	111	25.38 ± 0.03	25.37 ± 0.03	25.47 ± 0.04
	179	25.50 ± 0.03	25.52 ± 0.04	25.84 ± 0.03
	182	25.57 ± 0.02	25.61 ± 0.03	25.73 ± 0.05
	191	25.82 ± 0.02	25.94 ± 0.04	26.06 ± 0.03
	198	25.65 ± 0.02	25.60 ± 0.03	25.63 ± 0.03

https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/raw_data.html MNRAS 000,1-10 (2017)

MNRAS 000, 1-10 (2017)

This parametrization translates into a volume-averaged ionization parameter U = 9/4 U S (see equation 1 of Hirschmann et al. 2017).MNRAS 000,1-10 (2017)

We do not fit the reddest SDSS band (z) since the adopted (constant) star formation histories may introduce biases when fitting a band potentially dominated by evolved stellar populations.

Note that we do not use GALEX photometry, to avoid potential biases introduced by the large GALEX PSF.

https://docs.scipy.org/doc/scipy/reference/odr.html MNRAS 000, 1-10 (2017)

This value assumes a minimum UV magnitude M 1500 = -13, while adopting M 1500 = -17 would imply a higher log( f esc ξ ion /erg -1 Hz) = 24.9.MNRAS 000,1-10 (2017)

ACKNOWLEDGEMENTS

The authors are grateful to the referee, G. Stasińska, for several comments that significantly improved the paper. JC, SC, and AVG acknowledge support from the European Research Council (ERC) via an Advanced Grant under grant agreement no. 321323-NEOGAL. DPS acknowledges support from the National Science Foundation through the grant AST-1410155. TJ acknowledges support from NASA through Hubble Fellowship grant HST-HF2-51359.001-A awarded by the Space Telescope Science Institute. AF acknowledges support from the ERC via an Advanced Grant under grant agreement no. 339659-MUSICOS. The data described in this paper were obtained by the MMT Observatory, a joint facility of the University of Arizona and the Smithsonian Institution, and the NASA/ESA Hubble Space Telescope. Support for HST GO program #14168 was provided by NASA through a grant from the

Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., diunder NASA contract NAS 5-26555.


Mahadi Ddamulira
email: mddamulira@tugraz.at

Keywords: Padovan number, Repdigit, Linear form in logarithms, Baker's method Mathematics Subject Classification 11B39, 11D45, 11D61, 11J86

Let {P n } n≥0 be the sequence of Padovan numbers defined by P 0 = 0, P 1 = 1 = P 2 , and P n+3 = P n+1 + P n for all n ≥ 0. In this paper, we find all repdigits in base 10 which can be written as a sum of three Padovan numbers.

Introduction

Let {P n } n≥0 be the sequence of Padovan numbers given by P 0 = 0, P 1 = 1, P 2 = 1, and P n+3 = P n+1 + P n for all n ≥ 0. This is sequence A000931 on the On-Line Encyclopedia of Integer Sequences (OEIS) [START_REF]The On-Line Encyclopedia of Integer Sequences[END_REF]. The first few terms of this sequence are {P n } n≥0 = 0, 1, 1, 1, 2, 2, 3, 4, 5, 7, 9, 12, 16, 21, 28, 37, 49, 65, 86, 114, 151, . . . . A repdigit is a positive integer N that has only one distinct digit when written in base 10. That is, N is of the form

N = d 10 -1 9 ( 1
)
for some positive integers d and with 1 ≤ d ≤ 9 and ≥ 2. The sequence of repdigits is sequence A010785 on the OEIS.

Main Result

In this paper, we study the problem of writing repdigits as sums of three Padovan numbers. More precisely, we completely solve the Diophantine equation:

N = P n 1 + P n 2 + P n 3 = d 10 -1 9 , ( 2
)
in non-negative integers (N , n 1 , n 2 , n 3 , d, ) with n 1 ≥ n 2 ≥ n 3 ≥ 0, ≥ 2, and 1 ≤ d ≤ 9.

We discard the situations when n 1 = 1 and n 1 = 2 and just count the solutions for n 1 = 3 since P 1 = P 2 = P 3 = 1. For the same reasons, we discard the situation when n 1 = 4 and just count the solutions for n 1 = 5 since P 4 = P 5 = 2. Thus, we always assume that n 1 , n 2 , n 3 / ∈ {1, 2, 4}. Our main result is the following. This paper serves as a continuation of the results in [START_REF] Ddamulira | Repdigits as sums of three balancing numbers[END_REF][START_REF] Lomelí | Repdigits as sums of two Padovan numbers[END_REF][START_REF] Luca | Repdigits as sums of three Fibonacci numbers[END_REF][START_REF] Luca | Repdigits as sums of three Lucas numbers[END_REF][START_REF] Luca | Repdigits as sums of four Pell numbers[END_REF][START_REF] Normenyo | Repdigits as sums of three Pell numbers[END_REF]. The method of proof involves the application of Baker's theory for linear forms in logarithms of algebraic numbers, and the Baker-Davenport reduction procedure. Computations are done with the help of a simple computer program in Mathematica.

Preliminary results

The Padovan sequence

Here, we recall some important properties of the Padovan sequence {P n } n≥0 . The characteristic equation

x 3 -x -1 = 0
has roots α, β, and γ = β, where

α = r 1 + r 2 6 , β = -(r 1 + r 2 ) + √ -3(r 1 -r 2 ) 12 , (3)
and

r 1 = 3 108 + 12 √ 69 and r 2 = 3 108 -12 √ 69. (4)
Furthermore, the Binet formula is given by

P n = aα n + bβ n + cγ n for all n ≥ 0, (5)
where

a = α + 1 (α -β)(α -γ ) , b = β + 1 (β -α)(β -γ ) , c = γ + 1 (γ -α)(γ -β) = b. (6)
Numerically, the following estimates hold:

1.32 < α < 1.33, 0.86 < |β| = |γ | = α -1 2 < 0.87, 0.72 < a < 0.73, 0.24 < |b| = |c| < 0.25. ( 7
)
From ( 3), (4), and ( 7), it is easy to see that the contribution the complex conjugate roots β and γ , to the right-hand side of (5), is very small. In particular, setting

e(n) := P n -aα n = bβ n + cγ n it is seen that |e(n)| < 1 α n/2 , ( 8
)
holds for all n ≥ 1. Furthermore, by induction, it can be proved that

α n-3 ≤ P n ≤ α n-1 holds for all n ≥ 1. ( 9
)

Linear forms in logarithms

Let η be an algebraic number of degree d with minimal primitive polynomial over the integers:

a 0 x d + a 1 x d-1 + • • • + a d = a 0 d i=1 (x -η (i) ),
where the leading coefficient a 0 is positive and the η (i) 's are the conjugates of η. Then the logarithmic height of η is given by

h(η) := 1 d log a 0 + d i=1 log max{|η (i) |, 1} .
In particular, if η = p/q is a rational number with gcd( p, q) = 1 and q > 0, then h(η) = log max{| p|, q}. The following are some of the properties of the logarithmic height function h(•), which will be used in the next sections of this paper without reference:

h(η 1 ± η 2 ) ≤ h(η 1 ) + h(η 2 ) + log 2, h(η 1 η ±1 2 ) ≤ h(η 1 ) + h(η 2 ), h(η s ) = |s|h(η) (s ∈ Z). ( 10
)
For the proofs of [START_REF] Matveev | An explicit lower bound for a homogeneous rational linear form in the logarithms of algebraic numbers II[END_REF] and further details, we refer the reader to the book of Bombieri and Gubler [START_REF] Bombieri | Heights in Diophantine Geometry[END_REF].

We recall the result of Bugeaud, Mignotte, and Siksek [2, Theorem 9.4, p. 989], which is a modified version of the result of Matveev [START_REF] Matveev | An explicit lower bound for a homogeneous rational linear form in the logarithms of algebraic numbers II[END_REF], which is one of our main tools in this paper.

Theorem 2 Let η 1 , . . . , η t be positive real algebraic numbers in a real algebraic number field K ⊂ R of degree D K , b 1 , . . . , b t be nonzero integers, and assume that

Λ := η b 1 1 • • • η b t t -1, ( 11
)
is nonzero. Then log |Λ| > -1.4 × 30 t+3 × t 4.5 × D 2 K (1 + log D K )(1 + log B)A 1 • • • A t , where B ≥ max{|b 1 |, . . . , |b t |}, and
A i ≥ max{D K h(η i ), | log η i |, 0.
16}, for all i = 1, . . . , t.

Reduction procedure

During the calculations, we get upper bounds on our variables which are too large, thus we need to reduce them. To do so, we use some results from the theory of continued fractions.

For the treatment of linear forms homogeneous in two integer variables, we use the well-known classical result in the theory of Diophantine approximation. The following lemma is the criterion of Legendre.

Lemma 1 Let τ be an irrational number, p 0 q 0 , p 1 q 1 , p 2 q 2 , . . . be all the convergents of the continued fraction expansion of τ and M be a positive integer. Let N be a nonnegative integer such that q N > M. Then, putting a(M) := max{a i : i = 0, 1, 2, . . . , N }, the inequality:

τ - r s > 1 (a(M) + 2)s 2 ,
holds for all pairs (r , s) of positive integers with 0 < s < M.

For a nonhomogeneous linear form in two integer variables, we use a slight variation of a result due to Dujella and Pethő (see [START_REF] Dujella | A generalization of a theorem of Baker and Davenport[END_REF]Lemma 5a]). For a real number X , we write X := min{|X -n| : n ∈ Z} for the distance from X to the nearest integer.

Lemma 2 Let M be a positive integer, p q be a convergent of the continued fraction expansion of the irrational number τ such that q > 6M, and A, B, μ be some real numbers with A > 0 and B > 1. Furthermore, let ε := μq -M τ q . If ε > 0, then there is no solution to the inequality:

0 < |uτ -v + μ| < AB -w , in positive integers u, v, and w with u ≤ M and w ≥ log(Aq/ε) log B .
Finally, the following Lemma is also useful. It is Lemma 7 in [START_REF] Sánchez | Linear combinations of factorials and s-units in a binary recurrence sequence[END_REF].

Lemma 3 If r 1, H > (4r 2 ) r , and H > L/(log L) r , then L < 2 r H (log H ) r .

Bounding the variables

We assume that n 1 ≥ n 2 ≥ n 3 . From ( 2) and ( 9), we have

α n 1 -3 ≤ P n 1 ≤ P n 1 + P n 2 + P n 3 = d 10 -1 9 ≤ 10 and 10 -1 ≤ d 10 -1 9 = P n 1 + P n 2 + P n 3 ≤ 3P n 1 < α n 1 +3 ,
where we use α 4 > 3. Thus

(n 1 -3) log α log 10 ≤ and -1 ≤ (n 1 + 3) log α log 10 .
Since log α/ log 10 = 0.122123 . . . < 1/5, we can conclude from the above that

(n 1 -3)/10 < < (n 1 + 8)/5. ( 12
)
Running a Mathematica program in the range 0 ≤ n 3 ≤ n 2 ≤ n 1 ≤ 500, 1 ≤ d ≤ 9, and 1 ≤ ≤ 100 we obtain only the solutions listed in Theorem 1. From now on, we assume that n 1 > 500.

By using [START_REF] Luca | Repdigits as sums of three Lucas numbers[END_REF], equation ( 2) can be written as

aα n 1 + e(n 1 ) + aα n 2 + e(n 2 ) + aα n 3 + e(n 3 ) = d 10 -1 9 . ( 13
)
We then consider (13) in three different cases as follows.

Case 1

We have that

aα n 1 + e(n 1 ) + aα n 2 + e(n 2 ) + aα n 3 + e(n 3 ) - d • 10 9 = - d 9 .
This is equivalent to

aα n 1 - d • 10 9 = - d 9 -a(α n 2 + α n 3 ) -e(n 1 ) -e(n 2 ) -e(n 3 ).
This shows that

aα n 1 - d • 10 9 ≤ d 9 + a(α n 2 + α n 3 ) + |e(n 1 )| + |e(n 2 )| + |e(n 3 )| < 1 + 2aα n 2 + 3α -n 3 /2 < 5aα n 2 ,
and so

aα n 1 - d • 10 9 < 5aα n 2 . ( 14
)
We divide through (14) by aα n 1 to get

10 α -n 1 d 9a -1 < 5α n 2 -n 1 .
Thus, we have

10 α -n 1 d 9a -1 < 5 α n 1 -n 2 . ( 15
)
We put

Λ 1 := 10 α -n 1 d 9a -1.
To apply Theorem 2, we need to check that Λ 1 = 0. Suppose that Λ 1 = 0, then we have

aα n 1 = 10 • d 9 . ( 16
)
To see that this is not true, we consider the Q-automorphism σ of the Galois extension Q(α, β) over Q given by σ (α) := β and σ (β) := α. Now, since Λ 1 = 0, we get σ (Λ 1 ) = 0. Thus, conjugating the relation (16) under σ , and taking absolute values on both sides, we get 10

• d 9 = |σ (aα n 1 )| = |b||β| n 1 < |b| < 1 3 ,
which is false for ≥ 2 and d ≥ 1. Therefore, Λ 1 = 0. Therefore, we apply Theorem 2 with the data

t := 3, η 1 := 10, η 2 := α, η 3 := d 9a , b 1 := , b 2 := -n 1 , b 3 := 1.
It is a well-known fact that

a = α(α + 1) 3α 2 -1 ,
the minimal polynomial of a is 23x 3 -23x 2 + 6x -1 and has roots a, b, c. Since |b| = |c| < |a| = a < 1 (by ( 7)), we get

h(a) = 1 3 log 23.
Since η 1 , η 2 , η 2 ∈ Q(α), we take the field K := Q(α) with degree D K := 3. Since max{1, , n 1 } ≤ n 1 , we take B := n 1 . Further, the minimal polynomial of α over Z is x 3x -1 has roots α, β, γ with 1.32 < α < 1.33 and |β| = |γ | < 1. Thus, we can take h(α) = 1 3 log α. Similarly, h(10) = log 10 . Also, By comparing the above inequality with the right-hand side of (15) we get that n 1n 2 ≤ 6.18 × 10 14 log n 1 .

h(η 3 ) ≤ h(d) + h(9) + h(a) ≤ 4 log 3 + 1
(17)

Case 2

We have that

aα n 1 + e(n 1 ) + aα n 2 + e(n 2 ) - d • 10 9 = - d 9 -aα n 3 -e(n 3 ).
This is equivalent to

a(α n 1 + α n 2 ) - d • 10 9 = - d 9 -aα n 3 -e(n 1 ) -e(n 2
)e(n 3 ).

Thus, it follows that

a(α n 1 + α n 2 ) - d • 10 9 ≤ d 9 + aα n 3 + |e(n 1 )| + |e(n 2 )| + |e(n 3 )| < 1 + aα n 3 + 3α -n 3 /2 < 3aα n 3 ,
and so

a(α n 1 + α n 2 ) - d • 10 9 < 3aα n 3 . ( 18
)
We divide through (14) by a(α n 1 + α n 2 ) to get

10 α -n 2 d 9a(1 + α n 1 -n 2 ) -1 < 3α n 3 -n 2 1 + α n 1 -n 2 .
This means that

10 α -n 2 d 9a(1 + α n 1 -n 2 ) -1 < 3 α n 2 -n 3 . ( 19
)
We put

Λ 2 := 10 α -n 2 d 9a(1 + α n 1 -n 2 ) -1.
As before, to apply Theorem 2, we need to check that Λ 2 = 0. Suppose that Λ 2 = 0, then we have

a(α n 1 + α n 2 ) = 10 • d 9 . ( 20
)
To see that this is not true, we again consider the Q-automorphism σ of the Galois extension Q(α, β) over Q given by σ (α) := β and σ (β) := α. Now, since Λ 2 = 0, we get σ (Λ 2 ) = 0. Thus, conjugating the relation (20) under σ , and taking absolute values on both sides, we get 10

• d 9 = |σ (a(α n 1 + α n 2 ))| = |b|(|β| n 1 + |β| n 2 ) < 2|b| < 2 3 ,
which is false for ≥ 2 and d ≥ 1. Therefore, Λ 2 = 0. Therefore, we apply Theorem 2 with the data

t := 3, η 1 := 10, η 2 := α, η 3 := d 9a(1 + α n 1 -n 2 ) , b 1 := , b 2 := -n 2 , b 3 := 1.
Since η 1 , η 2 , η 2 ∈ Q(α), we take the field K := Q(α) with degree D K := 3. Since max{1, , n 2 } ≤ n 1 , we take B := n 1 . Further,

h(η 3 ) ≤ h(d) + h(9) + h(a) + h(1 + α n 1 -n 2 ) ≤ 4 log 3 + 1 3 log 23 + (n 1 -n 2 ) log α + log 2
< 1.77 × 10 14 log n 1 (by (17)). > -1.98 × 10 28 (log n 1 ) 2 log α.

Thus, we can take

By comparing the above inequality with the right-hand side of (19), we get that

n 2 -n 3 ≤ 2 × 10 28 (log n 1 ) 2 . ( 21
)

Case 3

We have that

aα n 1 + e(n 1 ) + aα n 2 + e(n 2 ) + aα n 3 + e(n 3 ) - d • 10 9 = - d 9 .
This is equivalent to

a(α n 1 + α n 2 + α n 3 ) - d • 10 9 = - d 9 -e(n 1
)e(n 2 )e(n 3 ).

Thus, we have

a(α n 1 + α n 2 + α n 3 ) - d • 10 9 ≤ d 9 + |e(n 1 )| + |e(n 2 )| + |e(n 3 )| < 1 + 3α -n 3 /2 < 3,
and so

a(α n 1 + α n 2 + α n 3 ) - d • 10 9 < 3. ( 22
)
We divide through (14) by a(α

n 1 + α n 2 + α n 3 ) to get 10 α -n 3 d 9a(1 + α n 1 -n 3 + α n 2 -n 3 ) -1 < 3α -n 1 (1 + α n 2 -n 1 + α n 3 -n 1 )
.

Thus, it follows that

10 α -n 3 d 9a(1 + α n 1 -n 3 + α n 2 -n 3 ) -1 < 5 α n 1 . ( 23
)
We put

Λ 3 := 10 α -n 3 d 9a(1 + α n 1 -n 3 + α n 2 -n 3 ) -1.
As before, in order to apply Theorem 2 we need to check that Λ 3 = 0. Suppose that Λ 3 = 0, then we have

a(α n 1 + α n 2 + α n 3 ) = 10 • d 9 . ( 24
)
To see that this is not true, we again consider the Q-automorphism σ of the Galois extension Q(α, β) over Q given by σ (α) := β and σ (β) := α. Now, since Λ 3 = 0, we get σ (Λ 3 ) = 0. Thus, conjugating the relation (24) under σ , and taking absolute values on both sides, we get 10

• d 9 = |σ (a(α n 1 + α n 2 + α n 3 ))| = |b|(|β| n 1 + |β| n 2 + |β| n 3 ) < 3|b| < 1,
which is false for ≥ 2 and d ≥ 1. Therefore, Λ 2 = 0. Therefore, we apply Theorem 2 with the data:

t := 3, η 1 := 10, η 2 := α, η 3 := d 9a(1 + α n 1 -n 3 + α n 2 -n 3 ) , b 1 := , b 2 := -n 3 , b 3 := 1.
Since η 1 , η 2 , η 2 ∈ Q(α), we take the field K := Q(α) with degree D K := 3. Since max{1, , n 3 } ≤ n 1 , we take B := n 1 . Furthermore

h(η 3 ) ≤ h(d) + h(9) + h(a) + h(1 + α n 1 -n 3 + α n 2 -n 3 ) ≤ 4 log 3 + 1 3 log 23 + ((n 1 -n 3 ) + (n 2 -n 3 )) log α + 2 log 2 < 6.6 log 3 + ((n 1 -n 2 ) + 2(n 2 -n 3 )) log α < 1.
72 × 10 28 (log n 1 ) 2 (by (17) and ( 21)).

Thus, we can take A 1 := 3 log 10, A 2 := log α, and A 3 := 5.16 × 10 28 (log n 1 ) 2 . So, Theorem 2 tells us that the left-hand side of ( 23) is bounded below by

log |Λ 3 | > -1.4 × 30 6 × 3 4.5 × 3 2 (1 + log 3)(1 + log n 1 )(3 log 10) × (log α)(5.16 × 10 28 (log n 1 ) 2 ) > -1.92 × 10 42 (log n 1 ) 3 log α.
By comparing the above inequality with the right-hand side of (19), we get that

n 1 ≤ 1.94 × 10 42 (log n 1 ) 3 . ( 25
)
Now, we apply Lemma 3 on the above inequality (25) with the data: r := 3, H := 1.94 × 10 42 , and L := n 1 . We obtain that n 1 < 2.7 × 10 48 . We now record what we have proved.

Lemma 4 Let (N , n 1 , n 2 , n 3 , d, ) be the nonnegative integer solutions to the Diophantine equation (2) with n 1 ≥ n 2 ≥ n 3 ≥ 0, 1 ≤ d ≤ 9, and ≥ 2. Then, we have

< n 1 < 3 × 10 48 .

Reducing the bounds

The bounds obtained in Lemma 4 are too large to carry out meaningful computations on the computer. Thus, we need to reduce these bounds. To do so, we return to (15), (19), and (23) and apply Lemma 2 via the following procedure. First, we put

Γ 1 := log 10 -n 1 log α + log d 9a , 1 ≤ d ≤ 9.
For technical reasons, we assume that n 1 -n 2 ≥ 20 and go to (15). Note that e

Γ 1 -1 = Λ 1 = 0. Thus, Γ 1 = 0. If Γ 1 < 0, then 0 < |Γ 1 | < e |Γ 1 | -1 = |Λ 1 | < 5 α n 1 -n 2 .
If Γ 1 > 0, then we have that |e Γ 1 -1| < 1/2. Hence e Γ 1 < 2. Thus, we get that

0 < Γ 1 < e Γ 1 -1 = e Γ 1 |Λ 1 | < 10 α n 1 -n 2 .
Therefore, in both cases, we have that

0 < |Γ 1 | = log 10 -n 1 log α + log d 9a < 10 α n 1 -n 2 .
Dividing through the above inequality by log α, we get

0 < log 10 log α -n 1 + log(d/(9a)) log α < 10 α n 1 -n 2 log α . ( 26
)
If we put τ := log 10 log α and μ(d

) := log(d/(9a)) log α , 1 ≤ d ≤ 9,
we can rewrite (26) as

0 < | τ -n 1 + μ(d)| < 36 • α -(n 1 -n 2 ) . ( 27
)
We now apply Lemma 2 on (27). We put M := 3 × 10 48 . A quick computer search in Mathematica reveals that the convergent

p 106 q 106 = 177652856036642165557187989663314255133456297895465 21695574963444524513646677911090250505443859600601
of τ is such that q 106 > 6M and ε(d) ≥ 0.0129487 > 0. Therefore, with A := 36 and B := α, we calculated each value of log(36q 106 /ε(d))/ log α and found that all of them are at most 432. Thus, we have that n 1n 2 ≤ 432. In the case n 1n 2 < 20, we have n 1n 2 < 20 < 432. Thus, n 1n 2 ≤ 432 holds in both cases. Next, we put

Γ 2 := log 10 -n 2 log α + log d 9a(1 + α n 1 -n 2 ) , 1 ≤ d ≤ 9.
For technical reasons, as before we assume that n 2n 3 ≥ 20 and go to (19). Note that e Γ 2 -1 = Λ 2 = 0. Thus,

Γ 2 = 0. If Γ 2 < 0, then 0 < |Γ 2 | < e |Γ 2 | -1 = |Λ 2 | < 3 α n 2 -n 3 .
If Γ 2 > 0, then we have that |e Γ 2 -1| < 1/2. Hence e Γ 2 < 2. Thus, we get that

0 < Γ 2 < e Γ 2 -1 = e Γ 2 |Λ 2 | < 6 α n 2 -n 3 .
Therefore, in both cases, we have that

0 < |Γ 2 | = log 10 -n 1 log α + log d 9a(1 + α n 1 -n 2 ) < 6 α n 2 -n 3 .
Dividing through the above inequality by log α, we get

0 < log 10 log α -n 2 + log(d/(9a(1 + α n 1 -n 2 ))) log α < 6 α n 2 -n 3 log α . ( 28
)
We put

τ := log 10 log α and μ(d, k) := log(d/(9a(1 + α k )) log α , 1 ≤ d ≤ 9, 1 ≤ k ≤ 432,
where k := n 1n 2 . We can rewrite (28) as

0 < | τ -n 2 + μ(d, k)| < 22 • α -(n 2 -n 3 ) . ( 29
)
We now apply Lemma 2 on (29). We put M := 3 × 10 48 . A quick computer search in Mathematica reveals that the 106-th convergent of τ is such that q 106 > 6M and ε(d, k) ≥ 0.000134829 > 0 for all 1 ≤ d ≤ 9 and 1 ≤ k ≤ 432 except for the case ε [START_REF] Luca | Repdigits as sums of four Pell numbers[END_REF][START_REF] Normenyo | Repdigits as sums of three Pell numbers[END_REF], which is always negative. Therefore, with A := 22 and B := α we calculated each value of log(22q 106 /ε(d, k))/ log α and found that all of them are at most 446. Thus, we have that n 2n 3 ≤ 446.

The problem in the case of ε [START_REF] Luca | Repdigits as sums of four Pell numbers[END_REF][START_REF] Normenyo | Repdigits as sums of three Pell numbers[END_REF] is due to the fact that

1 α 9 = 2α 3 + 1 α(α + 1)(α 11 + 1) . ( 30
)
Thus, if we consider the identity (30), the inequality (28) becomes

0 < τ - (n 2 + 9) < 6 α n 2 -n 3 log α . ( 31
)
In this case, we apply the classical result from Diophantine approximation given in Lemma 1. We assume that n 2n 3 is so large that the right-hand side of the inequality in (31) is smaller than 1/(2 2 ). This certainly holds if

α n 2 -n 3 > 12 / log α. ( 32
)
Since < n 1 < 3 × 10 48 , it follows that the last inequality (32) holds provided that n 2n 3 ≥ 415, which we now assume. In this case r /s := (n 2 + 9)/ is a convergent of the continued fraction of τ := log 10/ log α and < 3 × 10 48 . We are now set to apply Lemma 1.

We write τ := [a 0 ; a 1 , a 2 , a 3 , . . .] = [8; 5, 3, 3, 1, 5, 1, 8, 4, 6, 1, 4, 1, 1, 1, 9, 1, 4, 4, 9, . . .] for the continued fraction expansion of τ and p k /q k for the k-th convergent. We get that r /s = p j /q j for some j ≤ 106. Furthermore, putting a(M) := max{a j : j = 0, 1, . . . , 106}, we get a(M) := 564. By Lemma 1, we get Finally, we put

1 566 2 = 1 (a(M) + 2) 2 ≤ τ - r s < 6 α n 2 -n 3 log
Γ 3 := log 10 -n 3 log α + log d 9a(1 + α n 1 -n 3 + α n 2 -n 3 ) , 1 ≤ d ≤ 9.
We use the assumption that n 1 > 500 and go to (23). Note that e Γ 3 -1 = Λ 3 = 0. Thus, Γ 3 = 0. If Γ 3 < 0, then

0 < |Γ 3 | < e |Γ 3 | -1 = |Λ 3 | < 5 α n 1 .
If Γ 3 > 0, then we have that |e Γ 3 -1| < 1/2. Hence e Γ 3 < 2. Thus, we get that

0 < Γ 3 < e Γ 3 -1 = e Γ 3 |Λ 3 | < 10 α n 1 .
Therefore, in both cases, we have that (34)

We now apply Lemma 2 on (34). We put M := 3 × 10 48 . A quick computer search in Mathematica reveals that the 106-th convergent of τ is such that q 106 > 6M and ε(d, k, s) ≥ 0.000125 > 0. Therefore, with A := 36 and B := α we calculated each value of log(36q 106 /ε(d, k, s))/ log α and found that all of them are at most 485. Thus, n 1 ≤ 485. This contradicts our assumption that n 1 > 500. Hence, Theorem 1 holds.

3 log 23 < 5 log 3 Thus, we can take A 1 :

 31 = 3 log 10, A 2 := log α and A 3 := 15 log 3. So, Theorem 2 tells us that the left-hand side of (15) is bounded below by log |Λ 1 | > -1.4 × 30 6 × 3 4.5 × 3 2 (1 + log 3)(1 + log n 1 )(3 log 10)(log α)(15 log 3) > -6.16 × 10 14 log n 1 log α.

A 1 :

 1 = 3 log 10, A 2 := log α and A 3 := 5.31 × 10 14 log n 1 . So, Theorem 2 tells us that the left-hand side of (19) is bounded below by log |Λ 2 | > -1.4 × 30 6 × 3 4.5 × 3 2 (1 + log 3)(1 + log n 1 )(3 log 10)(log α)(5.31 × 10 14 log n 1 )

0 < |Γ 3 |< 10 α n 1 .

 31 = log 10n 3 log α + log d 9a(1 + α n 1 -n 2 + α n 2 -n 3 )Dividing through the above inequality by log α, we get0 < log 10 log α n 3 + log(d/(9a(1 + α n 1 -n 3 + α n 2 -n 3 )) α and μ(d, k, s) := log d/(9a(1 + α k + α s )) log α , 1 ≤ d ≤ 9, where 1 ≤ k := n 1n 3 = (n 1n 2 ) + (n 2n 3 ) ≤878 and 1 ≤ s := n 2n 3 ≤ 446. We can rewrite (33) as 0 < | τn 3 + μ(d, k, s)| < 36 • α -n 1 .

  This implies that n 2n 3 ≤ 435. Thus, in both cases we have that n 2n 3 ≤ 446. In the case n 2n 3 < 20, we get that n 2n 3 < 20 < 446. Thus, n 2n 3 ≤ 446 holds in all cases.

					α	,
	which gives
	α n 2 -n 3 <	566 × 6 log α	<	566 × 6 × 3 × 10 48 log α	.

Acknowledgements Open access funding provided by Austrian Science Fund (FWF). The author thanks the referee for the careful reading of the manuscript and the useful comments and suggestions that improved on the quality of the presentation of this paper. The author was supported by the Austrian Science Fund (FWF) projects: F5510-N26-Part of the special research program (SFB), "Quasi-Monte Carlo Methods: Theory and Applications" and W1230-'Doctoral Program Discrete Mathematics". Part of the work was done during the research stay of the author at the Max Planck Institute for Mathematics in Bonn, from September to November 2019. He thanks this institution for the hospitality and the fruitful working environment.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Sarah Wassermann

Michael Seufert

Pedro Casas

Li Gang

Kuang Li

ViCrypt: Real-time, Fine-grained Prediction of Video Quality from Encrypted Streaming Traffic

published or not. The documents may come   L'archive ouverte pluridisciplinaire

▪▪▪▪▪

  Training multiple ML models over more then 4.6M individual, 1 sec. slots (5-fold cross validation) -here using all 207 inputs ▪ Classification task: per second video resolution, 6-classes: 144p, 240p, 360p, 480p, 720p, 1080p VICRYPT: Real-time, Fine-grained Prediction of Video Quality from Encrypted Streaming Traffic S. Wassermann (1), M. Seufert (2), P. Casas (1), L. Gang (3), K. Li (3) Real-time (1-sec resolution) prediction of KQIs for video streaming ▪ Video chunk detection NOT NEEDED features are packet size/time based ▪ ML models for prediction of instant, per-sec: ▪ re-buffering events ▪ video resolution ▪ video bitrate COMPUTATIONAL TIME & IMPACT OF FEATURE SELECTION ▪ per-slot re-buffering estimation errors are high, stalling slots under-estimated… ▪ …but estimation of re-buffering ratio is perfect for almost 90% of the videos ▪ Data generation through semi-controlled testbeds ▪ 15.000+ YouTube video sessions streamed and recorded in late 2018/early 2019 ▪ Different ISPs, different geographic locations (Austria, Italy, Germany, China) ▪ Home/corporate WiFi networks, LTE networks ▪ QUIC and TCP sessions ▪ Bandwidth limitations: 20Mbps, 5Mbps, 3Mbps, 1Mbps, 300kbps + fluctuations ▪ JavaScript-based monitoring script to measure ground truth at the player ▪ Video stream-based analysis, using multiple sliding windows, capturing different temporal phenomena (current time, short-term trend, session-aggregated) ▪ Analysis is done in real time: for every video session and for every new time slot of 1 second, we consider the following set of 207 features: ▪ Features extracted from current time slot (C) -69 features ▪ Short-memory (trend) based features, extracted from last T (3) slots (CT) -69 features ▪ Cumulative based features, extracted from all past traffic for this video session (CS) -69 features ▪ Feature computation is done continually, in constant-memory boundaries, using sketches ▪ Automatic Feature Selection -CS features (F S ) alone provide the best results (69 features), improving overall performance. Top 20 features (F TOP20 ) provide similar improvement with much less features Regression task: estimation of per second average video bitrate ▪ ERT10 & BAGGING realize MAE below 100kbps, and RMSE below 190kbps (penalizes larger errors) ▪ 80% of the slots are estimated with errors below 100kbps Predictions are highly correlated with the target (PLCC = 0.93) Binary classification task: playback stalled/not-stalled at every new slot re-buffering ratio reconstruction error PDF ▪ laptop (i5 CPU, 8GB RAM) vs. server (Xeon Silver, 48 cores, 128GB RAM) ▪ server: avg. duration of full feature set update is 13 μs, prediction time below 1.4ms ▪ Laptop: avg. feature update takes 37 μs, prediction time below 16ms ▪ Evaluation of full feature set update time (done at every new incoming packet) and prediction time (done for every 1s slot), using an upper bound with all 207 features