First Light And Reionisation Epoch Simulations (FLARES) XII: The consequences of star-dust geometry on galaxies in the EoR

Using the First Light And Reionisation Epoch Simulations (${\rm F{\small LARES}}$), a suite of hydrodynamical simulations we explore the consequences of a realistic model for star--dust geometry on the observed properties of galaxies. We find that the UV attenuation declines rapidly from the central regions of galaxies, and bright galaxies have spatially extended star formation that suffers less obscuration than their fainter counterparts, demonstrating a non-linear relationship between the UV luminosity and the UV attenuation, giving a double power-law shape to the UVLF. Spatially distinct stellar populations within galaxies experience a wide range of dust attenuation due to variations in the dust optical depth along their line-of-sight; which can range from completely dust obscured to being fully unobscured. The overall attenuation curve of a galaxy is then a complex combination of various lines-of-sight within the galaxy. We explore the manifestation of this effect to study the reliability of line ratios to infer galaxy properties, in particular the Balmer decrement and the BPT diagram. We find the Balmer decrement predicted Balmer line attenuation to be higher (factor of $1$ to $\gtrsim10$) than expected from commonly used attenuation curves. The observed BPT line ratios deviate from their intrinsic values (median difference of 0.08 (0.02) and standard deviation of 0.2 (0.05) for log$_{10}$([N${\rm {\small II}}$]$\lambda 6585/$H$_{\alpha}$) (log$_{10}$([O${\rm {\small III}}$]$\lambda 5008/$H$_{\beta}$)). Finally, we explore the variation in observed properties (UV attenuation, UV slope and Balmer decrement) with viewing angle, finding average differences of $\sim0.3$ magnitudes in the UV attenuation.


INTRODUCTION
Attenuation by dust is ubiquitous when performing any astronomical observations in the rest-frame ultraviolet (UV) to the near-infrared (NIR), since dust is one of the key components of the interstellar medium (ISM) of galaxies.Dust modifies the intrinsic emission of galaxies, by partly absorbing light from the UV to the NIR, and reemitting at longer wavelengths.Even though dust only constitutes a small fraction of the baryonic mass of a galaxy, it is estimated that nearly half of the photons in the Universe have been reprocessed by dust (Dwek et al. 1998;Bernstein et al. 2002;Driver et al. 2016).In order to derive the physical properties of galaxies, such as their stellar mass, gas content, star formation rate, age and metallicity, it is important to separate the impact of dust from the intrinsic stellar and nebular emission (see reviews by Walcher et al. 2011;Conroy 2013).
In the current observational astronomy paradigm this decoupling ⋆ E-mail: apavi@space.dtu.dk is usually achieved through the use of dust attenuation curves/laws1 , which parameterise the wavelength dependence of the dust optical depth.These curves or laws are a function of the dust grain properties, such as their size, shape and composition, as well as the relative geometry of stars and dust (star-dust geometry) in the galaxy (see review by Salim & Narayanan 2020).Dust attenuation includes the aggregate effect of dust absorbing stellar and nebular emission along the line-of-sight, as well as the influence of unobstructed starlight that has not undergone any processing.Additionally, it encompasses the scattering of light into and out of the line-of-sight.Various works have measured dust attenuation curves or laws in the local Universe, such as those observed in the Milky Way (MW, Fitzpatrick & Massa 1990) or the Magellanic Clouds (LMC and SMC, Pei 1992).Generally, these are observed extinction laws, which corresponds to the attenuation law through a distant uniform screen geometry.This geometry is unrealistic for real galaxies; despite this, these curves are still regularly used to correct for dust attenuation of the spectral energy distribution (SED) of distant galaxies.
Real galaxies exhibit complex geometries, where the assumption of a uniform dust screen is insufficient.Thus, there is a concerted effort to obtain average / effective attenuation laws, that encompass the impact of different dust grain properties and compositions (silicates, carbonates, etc), along with the relative distribution of stars and dust within galaxies.These two aspects are hard to disentangle when dealing with unresolved observations, or when lacking multiwavelength coverage.One such average dust attenuation law, widely used in the literature, is known as the Calzetti-starburst attenuation curve.This curve has been measured for bright star forming galaxies (Calzetti et al. 1994;Calzetti 1997), and it is flatter than the SMC dust law.
Theoretical studies have also focused on geometries other than a simple uniform screen by taking into account the clumpiness of the dust and stellar distribution (e.g.Witt et al. 1992;Witt & Gordon 2000;Inoue 2005;Chevallard et al. 2013).Such a scenario leads to a more transparent medium than a uniform screen, due to stellar light escaping through dust free paths.These studies have also demonstrated that the attenuation law through a clumpy medium with a SMC type dust composition can be very similar to the Calzetti type dust law, which is significantly greyer at shorter wavelengths.Thus, galaxies which have intrinsically different dust properties can exhibit attenuation curves which are similar.
When modelling galaxy SEDs using various fitting methods, attenuation or extinction curves are commonly used to correct for the effect of dust.To simplify the fitting process, these curves are generally used in conjunction with the assumption that stars sit behind a uniform screen of dust (with the choice of the attenuation curve describing the effective attenuation through the screen, or how porous it is) with higher dust optical depths and different grain properties towards young star forming regions (e.g.Charlot & Fall 2000).Therefore, SED fitting already assumes a particular geometry between stars and dust, as well as dust properties, a priori when any particular dust curve is chosen.
Recent years have seen a significant increase in the number of detected sources in the high-redshift Universe (z ≥ 5).This has been mainly due to multi-wavelength observations taken through numerous ground and space based observatories, which has provided significant insight into the UV-optical emission of these galaxies.However, there are many obstacles when constructing a comprehensive picture of galaxy formation and evolution in the early Universe.This has been primarily due to the limitations of ground based observatories in terms of their accessible wavelength range for observing high-z galaxy spectral features (e.g.Lyman or Balmer break redshifting out of the observed wavelength, or not having rest-frame UV, optical and near-IR coverage to simultaneously constrain stellar masses and star formation rates).In addition, space based observatories are, generally, relatively smaller compared to their ground based counterparts, leading to poor spatial resolution as well as sensitivity.
As discussed, another obstacle is the presence of dust in galaxies.Even though the dust content of high-redshift galaxies is expected to be smaller than in the local Universe or at cosmic noon (z ∼ 2), recent studies of the high-redshift Universe have revealed a number of sources with highly dust-obscured star formation (e.g.Fudamoto et al. 2021;Ferrara et al. 2022;Algera et al. 2023).Understanding the physical properties of high-z galaxies is also complicated by the lack of spatially and spectrally resolved multi-wavelength observations, that can provide insight into the star-dust geometry within these galaxies.Thus, many works rely on the use of scaling relations (e.g. the [CII]λ 158 µm luminosity -SFR relation, De Looze et al. 2011) or proxies (e.g. the infrared excess -UV continuum slope, referred to as the IRX-β relation, Smit et al. 2012) to derive the phys-ical properties of galaxies.However, many of these relations are obtained from observations of the local Universe; one would expect that the physical conditions within high-redshift galaxies will differ substantially in the hot early Universe.Finally, observations of these high-z systems have shown that they are extremely clumpy (e.g.Chen et al. 2023), and hence generalising their properties from unresolved photometry will lead to incorrect inferences of their physical properties, such as stellar mass, SFR or dust attenuation (e.g.Giménez-Arteaga et al. 2023).
A comprehensive way to explore the impact of star-dust geometry is through the use of hydrodynamical simulations of galaxy formation and evolution.These simulations typically associate dust with the distribution of metals in the gaseous medium.The last decade has seen the development of many calibrated physical models run over large periodic volumes to create representative galaxy samples [e.g.EAGLE (Schaye et al. 2015;Crain et al. 2015), ILLUSTRIS-TNG (Marinacci et al. 2018;Naiman et al. 2018;Nelson et al. 2018;Springel et al. 2018;Pillepich et al. 2018), SIMBA (Davé et al. 2019), etc].With sophisticated forward modelling techniques these simulations can now be compared to observed galaxies in the same observational space (e.g.fluxes, luminosities or line equivalent widths).These advancements have led to many theoretical studies using simulations to understand the effect of dust on galaxy properties as well as the star-dust geometry (e.g.Narayanan et al. 2018;Wilkins et al. 2018;Trayford et al. 2020;Shen et al. 2020;Lovell et al. 2022;Lower et al. 2022;Hsu et al. 2023).They can help us explore and understand how galaxies -which are complex objects with multiple spatially separated stellar populations of different ages (star formation histories) and metallicities (metal enrichment histories) -contribute differently to the observed fluxes.
The recent launch of JWST has ameliorated many of the drawbacks of space based telescopes for high-z extragalactic astronomy, due to its size, as well as its access to the observed frame near-and mid-infrared.The analysis of the data has already revealed interesting features about the high-redshift Universe (e.g.Castellano et al. 2022;Ferreira et al. 2022;Naidu et al. 2022;Adams et al. 2023;Donnan et al. 2023;Harikane et al. 2023), with some studies detecting massive galaxies with very high inferred values of dust obscuration as early as z ∼ 10 (e.g.Rodighiero et al. 2023).The spatial resolution of JWST in the 1 − 5 µm range also provides a new window for exploring resolved properties of high-z galaxies (z ≥ 5) in the rest-frame UV-optical for the first time without gravitational lensing (e.g.Pérez-González et al. 2023).This will allow us to probe the distribution of different stellar populations, as well as the intervening dust content.In order to fully comprehend these observations, it is imperative to construct suitable theoretical frameworks that facilitate their interpretation.In this context, it is timely to examine the variations in the observed UV-optical fluxes / luminosities resulting from the interplay between dust and the associated star-dust geometry.This work will complement other theoretical works in this space, with a focus mainly on the brightest and most massive galaxies in the high-redshift Universe.
In this work, we explore the consequences that a realistic model for the star-dust geometry has on the observed properties of galaxies in the Epoch of Reionisation (EoR, z ≥ 5) using the First Light And Reionisation Epoch Simulations (FLARES, Lovell et al. 2021;Vijayan et al. 2021).FLARES re-simulates a wide range of overdensities in the EoR, utilising the EAGLE (Schaye et al. 2015;Crain et al. 2015) simulation physics.By simulating a very wide distribution of overdensities in the EoR, FLARES allows us to study the rare and bright objects, normally unavailable in periodic box simulations, due to their limited volume.This strategy also gives us access to a statistical sample of EoR galaxies (see Lovell et al. 2021).
This paper is structured as follows.In Section 2 we introduce the suite of simulations that we are using, our galaxy sample and explain briefly our spectral energy distribution (SED) modelling.We also introduce our methods to calculate the half-mass radii and observed SED for various viewing angles.In § 3 we explore the variation and consequence of the dust distribution within the simulated galaxies.In § 4, we also create toy galaxies to guide the interpretation of the impact of the variation of dust distribution in FLARES galaxies, exploring its effect on the attenuation curve in § 4 as well as on observed line ratios such as the Balmer decrement in § 5.1 and the BPT diagram in § 5.2.In § 5.3 we look at the effect of various viewing angle on a few observed properties.We finally summarise our findings and present our conclusions in § 6.

SIMULATIONS AND MODELLING
FLARES is a zoom-in hydrodynamical simulation suite of 40 regions, designed to study the evolution of galaxies hosted in different environments in the early Universe, ranging from extreme overdensities to extreme underdensities.These regions were chosen from a (3.2 comoving Gpc) 3 dark matter only (DMO) box, resimulated with full hydrodynamics using the AGNdT9 configuration of the EAGLE (Schaye et al. 2015;Crain et al. 2015) simulation physics to produce the galaxy population at z ≥ 5 with same resolution and cosmology as the reference EAGLE simulation (the physical model used to run the largest EAGLE volume (100 cMpc) 3 , see Table 3 in Schaye et al. 2015, for the differences in the used parameters of the two models).The AGNdT9 configuration produces similar mass functions to the Reference model but better reproduces the hot gas properties of groups and clusters (Barnes et al. 2017), with modifications that give less frequent, more energetic AGN feedback events.The model has been shown to be very successful in reproducing a number of observables in the low-redshift Universe not used in the calibration.The FLARES philosophy is to test the EAGLE physics model at high-redshift as well as forward model the simulated galaxies into the observed space.For this, we have deliberately selected extreme overdense regions (16), to obtain large sample of massive galaxies observable by telescopes in the high-redshift Universe.The properties of the high-redshift galaxies in FLARES such as the UVcontinuum slope evolution, galaxy stellar ages, metallicity have been shown to match the current observational constraints in a number of works (e.g.Vijayan et al. 2021;Tacchella et al. 2022;Wilkins et al. 2022).The simulations have also been used to provide predictions for future observations in the high-redshift Universe (e.g.Roper et al. 2022;D'Silva et al. 2023;Lovell et al. 2023;Wilkins et al. 2023a,b).
To obtain representative distribution functions, we apply an averaging scheme that weights the different regions based on their overall abundance in the parent box, with the very overdense and underdense regions contributing the least.For more details we refer the interested readers to Lovell et al. (2021) where the region selection and weighting strategy are described in Sections 2.2 and 2.4, respectively.The main driver in using this suite of simulations is that it provides a statistically significant sample of galaxies in the EoR and thus enables us to better explore the evolution and properties of high-z galaxies.

Galaxy Identification
Galaxies in FLARES, similar to the fiducial EAGLE model, are identified with the SUBFIND (Springel et al. 2001;Dolag et al. 2009) algorithm, which runs on bound groups found by a Friends-Of-Friends (FOF, Davis et al. 1985) algorithm.The integrated galaxy properties like stellar mass, SFR or luminosities are calculated within a 30 physical kpc (pkpc) 3D aperture centered on the most bound particle of the self-bound structures, similar to previous EAGLE and FLARES works.For the redshifts of the FLARES galaxies, this aperture ensures that we capture the majority of the mass/light belonging to each galaxy and remove any contribution from spurious distributions at larger radii.In this work, the quoted star formation rate (SFR) are measured using the mass of stars formed in the last 100 Myr.For the purpose of properly resolving and analysing the variation of different properties within galaxies, we assign a minimum threshold of 500 star particles for the galaxies in this study (corresponding to ∼ 10 8.8 M ⊙ in stellar mass).
The left panel of Figure 1 shows the sample of galaxies used in this study in the redshift range of 5 − 10.The figure plots the stellar mass as a function of the SFR of the galaxy.It can be seen that the sample extends over a stellar mass range of roughly 10 8.8 −10 11 M ⊙ , with the SFRs ranging from 0.1 − 1000 M ⊙ yr −1 .

Spectral Energy Distribution modelling
We follow the same prescription adopted in Vijayan et al. (2021), where we modelled the SED of each star particle in a galaxy, by treating them as a simple stellar population (SSP) based on their age and metallicity (the smoothed metallicity of the star i.e.SPH kernel weighted metallicity, see The EAGLE team 2017).The emission corresponding to a stellar population is then modelled using the the Binary Population and Spectral Synthesis (BPASS, V2.2.1) SPS (stellar population synthesis) models (Stanway & Eldridge 2018).Throughout this work, we assume a Chabrier initial mass function (IMF, Chabrier 2003).We also associate each young stellar particle (≲ 10 Myr) with a surrounding HII region powered by its LyC emission using the CLOUDY (V17.02)photoionisation code (Ferland et al. 2017); same as the photo-ionisation model described in Wilkins et al. (2020).
Since our simulations do not model the formation and destruction of dust, the dust attenuation is modelled by converting the line-ofsight metal column density along the z-axis of the simulation (the adopted viewing angle, thus giving a random orientation for each galaxy) to a dust column density.This conversion is done by adopting a dust-to-metal ratio (using the fitting function for the dust-tometal (DTM) ratio from Vijayan et al. 2019, equation 15) for the galaxy, implemented as follows: where τ ISM,V (x, y) and Σ (x, y) are the V-band (550nm) optical-depth and integrated metal column density respectively of the intervening diffuse ISM along the line-of-sight at position (x,y) of the star particle.κ ISM = 0.0795, is a proportionality constant, chosen to match the rest-frame far-UV (1500Å) luminosity function to the observed UV luminosity function from Bouwens et al. (2015) at z = 5.We also assume a contribution to the dust attenuation from the stellar birth clouds, which dissipate over a timescale of 10 Myr (Charlot & Fall 2000).For the stellar birth clouds, we scale their dust content with the smoothed star particle metallicity, described as The coloured solid lines represent the weighted median of the relation at different redshifts, along with the 84 th − 16 th percentile spread.Also shown is the histogram of the distribution of stellar mass and SFR in different bins for these redshifts.The total number of galaxies at these redshifts are indicated within brackets alongside the legend, with the colourbar indicating the number of galaxies within the hexbins.
follows, τ BC,V (x, y) = κ BC (Z ⋆ /0.01) t ≤ 10 7 yr 0 t > 10 7 yr , where τ BC,V (x, y) is the V-band optical-depth due to the stellar birthcloud and Z ⋆ is the smoothed metallicity of the young stellar particle.κ BC = 1.0, similar to κ ISM , is a normalisation factor, chosen to to match the UV-continuum slope observations from Bouwens et al. (2012Bouwens et al. ( , 2014) ) at z = 5 and the [OIII]λ 4959, 5007 + Hβ equivalent width distribution at z = 8 from De Barros et al. (2019).The value also encapsulates the dust-to-metal ratio in the stellar birth clouds.
By fixing the values of κ ISM and κ BC for all redshifts, we assume there is no evolution in the general properties of the dust grains in galaxies such as the average grain size, shape, and composition.We link the optical-depth in the V-band to other wavelengths using a simple powerlaw relation: (3) This produces an extinction curve for a stellar particle flatter in the UV than the Small Magellanic Cloud curve (Pei 1992), but not as flat as the Calzetti et al. (2000) curve.We refer the interested readers to Section 2.3, 2.4 and Appendix A of Vijayan et al. (2021) for an in-depth discussion of the photometry generation in FLARES.
Our modelling differs from using a uniform screen of dust across a galaxy (usually assumed in SED fitting of observed galaxies), with spatially distinct stellar populations experiencing differing amounts of dust attenuation, thus a more realistic model for dust.
The right panel of Figure 1 shows the UV luminosity of the galaxies obtained using the described forward modelling, plotted against the galaxy SFR for z ∈ [5, 10].From the histogram of the distribution of the UV luminosity, it can be seen that the sample contains more number of brighter (< −22 Mag) and fainter (> −19 Mag) galaxies as one goes towards lower redshifts, similar to the SFR distribution.The lower SFR cut at z = 10 is due to our star particle selection limit.The galaxies with similar SFR have different UV luminosities due to the variation in the amount of dust.Moving down in redshift, the distribution extends towards lower luminosities (> −19 Mag) due to the build-up of dust as well as galaxies with older stellar populations, while the presence of very massive young star forming galaxies extends the distribution at the bright end (< −22 Mag).

Half-mass radii
Hydrodynamical simulations like FLARES provide information on the spatial distribution of stars and gas within a galaxy.We can utilise this to calculate the mass distribution within the galaxy and thus derive the half-mass radii for either of these components.To calculate the half-mass radii we compute the distance of each star/gas particle within a galaxy from its centre of potential.By ordering them in increasing distance one can find the radius that encloses half the mass in that component.The half-mass radii will be used to normalise the radial variation in observed galaxy properties in § 3.
Galaxy sizes in FLARES have been extensively discussed in Roper et al. (2022Roper et al. ( , 2023)); who find that galaxies at the massive end in the early Universe have most of their star formation localised to extremely compact central cores.A main reason for these compact cores are due to the gas densities required for star formation in the early Universe being very high, because of the very low gas phase metallicities.The EAGLE model also defines a gas critical density at which a star particle can be formed, which depends inversely on the gas phase metallicity, thus driving the critical density higher in the early Universe.And once star formation is triggered in these regions, they become quickly enriched by metals, starting a process of runaway star formation.Thus these compact cores are very metal rich and suffer from extreme dust attenuation.This causes these galaxies to have observed (i.e.dust attenuated) UV sizes that are significantly larger than the intrinsic (i.e.no obscuration by dust) UV sizes.These compact systems evolve to larger intrinsic as well as observed UV sizes due to the surrounding gas becoming enriched at later times.For more details about the size evolution of the FLARES galaxies, interested readers are directed to the referred works.

Re-projection
We use the Hierarchical Equal Area isoLatitude Pixelisation (HEALPix 2 ) scheme as implemented in healpy 3 to pixelize the lineof-sight angles.We assume N side = 8 yielding 768 individual linesof-sight.We re-project the star and gas particles of the galaxy along these lines-of-sight, with the viewing angle along the z-axis as before.We then calculate the dust-attenuated SED.This allows us to explore the observed properties of galaxies as a function of viewing angle (e.g.see Figure 11).In § 5.3 we use these re-projected values to quantify the effect of different viewing angles on the observed galaxy properties.

VARIATION OF UV ATTENUATION WITHIN GALAXIES
In this section we explore how dust attenuation varies within individual galaxies.While galaxies can be characterised by a global attenuation, it is clear that they are complex objects with different regions experiencing varying degrees of dust attenuation.The lineof-sight dust attenuation model implemented in FLARES allows us to explore this variation.We begin in § 3.1 by exploring how dust attenuation and the star formation traced by the UV emission (i.e. the unobscured SFR, denoted by SFR UV ) correlates with distance from the galaxy centre and its implications.From here onwards, to make the distinction from SFR UV , we will denote the total SFR of the galaxy as SFR Total , which includes contribution from both obscured and unobscured star formation.We then explore in § 3.2 the spread in the degree of attenuation within individual galaxies of different star particles.

Radial variation and the luminosity function
Using the half-mass radii calculated in §2.3 we describe how the average UV attenuation within galaxies varies with distance from the centre.In Figure 2, we plot the median variation with radius of the UV star formation rate surface density (Σ SFR,UV ), and the attenuation in the far-UV (in magnitudes, A FUV ), for galaxies at z ∈ [5, 10].We compute SFR UV in a galaxy by defining the observedto-intrinsic fraction of UV emission, f esc , using 2 http://healpix.sourceforge.net/ 3https://github.com/healpy/healpy/ We can then define where r o and r i are the radii of the outer and inner bins (in kpc) respectively, considered for calculating the surface area.We reiterate here that f esc , in the context of this work, does not refer to the escape fraction of ionising photons from HII regions, but it refers to the fraction of UV emission that is observed or unobscured.In Figure 2, we have normalised the radial bins with respect to the stellar half mass radii (r 1/2,⋆ ).
From Figure 2, it can be seen that the attenuation of the stellar light is highest in the centre of the galaxies, rapidly decaying towards the outskirts.This is not surprising considering that stars are formed first in the dense central cores, with the subsequent generation of stars forming from gas enriched by this population.This inside-out growth of star formation within the galaxy builds up higher dust surface densities within the inner parts of the galaxy, leading to higher dust attenuation in the centres.At z = 10, the median radial variation of Σ SFR,UV for different galaxies with different observed FUV luminosities are almost the same.This is due to the galaxies in the most luminous bin (≤ −20 Magnitude) having higher UV obscuration in the centre.The galaxies on the fainter side at z = 10 show almost uniform dust obscuration within the inner 2×r 1/2,⋆ .This is also accompanied by slightly lower Σ SFR,UV in the inner regions.Moving to lower redshifts, the brightest luminosity bins (≤ −22 Magnitude) show significant amount of dust obscuration within the central region4 .At z = 5 the brightest galaxies (> −21 Magnitude) are still intensely star-forming with most of the activity being highly-obscured in the centre.This is due to highly dust enriched central regions in these galaxies (see Roper et al. 2022).
However, the radial decline in the amount of UV star formation is not as dramatic as the UV attenuation, with the figure clearly indicating there is a significant amount of UV luminosity that is coming from the outskirts, especially for the brightest galaxies.This effect is also reflected in the UV size evolution of galaxies across redshift, with the intrinsic sizes as much as 50 times smaller than the observed sizes (see Roper et al. 2022Roper et al. , 2023)).In short, galaxies at the bright end of the UVLF have started forming stars away from the central dust rich core towards the outskirts.However, since the gas in the outskirts is less enriched, there is thus less dust attenuation.This phenomenon makes these galaxies brighter than expected compared to the galaxies in the lower luminosity bins.Thus, in FLARES, for the bright galaxies we are seeing a non-linear relationship (in the magnitude space) between the UV luminosity and the UV attenuation.
To probe the effect of this extended star formation and the resultant lower attenuation in these galaxies have on the shape of the UV LF at the bright end, we perform the following simple exercise.We convert the intrinsic UV LF of the FLARES galaxies to the dust attenuated value using a linear scaling relation between the two, which is commonly used in literature (e.g.Smit et al. 2012;Tacchella et al. 2013).The scaling relation assumes that there is a one-to-one relation between the intrinsic UV luminosity and dust attenuation, which means that it would fail to properly capture the effect of unobscured star formation in the outskirts.This is done by following the parameterisations for the mean UV-continuum slope and the UV attenua- . The intrinsic (dotted blue) UV LF and the observed UV LF for the line-of-sight model (solid green) and the one obtained using equation 9 (dashed orange) for the FLARES galaxies in z ∈ [5, 10] are shown.The shaded region denotes the 1 − σ Poisson errors on the numbers, and is only shown for bins with ≥ 5 galaxies.Note that the star particle cut mentioned in § 2.1 is not enforced here to preserve the full shape of the UV LF. tion: which feeds into and where M UV,int is the intrinsic UV luminosity, M UV,att the dust attenuated UV luminosity.M 0 = −19.5,σ β = 0.34 and the values of dβ dM UV (z) and β M 0 (z) are taken from Table 1 of Tacchella et al. ( 2013) 5 .These are in turn from Bouwens et al. (2012); Bradley et al. (2012).
Figure 3 shows the intrinsic (dotted blue line) UV LF, with the dashed orange line showing the attenuated UV LF using Equation 9(linear model, in figure).We also plot the observed UV luminosity using the line-of-sight dust attenuation model (the default model for the FLARES galaxies described in §2.2; LOS model in figure) as the solid green line.As can be seen from the plot, the two models produce similar number densities for the fainter galaxies, while from z ≤ 9, the galaxy number densities are extended to higher luminosities compared to the scaled model.This extension to higher luminosities is not strongly visible as we go to higher redshifts due to FLARES missing the extreme bright galaxies.However, it is clear that across redshifts, the brightest galaxies have less global dust attenuation in the LOS model, and thus the number densities have been pushed to higher values6 than the scaled model.This is due to the latter's simpler parameterisation of the relationship between intrinsic UV luminosity and dust attenuation, i.e. a linear relation between attenuation, β slope and UV luminosity.For the FLARES galaxies, the relationship between these values deviate from a linear relationship and becomes flat at the bright end.This phenomenon is partly responsible for the FLARES UV LF exhibiting the double power-law (DPL) shape measured in (Vijayan et al. 2021): the double Schechter shape of the intrinsic UV luminosity function in FLARES, when attenuated (or observed), becomes a DPL.
An implication of the above is that the use of a linear scaling between the UV attenuation and luminosity will over-predict (since there is agreement at the faint-end) the UV attenuation at the bright end of the UV luminosity function.Thus the estimated intrinsic UV luminosity will be higher.Also the predicted stellar mass and star formation rates will be wrong due to widely varying mass-to-light ratios across the galaxy.Thus it is important to take into account such biases when interpreting derived intrinsic values from observed UV luminosity, when the derived values are estimated with the use of UV-to-SFR conversion factors or constant mass-to-light ratios.

Variation across stellar populations
In this subsection we will explore how much scatter there is in the dust attenuation seen by different star-forming regions (or particles in the simulation).Figure 4  Figure 4 shows this scatter for different galaxies in the redshift range [5, 10] as a function of the observed UV luminosity, with the hexbins coloured by the median SFR Total within that bin.We also plot alongside the probability distribution function (PDF, normalised by the bin width) of f esc,84 − f esc,16 for our galaxy sample.The figure shows that there is a huge spread within most galaxies at z = 10 (see PDF), with the number of galaxies with small spread in the the f esc values increasing with decreasing redshift.At all redshifts, in FLARES there is a large spread in f esc,84 − f esc,16 that is correlated neither with the galaxy UV luminosity (as can be seen directly from Figure 4) nor the stellar mass.We have also looked at the contribution to this scatter due to young stars still embedded in their birth clouds, and found it to be negligible.From the hexbin colours, it can be seen that there is a weak correlation with the SFR Total of the galaxy at z = 10, with the correlation increasing with decreasing redshift.I f we concentrate on the highly star forming galaxies in FLARES, by splitting them using the median SFR Total , it can be seen that they predominantly occupy the upper part of f esc,84 − f esc,16 , while the the lower half consists of less star forming galaxies.A similar pattern is also seen if one were to replace SFR Total with UV attenuation, however for the rest of this work, we use SFR Total to quantify this spread.
This observed spread is mostly a result of the complex star-dust geometry within the galaxy, thus having mix of stars with highly obscured and unobscured sightlines.The observed spread is seen to be correlated with the galaxy SFR in FLARES.The correlation implies that highly star forming galaxies have more regions that have sightlines suffering different levels of dust attenuation.This is not unexpected, since we also saw in §3.1 that extreme UV-bright galaxies have star formation that is less obscured, due to stars forming outside the dust-rich core.This result is robust in the mass range we are investigating, due to the conservative particle number cut employed in this study, which implies that it is highly unlikely that the dust distribution is poorly sampled (creating holes within the distribution).
Another interesting feature we observe is that with decreasing redshift there are more galaxies that start to populate the region with f esc,84 − f esc,16 ∼ 0. This is indicative of galaxies with a more homogeneous dust distribution, or with very low levels of dust attenuation (the spread can still be significant if the attenuation is very low).In these cases, the application of a homogeneous screen model will be suitable.There are two reasons for the surge in the number of such galaxies.One of the reason for more galaxies is that in general as one goes to lower redshift, we have more galaxies in the simulation fulfilling our selection criteria.Secondly, galaxies at higher redshift are expected to have undergone a more bursty phase of star formation (Wilkins et al. 2023a), and thus more likely to exhibit significant variation in the dust attenuation across the galaxy (due to the seen corelation with SFR).However, as these galaxies move to lower redshifts, their star formation rates are expected to gradually decline and undergo fewer frequent bursts.Hence they are expected to have more homogeneous distribution of dust along the line-of-sight (also see § 4).
Recently, Lower et al. (2022) implemented an additional parameter in SED fitting, f unobscured , to take into account contribution from stars which lie across sightlines which suffer from negligible dustattenuation.The model was implemented in the SED fitting code Prospector (Johnson et al. 2021) and was used to fit the observed Figure 4. Shows the variation of the dust obscuration in the UV within galaxies as a function of the observed UV luminosity.The spread is described as the difference between the 84th and 16th percentile of the UV attenuation of stars within the galaxy.The hexbins have been coloured by the median SFR Total of the galaxies in the bin.The plotted distribution function alongside shows the normalised counts of galaxies in different bins.The dashed line is the distribution of fraction of galaxies with their star formation rate higher than the median value at that redshift, with the dotted line showing the fraction of galaxies below the median.
fluxes (19 bands, including GALEX, HST, Spitzer, and Herschel) of mock galaxies at z = 0 from the SIMBA simulations generated using the radiative transfer code Powderday (Narayanan et al. 2021).The results from the study showed that incorporating an extra parameter to allow for unobscured sightlines provide better constraints on derived values such as stellar mass and star formation rate.Their results are in agreement with ours, showing the value of incorporating an unobscured parameter in SED fitting.

ATTENUATION CURVE
The observed spectra of galaxies are a result of different contributions from the stars of different ages and metallicities that are inhomogeneously distributed in a clumpy soup of dust.A number of steps are required to model the intrinsic spectrum of galaxies (see Walcher et al. 2011;Conroy 2013, for a detailed review).Dust attenuation curves describe how the intrinsic emission (from stellar and nebular emission) from a galaxy is modified by dust at different wavelengths.An understanding of these attenuation curves is essential for accurately determining physical properties of galaxies, through SED fitting of the observed galaxy photometry or spectrum.Therefore, they play a pivotal role in galaxy formation and evolution studies.
We already saw in the previous section that the galaxies in FLARES have a wide distribution of dust attenuation across different star particles.In order to understand the effect this variation has on the attenuation curve of galaxies as well as on the conclusions drawn from the observed SEDs, we create a controlled experiment.This will help us better understand the observed behaviour exhibited by the FLARES galaxies in this and the next section.For this we build toy galaxies at z = 5 in the SED fitting code Bagpipes (Carnall et al. 2018) with the 2016 updated version of BC03 SPS (Bruzual & Charlot 2003) library, modelling it as comprising of 50 star forming clumps.The resulting galaxy has a composite SED obtained by summing up the SEDs of the individual clumps.To generate the SEDs for these star forming regions we assume a single burst star formation history (SFH) with the same stellar age (5 Myr), metallicity (0.1 Z ⊙ ), mass of stars formed (10 7 M ⊙ ) and logU (ionisation parameter, -3).We assume a SMC type dust law, and change the line-of-sight dust attenuation by varying the attenuation in the V-band (A V magnitude) for the clumps.Thus the intrinsic spectrum of the galaxies are the same, however, the observed spectrum is different.It should be noted that in general SED fitting codes assume the same extinction or attenuation (commonly referred to as a dust screen, with the choice of curve determining how porous the screen is) for all the different stellar components (with multiple ages and metallicities), with extra attenuation only for young stars which are still embedded within their birth clouds.We also produce these clumps with no added contribution from nebular emission or extra dust from young stars and find similar conclusions to the results presented later.
We assess the impact of the variation of dust attenuation across the galaxy on the resultant attenuation curve by randomly drawing the A V from a normal distribution with a mean (µ) of 2 and standard deviation (σ 0 ) of 0.3.We quantify the spread within the galaxy by varying the standard deviation from 2σ 0 − 14σ 0 i.e.where A V,i is the A V value of the i th clump and n ∈ [2, 14], quantifying the spread.Any negative A V,i value is changed to 0.01.The reason for choosing this particular mean A V along with wide values of the spread in attenuation, is to replicate the variations in the dust attenuation that are similar to the galaxies in the FLARES simulation.This will become increasingly clear from the following discussions.
In Figure 5 we show the result of our experiment.The left subplot shows the intrinsic SED of the galaxy in dotted blue line, with the various coloured solid lines showing the observed SED after assigning the different values of dust attenuation for the star forming clumps forming the galaxy.The right subplot shows resultant attenuation curves for the various n σ 0 choices in solid coloured lines.The input SMC curve is also shown as the green dotted line.Emission line features are not seen in the resultant curves since all the star forming clumps have the same age and metallicity.It can be clearly seen that with increasing value of n (i.e. more spread in dust attenuation), the resultant attenuation curve starts to deviate from the input SMC curve (one recovers the input SMC curve, if the different clumps in the toy galaxy were provided with similar values of A V magnitude).Higher values of n make the curve significantly flatter than the input curve, and at the longer wavelength end the normalisation is higher as well.Previous studies have also shown that in a clumpy medium, one can go from SMC type dust law or grain composition to produce a Calzetti like law which is significantly greyer than the SMC dust law (e.g.Witt & Gordon 2000;Inoue 2005;Narayanan et al. 2018).
In forward modelling the FLARES galaxies, we adopted an attenuation curve which was similar to the SMC curve, but less steep, as described in Equation 3. In Figure 6, we show the shape of the attenuation curve for the FLARES galaxies at z ∈ [5, 10] in different SFR bins.One immediate insight is that the attenuation curve gets steeper with decreasing SFR.The galaxies with the lowest star formation rate at z > 8 have a similar shape to our input extinction curve (Equation 3).In case of the lower redshift galaxies, the shape becomes steeper than our input curve as well as the SMC curve.This has also been shown to be the case by other simulation works (see Narayanan et al. 2018) in terms of the V-band optical depth.This is mainly due to the optical light of the galaxy being dominated by old stars with less dust obscuration.The young stars in these galaxies are formed in regions of higher metallicity.Since our dust model links the metal column density to the dust optical depth7 , this will cause these stars to exhibit higher dust obscuration.This increases the A λ /A V value in the UV, thus giving rise to steeper attenuation curves.
For highly star-forming galaxies the opposite effect is apparent, with the shape lying very close to the Calzetti starburst curve.This is not surprising considering our discussion based on our toy galaxies, as well as previous works (e.g.Városi & Dwek 1999;Witt & Gordon 2000;Inoue 2005;Lin et al. 2021) that have already shown that the attenuation law through a clumpy medium or with increasing density along the line-of-sight with SMC dust type can look similar to the Calzetti law.As we saw in the previous section, there is a direct correlation between clumpiness and higher SFR for galaxies in FLARES.The effect of this correlation can be seen from the fact that with the increasing star formation rate, the resultant attenuation curve becomes flatter.Another feature seen is the higher normalisation of the resultant attenuation curve at longer wavelengths for highly star-forming galaxies, similar to what we saw for our toy galaxies, which had higher spread in the dust attenuation.The change in the steepness of the attenuation curve between starforming and passive galaxies was also shown in Wuyts et al. (2011), who attributed it to changes in the dust grain size distribution, while other works (e.g.Inoue 2005;Chevallard et al. 2013) have shown that different dust optical depths can also contribute to that effect.Here we demonstrate the effect of the latter, with changes in the attenuation curve arising out of geometrical effects rather than changes in the dust grain properties.
The above results imply that by using a single extinction or effective attenuation curve for galaxies with varying SFR or UV attenuation can lead to incorrect derivation of physical properties (discussed further in §5).Fig. 7 is the same figure as 6, but with a selection criteria imposed to mimic the galaxy sample in Scoville et al. (2015) (galaxies brighter than I AB =25 magnitude and in the redshift range [4.77, 6]), to directly compare to their observed dust attenuation curve.We can see that there is a good match to the shape of the curve at wavelengths short of the rest-frame near-UV, where the curves agree very well within the spread.At longer wavelengths the normalisation is slightly higher compared to the plotted attenuation curves.This implies that in the bright FLARES galaxies there is either more dust attenuation at longer wavelengths, or the attenuation at 1300Å is lower.This is a direct consequence of the clumpy dust geometry in the star forming galaxies in FLARES, as shown also in our toy galaxies, due to the resultant A V values moving towards lower values.We also tested using our toy galaxy sample that this persists even when excluding contribution from nebular emission.
It should be noted that the input extinction curve we used did not have the UV-absorption feature (commonly attributed to absorption by PAHs) to produce the bump feature, however, the output curve does indeed show a very weak UV-bump similar to the Scoville et al. (2015) curve.The presence of the bump feature is also visually magnified due to our normalisation being higher at longer wavelengths.Alternative explanations of the bump also include graphite grains, which are included in our nebular emission model.Even though the graphite absorption is not exactly centered at 2175Å, the presence of graphite can thus introduce a weak UV-bump in our high-redshift sample of bright galaxies.

SOME OBSERVATIONAL CONSEQUENCES
In this section we will discuss a few of the observational consequences of the effects described in the previous sections.The spread in age and chemical composition of star-forming regions, as well as the variations in their dust attenuation, have an impact on emission line ratios and the inferences drawn from analysing them.To explore this we will concentrate on two well used line ratio methods, the Balmer decrement to correct for dust attenuation, and the BPT diagram to classify galaxies into star-forming and AGN dominated.Finally, we also look at the effect that viewing angle has on the photometry by analysing the spread in UV attenuation, UV-continuum slope and the Balmer decrement for different galaxy viewing angles.

Balmer Decrement
The Balmer decrement, is calculated as the ratio of the Balmer transitions (transitions to the second energy level of hydrogen), usually using H and H γ (level 5 − → 2, 4340.46Å), with each becoming progressively weaker.These transition strengths can be calculated as they are determined by atomic constants.Thus the ratio of these lines (flux or luminosity value) are usually used as a good indicator of the amount of dust attenuation within a galaxy.The dust-free or intrinsic line ratios are for electron density 10 2 cm −2 and temperature 15, 000 K, conditions typically expected for high-redshift galaxies (Sanders et al. 2023).These ratios have a weak dependence on the electron density and are more sensitive to the temperature.However, these variations are small compared to those that can be introduced by the presence of dust (Groves et al. 2012).
The colour excess due the effect of dust on the observed ratios can be written as, where i and j stands for the different Balmer lines (α, β and γ), with the subscripts 'int' and 'obs' corresponding to the intrinsic (dustfree) and observed ratios of the Balmer lines, respectively.Using, where k(λ ) is the shape of the attenuation curve, one can then write, By using an assumed attenuation curve, the values of k(H i ) and k(H j ) can be derived to infer the colour excess.
Figure 8 shows the relationship between the Balmer decrement (here quantified based on the negative of Equation 12) and the attenuation (calculated by integrating the dust corrected luminosity from each star particle as described in § 2.2) suffered by the H α (left) and H β (right) line for the FLARES galaxies in z ∈ [5, 10].The points are coloured by their redshift.Also shown is the expected relations arising from the Calzetti and SMC attenuation curves, as well as our input curve for values of the balmer decrement using Equations 13 and 14.In case of the Calzetti attenuation curve, we also plot 2 modifications of it, by changing the slope of the curve by multiplying it with a wavelength dependent power law, with index δ i.e. λ δ .Now one can modify the Calzetti curve (δ = 0), to a steeper (δ = −1) or shallower (δ = 0.5) one.Also plotted on the figure (coloured stars) are the values for the toy galaxies presented in §4 as a guide to understand the relationship between the variation in attenuation within a galaxy and the deviation from expected relation.We now add a few modifications to the intrinsic spectrum of the galaxies, by introducing a spread in the stellar metallicities (hexagons) and stellar ages (squares) of the star forming clumps.This means for one sample we vary the metallicities, keeping the ages fixed, while in the other sample, we vary the ages, keeping the metallicities fixed.The metallicities and ages are randomly sampled from a uniform distribution spanning [0.01Z ⊙ , 1Z ⊙ ] and [1Myr, 10Myr] respectively.We then perform the same activity as in § 4; obtain the observed spectrum of the galaxies by varying the A V magnitude of the different clumps.
We see that H α and H β attenuation is well correlated with the Balmer decrement, with Pearson correlation coefficients of ∼ 0.85.However, as we can see from the figure, the correlation is in-fact accompanied by large scatter, even at very low values (≲ 0.05, where we expect to have very low dust attenuation), from the expected relation under assumption of any particular attenuation/extinction curve.In case of the H β /H γ ratio, the scatter is higher at low values compared to the H α /H β ratio.When the Balmer decrement approaches extremely low values (∼ 0), the median attenuation of the lines is similar to our input curve.However, as the Balmer decrement increases, the inferred line attenuation deviates from both the expected values on the input curve as well as from all the plotted theoretical curves, trending towards higher values.It can also be seen that, even when the observed Balmer decrement is small, the discrepancy from the expected line attenuation can be a factor of ≳ 10.This is not unexpected, since the relation assumes that all the line emitting regions are equally attenuated by dust.We have already seen that differential attenuation arising from star-dust geometry, can indeed change the shape of the effective attenuation curve.This is indeed seen in Figure 8, with the shallower Calzetti curve (δ = 0.5) tracing the galaxies better.The steepening of the attenuation curve with increasing dust attenuation has also been shown for the high-redshift (z ∈ [4.4,5.5]) ALPINE galaxy sample in Boquien et al. (2022).However, we note here that, in many cases the high value of δ = 0.5 too fails to capture the scatter fully.Therefore, in many cases, using these line ratios can be somewhat circular, since we do not know a priori the effective attenuation curve of the galaxy.It should also be noted that young star-forming regions (nebular emission regions) are expected to suffer higher dust attenuation than the older stars, hence the estimates from simple screen type attenuation curves can be underestimated.
Another reason for the scatter is the distribution of stars with different stellar ages and metallicities, due to variations in star formation histories.This can also be seen from the figure, where we have plotted our toy galaxies where we varied the stellar metallicity and age, adding dust distribution variations on top of that.However, it is hard to draw conclusions on how exactly such variations skew our inferences.This has also been seen from SDSS-IV MaNGA (Ji et al.For reference also plotted is the median relation of the FLARES galaxies.The expected relation for the Calzetti-starburst (for slopes of -1, 0 and 0.5), SMC and input (for the FLARES galaxies) attenuation curves are plotted.The star, square and hexagon shaped coloured points are the relationship for the generated toy galaxies by only varying the dust attenuation towards the stars, by varying the metallicity and dust attenuation of stars, and varying the ages and dust attenuation of the stars respectively (see text for more details).
2023) analysis of dust attenuation from using different optical-to-NIR emission lines.They also infer that no single attenuation curve from literature can explain the attenuation of all the nebular lines within their different galaxies.Chen et al. (2023) using JWST NIR-Cam data have already provided some insight regarding the use of line ratios to trace reionisation era galaxy properties.In their work they study resolved galaxy images, and probe the ages and optical depth of the observed clumps within the galaxy to understand the variation of the [OIII]λ 4959, 5007+Hβ line strength across the galaxy, which is a probe of young star forming regions.They conclude that the variation in ages (and thus metallicities) can have a large impact on the inferred average properties of galaxies.
In Figure 9, we plot the intrinsic and observed H α luminosity function for FLARES galaxies in z ∈ [5, 10].We plot alongside the Balmer decrement derived H α luminosity function using our input attenuation curve as well as Calzetti curve with powerlaw slopes of δ = 0 and 0.5, due to better match with the bulk of the data points in Figure 8.As can be seen from the figure, the input as well as the the Calzetti δ = 0 attenuation curve produce very identical dustcorrected luminosity functions, while the Calzetti δ = 0.5 curve produces a higher correction than the both of them.The effect of this can be clearly seen as we progress down in redshift.At z = 10, 9 and 8, all the used curves manage to reproduce the intrinsic luminosity function well, owing to the fact that there is less dust obscuration in most galaxies at these redshifts.At z = 7, 6 and 5, all the curves agree at the fainter end (≤ 10 42.5 erg/s), while at the brighter end, the FLARES input curve as well as the Calzetti δ = 0 curve underestimate the dust correction, while the Calzetti δ = 0.5 curve provides a reasonable match.This implies that given different attenuation curve choices, one can get reasonable dust corrections by employing a very shallow curve.

BPT diagram
The BPT (Baldwin, Phillips & Terlevich) diagram (Baldwin et al. 1981), is a set of nebular emission line ratio diagrams used to distinguish the most dominating ionising source within a galaxy.The most famous version, uses the nebular line ratios of [OIII]λ 5008/H β vs [NII]λ 6585/H α (referred to as N2-BPT) to classify galaxies into star-forming and AGN hosts.Several studies (e.g.Brinchmann et al. 2008;Shapley et al. 2015;Strom et al. 2018;Garg et al. 2022) have looked at the loci of the galaxies inhabiting this space to understand how differences in the ionisation parameter, abundance ratios and metallicities contribute to galaxies moving across this space.Here we will look at how the variation in dust attenuation across different stars in a galaxy affect the BPT diagram.
Figure 10 shows the BPT diagram of the intrinsic line ratios (dust-free, left) and the dust attenuated line ratios (right).It can be seen that the locus of the points are very similar, however the density of the points has now shifted towards the left side in the dust attenuated plot compared to the intrinsic one.We can quantify the shift in terms of the median difference in the absolute values of the intrinsic and observed line ratios.The median shift in log 10 ([NII]λ 6585/H α ) is 0.08 with a standard deviation 0.2, while in the case of log 10 ([OIII]λ 5008/H β ) the median shift is 0.02 with a standard deviation of 0.05 across the plotted redshift range.It is )) Figure 9. H α luminosity functions of galaxies in the FLARES for z ∈ [5, 10].The intrinsic and the dust attenuated value from the FLARES galaxies are plotted in blue and orange colours respectively, with the shaded region denoting the Poisson 1-σ errorbar (for only bins with more than 5 data points).Plotted alongside is the H α luminosity function obtained from using the Balmer decrement (using the observed H α and H β ratios from FLARES) corrected H α applying the FLARES input attenuation curve and the Calzetti-starburst (Calzetti et al. 2000) curve (for δ = 0 and 0.5).Note that the star particle cut mentioned in § 2.1 is not enforced here to preserve the full shape of the luminosity function.Also the Poisson 1-σ errorbar has only been shown for the Calzetti δ = 0.5 curve.
quite reasonable to wonder why there is such a difference in the intrinsic and observed line ratios, when the ratios are calculated for emission lines that are very close in wavelengths.The variation occurs because galaxies are a composite mixture of stellar populations having different ages and metallicities, which in turn determine the intrinsic line strengths of [NII]λ 6585 and H α or OIII]λ 5008 and H β .The observed (i.e.attenuated) line strength depends on the dust along the line-of-sight to these different stellar populations.The total line strength is then the sum of these individual values.However, in the BPT diagram, the x and y axes represent the ratio of the sums of these individual line strengths, not an average of these individual ratios.The variation in line ratios occurs because the sum of observed ratios is not exactly equal to the sum of the intrinsic ratios multiplied by the effective attenuation i.e.
where F i,k and F j,k are the line fluxes of the respective lines used to compute the ratios from the k th component in the galaxy.T i,k is the attenuation experienced by the line 'i' for the k th component, whereas T eff is the assumed effective attenuation for the galaxy.
Since the ratio T eff,i /T eff, j is calculated for wavelengths that are very similar, it can be roughly approximately as 1.Therefore significant deviations from this assumption of effective attenuation occurs when there is a varying amount of dust along the line-of-sight of stellar populations with differing physical properties.The seen shift in location can be explored with the toy galaxies, by plotting them on Figure 10 (same marker types as in § 5.1, Figure 8).If we look at the simplest case of varying dust attenuation and keeping all the other properties same, we can see that a spread in dust attenuation moves the point vertically (up and down).Once you add to this, a spread in the metallicity of the stars, this moves the point along the right-to-left diagonal while a spread in ages moves the point along the left-to-right diagonal.Thus we can see that a combination of these effects can actually move a point from the intrinsic space towards the left in the observed space.Therefore, when interpreting physical conditions based on line ratios in dusty galaxies, it is advisable to exercise caution, even when dealing with line ratios from wavelengths that are closer to each other.

Line-of-sight variation
The viewing angle can also change the observed properties of the galaxies that we have discussed in the previous section.This is because the dust distribution seen by spatially distinct stellar population is highly dependent on the viewing angle.This effect of the viewing angle can be explored with hydrodynamical simulations on various observed galaxy properties.In the previous sections, all the discussion assumed a single line of sight (i.e.along the z-axis) to the galaxy.In this section we investigate the variation in a few galaxy observables, such as the UV-attenuation, the UV-continuum slope and the balmer decrement, when our sample galaxies are observed along different lines-of-sight (as described in § 2.4).
In §2.4 we explained how by using HEALPY we calculate different line-of-sight re-projections.In Figure 11 we show examples of UV-attenuation maps for the different lines-of-sight for 6 galaxies in FLARES at z = 5.The sightlines, represented by the healpix bins, have been coloured by the galaxy UV attenuation.The blue and yellow crosses denote the sightlines looking along the angular momentum unit vector corresponding to the gas and stellar particles, respectively, of the galaxy disc.The blue and yellow lines denote the plane perpendicular to the corresponding angular momentum unit vectors, aligning with sightlines going through the gas and stellar disc of the galaxy.

Population Trends
In Figure 12, we plot the 1σ spread (the difference between the 84 th − 16 th percentile), across different lines-of-sight, of the UVattenuation (A FUV ), the UV-continnum slope (β ) and the balmer decrement 2.5 log 10 against their median values.
The different coloured lines are the median spread for bins of different SFRs.From the figure it can be seen that the spread in the value of the observed A FUV from its median value can be quite significant, reaching a maximum of ∼ 0.4 at z = 10 to ∼ 0.7 at z = 5.The variation is extremely low (0 − 0.2) at low values of the median UV attenuation (< 1 Mag), while it can range from 0.2 − 0.8 in heavily attenuated (A FUV > 2) galaxies.It is a similar case for β as well, reaching a maximum spread from the median β of ∼ 0.1 at z = 10 to ∼ 0.2 at z = 5.Here, we also see a weak dependence of the spread on the star formation rate of the galaxy, with galaxies in higher SFR bins exhibiting a greater spread from their median value.This is another indication that the UV slope of low-SFR galaxies, or those in the passive regime, are not affected by dust to the same extent as the others.In case of the Balmer decrement, the variation is extremely low at z = 10 (maximum of ∼ 0.05), with the deviation from the median value increasing with decreasing redshift (reaching a maximum of ∼ 0.15 at z = 5).However, from Figure 8 we know that there is considerable spread in the expected attenuation even at a single value of the Balmer decrement, hence small shifts in the value due to inclination effects can further add to the uncertainty of the derived parameter.The increased dispersion around the median relation towards lower redshifts can be attributed to the buildup of dust in galaxies, which means that the star-dust geometry has a more significant impact.Additionally, there is a higher prevalence of galaxies with well-formed discs.And in such galaxies, the difference between edge-on and face-on orientations will be more pronounced.

Correspondence with axis of rotation
Most flagship periodic cosmological simulations (e.g.EAGLE, ILLUSTRIS-TNG, SIMBA, etc), at the redshifts we are looking at have a dearth of disc galaxies due to the smaller physical volume they probe.Therefore, they sample fewer instances of the highly overdense regions, which are rare and expected to host the massive galaxies in the early Universe, including the formation of disc galaxies.This is not a problem in FLARES, due to the novel resimulation strategy employed, with regions ranging from extremely underdense to extremely overdense regions.And with the resimulated regions in FLARES picked from a (3.2cGpc) 3 (see 2021).However, for this work, we quantify a galaxies' disciness by computing the energy invested in ordered rotation (see Sales et al. 2012;Correa et al. 2017, for more details).The following equation is used to compute the value: where the sum is over all the particles of a particular type (gas or stars) within our aperture centred on the potential minimum, m i is the mass of the particle and K is the total kinetic energy (∑ i 1 2 m i v i ) 2 , with v i the velocity of the particle).L z,i is the particle angular momentum along the direction of the total angular momentum of the component type of the galaxy and R i the projected distance to the axis of rotation.We compute κ co for the gaseous and stellar component of the galaxy, denoting them as κ co,gas and κ co,⋆ respectively.
To understand the relationship between the spread in the UV attenuation and β , we plot a hexbin distribution of these values against the dust surface density in Figure 13; the dust surface is calculated within 2 × R 1/2,dust , where R 1/2,dust is the half-mass dust radii.We colour the hexbin distribution by the median κ co,⋆ value within the bin.It can be clearly seen that there is a weak evolution in the spread of both the UV attenuation and β with the dust surface density, with the extreme variation in the spread taken up by galaxies that are very discy.There is a clear gradient in κ co,⋆ values moving from low to high spread.This can also be seen in Figure 11 where the maximum of the attenuation traces the line-of-sight corresponding to the plane of the disc in these galaxies i.e. viewing through the galaxy gaseous/stellar disc.

CONCLUSIONS
In this work we have explored how variations in dust attenuation across a galaxy, affect their observed properties using the First Light And Reionisation Epoch Simulations (FLARES).Our conclusions are as follows: • We explore the radial profile of the UV traced star formation rate surface density as well as the UV attenuation, finding that the star formation in the brightest galaxies in FLARES have extended towards the outskirts.However these stars experience less dust attenuation than the stars formed in the galaxy cores.This results in a non linear relationship between the intrinsic UV luminosity and the UV attenuation at the bright end, with the FLARES UVLF exhibiting a double power-law shape.
• We also look at the spread in the UV attenuation of different star particles within a galaxy, finding huge spread, with many galaxies having stars that are fully obscured to fully unobscured.We found that the magnitude of this spread is correlated with the star formation rate of galaxies in FLARES, with the extreme star forming galaxies exhibiting the largest spread, while less star-forming galaxies have more homogeneous dust distribution or very low dust attenuation.
• We explore the consequences of this spread by looking at the resultant attenuation curve of galaxies in FLARES.We find that low SFR galaxies exhibit very steep attenuation curve while galaxies with extreme SFRs show flatter attenuation curves.
• We explore some of the observational consequences of this spread in dust attenuation across the galaxy, by probing its effects on the Balmer decrement and the BPT diagram.Dust complicates the interpretation of these line ratios.The expected attenuation calculated from the Balmer decrement does not follow the widely used attenuation/extinction curves from literature.The location of points in the BPT diagram compared to the points before dust attenuation are also changed.Inhomogeneous dust distribution coupled with the variation in stellar ages and metallicities across the galaxy is the main driver of these changes.
• We also explore the effect of variations in the viewing angle towards the galaxy.Different viewing angles results in different stardust geometry, which also introduces spread in the observed values of the UV attenuation, UV-slope and the Balmer decrement.The extreme variations in this space are observed for galaxies that are very discy.
With recent influx of data from JWST of the high-redshift Universe, it is now possible to explore galaxy ISM conditions using nebular emission lines (e.g.Arellano-Córdova et al. 2022;Endsley et al. 2023;Cameron et al. 2023;Sanders et al. 2023).These have already provided some insights into the ISM conditions such as the ionisation parameter, metallicity and elemental abundances.There is promise of even more spectroscopic data, with programmes covering large areas of the sky (e.g.Cosmos-Web (Casey et al. 2023), UNCOVER (Bezanson et al. 2022)), that will be able to provide statistically significant samples of galaxies.These will enable to understand the diverse conditions within various galaxy populations across redshift.However, our results suggest that we should be careful in the interpretation of these results.Galaxies are com-plex structures consisting of distinct regions with varying physical conditions and properties.Consequently, many commonly used diagnostic measures, such as the Balmer decrement and BPT diagram, are a composite measure of these diverse physical conditions.This diversity in properties can conspire together to yield diagnostic measurements that differ from the actual conditions present within the galaxy, ultimately impacting our conclusions.

Figure 1 .
Figure 1.Shows the relationship of the galaxy stellar mass (left) and the galaxy UV luminosity (right) against their star formation rate (SFR) for z ∈ [5, 10].The coloured solid lines represent the weighted median of the relation at different redshifts, along with the 84 th − 16 th percentile spread.Also shown is the histogram of the distribution of stellar mass and SFR in different bins for these redshifts.The total number of galaxies at these redshifts are indicated within brackets alongside the legend, with the colourbar indicating the number of galaxies within the hexbins.

FUVFigure 2 .
Figure 2. Variation in the UV star formation surface rate density (solid lines) and the UV-attenuation (dotted lines) with distance from the galaxy centre, in units of the half-mass radius.Only M FUV bins with ≥ 10 galaxies are plotted.There is only negligible change if we use an SFR averaging timescale of 10 Myr instead of 100 Myr.
shows the spread in the UV attenuation (or its proxy, f esc in equation 4, now for individual star particles; where f esc − → 0 indicates a fully dust obscured star form-ing region, and vice versa) in the UV for z ∈ [5, 10].Here we describe the scatter within a galaxy, by the difference between the 84th and 16th percentile of the different f esc values of star particles, i.e. f esc,84 − f esc,16 .A high value of f esc,84 − f esc,16 (− → 1) is indicative of extremely inhomogeneous dust distribution within the galaxy.

Figure 5 .
Figure5.Left: The intrinsic (blue, dotted) and observed (for different spread in attenuation, in coloured solid) spectrum of the galaxy by summing up the contribution from the different star-forming clumps made with SED fitting code Bagpipes.The regions were all assigned the same stellar ages with a single burst SFH, metallicity, ionisation parameter and mass (see §4).Each clump follows an SMC dust attenuation law, where we vary the line-of-sight attenuation by varying the A V magnitude.Right: Attenuation curves of the resultant galaxy, which is shown in different colours, obtained varying the dust attenuation of the individual clumps.The green dotted line is the input SMC attenuation curve.

Figure 6 .Figure 7 .
Figure6.Galaxy attenuation curve for the FLARES galaxies in different bins of SFR Total .The different panels shows the attenuation curve for z ∈ [5, 10].The curves have been smoothed and the emission lines removed for better visual clarity.We also plot alongside attenuation curves from literature: SMC(Pei 1992),Calzetti-Starburst (Calzetti et al. 2000)  and the Milky Way curve fromNarayanan et al. (2018, MW N18)  as well as the FLARES input attenuation curve.

Figure 8 .
Figure 8. Balmer decrement inferred H α (left) and H β (right) attenuation for the FLARES galaxies in z ∈ [5, 10].The galaxies are coloured by their redshift.For reference also plotted is the median relation of the FLARES galaxies.The expected relation for the Calzetti-starburst (for slopes of -1, 0 and 0.5), SMC and input (for the FLARES galaxies) attenuation curves are plotted.The star, square and hexagon shaped coloured points are the relationship for the generated toy galaxies by only varying the dust attenuation towards the stars, by varying the metallicity and dust attenuation of stars, and varying the ages and dust attenuation of the stars respectively (see text for more details).

Figure 10 .
Figure 10.BPT diagram for galaxies in z ∈ [5, 10].The galaxies are coloured by their redshift.For reference also plotted is the median relation from the Keck Baryonic Structure Survey (KBSS) at z ∼ 2 (Strom et al. 2017).The star, square and hexagon shaped coloured points are the relationship for the generated toy galaxies by only varying the dust attenuation towards the stars, by varying the metallicity and dust attenuation of stars, and varying the ages and dust attenuation of the stars respectively.

Figure 11 .
Figure11.Healpix maps for different sightlines for 6 galaxies at z = 5 in our sample.The sightlines have been coloured by the galaxy UV attenuation.The blue and yellow crosses denote the sightlines looking along the angular momentum unit vector corresponding to the gas and stellar particles within the galaxy, with the lines denoting the plane perpendicular to the unit vectors.Also shown are the SFR, κ co,⋆ and κ co,gas values.