# Constraining the Milky Way’s Ultraviolet to Infrared SED with Gaussian
Process Regression

Catherine E. Fielder,1,2 Jeffrey A. Newman,1,2 Brett H. Andrews,1,2 Gail
Zasowski,3 Nicholas F. Boardman,3Tim Licquia,1,4 Karen L. Masters,5 Samir
Salim,6
1Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh, PA
15260, USA
2Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC),
University of Pittsburgh, Pittsburgh, PA 15260, USA
3Department of Physics and Astronomy, University of Utah, Salt Lake City, UT,
84112, USA
4Core R&D: Information Research, Dow Inc., 633 Washington St, 1776 Building,
Midland, MI, 48667, USA
5Departments of Physics and Astronomy, Haverford College, 370 Lancaster
Avenue, Haverford, PA 19041, USA
6Department of Astronomy, Indiana University, Bloomington, IN, 47405, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Improving our knowledge of global Milky Way (MW) properties is critical for
connecting the detailed measurements only possible from within our Galaxy to
our understanding of the broader galaxy population. We train Gaussian Process
Regression (GPR) models on SDSS galaxies to map from galaxy properties
(stellar mass, apparent axis ratio, star formation rate, bulge-to-total ratio,
disk scale length, and bar vote fraction) to UV (GALEX $FUV/NUV$), optical
(SDSS $ugriz$) and IR (2MASS $JHKs$ and WISE $W1/W2/W3/W4$) fluxes and
uncertainties. With these models we estimate the photometric properties of the
MW, resulting in a full UV-to-IR spectral energy distribution (SED) as it
would be measured externally, viewed face-on. We confirm that the Milky Way
lies in the green valley in optical diagnostic diagrams, but show for the
first time that the MW is in the star-forming region in standard UV and IR
diagnostics—characteristic of the population of red spiral galaxies. Although
our GPR method predicts one band at a time, the resulting MW UV–IR SED is
consistent with SEDs of local spirals with characteristics broadly similar to
the MW, suggesting that these independent predictions can be combined
reliably. Our UV–IR SED will be invaluable for reconstructing the MW’s star
formation history using the same tools employed for external galaxies,
allowing comparisons of results from in situ measurements to those from the
methods used for extra-galactic objects.

###### keywords:

Galaxy:general – Galaxy:fundamental parameters – Galaxy:structure

††pubyear: 2021††pagerange: Constraining the Milky Way’s Ultraviolet to
Infrared SED with Gaussian Process Regression–D

##  1 Introduction

Within the Milky Way we have a unique opportunity to study the nuances of
galactic properties, allowing us to test galaxy formation and evolution models
at an unrivalled level of detail. For example, chemical abundances for
hundreds of thousands of stars have been obtained from spectroscopic surveys
(Majewski et al., 2017; Martell et al., 2017), and stellar surveys that
catalogue distance and dynamical measurements on millions of stars have been
performed (Gaia Collaboration et al., 2018). This exhaustive stellar
information has helped constrain the Milky Way’s evolutionary history. In
turn, high-resolution dynamical simulations have been able to produce galaxies
of increased similarity to the Milky Way (Guedes et al., 2011; Sawala et al.,
2016; Wetzel et al., 2016), matching fundamental galaxy properties such as
dwarf satellite populations and reproducing characteristics of our Galaxy’s
gas, dust, and stellar components. Comparisons between Milky Way stellar data
and high-resolution hydrodynamical simulations of Milky Way-like galaxies are
an increasingly useful way to improve our understanding of galaxy formation.

However, the marriage between observations and models is delicate: incorrect
assumptions on one side can propagate into the other. Simulators must make
choices about how to implement crucial parameters that affect the galaxy
evolution process, such as the gas density threshold for star formation to
occur and the efficiency with which it proceeds; this is sometimes done by
attempting to match observed properties of the Milky Way. But without knowing
how our Galaxy fits in amongst the broader galaxy population, it is difficult
to determine whether simulations match the Milky Way because they have the
correct physics or because they have incorrectly tuned parameters that match
by coincidence or design.

This is complicated by the exceptional difficulty of obtaining a global
picture of the Milky Way, given our location in the disk and the obscuration
caused by interstellar dust. As a result, there are properties that we can
easily measure in external galaxies that are impossible to measure directly
within our own, making it difficult to determine where we fit within the
broader galaxy population.

Creating an outside-in picture of the Milky Way that spans a multitude of
broad-band wavelengths will enable simulators to more accurately tune their
physics assumptions, as it will then be possible to test whether quantities
that can only be determined from large-scale stellar surveys and those that
can be measured directly only for extragalactic objects are reproduced. The
most basic, easiest-to-measure intrinsic quantities we can use to study
galaxies are their luminosities and colours; once redshift is known, these can
be inferred from broad-band photometry. Hence the focus of this paper will be
in determining these properties for the Milky Way. This will enable our Galaxy
to be placed on standard colour-magnitude and colour-colour diagrams and
result in a multi-wavelength spectral energy distribution (SED) for the Milky
Way.

Astronomers have found that galaxies in the local Universe predominantly fall
into two populations: passively-evolving red galaxies with older stellar
populations; or blue galaxies that are still forming stars. In the optical
colour-magnitude diagram (CMD), these two galaxy populations are commonly
referred to as the “red sequence” and the “blue cloud”. The colour bimodality
of galaxies has been observed at both low redshift, ($z\sim 0.1$, e.g.,
Strateva et al., 2001; Baldry et al., 2004) and up to redshift of $z\sim 1$
(e.g., Bell et al., 2004; Weiner et al., 2005). The region of the CMD between
these two distinct galaxy populations is often referred to as the “green
valley”. This locus is thought to contain a transitional population of
galaxies that are “passively” evolving in the sense that no new star formation
is occurring (e.g., Bell et al., 2004; Faber et al., 2007), though they still
can contain some younger stars. The increase in the fraction of red galaxies
over time has led many astronomers to conclude that a galaxy first lives in
the star forming blue cloud and then transitions into the green valley and
ultimately into the red sequence in complete quiescence, with the galaxy
growing more and more red over time due to the ageing stellar population.
Green valley galaxies are presumed to be undergoing some form of quenching of
their star formation (Salim et al., 2007; Schawinski et al., 2014; Smethurst
et al., 2015) \- either late-type galaxies that are gradually running out of
their cold gas reservoir or having their star formation suppressed, or early-
type galaxies that had their gas reservoirs rapidly destroyed. Objects found
within the green valley may also simply be in the tails of the blue cloud or
red sequence, rather than being in a transitional state (Taylor et al., 2015).
While the precise details of quenching processes and the origin of the galaxy
colour bimodality have yet to be determined, a galaxy’s location on the CMD
remains a very useful tool for determining how a galaxy fits into the broader
population.

The radiation emitted by a galaxy is characterised by its spectral energy
distribution (SED), or flux as a function of wavelength. Galaxy SEDs contain
the imprints of the physical processes occurring within - the stellar
population’s ages and abundances (i.e., the star formation history and
metallicity of the galaxy), the dust and gas content and the chemistry and
physical state of the interstellar medium (ISM). Because different sources
dominate the emission at different wavelengths, long-wavelength-baseline SEDs
allow one to disentangle the contributing effects. This makes SEDs one of the
best direct probes for studying galaxy formation and evolution from both an
observational and theoretical modelling perspective.

However, comparing colours and luminosities of the Milky Way to external
galaxies is not trivial, regardless of whether we compare to observed galaxies
or to mock images from high-resolution hydrodynamical simulations (such as
Eris Guedes et al. 2011, APOSTLE Sawala et al. 2016, and Latte Wetzel et al.
2016). Much of our view of the Galaxy is obscured by interstellar dust,
especially at UV and optical wavelengths (e.g., Cardelli et al., 1989;
Schlegel et al., 1998). Stars outside of the local solar region are reddened
as a result of the dust obscuration. Determining the integrated light of
stellar populations in the Milky Way is challenging due to the spread of stars
over large and varying distances, with correspondingly large and varying dust
extinction along lines of sight to the Earth. This makes the study of any
portions of the Galactic disk beyond the solar neighbourhood exceptionally
difficult, and results in a fragmented picture of the Milky Way. Integrated
properties that are relatively painless to obtain in external galaxies (though
dust obscuration can affect these observations as well, see e.g., Masters et
al. (2003); Masters et al. (2010a)), are impossible to obtain directly within
our own Galaxy (e.g., Mutch et al., 2011). As a result, simulators often must
resort to comparing their simulated Milky Ways to very general galaxy
populations (such as sets of Sbc or late type galaxies) that, while
superficially resembling the Milky Way, have a wide range of other global
properties (e.g., Guedes et al., 2011).

In an effort to circumvent our limited view of the Milky Way, we can study
galaxies that mimic the properties of our Galaxy but can be observed from
outside, which we label Milky Way Analogues (MWAs). This method hinges on the
Copernican assumption that the Milky Way should not be extraordinary among a
galaxy population that shares some key properties with it. These comparisons
are enabled by working within volume-limited subsets of large surveys, which
ensure that the observed objects constitute a representative population.
Previous work suggests that galaxies with similar stellar mass and star
formation rates are also similar in other properties, as the observed galaxy
population is well-matched by models that parameterise galactic star formation
histories with a limited collection of curves (Behroozi et al., 2013; Gladders
et al., 2013; Abramson et al., 2014; Kelson et al., 2014). Even further, Bell
& de Jong (2001) showed mass and star formation rate are strongly correlated
with the photometric properties of a galaxy. Therefore we can exploit the fact
that two galaxies of identical mass and star formation rate should have
similar luminosities and colours, with some scatter given the range of galaxy
photometric properties at fixed physical parameters. Licquia et al. (2015),
hereafter LNB15, utilised this to constrain the Milky Way’s optical colours
and magnitudes based on the range of observed properties of MWAs that were
matched in stellar mass and star formation rate. MWAs also allow direct
comparison of properties of our Galaxy to its closest peers (e.g., Licquia &
Newman, 2016; Licquia et al., 2016; Fraser-McKelvie et al., 2019; Boardman et
al., 2020a; Krishnarao et al., 2020) and have been a successful tool for
improving our understanding of the Milky Way in an extra-galactic context.

The Milky Way, however, has some characteristics that are atypical (at the
$<2\sigma$ level) amongst its peers - e.g., the Milky Way has an unusually
compact disk (i.e., a small disk scale-length) (Bovy & Rix, 2013; Bland-
Hawthorn & Gerhard, 2016; Licquia & Newman, 2016), and an unusually quiescent
merger history (from observation Unavane et al., 1996; Ruchti et al., 2015),
(and simulation, e.g., Fielder et al., 2019; Carlesi et al., 2020). The
deviations of the Milky Way from the average suggest that we should consider
parameters beyond just stellar mass and star formation rate in order to
identify samples of objects that more closely resemble the Milky Way.

Galaxy morphological characteristics such as disk scale length ($R_{d}$) and
bulge-to-total ratio ($B/T$) are tied to a galaxy’s evolutionary history and
therefore should connect to its photometric properties (Cappellari, 2016; Saha
& Cortesi, 2018) as well as to the ways in which the Milky Way is atypical. We
would therefore wish to incorporate these properties in addition to stellar
mass and star formation rate in defining an MWA. However, as the number of
parameters required to match the Milky Way increases, the number of MWAs
correspondingly reduces dramatically. For example Fraser-McKelvie et al.
(2019) only found 179 analogues when selecting on stellar mass, bulge-to-total
mass ratio, and morphology; Boardman et al. (2020b) and Boardman et al.
(2020a) found no MWAs within $1\sigma$ of the MW when selecting on stellar
mass, star formation rate, bulge-to-total ratio, and disk scale length in
either the SDSS-IV MaNGA survey (Bundy et al., 2015) or a larger photometric
sample drawn from the GALEX-SDSS-WISE Legacy catalogue (GSWLC; Salim et al.
(2016)) respectively.

LNB15 found that the colour of the Milky Way is consistent the green valley
region of the CMD as it has been defined using purely optical passbands.
Characterising the UV and IR colours of the Milky Way can provide more
sensitive probes of whether it would be classified as in the process of
quenching if seen from outside, allowing us to better understand what type of
population the Milky Way may belong to. LNB15 speculated that the Milky Way
might belong to the population of massive “red spiral” galaxies, which are
characterised by their red optical colours despite ongoing star formation
(Masters et al., 2010b; Cortese, 2012). Galaxies within this population may be
moving into the green valley due to slow quenching (cf. Schawinski et al.
(2014)). This conjecture can only be fully tested by examining wavelengths
outside of the optical range; in $g-r$, the colours of massive spiral galaxies
on the star-forming main sequence (a population that should include the
Galaxy) overlap with both the red sequence and the blue cloud (Cortese, 2012;
Salim, 2014). However, samples of MWAs that have high-quality photometry over
a broader wavelength range will have reduced numbers due to the limited
coverage of sufficiently deep photometry in GALEX.

To address the lack of analogues when multi-dimensional parameter spaces are
used, and the smaller overall sample size resulting from the increase in
wavelength coverage, we introduce a Gaussian process regression (GPR) approach
in this work. GPR is an emergent tool in astrophysics. For example, Bocquet et
al. (2020) employed GPR to emulate results of cosmological simulations, while
Gordon et al. (2020) used GPR to detect and classify exoplanets. We can use
GPR to leverage information from a wider variety of galaxies, instead of just
the closest Milky Way analogues, in order to extract information from large-
scale trends between galaxy physical and photometric properties. Thanks to the
probabilistic framework that underlies GPR, we obtain uncertainty estimates
for all predicted quantities for free. The primary result from this paper will
be an ultraviolet to infrared SED of the Milky Way as viewed face-on,
determined via GPR based on star formation history and structural parameters
(i.e., galaxy physical parameters) that have been measured well for both the
Milky Way and galaxies from the Sloan Digital Sky Survey (Aihara et al.,
2011).

The paper is organised in the following manner. In Section 2 we describe the
observational data used, including the external galaxy data in Section 2.1 and
Section 2.2, and estimated properties of the Milky Way in Section 2.3. Section
3 details our new Gaussian process regression-based methodology. In Section 4
we compare the luminosity and colours of the Milky Way at multiple wavelengths
and the Milky Way’s predicted SED to properties of other galaxies. Finally we
summarise our results, and discuss implications and future work in Section 5.
Appendix A provides a summary of the galaxy parameters and tables that list
predicted photometry for the Milky Way. Appendix B describes tests of the
accuracy of the GPR procedures used here, and Appendix C describes how we
address the systematic corrections needed for $k$-corrections and Eddington
bias.

In this paper all magnitudes are reported in the AB system, except for the
Johnson-Cousin $UBVRI$ magnitudes which are presented in the Vega system.
Absolute magnitudes are derived using a Hubble constant
$H_{0}=100\;\text{km}\,\text{s}^{-1}\rm{Mpc}^{-1}$, so they are equivalent to
$\rm{M}_{y}-5\log{h}$ (where $\rm{M}_{y}$ is the $y$-band absolute magnitude
and $h=H_{0}/100$) for other values of $h$. For other properties in which
measurements for the Milky Way are compared to extra-galactic galaxy
measurements we assume $H_{0}=70\;\text{km}\,\text{s}^{-1}\rm{Mpc}^{-1}$
($h=0.7$) in accordance with Licquia et al. (2015), for a standard flat
$\Lambda$CDM cosmology with $\Omega_{m}=0.3$. Parameters such as log stellar
mass and log star formation rate can be modified for different $h$ values by
subtracting $2\log{h/0.7}$. We do this to avoid confusion and to allow for
potential updates to future $h$ measurements.

##  2 Observational Data

In this section we describe the many galaxy catalogues utilised in this work.
We break this up by photometry (Section 2.1) and inferred galaxy properties
(Section 2.2), with the Milky Way measurements included in the final
subsection (Section 2.3).

###  2.1 Photometry

####  2.1.1 SDSS Galaxies

The sample of galaxies that we use as a starting point originates from the
eighth data release (DR8; Aihara et al. (2011)) of the Sloan Digital Sky
Survey III (SDSS-III; York et al. (2000)). DR8 provides both images and
photometry of thousands (almost $10^{6}$) of local galaxies. The optical
broadband passbands, $u$, $g$, $r$, $i$, and $z$ were the subjects of previous
Milky Way analogue work by LNB15 and are used in this study in addition to
bands outside of the optical range.

We make use of both the ‘‘model’’ and ‘‘cmodel’’ magnitudes from SDSS. The
former refers to magnitudes derived from the better of either a de Vaucouleurs
or an exponential profile fit to the galaxy surface brightness distribution.
These types of magnitudes are expected to produce the highest signal-to-noise
estimate of galaxy colours; thus when we refer to galaxy colours derived from
SDSS we will be using model magnitudes for the calculations. Alternatively,
cmodel magnitudes are derived from the best fit to a linear combination of a
de Vaucouleurs and exponential profile. These magnitudes provide the best
estimate of the total flux of a galaxy in each passband. When we refer to
galaxy absolute magnitudes for SDSS bands we will use ‘‘cmodel’’ magnitudes.

$k$-corrections on these passbands to rest-frame $z=0$ were calculated via the
kcorrect v4.2 software (Blanton & Roweis, 2007), as described in LNB15. This
provided AB absolute magnitudes for the SDSS $ugriz$ photometry. Additionally,
kcorrect was used to convert the SDSS $ugriz$ photometry to restframe
Johnsons-Cousins UVBRI Vega magnitudes in order to make easy comparisons to
literature values. Results are presented with the adoption of the Blanton &
Roweis (2007) and LNB15 notation, where an absolute magnitude of passband $y$
at redshift $z$ is denoted as ${}^{z}M_{y}$.

Our main galaxy sample is derived from the volume-limited sample presented in
LNB15. A volume-limited sample is required for accurate results from Milky Way
analogues in order to alleviate a radial selection effect known as Malmquist
bias, i.e., the preferential inclusion of intrinsically bright galaxies. At
higher redshifts within the main SDSS sample (Strauss et al., 2002), only the
most luminous galaxies will be brighter than the sample magnitude limit and
followed-up spectroscopically. By using a volume-limited sample we ensure that
galaxies within the range of the Milky Way’s parameters are included equally
at all distances considered. LNB15 determined the limits for their volume-
limited sample from an initial draw of Milky Way analogues from the full SDSS
DR8 parent catalogue without any redshift cuts. Then in ${}^{0}(g-r)$ vs.
${}^{0}M_{r}$ (i.e., restframe $g$-$r$ colour derived using $z=0$ passbands
versus $r$-band absolute magnitude, again evaluated with the $z=0$ passband)
colour-magnitude space a maximum redshift was chosen such that all objects as
low in luminosity as the faintest Milky Way analogues would still be included
at that $z$. A minimum redshift was also applied to limit the impact of the
finite SDSS fiber aperture on measured galaxy properties. The resulting
volume-limited sample contains a total of 124,232 target galaxies within the
redshift range of $0.03<z<0.09$. Some initial cuts on SDSS quality flags were
also employed; for further details on the construction of this volume-limited
sample refer to LNB15 Section 3.1. All cross-matches from SDSS to other
catalogues presented here were constructed only using the volume-limited
sample. Both the SDSS sample used here and the cross-matched catalogues
presented in this paper are available at our catalogue GitHub
repository111https://github.com/cfielder/Catalogs.

####  2.1.2 GALEX–SDSS–WISE Legacy Catalogue

Photometry in ultraviolet and infrared wavelengths used in this work comes
from the GSWLC (Salim et al., 2016, 2018). We use GSWLC-M2, the medium-deep
catalogue of GSWLC-2, which covers $49\%$ of the SDSS DR10 footprint. While
this reduces the number of targets for study, the improved signal-to-noise in
the UV-imaging over the shallow catalogue enables tighter results. SDSS
Photometry between DR7/DR8 and DR10 is the same, so no issue arises between
cross matches of the SDSS volume-limited sample and the GSWLC-M2 sample. In
order to account for any differences in astrometry we consider matches to be
separations within 1.5 arcseconds.

The ultraviolet sample in the GSWLC catalogue originates from the GALEX survey
($FUV,NUV$), and has been corrected for galactic reddening and calibration
errors (Martin & GALEX Team, 2005). UV detections are available for $74\%$ of
SDSS targets in the GSWLC-M2 catalogue. The infrared sample in the GSWLC
catalogue originated from the 2MASS and WISE surveys. The 2MASS photometry
($JHKs$) is available for $48\%$ of SDSS targets and WISE photometry
($W1,W2,W3,W4$) is available for $41\%$ of SDSS targets.

####  2.1.3 DESI Legacy Imaging Surveys

Data release 8 of the DESI (Dark Energy Spectroscopic Instrument) Legacy
imaging surveys includes $g$, $r$, $z$, and WISE photometry of sources
detected in DECam or BASS/MzLS imaging (Dey et al., 2019). We use these
catalogues for all WISE-band photometry presented here, as the matched-model
measurements from the Legacy Surveys Tractor catalogues go substantially
deeper and have lower errors than other public WISE data products. We also
calculate $k$-corrections for the WISE bands using this photometry; we discuss
our procedures for this, which allow accurate $k$ corrections in the IR
without making any assumptions about SED templates at those wavelengths,
briefly in Section C.1 and more extensively in a forthcoming paper (Fielder,
Newman, & Andrews, in prep.).

DESI Legacy and SDSS contain photometry that comes from differing filter sets
and detectors. DESI Legacy North uses BASS and MOSAIC filters while DESI
Legacy South used DECam filters (Dey et al., 2019). Our corrections for
offsets between the DESI Legacy Survey and SDSS photometry is described in
Section C.2.

###  2.2 SDSS-based Properties

Only one of the galaxy properties which we use to create SEDs is recorded in
the standard SDSS photometric catalogues and required no further analysis;
specifically, the best fit axis ratio ($b/a$) of a galaxy. We use the $b/a$
value determined from an exponential fit to the galaxy’s surface brightness
density in the $r$-band throughout this work.

Axis ratio is used here as a proxy for galaxy inclination. Disk galaxies with
very small axis ratios are likely heavily tilted from our perspective, in
contrast to face-on disks which would typically have axis ratios of $b/a\sim
0.9$ (Cho & Park, 2009; Maller et al., 2009). Highly-inclined star forming
galaxies suffer greatly from dust obscuration which can effect their colours
and magnitudes by making them appear redder and fainter than their intrinsic
properties (e.g., Conroy et al., 2010; Masters et al., 2010a; Morselli et al.,
2016; Kourkchi et al., 2019). This effect is important to consider when we are
making predictions for a disk galaxy like the Milky Way.

####  2.2.1 MPA-JHU Masses and Star Formation Rates

The MPA-JHU galaxy property catalogue provides total mass and star formation
rate estimates for galaxies within SDSS DR7. LNB15 developed an updated
version of this catalogue using the SDSS DR8 photometry. To summarise, the
masses are calculated by fitting stellar population synthesis model to the
galaxy’s photometry (instead of spectral features, as in Kauffmann et al.
(2003); Gallazzi et al. (2005)), similar to Salim et al. (2007). The star
formation rates are calculated by emission-line modelling based upon
Brinchmann et al. (2004) with some updates. In order to account for aperture
bias LNB15 follows the method of Salim et al. (2007) in calculating photometry
for the light that falls outside of the fibre and fitting stochastic stellar
population synthesis models to it. Thus each galaxy SFR measurement consists
of a combination of the SFR measured from inside and outside of the SDSS
fibre. For more details on these calculations refer to LNB15 Section 2.2.2 and
Brinchmann et al. (2004). In both cases a Kroupa initial mass function (broken
power law) is assumed (Kroupa & Weidner, 2003).

Each galaxy in our volume-limited sample is assigned a posterior probability
distribution function (PDF) of log stellar mass and log star formation rate,
as well as a corresponding cumulative distribution function (CDF) determined
by the posterior. Nominal values used for the $\log{\dot{\rm{M}_{*}}}$ and
$\log{\rm{SFR}}$ of each object are taken to be the mean value of the
posterior. The standard deviation ($\sigma$) on these values are calculated
from percentiles ($P$) in the CDF provided in the MPA-JHU measurements; we
take $\sigma=(P_{84}-P_{16})/2$, i.e., half the difference between the
$16^{\rm{th}}$ and $84^{\rm{th}}$ percentile values. This gives us estimates
of the stellar mass and star formation rates as well as their associated
errors for each galaxy in the sample.

The GSWLC-2 catalogues which we use for photometry also provide computations
of stellar masses and star formation rates. We chose to use the MPA-JHU
results instead in order to avoid any systematic effects that may result from
using the same measurements for our galaxy properties that were used to
determine these derived quantities. We have tested how our results change if
we use GSWLC-2 masses and star formation rates and find minimal differences,
as summarised in Section 4.2.4.

####  2.2.2 Simard et al. Bulge and Disk Decompositions

In order to characterise the light-weighted bulge-to-total ratios ($B/T$) and
disk scale lengths (Rd) for SDSS galaxies, we use the catalogue of Simard et
al. (2011). This work performed galaxy image decompositions for objects within
the Legacy area in SDSS via the GIM2D software package (Simard et al., 2002).

Specifically, we use the catalogue of fits based on composites of a Sersic
$n=4$ bulge and a pure exponential disk (Sersic $n=1$). For our analysis we
use the $B/T$ computed in $r$ band, as it is expected to be more stable than
$g$-based measurements. As a check we have performed our analyses using
$B/T_{g}$ instead of $B/T_{r}$ and found minimal differences (cf. Section
4.2.4). It is worth noting that Kruk et al. (2018) finds that these bulge+disk
decompositions can be somewhat inaccurate when applied to strongly-barred
galaxies, which may lead to some biases within our galaxy sample. However,
fits optimised for barred galaxies have been performed for only samples of a
few thousand objects, inadequate for our purposes, so we rely on the Simard et
al. (2011) results here. The Simard et al. (2002) bulge and disk
decompositions are derived using a $H_{0}=70$ km$\rm{s}^{-1}\rm{Mpc}^{-1}$.

####  2.2.3 Galaxy Zoo 2 Bar Presence

The presence of bars have been speculated to be correlated to the star
formation history within a galaxy. However, the sense of the effect is
unknown; a bar may be related to an increase in star formation (e.g., Hawarden
et al., 1986; Hawarden et al., 1996; Devereux, 1987; Hummel et al., 1990), no
effect (Pompea & Rieke, 1990; Martinet & Friedli, 1997; Chapelon et al.,
1999), or decreased star formation (e.g., Vera et al., 2016; Díaz-García et
al., 2020; Fraser-McKelvie et al., 2020), along with other potential effects
on colour and metallicity (see e.g., Masters et al., 2011). Regardless of the
source of any correlations, we wish to incorporate any possible differences
having a bar may cause when we determine the Milky Way’s SED, since our Galaxy
does exhibit clear evidence of a bar (e.g., Blitz & Spergel, 1991). The Galaxy
Zoo 2 (GZ2) catalogue (Willett et al., 2013; Hart et al., 2016) contains
identifications of detailed morphological features in disk galaxies, such as
spiral arms and bar presence. The galaxies classified in Galaxy Zoo 2 are a
sub-sample of the brightest and largest galaxies in the SDSS Main Galaxy
Sample.

The Galaxy Zoo projects are open public projects in which members of the
community identify whether they found a variety of galaxy features in the
images provided. In the Galaxy Zoo 2 catalogue the number of raw votes for
each morphological feature is weighted (to account for user consistency) and
adjusted to mitigate the impact of redshift-dependent biases, yielding a
corrected fraction of volunteers who identified a given morphological feature
for each galaxy. In our case we focus on the debiased fraction of volunteers
who identified a galaxy as having a bar, which we denote by $p_{bar}$
(following Galaxy Zoo labelling conventions; however, this is a fraction of
votes, not a probability). Above some threshold in $p_{bar}$ (after cuts in
related parameters) we can be confident that a given galaxy indeed hosts a
bar. Willett et al. (2013) developed the initial version of Galaxy Zoo 2 bar
thresholds; more conservative thresholds were later defined by Galloway et al.
(2015).

Often, when one uses Galaxy Zoo results it is necessary to consider responses
to previous questions that influence whether the question of interest is even
presented to the volunteers. For example, as described in Willett et al.
(2013), a voter is only asked “Is there a sign of a bar feature through the
centre of the galaxy?” if they first selected that the object has “features or
disk” when asked “Is the galaxy simply smooth and rounded, with no sign of a
disk?”, and then responded in the negative when asked “Could this be a disk
viewed edge-on?”. One would then identify galaxies that have received a large
fraction of “yes” votes as containing bars. It is worth noting that Willett et
al. (2013) states galaxies that receive fewer than $10$ net votes for a given
question may not have a reliable classification. Therefore to construct a bar
sample using the minimum allowances of Willett et al. (2013), using the
debiased vote fractions from Hart et al. (2016) (which was an improvement to
the debiasing of Willett et al. (2013)), one would use the cuts in the second
and third row of the third column of Table 3 in Willett et al. (2013) with
$p_{\rm{bar}}\geq 0.3$ as recommended by Galloway et al. (2015) and additional
voting number thresholds $N_{\rm{not\;edgeon}}\geq 10$, and $N_{\rm{bar}}\geq
10$.

However, here we do not want to consider only information on the relationship
between galaxy star formation history, structural properties, and photometric
properties that comes from those galaxies that most definitively host bars.
Rather, it is desirable for the training sample for our Gaussian process
regression model to include objects spanning a broad range of parameter space:
this yielded smaller net prediction errors in Milky Way photometric parameters
than when small, restricted samples are used for training. For additionally
discussion on this choice refer to Section 4.2.4.

To summarise, the primary set of data that we employ in this paper consists of
a cross-match between the original SDSS DR8 volume-limited sample reported in
LNB15, an updated version of the MPA-JHU catalogue (Brinchmann et al., 2004)
of masses and star formation rates in SDSS, the Simard et al. (2011)
morphological bulge-disk decomposition catalogue, the GSWLC-2 medium-deep
photometry catalogue (Salim et al., 2016, 2018), the Galaxy Zoo 2 catalogue
(Willett et al., 2013; Hart et al., 2016), and the DESI Legacy survey DR8
catalogue (Dey et al., 2019). These cross-matches were performed using the
$\tt{astropy}$ coordinates package. We required that objects be separated by
less than 1.5 arcseconds from a counterpart in the SDSS volume-limited sample
to be considered a match, and discard any galaxies that are not included in
all catalogues considered.

Our sample is smaller than the original volume-limited sample as a result. The
MPA-JHU catalogue contains all objects from the volume-limited sample, so we
begin with the same set of 124,232 galaxies utilised by LNB15. After matching
to Simard et al. (2011), 123,167 galaxies remain; some objects are lost due to
minor differences between the DR7 and DR8 SDSS catalogues. When we require
GSWLC-M2 measurements, we are left with only 60,857 galaxies, as roughly half
of the SDSS footprint is covered by medium-deep GALEX, 2MASS, and WISE
photometry (see Section 2.1.2). After matching to GZ2 29,836 galaxies remain,
a consequence of the brighter magnitude limit used to select GZ2 objects
compared to the original SDSS Main Galaxy Sample (see Section 2.2.3). We note
that due to the brighter magnitude limit of GZ2, our final sample is no longer
volume-limited. This would yield biased results if we were measuring aggregate
properties of Milky Way analogues. However, for our GPR methodology this only
modulates the density of our training set within parameter space, causing
larger prediction errors due to the sparser sampling, but not leading to a
bias.

Finally, matching to DESI Legacy leads to only a minor reduction in galaxy
number as it covers a super-set of the SDSS area with deeper photometry (but
different reduction algorithms). After these cross matches we then remove
objects with photometric values of “NaN”, infinity, or $-99$, which all
indicate missing photometry in one catalogue or another. This is only done
when necessary for the evaluation of the GPR. For our WISE $k$-correction
calculations (see Appendix C.1) we will exclude objects with large WISE
photometry errors in a given band, which we do not propagate into our main
galaxy sample.

The final data sample consists of 29,588 galaxies in total from redshift
$0.03<z<0.09$, which is publicly available at our catalogue GitHub. The
parameters which we will use to predict photometric properties are stellar
mass ($\rm{M}_{*}$), star formation rate (SFR), galaxy axis ratio as a proxy
for inclination ($b/a$), bulge-to-total ratio ($B/T$), disk scale length
($R_{d}$), and corrected bar vote fraction ($p_{\rm bar}$). Covariances
amongst these parameters are minimal, as discussed further in Appendix A.1; we
show joint distributions for these parameters in Fig. 12.

###  2.3 Milky Way Properties

A number of the properties of the Milky Way used in this study have been
derived using the Bayesian mixture model meta-analysis method first presented
in Licquia & Newman (2015). This technique combines information from multiple
measurements in order to obtain aggregate constraints on the Milky Way’s
properties, taking into account the possibility that individual measurements
could be incorrect or have their errors miss-estimated. This Bayesian method
is combined with Monte Carlo simulations in order to account for uncertainties
on the Sun’s measured Galactocentric radius, $R_{0}$; the Galactic exponential
disk scale length; and uncertainties in the local surface density of stellar
mass. The inferred mass of the Milky Way’s stellar disk depends upon all three
of these parameters.

Gravity Collaboration et al. (2019) recently obtained a greatly improved
geometric measurement of our distance from the centre of the Milky Way,
$R_{0}=8.178\pm 0.026$ kpc (Licquia & Newman (2015) used the value $8.3\pm
0.35$ kpc from Gillessen et al. (2009)). We have rerun the Bayesian mixture
model inference from Licquia et al. (2015) and Licquia & Newman (2016) with
this updated measurement of the Galactocentric distance to obtain updated
measurements of the Milky Way’s mass, bulge-to-total mass ratio, and disk
scale length; see those papers for all details of the data sets used and the
calculations (the star formation rate estimate from Licquia & Newman (2015) is
not affected by the value of $R_{0}$, so we adopt it unchanged). The resulting
Milky Way parameters used in this study are as follows:

  * •

Bulge Mass $\rm{M}_{*}^{B}=0.90\pm 0.06\times 10^{10}\;M_{\sun}$

  * •

Disk Mass $\rm{M}_{*}^{D}=4.58^{+1.18}_{-0.94}\times 10^{10}\;M_{\sun}$

  * •

Total Stellar Mass $\rm{M}_{*}=5.48^{+1.18}_{-0.94}\times 10^{10}\;M_{\sun}$

  * •

Star Formation Rate SFR$=1.65\pm 0.19\;\rm{M}_{\sun}\rm{yr}^{-1}$ (Licquia &
Newman, 2015)

  * •

Log Specific Star Formation Rate $\log{\frac{\rm{SFR}}{\rm{M}_{*}}}=-10.52\pm
0.10$

  * •

Bulge-to-total Mass Ratio $B/T$ $=0.16\pm 0.03$

  * •

Disk Scale Length $R_{d}=2.48^{+0.14}_{-0.15}$ kpc

The revised stellar mass estimate for the Milky Way is smaller than the
previous estimate, but by an amount that is significantly smaller than the
previously-estimated uncertainties.

By design these physical parameters are constructed such that they can be
directly compared to the extragalactic catalogues used to predict photometric
properties. For example, the stellar mass and star-formation rate estimates
for the Milky Way are constructed assuming a Kroupa initial mass function
(Kroupa & Weidner, 2003) and an exponential disk model, which is the same way
in which mass and star formation rates were calculated in the MPA-JHU
catalogue (Brinchmann et al., 2004) that we use in this analysis. The only
parameter where this is not fully the case is the bulge-to-total ratio, $B/T$.
For the Milky Way we can securely estimate only a mass-weighted value for this
quantity. In external galaxies, however, the mass-weighted $B/T$ is much more
difficult to obtain, and light-weighted measurements tend to be more reliable.
However, our predicted Milky Way properties do not change significantly when
we switch from $r$-band-based to $g$-band-based $B/T$ values, even though
mass-to-light ratios differ significantly between these bands, suggesting that
this is not a major issue.

Our analysis also depends on two additional parameters that are determined
independently of any Bayesian mixture model meta-analysis. These quantities
have values assigned to them based on their meaning in their respective
catalogues and our understanding of the Milky Way:

  * •

Axis ratio (inclination proxy) $b/a=0.9\pm 0.\textbf{}1$

  * •

Bar vote fraction $p_{\rm{bar}}=0.45\pm 0.15$

Galaxy inclination has a strong effect on colour and luminosity measurements
for disk galaxies. As mentioned in Section 2.2, dust alters the observed
colours and magnitudes of star-forming galaxies that are highly inclined or
edge-on. Our perspective within the Milky Way makes it somewhat equivalent to
being edge-on to us (though our position within the Galaxy, rather than
outside, does cause some differences). However, photometric properties are
most cleanly determined for those objects which are observed face-on.
Therefore, we predict the SED that would be observed for the Milky Way for
axis ratio values drawn from a uniform distribution spanning from $b/a=0.8$ to
$1.0$, consistent with the intrinsic axis ratios of spiral galaxy disks as
described by Maller et al. (2009) and in Section 2.2. Our results should
therefore correspond to the properties of our Galaxy if it were observed face-
on. While axis ratio is a good proxy for inclination, it is not a perfect
substitute. For example, van de Sande et al. (2018) found that disk galaxies
that are rounder ($b/a\sim 1$) tend to be older and therefore intrinsically
redder. This means that small biases could result from treating the Milky Way
as having a face-on axis ratio. However, there is no straightforward way to
avoid this, and the effect should be small compared to other sources of error.

The Milky Way exhibits clear evidence that it contains a bar (e.g., Blitz &
Spergel, 1991; Shen & Zheng, 2020). However, very few Galaxy Zoo 2 galaxies
have $p_{\rm{bar}}=1.0$, and those galaxies with the highest bar vote
fractions are expected to have very strong bars, which may not match our
Galaxy. Therefore we assume that in Galaxy Zoo 2 the Milky Way would have a
vote fraction above the threshold for defining a bar, but not a value higher
than the bulk of barred galaxies. Galloway et al. (2015) and Willett et al.
(2013) find that $p_{\rm{bar}}\geq 0.3$ serves as a reliable threshold between
bar presence and lack thereof. Galaxies with $0.3\leq p_{\rm{bar}}<0.5$ likely
have weaker bars while galaxies with $p_{\rm{bar}}>0.5$ likely have stronger
bars. Because the bar strength of the Milky Way’s bar as it would be
determined from outside our Galaxy is not well-constrained, we treat the bar
vote fraction for the Milky Way as uniformly distributed between
$p_{\rm{bar}}=0.3$ and $p_{\rm{bar}}=0.6$. Choosing a larger mean vote
fraction has small effect on our results, given the large range of fractions
considered (compared to the distribution of $p_{\rm{bar}}$ in GZ2).

##  3 Gaussian Process Regression for Predicting Milky Way Photometry

In this subsection we describe how Gaussian process regression can be used to
estimate photometric properties for the Milky Way. First we explain the need
to transition from Milky Way analogue-based methods to GPR when we consider
higher-dimensionality parameter spaces in Section 3.1. In Section 3.2 we
explain the basic concepts behind GPR, and in Section 3.2.3 we highlight the
fundamental differences between GPR and analogue galaxy methods. Section 3.2.1
describes the kernel used to set up our GPR, which guides how information is
propagated from training objects to predictions. We briefly describe the
computational limitations of the GPR implementation we are using in Section
3.2.2. Lastly, Section 3.2.4 investigates the contributions of various sources
of uncertainty to our GPR predictions.

###  3.1 Limitations of Using Analogue Galaxies

Using Milky Way analogues to predict the photometric properties of the Milky
Way, as was done in Licquia et al. (2015), has been a very useful methodology
but also has limitations. Of particular concern is the dramatic reduction in
MWA sample size that occurs as the number of parameters that must be matched
increases; we would like to move from the two parameters considered by Licquia
et al. (2015) ($\rm{M}_{*}$and SFR), to a total of six, adding $b/a$, $B/T$,
Rd, and $p_{\rm{bar}}$. Requiring that analogues be Milky Way-like in more
ways should reduce the spread in photometric properties of the resulting
sample, potentially enabling stronger constraints. There are only limited
correlations between these six parameters (cf. Fig. 12 in Appendix A.1), so
degeneracies between them are minimal: they each add new information. For
instance, galaxies with stellar mass and star formation rate matching the
Milky Way exhibit bulge-to-total ratios ranging from zero to one: structural
and star-formation history parameters carry distinct information. However,
while matching on additional galaxy parameters produces a population of
analogues that must each be closer in properties to the Milky Way, the
resulting MWA sample becomes much smaller (with as few as $\sim 5$ analogues
in the sample that are within $3\sigma$ of the Milky Way in all of the
properties considered, and none within $2\sigma$ for every aspect).

The reduction in the size of analogue galaxy samples as more parameters are
considered is illustrated in Fig. 1. We plot the total number of galaxies that
are within a given number of $\sigma$ for every Milky Way parameter considered
(where $\sigma$ represents the uncertainty in the Milky Way value for a given
property) as a function of the number of $\sigma$ used as a threshold. The
lightest shade of blue denotes the number of analogues within a given
threshold when only considering stellar mass ($\rm{M}_{*}$) and axis ratio
($b/a$). The consecutive additions indicated in the legend represent the
inclusion of the listed parameter in addition to all previous ones; i.e., we
consecutively incorporate star formation rate, disk scale length, bulge-to-
total ratio, and finally bar vote fraction. Hence, the darkest purple line
shows the number of analogues when using all six parameters (which we have
used in order of decreasing constraining power on Milky Way colours). The
$\sigma$ tolerances used for each parameter are the same Milky Way measurement
errors defined in Section 2.3, except for $p_{\rm{bar}}$. The significant
uncertainty we fiducially ascribe to bar vote fraction would cause the 5
parameter and 6 parameter lines to be degenerate with one another; to avoid
confusion we use $\sigma_{p_{\rm{bar}}}=0.05$ instead of $0.15$ when
constructing this plot.

Figure 1: Total number of galaxies within an allowed tolerance
($N_{\rm{gal}}<N_{\sigma}$) as a function of the maximum deviation away from
the Galactic values allowed for each parameter, in units of the uncertainty in
the Milky Way value of that parameter, $N_{\sigma}$. The lightest shade of
blue denotes galaxy counts when analogues are selected using only stellar mass
($\rm{M}_{*}$) and axis ratio ($b/a$). Consecutive parameter additions are
cumulative. Thus the darkest purple line incorporates all six parameters of
stellar mass, axis ratio, star formation rate (SFR), disk scale length
($R_{d}$), bulge-to-total ratio ($B/T$), and bar vote fraction
($p_{\rm{bar}}$). The $\sigma$ values used for each parameters are defined in
Section 2.3, except in the case of bar vote fraction, for which we use a
smaller error ($\sigma_{p_{\rm{bar}}}=0.05$ instead of $0.15$) for
illustrative purposes. Note that there are very few galaxies even within large
tolerance windows when using several parameters, e.g., only $\sim 10$ MWAs
within $3\sigma$ for 5 parameters simultaneously, and no objects that are
close to the Milky Way in all aspects, since the MWA method scales poorly with
increased dimensionality.

The previous work by LNB15 would most closely correspond to the +SFR (3
parameters), medium blue line. If one were to restrict to objects within $\pm
2\sigma$ of the Milky Way value for all parameters employed, there would be a
total of zero Milky Way analogue galaxies when using five or more parameters.
In the six-dimensional space that we employ below, there are only $~{}200$
galaxies that are within even $\pm 6\sigma$ of the Milky Way value in all six
properties; at that extreme, analogue samples would be selecting objects that
are not very close to the Milky Way at all. The lack of close analogues in
high-dimensional parameter spaces makes constraints on Milky Way properties
from the MWA method weak in that limit, with correspondingly large
uncertainties.

###  3.2 Gaussian Process Regression: A Powerful Method for Interpolation and
Prediction

To address the lack of Milky Way analogue galaxies in our multi-dimensional
parameter space, we have developed alternative methods for predicting Milky
Way properties based upon Gaussian process regression. In this sub-section we
summarise the basic properties of GPR relevant for this work. For in-depth
discussion, we refer the reader to Rasmussen & Williams (2006) and Görtler et
al. (2019).

Gaussian process regression (sometimes called kriging) is effectively a method
of interpolation where information from training data is accounted for by a
smooth and continuous weighting function, called a “kernel” or covariance
function. The joint probability distribution of the values of a Gaussian
process at any finite set of points in parameter space will be a multi-variate
Gaussian (with a number of dimensions set by the number of points in the set);
the kernel specifies how the covariance between points depends upon their
separation. The kernel should be a smooth and continuous function, with a
length scale (which governs how far information propagates from a given point)
that is optimised by training on the observed data. We can then predict what
the value and the uncertainty of the desired quantity would be at any
arbitrary point in space by applying this kernel to the training data. This is
in contrast to other supervised learning algorithms which typically make
single-valued, “point” predictions rather than predicting PDFs. It can be
shown that GPR yields the minimum variance out of any unbiased interpolation
method that depends only linearly on the training data; this makes GPR an
optimal interpolation algorithm.

For our application of GPR, the galaxy sample described in Section 2 will
serve as the training set. The six parameters we defined in Section 2.2
(stellar mass, star formation rate, etc.) serve as the features we will use
for prediction. Our goal is to determine an optimised mapping from these
physical parameters to a single output photometric parameter in our catalogue,
e.g., the $r$-band absolute magnitude, ${}^{0}M_{r}$.

Once our training data is selected we then go into the model-selection phase
of GPR, during which the mean function and covariance function (or kernel)
used for GPR are selected and tuned. We detail our selection of the covariance
function in Section 3.2.1. Effectively, the kernel determines how information
from a given training point will be propagated to make predictions at other
points in parameter space. Hyper-parameters describing the kernel are tuned at
this step to maximise the log-marginal likelihood of the training data. After
this step we consider the model to be “fit”.

Finally we enter the inference phase of GPR. At any point in our six-
dimensional parameter space, we can now determine the posterior probability
distribution for the parameter of interest by applying the kernel to the
training data. When evaluated at a single point in parameter space, a Gaussian
process corresponds to a 1-D Gaussian; we thus obtain both a predicted mean
for the property of interest (e.g., ${}^{0}M_{r}$) and the standard deviation
of the Gaussian describing its uncertainty. In our example we would pass in a
set of physical parameters measured for the Milky Way and obtain a predicted
value for the ${}^{0}M_{r}$ of the Milky Way, as well as the uncertainty in
that value.

In reality we do not only query the GPR at the mean measured Milky Way
properties presented in Section 2.3; rather, we perform random draws from the
PDFs describing $\rm{M}_{*}$, SFR, $b/a$, $B/T$, $R_{d}$, and $p_{\rm bar}$ in
order to incorporate the uncertainties in the Milky Way’s measured properties
into our analysis. For $\log{\rm{M}_{*}}$, $\log{\rm{SFR}}$, $B/T$, and
$R_{d}$ we assume a normal distribution. For $b/a$ and $p_{\rm{bar}}$ we draw
from uniform distributions, as described in Section 2.3.

We perform these random draws 1000 times (so that we have 1000 full sets of
Milky Way parameters). We then evaluate the GPR predictions for Milky Way
photometric properties at each of these points in parameter space. This gives
us a prediction and error estimate corresponding to each draw from the PDFs of
Milky Way characteristics. Thus we end up with 1,000 total predictions for
each Milky Way photometric property. Our mean prediction for the Milky Way in
a given photometric band corresponds to the arithmetic mean of all these
predictions.

In this work all Gaussian Process regression calculations have been done with
the Python scikit-learn Gaussian Process module,
$\tt{sklearn.gaussian\\_process}$ (Pedregosa et al., 2011). For details on the
implementation of this module, refer to Pedregosa et al. (2011) and the
scikit-learn documentation.

The following sub-sections provide more details on some aspects of our GPR
methods and their advantages.

####  3.2.1 Choice of Kernel for GPR

In this work we use a combination of two kernels for Gaussian Process
Regression: a Radial Basis Function kernel and a white noise kernel. The
Radial Basis Function (RBF) kernel decreases proportionally to $\exp{-\gamma
D^{2}}$, where $\gamma$ is a free parameter and $D$ is the Euclidean distance
between points; this kernel will cause the covariance between the predicted
values from GPR at different points in parameter space to decrease as a
Gaussian in distance as the separation between those points increases.

However, there is also scatter in galaxy photometric properties even for
objects measured to have the same physical properties. In order to capture
that, we also incorporate a white noise kernel, which models the spread in
values for the predicted property at a fixed point in parameter space with
normally distributed noise (Rasmussen & Williams, 2006). The net covariance
used for the Gaussian process regression is then the sum of the distance-
dependent covariance from the RBF kernel and the (diagonal) covariance matrix
corresponding to the white noise kernel.

For a given set of training data, there is a nearly endless number of
functions that can fit the given data points, each one a realisation of the
Gaussian process. The kernel creates a prior on the GP to constrain which
functions from that set are most likely to describe the parameter space. The
posterior is then determined using the training data values. Due to this
probabilistic approach, the Gaussian process provides both predicted values
and uncertainties at any points within the parameter space. Uncertainties due
to the finite training sample size and its distribution in parameter space and
those corresponding to intrinsic variation between training objects that have
the same physical parameters are both captured. In regions of parameter space
that are poorly constrained by the training data, the prediction uncertainties
are correspondingly larger.

The kernels we use in this work are available in the
$\tt{sklearn.gaussian\\_process.kernels}$ base class. For the white noise
kernel (WhiteKernel in scikit-learn) we initialise the noise level to be $1$;
similarly, we initialise the length scale for the RBF kernel to be $1$ for
each parameter. We opt to normalize the output photometric property to have
mean zero and variance one across the training set, which helps to ensure that
these initial guesses will have the right order of magnitude. The noise level
and length scales are then optimised and the regression model is built via the
$\tt{sklearn.gaussian\\_process.GaussianProcessRegressor}$ class. In our model
we allow the optimiser to restart 10 times in order to find the kernel
parameters that maximise the likelihood without being trapped in a false
maximum.

####  3.2.2 Optimizing Training Samples

The computation time and required memory for the $\tt{scikit-learn}$
implementation of GPR scales as the number of data points used to train the
model squared and cubed, respectively. As a result, we find that the maximum
training sample size we can use without running out of memory on the computers
used for this work is $\sim 6,000$; it is infeasible to train from our entire
catalogue when using this GPR implementation.

We have therefore tested the effects of either restricting to objects with
physical parameters within some tolerance of the MW fiducial values or
randomly selecting a subset of objects in order to reduce the training set
size. We have focused on the root mean squared error (RMSE) of predicted Milky
Way photometry for the $NUV$, $r$, and $J$ bands for this optimisation. We use
five-fold cross-validation for all the tests; i.e., we always train with 80%
of the data and test with 20%, but rotate what objects are used for training
and testing through the whole dataset, and only retain the values for an
object when it was in the test set. This provides unbiased estimates of the
RMSE for a training set 80% as large as the one we actually have. We find that
the combination that offered the lowest RMSE across all bands while keeping
computational time manageable was to randomly select 2,000 galaxies out of the
set of objects that are within $12\sigma$ of the Milky Way for every
parameters of interest. We therefore adopted this training strategy for all
results below.

####  3.2.3 Comparison to Results from Analogue Samples

GPR can provide more accurate predictions than many other techniques thanks to
its ability to leverage information from both nearby objects in the training
set as well as from more distant objects that characterise larger-scale
trends. In our application, this allows the GPR to map from the Milky Way’s
physical properties to its photometric properties much more accurately than if
we had only used the few objects that are similar to the Milky Way in all
respects (i.e., those which would be classified as MWAs) to inform the
mapping.

Fig. 2 illustrates the fundamental difference between how properties are
constrained by the Milky Way analogue selection method versus Gaussian process
regression. For simplicity’s sake we perform this comparison based on only
three parameters (stellar mass, star formation rate, and axis ratio, the same
ones utilised in LNB15), as the analogue method starts to break down when more
parameters are included. We also do not correct for Eddington bias in either
measurement (q.v. Appendix D) for simplicity. Objects included in a set of
MWAs based on 5,000 samples from the distribution of possible Milky Way
properties via methods equivalent to those from Licquia et al. (2015) are
depicted by the orange points. The analogues fall within a narrow range of
stellar mass, limiting the set of objects that contribute information. The
orange star represents the mean prediction for the Milky Way’s ${}^{0}(g-r)$
colour (${}^{0}(g-r)=0.682$) resulting from this set (derived via the
Hodges–Lehmann robust estimator, Hodges & Lehmann (1963)). The orange ellipse
depicts the $1\sigma$ confidence region for ${}^{0}(g-r)$ colour and stellar
mass. In contrast to the Milky Way analogues, the sample of galaxies used to
train the GPR are shown by purple points. These cover a much broader range of
parameter space than the MWAs. Similarly, the prediction for the Milky Way’s
${}^{0}(g-r)$ colour using GPR is shown by a purple star (with
${}^{0}(g-r)=0.668$), along with the corresponding $1\sigma$ confidence region
which is shaded in purple. When we have many analogues, both techniques yield
very similar results, but unlike MWAs the GPR technique still provides strong
constraints when we consider many parameters at once.

Figure 2: A comparison of galaxy samples and results from the Milky Way
analogues method versus Gaussian Process regression. For this analysis
analogues were selected based upon stellar mass, star formation rate, and axis
ratio, and we use the same parameters to predict colour via GPR. Orange points
represent Milky Way analogues selected by methods equivalent to Licquia et al.
(2015), while purple points indicate the galaxies used to train the Gaussian
process regression. The predicted mean Milky Way ${}^{0}(g-r)$ colour is
denoted by a star in the corresponding colour, and ellipses depict a $1\sigma$
joint confidence region for colour and stellar mass. Both methods yield a
similar predicted ${}^{0}(g-r)$ colour for the Milky Way when we have large
numbers of analogues. In contrast to the small window that the Milky Way
analogues lie in, the GPR utilises a wider variety of galaxies to capture
larger-scale trends.

####  3.2.4 Characterising Sources of Uncertainty

We can quantify the contributions of different sources of uncertainty to our
GPR estimates by changing how we perform the regression. We illustrate our
methods by evaluating how the prediction from a six-parameter Gaussian process
regression fit changes as a single parameter varies. The left panel of Fig. 3
shows one of the physical parameters being regressed from, star formation
rate, on the x axis and the target value, ${}^{0}(g-r)$ colour, on the y axis.
We chose this pair as galaxy star formation rate is expected to correlate with
galaxy colour well at fixed stellar mass.

First we isolate the scatter in colour at fixed properties; this corresponds
to the contribution of the white noise kernel to the covariance function of
the GPR. To determine the magnitude of this scatter we query the GPR at the
Milky Way’s fiducial physical properties to obtain a predicted PDF of
${}^{0}(g-r)$ at this point in parameter space from which we can draw samples.
The standard deviation of the colour of these samples corresponds to the
scatter encoded in the white noise kernel.

In the panel at the right of Fig. 3 we plot a histogram of 10,000 possible
${}^{0}(g-r)$ colours drawn from the GP’s predicted PDF, evaluating it at the
fiducial value of the Milky Way’s SFR. In the left panel the lavender star
denotes the mean predicted ${}^{0}(g-r)$ for the Milky Way with this model,
and the error bar corresponds to the standard deviation of the sample values.
This error bar therefore corresponds to the $1\sigma$ scatter in colour at
fixed properties. The grey shaded area corresponds to the $\pm 1\sigma$ band
for the GPR prediction of colour for a range of $\log{(\rm{SFR})}$ values. By
construction the half-width of the band must match the standard deviation of
samples from the PDF at fixed properties.

To instead isolate the contribution to errors resulting from the uncertainty
in Milky Way properties, we determine the distribution of mean predicted
colours evaluated at varying values of SFR drawn from the PDF for the Milky
Way. We perform 1000 draws from the fiducial MW $\log{\rm{SFR}}$ PDF in total
and evaluate the GPR mean predicted ${}^{0}(g-r)$ for each. In Fig. 3 the dark
purple points in the left panel shows the resulting predictions, which all
fall on a continuous curve by construction. The panel to the right shows the
histogram of this set of predictions in purple. We quantify the scatter in
colour attributable to uncertainties in the Milky Way’s measurements via the
standard deviation of the ${}^{0}(g-r)$ values at these 1,000 points.

To illustrate the full range of values obtained via GPR we plot ten samples
from the distribution of predictions for each of the thousand Milky Way SFR
draws as faint blue points in Fig. 3. These samples vary in colour due to both
the scatter in colour at fixed properties and due to the uncertainty in Milky
Way properties.

We can extend these same ideas in 1-D to evaluate the relative contributions
of uncertainties in Milky Way properties and of the scatter in properties at
fixed colour to the error in Milky Way ${}^{0}(g-r)$ colour, for any GPR model
of interest. The key difference from Fig. 3 is that, in order to quantify the
full scatter due to uncertainties in Milky Way properties, we allow all the
parameters to vary, not only SFR. We present the results of this analysis for
GPR models based on two to six physical parameters in Fig. 4. We use a stacked
bar plot to display the contribution to the variance from each error source,
where the x-axis is labelled according to the physical parameters used to
train the GPR (where additions are all cumulative, so entries further to the
right incorporate more parameters). Since independent errors will add in
quadrature, the contribution to the net variance from each factor is
proportional to the height of its bar. The variance due to the scatter at
fixed properties is shown in a lighter violet shade, while that resulting from
the uncertainty in Milky Way properties is shown in dark purple. The scatter
at fixed properties contributes to the majority of the error in the Gaussian
process regression for ${}^{0}(g-r)$ colour.

We have also evaluated the contribution to uncertainties resulting from the
finite size of the training sets used. This scatter is isolated by varying the
randomly-selected training sample 100 times. For each training sample we
evaluate the mean predicted value from GPR at the Milky Way’s fiducial
physical parameters. The contribution to uncertainties from the finite size of
training samples is obtained by calculating the variance of the GPR
predictions across all of the training sets. These values are minute: the
variances are an order of magnitude smaller than error attributed to the
scatter due to Milky Way measurements. Thus any contribution to errors
resulting from the finite training sample size is negligible.

We have performed the same error budget test for colours in the UV, near-IR,
and mid-IR. The results mirror those presented in Fig. 4: the errors are
dominated by scatter at fixed properties, followed by scatter from the Milky
Way measurements. In all cases errors attributable to finite training set size
are negligible compared to other sources. While the cumulative variance
decreases for every parameter added when we predict colours, this is not the
case for absolute magnitudes. In that case, the cumulative variance decreases
as we add parameters until the sixth parameter, $p_{\rm{bar}}$, is
incorporated. At that point, the variance increases and the scatter due to
finite training becomes more important. For this reason all absolute magnitude
predictions within this paper are performed using only 5 parameters, excluding
$p_{\rm{bar}}$.

Figure 3: A breakdown of uncertainty due to scatter at fixed properties and
scatter due to Milky Way measurement uncertainty. We plot results from a six
parameter Gaussian process regression trained to predict ${}^{0}(g-r)$ colour.
For this example we vary only one of the input parameters, the star formation
rate (SFR). The light violet colour corresponds to results when we use a fixed
training set and evaluate the GPR at fixed MW properties; this isolates the
variation attributable to the white noise kernel, corresponding to the scatter
in colour at fixed galaxy properties. The dark purple colour corresponds to
results when the training set remains fixed but the SFR value is drawn
randomly from the fiducial MW PDF, so that the contribution to the colour
error attributable to the uncertainty in the star formation rate of the Milky
Way can be evaluated (the other five parameters used to train the GPR are held
constant for simplicity for this example). In the left panel the shaded region
depicts the $\pm 1\sigma$ range predicted by the GPR fit, plotted out to $\pm
3\sigma$ of the Milky Way’s SFR. The predicted ${}^{0}(g-r)$ colour and
$1\sigma$ error for the Milky Way’s fiducial SFR (and hence only including
error at fixed properties) is shown by the star-shaped point and error bar in
light purple. The predicted ${}^{0}(g-r)$ colour from each of 5,000 draws from
the SFR PDF is shown by the purple points, which accounts for scatter due to
the SFR measurement. These predictions perfectly trace the GPR fitand tend to
fall within $\sim 3\sigma$ of the Milky Way’s SFR, by construction. For
reference we also show 10 samples from the GPR-predicted PDF for each of the
5,000 random SFR values as faint blue points. The distribution of these points
reflects both the scatter of galaxy colours at fixed properties and the
uncertainty in the MW measurements. Histograms of the colours for each sample,
whose distributions correspond to the scatter at fixed properties and the
scatter due to Milky Way SFR measurement uncertainties, are shown in the right
panel. It is evident that the spread in GPR predictions at fixed properties is
much larger than the scatter that results from uncertainties in the Milky
Way’s SFR. Figure 4: Contributions to the variance in restframe ${}^{0}(g-r)$
colour for Gaussian process regression employing varying sets of galaxy
physical parameters. The x-axis shows the set of parameters used to predict
colour; they increase cumulatively as we go from left to right, in order from
the most constraining to the least constraining parameter. The variance
decreases monotonically as the number of parameters increases. In every case
the scatter in colour at fixed properties dominates errors; uncertainties in
the physical properties of the Milky Way are sub-dominant. Contributions to
uncertainties due to having a finite training set are small enough to be
considered negligible, enough so that they would not be visible on this plot
if we included them.

While uncertainties in Milky Way characteristics will contribute to the random
errors in the derived photometric properties of our Galaxy, uncertainties in
the physical parameters of the training galaxies can cause systematic errors.
If the density of objects in parameter space varies quickly (with non-
negligible second or higher derivatives), objects will more often scatter from
well-populated regions of parameter space into sparser regions than vice
versa. The resulting systematic shift in the measured distribution of
parameters compared to the underlying distribution with no scatter is known as
Eddington bias.

In the context of this work, Eddington bias will lead to shifts in the colour
and luminosity predicted for the Milky Way. We derive corrections for
Eddington bias using methods similar to those of LNB15; we detail our
procedures in Appendix D. The estimated Eddington bias is subtracted off from
the GPR-predicted colours and luminosities for the Milky Way to produce our
final estimates for the Galaxy’s photometric properties and likewise the
uncertainty on the Eddington bias calculations is propagated into our final
error estimates. In general Eddington bias has small but nonzero effects on
our results ($<1\sigma$ for almost all parameters, as listed in Appendix A).

In the following sections, the errors on GPR results presented include the
contributions from scatter at fixed properties, uncertainties in Milky Way
properties, and uncertainty from Eddington bias.

####  3.2.5 Summary of the GPR Algorithm for Determining Milky Way Photometric
Properties

Here we summarise the steps taken to predict Milky Way photometric properties
via GPR. Our method proceeds as follows:

  1. 1.

Construct the training sample by restricting to objects within $12\sigma$ of
the Milky Way in all physical parameters considered and then randomly down-
sampling to 2,000 objects.

  2. 2.

Adopt the combination of a Radial Basis Function (RBF) kernel and a white
noise kernel as the covariance function to be used for GPR.

  3. 3.

Train the GPR using a single photometric property (normalised to have mean
zero and variance one) as the output or “y” value and the physical galaxy
parameters as the “x” values. This training will tune the hyperparameters of
the kernel.

  4. 4.

Perform 1,000 random draws from the PDFs that describe the fiducial Milky
Way’s properties. This will allow us to incorporate uncertainties in the Milky
Way measurements into our results.

  5. 5.

Use GPR to apply the optimised kernel to the training set and predict the
photometric property of interest. For each randomly-drawn set of physical
properties for the Milky Way, we obtain the mean prediction, predicted
variance, and a set of 1,000 values drawn from the GPR-predicted PDF
corresponding to that position in physical parameter space (which we refer to
as a set of samples).

  6. 6.

The mean photometric prediction for the Milky Way is then calculated as the
mean of the set of GPR output means at the position of each MW draw. The error
on the prediction is calculated as the standard deviation of the values from
the complete set of samples generated, allowing us to incorporate both
uncertainties associated with the scatter at fixed properties and errors
resulting from the uncertainties in MW properties.

The code used to construct the GPR is provided on our GP GitHub page for
public use here222https://github.com/cfielder/GPR-for-Photometry. At this site
we provide sample code for determining photometry estimates, addressing
systematics, and constructing an SED.

##  4 Results

Via GPR predictions for Milky Way photometric properties across the spectrum,
we can produce a comprehensive outside-in portrait of the Milky Way SED,
allowing comparisons to the colors and luminosiities of other galaxies. In
this section we apply a variety of diagnostics from the literature, such as
colour-luminosity, colour-mass, and colour-colour diagrams, in order to assess
how the Milky Way compares to the broader population. We also construct a
multi-wavelength SED for the Milky Way and compare our results to templates
from the literature.

###  4.1 The Milky Way Compared to the Broader Galaxy Population

As discussed in Section 3, we have predicted the Milky Way colours and
luminosities based upon the six parameters of stellar mass ($\rm{M}_{*}$),
star formation rate (SFR), axis ratio ($b/a$), bulge-to-total ratio ($B/T$),
disk scale length ($R_{\rm{d}}$), and bar vote fraction ($p_{\rm bar}$). In
the following colour diagrams, all magnitudes and colours are presented as
rest-frame AB magnitudes (evaluating all passbands at redshift zero).

Our quantitative results are summarised in Table 2-Table 4. The values
provided correspond to the mean rest-frame predictions based upon the Gaussian
process regression derived via the methods presented in Section 3, and have
been corrected for Eddington bias as described in Appendix D. Colours and
magnitudes are all calculated independently of one another. For example, we
use GPR to predict ${}^{0}(g-r)$ galaxy colour directly, as opposed to
deriving this value by subtracting the predicted ${}^{0}M_{r}$ from the
predicted ${}^{0}M_{g}$. For SDSS photometry, our derived colours are based
upon model magnitudes, as these yield the most accurate colour estimates for
SDSS galaxies; however, the absolute magnitudes provided are based upon cmodel
magnitudes, as those most accurately represent the total brightness of an
object.

Log-spaced density contours corresponding to the cross-matched galaxy sample
of $29,836$ galaxies described in Section 2.1 and Section 2.2 are plotted in
grey-scale on all of the following colour-based diagrams. We also overlay red
and blue ellipses which denote the rough locus of the red sequence and blue
cloud, respectively, in each plot. These shadings are intended to guide the
eye and should not be interpreted in a quantitative manner. In a corner of
each plot we provide error bars that correspond to the mean uncertainties in
each galaxy property being plotted for the training set.

In each diagram we also show the locations of the 36 red spiral galaxies
selected in Masters et al. (2010b) that overlap with our cross-matched sample
(out of 294 in the original catalogue). This sample of objects was selected
based upon their colour, presence of spiral features, and shape/structural
parameters from SDSS. They are required to have colour
${}^{0}(g-r)>0.63-0.02(^{0}M_{r}+20)$, overlapping the blue edge of the red
sequence. They are also selected to have a spiral likelihood
$p_{\rm{spiral}}\geq 0.8$ in the prescription of Bamford et al. (2009), and
are required to have visible arms in Galaxy Zoo 1 $p_{\rm{CW}}>0.8$ or
$p_{\rm{ACW}}>0.8$ (Lintott et al., 2011) in order to ensure they have spiral
morphology. These objects are also selected to be approximately face-on
(equivalent to an axis ratio requirement $b/a>0.63$), as dust reddening is
expected to have a substantial impact on the apparent colours of spirals
(Masters et al., 2010a). However, in that paper the axis ratio values wre
calculated via $r$–band isophotal measurements, while ours are determined from
an exponential profile fit. Therefore we apply a profile-fit-based cut of
$b/a>0.6$ to this sample to enable a more direct comparison to our face-on
results for the Milky Way. Finally, Masters et al. (2010b) requires that the
red spiral sample contain galaxies with an SDSS $f_{\rm{deV}}\leq 0.5$, where
$f_{\rm{deV}}$ is defined as the weight of the de Vaucouleurs profile in the
best-fit linear combination with the exponential profile matched to the
object’s image. This ensures that S0 galaxies do not contaminate the sample,
although they are already only a small percentage of the GZ1 sample.

The resulting red spiral sample is represented by red points in our plot. We
overlay the positions of these objects in each parameter space to help assess
the consistency of the inferred properties of the Milky Way with this
population. Two objects whose ${}^{0}(g-r)$ colours in the cross-matched
catalogue differed by $>0.1$ mag from the photometry used in Masters et al.
(2010b) due to changes in SDSS pipelines were excluded.

####  4.1.1 Optical Colours

We first present results at optical wavelengths, as they allow us to compare
directly to previous work done with Milky Way analogues in LNB15. We focus on
the SDSS $ugriz$ bands (cf. Section 2.1.1). Fig. 5 presents predictions for
Milky Way optical optical colours as a function of stellar mass ($\rm{M}_{*}$)
in solar mass units.

The upper panel shows ${}^{0}(u-r)$ colour and the lower panel shows
${}^{0}(g-r)$ colour versus mass. Both panels have overlaid dashed reference
lines which can be used to distinguish general regions of the diagrams. The
top portion contains galaxies that are on the red sequence, while the middle
portion contains green valley galaxies, and the lower portion corresponds to
the blue cloud. The dashed lines bracketing the green valley in the upper
panel correspond to ${}^{0}(u-r)=-0.24+0.25\rm{M}_{*}$ and
${}^{0}(u-r)=-0.75+0.25\rm{M}_{*}$ (Schawinski et al., 2014). In the lower
panel the plotted lines correspond to ${}^{0}(g-r)=0.6+0.06(\rm{M}_{*}-10)$
and ${}^{0}(g-r)=(0.6+0.06(\rm{M}_{*}-10))+0.1$ (Mendel et al., 2013).

Previous results from LNB15 are plotted in orange. The star represents the
mean prediction and the ellipse encompasses the $1\sigma$ confidence region.
We remind the reader that these constraints were determined based only on
stellar mass and star formation rate, along with a cut on axis ratio. In
comparison, the results of the 6-parameter Gaussian process regression are
plotted in purple. In red we plot the GPR result evaluated with an axis ratio
of $b/a=0.3$ rather than 0.9, to illustrate the impact that the assumed
inclination has on the inferred SED. The GPR confidence regions are calculated
from the covariance between the the samples drawn from the regression
predictions; the distribution of these samples incorporates both uncertainties
in Milky Way properties and scatter in colours at fixed properties (cf.
Section 3.2.4). In the lower right corner of each panel we show the mean error
in optical colour and log stellar mass amongst the galaxies in our final
sample. Per-object uncertainties in the optical colours account for roughly
half of the total scatter in our GPR colour prediction for a face-on Milky
Way. In contrast, the average error in stellar mass in the training sample
does not affect the uncertainty in the stellar mass of the Milky Way (it will,
however, contribute to Eddington biases, as discussed in Appendix C. Note that
LNB15 did not use the same stellar mass as we do, reflecting our updated
estimate for the mass of the Milky Way (cf. Section 2.3).

In both ${}^{0}(u-r)$ and ${}^{0}(g-r)$ our results are consistent with,
though marginally redder than, those reported in LNB15. This is no surprise as
we do not expect the Milky Way to move far in optical colour space when
constraints tighten. Even when we make predictions for a much more inclined
Milky Way our results do not change dramatically in the optical, with shifts
well within the uncertainties in both our face-on results and those from
LNB15, although the colour does become marginally redder as expected. In the
optical the Milky Way appears to lie in the “saddle” of the the galaxy-colour
bimodality, implying that the Milky Way is redder than the average spiral
galaxy in the local Universe in the optical bands. That said, if one were to
only consider spiral galaxies of similar mass to our Galaxy, the Milky Way is
not as unusually red as it would be if compared to lower-mass spiral galaxies.

Figure 5: Rest-frame optical colour as a function of stellar mass. The grey-
scale, log-spaced contours depict the density of $0.03<z<0.09$ galaxies in our
cross-matched sample, without any limiting magnitude applied. The dashed grey
lines correspond to divisions of the galaxy population used in the literature.
We expect “red sequence” galaxies to be in the upper portion of each plot,
with the “blue cloud” corresponding to the bluest colours. The region between
the two lines corresponds to the “green valley” population. For ${}^{0}(u-r)$
we use the divisions of Schawinski et al. (2014) and in ${}^{0}(g-r)$ we
follow Mendel et al. (2013) (cf. Section 4.1.1). The results and $1\sigma$
confidence region from LNB15 are marked in orange. Our results from applying
GPR to all six galaxy physical parameters considered are marked in purple; the
$1\sigma$ region is determined by the covariance between Gaussian process
samples. For comparison, in lighter red we show the six-parameter GPR result
obtained when setting the axis ratio of the Milky Way to $b/a=0.3\pm 0.1$
rather than 0.9. The stellar masses differ between our prediction and LNB15
due to the updates to the mass estimate for the Milky Way described in Section
2.3. Red points correspond to members of the red spiral galaxy sample of
Masters et al. (2010b). In the lower right corner of each panel we depict
error bars representative of the mean uncertainties for galaxies in the
comparison sample. Our results are consistent with LNB15 and indicate that at
optical wavelengths the Milky Way is redder than the typical star-forming
spiral galaxy, in addition to being more massive.

The green valley (and by extension the galaxy colour-bimodality) has been used
as a basic tool to distinguish transitional galaxies from the general galaxy
population. Transitional galaxies have lower specific star formations rates
(sSFR = ${\rm{SFR}/\rm{M}_{*}}$) than a star-forming galaxy of the same mass;
specific star formation rate can be used as a proxy for the evolutionary state
of a galaxy and its star forming history. Salim (2014) defines the
transitional region in sSFR space to be below the sSFR of massive Sbc
galaxies, as these Sbc’s are the earliest galaxy type expected to proceed with
regular star formation free of quenching, but above the sSFR at which galaxies
appear to no longer be star forming in the UV. As described in that work, this
range corresponds to $-11.8<\log{(\rm{SFR}/\rm{M}_{*})}<-10.8$. Note that here
and throughout this paper $\log$ refers to the base 10 logarithm.

In Fig. 6 we plot the same rest-frame colours as in Fig. 5 but as a function
of log specific star formation rate. The vertical dashed lines denote the
transitional region in $\log{\rm{SFR}/\rm{M}_{*}}$, as defined by Salim
(2014). The region with log sSFR above -10.8 corresponds to galaxies that are
actively forming stars while objects with log specific star formation rate
below -11.8 are quiescent; transitional objects are between them. In the lower
panel the green horizontal lines correspond to the green valley definition of
Lackner & Gunn (2012) evaluated with the predicted $r$ band absolute magnitude
for the Milky Way (${}^{0}M_{r}=-20.65$). Galaxies residing within this range
in ${}^{0}(g-r)$ are expected to reside within the green valley. Galaxies
above this designation are expected to lie in or near the red sequence, and
galaxies below are expected to lie in or near the blue cloud.

Figure 6: Rest-frame optical colour as a function of log specific star
formation rate ($\log{(\rm{SFR}/\rm{M}_{*})}$) in units of $\rm{yr}^{-1}$. As
before the grey-scale log contours depict the density of the parent cross-
matched sample, red points correspond to the red spiral sample of Masters et
al. (2010b), the orange star and ellipse correspond to the results of Licquia
et al. (2015), and the purple star and ellipse are the results of our GPR
analysis. The vertical dashed grey lines designate divisions of galaxy
populations according to their specific star formation rates, following the
definitions of Salim (2014): quiescent galaxies are at left, transitional
objects are in the middle, and star-forming objects are at right. The Milky
Way lies on the star-forming side of these divisions. The green horizontal
lines in the bottom panel correspond to the green valley definition of Lackner
& Gunn (2012), evaluated at the $r$-band absolute magnitude of the Milky Way.
According to the prescription by Mendez et al. (2011) the green valley has a
width of 0.1 in ${}^{0}(g-r)$, leading to the limits shown here here. The
Milky Way mean value falls above this “green valley” region. While the optical
colour of the Milky Way is redder than most star-forming galaxies in the local
Universe, based on its specific star formation rate the Milky Way would not be
considered a transitional galaxy.

Based on its specific star formation rate, the Milky Way must lie within the
star-forming population, rather than in the transitional range. While one
might expect an object that meets optical definitions of the green valley to
have a transitional sSFR this is not necessarily true, as galaxies of
different evolutionary states can share the same optical colour (see e.g.,
Cortese, 2012; Salim, 2014).

In ${}^{0}(g-r)$ if we take the green valley to be 0.1 in width, as defined by
Mendez et al. (2011), the Galaxy would either fall within or be redder than
the green valley, consistent with the results shown in Fig. 5. Despite ongoing
star formation, more massive spiral galaxies tend to be redder in the optical
than their lower-mass counterparts (e.g., Masters et al., 2010b). However, our
estimated properties for the Milky Way lie in the middle of the distributions
of colours, masses, and sSFRs of the red spiral sample from Masters et al.
(2010b); these objects are redder in the optical than is typical for even the
most massive spirals.

As these plots exemplify, differentiating between galaxy populations based
only on optical photometry is challenging. For instance, in the lower panel in
Fig. 6 we can see that red sequence, transitional, and star-forming objects
can all have colours of ${}^{0}(g-r)\sim 0.7$. In Fig. 5 the blue cloud
becomes difficult to distinguish from the red sequence at high masses as the
most massive spirals have lower sSFRs and, therefore, redder colours. In the
following subsections we investigate constraints on the colour of the Milky
Way at UV and IR wavelengths where galaxy populations may separate more
clearly.

####  4.1.2 UV Colours

We utilise far-ultraviolet (FUV) and near-ultraviolet (NUV) photometry from
GALEX provided in the GSWLC-M2 catalogue (Martin & GALEX Team, 2005; Salim et
al., 2016, 2018), as discussed in Section 2.1.2. Thermal emission from massive
stars with lifetimes $<100$ Myr peaks at ultraviolet wavelengths, while lower-
mass, longer-lived stars play a larger role at optical wavelengths (Salim,
2014; Tuttle & Tonnesen, 2020). Because UV radiation is produced by short-
lived but high-luminosity stars it provides a sensitive indicator of _recent_
star formation. As a result, UV photometry can more clearly differentiate
star-forming from quiescent galaxies than optical measurements can.

Much as in Section 4.1.1, we can use our GPR results to place the Milky Way on
UV-based diagnostic diagrams from the literature. In Fig. 7 we plot
${}^{0}(\rm{FUV}-r)$ and ${}^{0}(\rm{NUV}-r)$ UV-optical colours versus
specific star formation rate. The contours and vertical reference lines shown
are defined in the same way as in Fig. 6. For ${}^{0}(\rm{NUV}-r)$ we show
horizontal lines corresponding to the “green valley” definition of Salim
(2014), bounded at $4<\;^{0}(\rm{NUV}-r)<5$. As before, our face-on MW
prediction and the corresponding $1\sigma$ confidence region are plotted in
purple and the inclined MW prediction is plotted in red; unlike in the
optical, there are no previous estimates of Milky Way properties in this space
that we could plot. Much as in the optical, uncertainties in the UV photometry
for individual objects account for roughly half of the total scatter ascribed
to our Milky Way UV predictions.

Figure 7: As Fig. 6, but for UV-optical colours. The vertical reference lines
come from Salim (2014) where objects with a $\log{\rm{sSFR}}>-10.8$ are
actively forming stars, objects with $\log{\rm{sSFR}}<-11.8$ are quiescent,
and objects in between are considered transitional. In ${}^{0}(\rm{NUV}-r)$ we
also plot the bounds of the UV-optical “green valley” as defined in Salim
(2014). Our predictions for the Milky Way show that it lies on the star
forming side of the transitional regions, and is most likely on the blue side
of the UV-optical green valley border. Our Galaxy lies in a similar region of
these colour spaces to the red spiral sample of Masters et al. (2010b).

Compared to typical star-forming galaxies in the local Universe, the Milky Way
has redder than average UV colours and lower than average sSFR. The Milky Way
appears to lie on the blue side of the ${}^{0}(NUV-r)$ green valley border, in
contrast to its location in the optical (cf. Fig. 5 and Fig. 6). This reflects
the limited discriminating power of optical colour; the green valley is only
0.1 mag wide in ${}^{0}(g-r)$, allowing objects to easily scatter over its
borders due to even small photometric errors or inclination effects, but it
spans an entire magnitude in ${}^{0}(NUV-r)$. A more inclined Milky Way is
predicted to be notably more red in the UV than in the optical, so much so
that it could be consistent with the UV green valley in colour.

As in the optical, our estimates for the Milky Way in the UV-sSFR plane lie in
the middle of the Masters et al. (2010b) red spiral population. Red spiral
galaxies tend to lie outside of the UV green valley as they have star
formation rates comparable to typical blue spirals of the same mass (Cortese,
2012), and UV colour is more sensitive to recent star formation rate than the
optical is.

####  4.1.3 Infrared/WISE Colours

The infrared data for our galaxy sample originates from the 2MASS and WISE
surveys, as included in the GSWLC-M2 (Salim et al., 2016, 2018) and DESI
Legacy catalogues (Dey et al., 2019), respectively (cf. Section 2.1). Similar
to in the ultraviolet, the infrared brightness of a galaxy is sensitive to
recent star formation due to re-emission of UV photons absorbed by dust. The
infrared colours of galaxies also exhibit a colour bi-modality, but star-
forming galaxies exhibit redder IR colours than the passively-evolving
population, rather than bluer. Instead of the “green valley,” the region
between the star-forming and quiescent populations in the IR is commonly
referred to as the infrared transition zone (IRTZ), following Alatalo et al.
(2014).

In 8(a) and 8(b) we plot WISE colour-colour diagrams for the cross-matched
galaxy sample in addition to the GPR prediction for the Milky Way. The Milky
Way colour is poorly constrained in some WISE bands due to the lower signal-
to-noise of these detections. If we compare our covariance ellipse in the
${}^{0}(W1-W2)$ direction to the average errors in the photometry (lower right
error bar) the photometric errors account for a modest fraction of the total
uncertainties in our GPR prediction. However, in ${}^{0}(W2-W3)$ and
${}^{0}(W3-W4)$ errors in the photometry for individual objects dominate the
estimated uncertainties in the Milky Way GPR predictions. The vertical lines
in these plots designate the IRTZ from Alatalo et al. (2014); objects with
${}^{0}(W2-W3)>0.565$ in AB magnitudes correspond to late-type galaxies, those
with ${}^{0}(W2-W3)<-1.035$ are early-type galaxies, and in between lies the
IRTZ. Magnitudes were converted from Vega to AB magnitudes via the
prescription of Jarrett et al. (2011).

(a)

(b)

Figure 8: WISE colour-colour diagrams for both our parent sample (log density
contours) and the predicted results for the Milky Way from Gaussian process
regression. Reference lines from Alatalo et al. (2014) designate the infrared
transition zone ($1.035<^{0}(W2-W3)<0.565$; IRTZ). The Milky Way appears to
lie on the star forming side of the IRTZ, much closer to the median colours of
typical spiral galaxies, in contrast to the UV and IR. Again our Galaxy lies
in a similar region of these colour spaces to the red spiral sample of Masters
et al. (2010b).

As before, we show predictions for the Milky Way if it were approximately
face-on or more steeply inclined by evaluating the GPR at different axis
ratios. An inclined Milky Way appears the be most notably different in the W3
band, which traces prominent dust emission features, particularly those
associated with polycyclic aromatic hydrocarbons (Wright et al., 2010). It
appears that an inclined spiral galaxy would be measured to be more IR bright
compared to a face-on counterpart matching it in all other ways; this may
represent a systematic effect related to data processing, since the dust
emission would be expected to be optically thin. We find similar results in
Section 4.2.3.

In both diagrams the prediction for the Milky Way lies on the star-forming
side of the infrared transition zone. If we compare 8(a) to the classification
scheme in Figure 12 of Wright et al. (2010) (note that this requires
converting our AB magnitudes to Vega magnitudes), the Milky Way lies within
the region of colour space they label as typical for spiral galaxies, as would
be expected. Similarly, the classification scheme of Figure 11.b of Jarrett et
al. (2017) would place the Milky Way in the intermediate disk region,
consistent with the expectation for a massive spiral galaxy. Intermediate disk
objects are thought to be in transition towards being quenched due to star
formation rates that are decreasing with time; our estimate for the Milky
Way’s star formation rate is slightly below average for a spiral galaxy of the
same mass, consistent with this picture.

The infrared results mirror what we find from the UV: the Milky Way is still
forming enough stars to appear bright in the IR due to re-emission from dust.
Galaxies are expected to transition in the optical before they do in the
infrared (Alatalo et al., 2014; Tuttle & Tonnesen, 2020). If we follow the
narrative of Alatalo et al. (2014) and Smethurst et al. (2015), the Milky Way
may be in the early transition phase from star-forming to quiescent. It is
brighter and redder in the optical compared to the typical star-forming
galaxy. In the UV the Milky Way is on the blue side of the green valley but
near it. In the IR the Milky Way’s inferred colour is more typical for a star-
forming galaxy, though uncertainties are substantial. This would track with
the expectation that a galaxy transitions in the optical before it does in the
IR.

As before, the inferred colours of the Milky Way in the mid-IR are consistent
with the range of values for the Masters et al. (2010b) red spiral sample,
though some WISE bands are not well constraining due to the low signal-to-
noise of the underlying measurements used for prediction. Based on this, it
remains plausible that the Milky Way is a part of the red spiral population.
We discuss how the Milky Way’s colours compare to other galaxy populations
further in Section 5.2.

###  4.2 The Multiwavelength Spectral Energy Distribution of the Milky Way

Thus far, we have focused on predictions for a single Milky Way colour at a
time. However, we can assemble colour information across all passbands to
construct a spectral energy distribution (SED) for the Galaxy. SEDs, which
quantify the the total energy of emitted photons as a function of wavelength
or frequency, are valuable tools in the study of galaxies. Many physical
characteristics of galaxies can alter their SEDs - the age of their stellar
population, stellar abundances, gas and dust content, inter-stellar medium
(ISM) chemistry, details of star formation history, and the presence of an AGN
can all leave distinct signposts that give observers insight into the
formation and evolution of a given galaxy (see e.g., Silva et al., 2011).
Because these effects each tend to alter the SED at specific portions of the
spectral range, with broad enough wavelength coverage and detailed enough
spectral information one can disentangle the dominant processes in a given
galaxy.

Detailed modelling of the Milky Way’s SED will be the focus of a follow-up
paper. In this work, we will present a proof-of-concept for a GPR-constructed
SED for the Milky Way and provide an initial analysis of its properties. In
the following sub-section we outline our GPR-based methods for determining the
SED of the Galaxy before presenting quantitative results and assessing the
effects the galaxy physical parameters used for prediction each has on the
SED.

####  4.2.1 Algorithm for Calculating the SED for the Milky Way

We work in frequency ($\nu$) space instead of wavelength ($\lambda$) space
when calculating the SED of the Milky Way, as spectral energy distributions
are most typically presented in units of energy per unit frequency. Our
algorithm for calculating the SED proceeds as follows:

  1. 1.

Estimate ${}^{0}M_{r}$ and colours for the Milky Way –
Our SED is calculated in reference to the $r$-band. Therefore using the GPR
(described in Section 3) we predict the $r$-band AB absolute magnitude
(${}^{0}M_{r}$) and all colours with restframe $r$ as the reference band;
i.e., ${}^{0}(x-r)$ where $x$ spans from FUV to W4 (e.g.,
${}^{0}(FUV-r),...,^{0}(W4-r)$).
Eddington bias is subtracted off separately from our predicted colours and
${}^{0}M_{r}$ before we combine them. Similarly, the uncertainty in the
Eddington bias is added in quadrature to the uncertainty in the GP
calculations (which incorporates both scatter at fixed properties and errors
due to the Milky Way property uncertainties; cf. Section 3.2.4). For details
on the Eddington bias calculations refer to Appendix D.

  2. 2.

Calculate flux ratios –
We calculate the flux in each band $f_{\nu,x}$ relative to the flux in the
$r$-band $f_{\nu,r}$ via the relation

$\log{\left(\frac{f_{\nu,x}}{f_{\nu,r}}\right)}=\frac{{}^{0}(x-r)}{-2.5},$ (1)

where ${}^{0}(x-r)$ is the Eddington bias-corrected colour in the $x$ band
compared to the $r$ band, which we have predicted via GPR.

  3. 3.

Calculate luminosity –
From the $r$-band absolute magnitude combined with the flux ratios it is
straightforward to convert to luminosity. We first calculate the $r$ band
luminosity as

$\log{(L_{\nu,r})}=\frac{(^{0}M_{r}-34.04)}{-2.5}\;\left[\log{\left(\frac{W}{Hz}\right)}\right],$
(2)

where ${}^{0}M_{r}$ is the Eddington bias-corrected $r$-band absolute
magnitude for the Milky Way obtained via GPR. This formula is derived via the
relation for converting flux to AB magnitudes in combination with the area of
a 10 pc radius sphere to convert flux to luminosity. The luminosity in any
other band can then be calculated via the relation

$\log{(L_{\nu,x})}=\log{\left(\frac{f_{\nu,x}}{f_{\nu,r}}\right)}+\log{(L_{\nu,r})}.$
(3)

We can then add log frequency ($\log{\nu}$) to obtain $\log{\nu L_{\nu,x}}$.

  4. 4.

Calculate errors –
To convert uncertainties in magnitudes and colours to uncertainties in
$\log{\nu L_{\nu,x}}$ we make use of propagation of errors. First to determine
$\sigma_{\log{\nu L_{\nu,r}}}$ we calculate the partial derivative of
$\log{\nu L_{\nu,r}}$ (which is Equation 2 $+\log{\nu}$) with respect to
${}^{0}M_{r}$. This yields

$\sigma_{\log{\nu L_{\nu,r}}}=0.4\sigma_{{}^{0}M_{r}}.$ (4)

Errors in the other bands are calculated in a similar manner, but they depend
on both the error in colour and the error in ${}^{0}M_{r}$:

$\sigma_{\log{\nu
L_{\nu,x}}}=0.4\times\sqrt{\sigma_{{}^{0}(x-r)}^{2}+\sigma_{{}^{0}M_{r}}^{2}}.$
(5)

In the plots that follow we do not plot the contribution to errors from
$\sigma_{{}^{0}M_{r}}$ as it is fully covariant across all bands; as a result,
when templates are normalised to match the observed SED, any error in
${}^{0}M_{r}$ would simply change the normalisation. We do provide its value
for reference. Thus the error bars presented in the Milky Way SED plots are
equivalent to $0.4\sigma_{{}^{0}(x-r)}$.

The colour predictions used to derive the luminosities used for our SED are
provided for reference in Table 2 and the value of ${}^{0}M_{r}$ is provided
in the absolute magnitude table Table 4. We present our estimated luminosities
and associated uncertainties, incorporating Eddington bias corrections in each
case, in Table 1.

Passband | $\lambda_{eff}\;[\mu m]$ |  $\log{\nu L_{\nu}}$ [$\log{\rm{W}}$] |  $\sigma_{\log{\nu L_{\nu}}}$ [$\log{\rm{W}}$]
---|---|---|---
FUV | 0.155 | 35.53 | 0.20
NUV | 0.2275 | 35.63 | 0.20
u | 0.354 | 36.01 | 0.06
g | 0.4750 | 36.44 | 0.02
r | 0.622 | 36.62 | 0.09
i | 0.763 | 36.66 | 0.01
z | 0.905 | 36.70 | 0.02
J | 1.25 | 36.67 | 0.04
H | 1.65 | 36.54 | 0.05
Ks | 2.15 | 36.39 | 0.06
W1 | 3.368 | 35.59 | 0.08
W2 | 4.618 | 35.58 | 0.10
W3 | 12.082 | 35.61 | 0.20
W4 | 22.194 | 35.49 | 0.30
Table 1: The passbands and corresponding power and uncertainties for the
predicted SED of the Milky Way, as plotted in Fig. 9. These values have
already had Eddington bias subtracted out.

####  4.2.2 Interpreting the SED of the Milky Way

In Fig. 9 we present our full predicted spectral energy distribution for the
Milky Way, along with a variety of empirical template galaxy spectra from the
literature. We plot $\log{\nu L_{\nu}}$ (the power emitted per log interval in
frequency) on the vertical axis, in units of log Watts. The horizontal axis
corresponds to restframe wavelength in units of $\mu m$; each photometric band
used is labelled at its effective wavelength along the top of the plot. The
black open circles represent the estimates for the Milky Way’s luminosity
along with their associated errors, calculated as described in Section 4.2.1.
We also show the error bar corresponding to the uncertainty in ${}^{0}M_{r}$
near the bottom of the 9(a) for reference. The numerical values for the Milky
Way SED corresponding to the plotted points are provided in Table 1.

(a)

(b)

(c)

(d)

Figure 9: The predicted SED for the Milky Way from a six-parameter Gaussian
process regression, depicted by black open circles in all panels; the method
of calculation is described in Section 4.2.1. (a) We compare to empirical
galaxy templates from Benítez et al. (2004) (re-calibrated from Coleman et al.
(1980)) normalised to match the Milky Way in the $r$-band. The Milky Way SED
is consistent with the ’Sbc’ template from this set, which is labelled as such
because it was originally based on the average of the SEDs of two blue
galaxies of morphological type Sbc. (b) We compare the SED of the Milky Way to
the most closely-matched templates from Brown et al. (2014), which are based
on spectra and model fits to bright nearby galaxies. We show all templates
whose $\chi^{2}$ values in comparison to the Milky Way SED are below the 95%
upper limit value of $\chi^{2}$ for 13 degrees of freedom; again we normalise
in the $r$-band. The two galaxies with the smallest $\chi^{2}$ (which fall
below the 68% upper limit) have their photometry plotted as round points,
which are been offset in the wavelength direction for clarity. Their
accompanying spectra from the SED atlas are also plotted, without an offset.
We also show as fainter curves the SEDs for the remaining two galaxies with
potentially matching SEDs, which have $\chi^{2}$ values below the 68% limit.
Images of the four best-fitting galaxies are provided in Fig. 10. Portions of
the spectra from the SED atlas that are based on models are depicted using
dotted lines, while those that are directly based on observations are shown
using solid lines. (c) The optical and near-IR portion of the Milky Way SED,
with Brown et al. (2014) templates, as in (b). (d) The mid-IR portion of (b).
The GPR method produces a Milky Way SED that is consistent in shape with both
composite SED templates and individual observed SEDs. Figure 10: Postage stamp
images corresponding to galaxies within the Brown et al. (2014) SED atlas with
$\chi^{2}$ values below the 95% upper limit when compared to the estimated SED
of the Milky Way. At the top of each image we indicate the galaxy NGC number
as well as the $\chi^{2}$ value for the comparison of its SED to the Galactic
one; we present them in order from smallest to largest $\chi^{2}$. At the
bottom left of each tile is a letter marked (A) or (B) which denote the source
of the given image. (A) images are from ESA/Hubble (2021) and (B) images are
from SDSS (2021). The border around each image matches the colour used for its
SED in 9(b). It is evident that galaxies that may have similar SEDs to the
Milky Way exhibit a broad range of visual morphologies, further emphasising
the lack of a unique mapping between morphology and galaxy SED.

While detailed fitting of the Milky Way’s SED using physical models lies
beyond the scope of this paper, we will compare to observed SEDs of individual
galaxies and composite galaxy templates from the literature as a sanity check
on the realism of our results. Using photometry for extragalactic samples to
constrain the SED of the Milky Way, and then comparing the results to observed
galaxy photometry (albeit for different objects) is somewhat circular.
However, given that our analysis has treated every band completely
independently, there were no guarantees that we should get a sensible SED when
combining GPR results across the spectrum.

In 9(a) we compare our predicted Milky Way SED to templates from Benítez et
al. (2004), which are refinements to the templates from Coleman et al. (1980)
and (for starburst galaxies) Kinney et al. (1996). The Coleman et al. (1980)
templates were based upon averaging the observed SEDs of relatively blue
galaxies of a given morphological type. Given the broad range of observed SEDs
for objects with similar morphological classification, these should not be
considered universally applicable for all galaxies of a given type; however,
we use the same labelling as Benítez et al. (2004) for consistency with the
literature.

In this and successive plots, we have normalised all templates to match the
estimated Milky Way SED in the $r$ band. We normalise in an optical band as
those bands have the smallest errors in the predicted SED; which particular
optical band we choose has minimal effect on our comparisons. It is apparent
that the Milky Way SED is generally consistent with the (Benítez et al., 2004)
“Sbc” galaxy template. The galaxies used to construct this template, M51 (NGC
5194) and NGC 2903 are undergoing moderate amounts of star formation, with log
specific star formation rates of -10.1 and -10.4 (versus -10.5 for the Milky
Way) (Muñoz-Mateos et al., 2007). The Milky Way is generally expected to have
an SBbc morphological type (see, e.g., Bland-Hawthorn & Gerhard (2016) for a
recent review on Milky Way structure, as well as Hodge (1983); Kennicutt
(2001); Efremov (2011)), though given how it was constructed, we should not
read too much into the agreement of our SED with the Benítez et al. (2004)
“Sbc” template in particular. It is clear, however, that the GPR method yields
results that resemble composite SED templates from the literature.

Benítez et al. (2004) provides only a sparse set of composite templates that
may not match the SEDs of an individual galaxy. Additionally, those templates
do not span the full wavelength range of the Milky Way SED we have produced.
Thus we also compare to the set of 129 observed galaxy spectral energy
distributions from Brown et al. (2014). This SED atlas encompasses a variety
of bright galaxies in the very local Universe. Unlike the templates shown in
9(a), these correspond to SEDs of individual galaxies (not averages) with
minimal modelling used to interpolate between photometric and spectroscopic
coverage.

The data tables from Brown et al. (2014) provide extinction-corrected
photometry as well as a variety of summary values such as luminosity distance.
Because the magnitudes are presented in the AB system we can use the relation:

$\log{f_{\nu}}=\frac{m_{AB}-8.9}{-2.5},$ (6)

to calculated flux, where $f_{\nu}$ is flux in units of Jansky and $m_{AB}$ is
the observed magnitude in each band. We neglect $k$-corrections as these
galaxies are very nearby ($3.1<D_{L}<249.2\rm{Mpc}$ at the most extreme), so
the corrections are generally negligible. We then can use the flux-luminosity
relation $L_{\nu}=4\pi D_{L}^{2}f_{\nu}$ to calculate the luminosity in each
band. Via propagation of errors the uncertainty in $\log{\nu L_{\nu}}$ for
these templates is equivalent to $\sigma_{\log{\nu
L_{\nu}}}=0.4\sigma_{m_{AB}}$.

We take a few further steps before comparing the Brown et al. (2014) SEDs to
our Galaxy’s. First, we have obtained the axis ratios (as a proxy for
inclination) for 89 out of the 129 galaxies in the Brown et al. (2014) sample
from the Siena Galaxy Atlas, which have been distributed as part of DESI
Legacy imaging surveys Data Release 9. All galaxies with $b/a$ below $0.5$ or
unknown axis ratios were excluded from comparisons. Lower axis ratios should
correspond to highly-inclined galaxies for which reddening will strongly
affect the SED (e.g., Unterborn & Ryden, 2008; Maller et al., 2009), making
them inappropriate comparisons to our face-on SED for the Milky Way. We then
normalise the observed SEDs to match the Galactic SED in the $r$-band, as was
done for the Benítez et al. (2004) templates before.

Finally, we calculate the $\chi^{2}$ difference between our predicted SED for
the Milky Way and each of the galaxy SEDs presented in the Brown et al. (2014)
atlas. We emphasise that no fitting is performed in this comparison other than
matching in the $r$ band. We then calculate $\chi^{2}$ using log quantities,
as that is the space in which we perform our predictions and for which errors
are (by construction) symmetric:

$\chi^{2}=\sum\left(\frac{\log{\nu L_{\nu}^{\rm{atlas}}}-\log{\nu
L_{\nu}^{\rm{MW}}}}{\sigma_{\log{\nu L_{\nu}}}}\right)^{2},$ (7)

where $\sigma_{\log{\nu L_{\nu}}}$ combines in quadrature the total error in
the MW SED for a given band, the uncertainties in the Brown et al. (2014)
photometry, and $\log_{10}$ (1.1), which corresponds to a 10% error in $\nu
L_{\nu}$. This extra error is added to account for systematic uncertainties in
the photometry for a given band relative to others; if this were not included,
optical bands would dominate the $\chi^{2}$ value due to their small nominal
uncertainties. We calculate $\chi^{2}$ using the 14 bandpasses in which we
have measured the Milky Way’s predicted SED, which yields 13 total degrees of
freedom (one is lost due to the $r$-band normalisation performed).

9(b) over-plots the Brown et al. (2014) SEDs for galaxies whose $\chi^{2}$
values fall within the $68\%$ upper limits for a $\chi^{2}$ distribution with
13 degrees of freedom (corresponding to $\chi^{2}=14.8$). We also examine
galaxies that fall below the 95% limit (with $\chi^{2}<22.4$), but exclude
them from plotting for brevity. Objects below the $68\%$ upper limit (four in
total) are clearly consistent with the SED of the Milky Way, while those
between the $68$ and $95$ percent limits (comprising three objects) are in
some tension with the SED of the Milky Way, but could still be a match. While
we calculate $\chi^{2}$ using the broadband photometry for each galaxy, we
also plot the full SED from Brown et al. (2014) for each of these galaxies for
reference.

In 9(b) the two galaxies with the smallest $\chi^{2}$ are labelled and plotted
with the highest opacity in teal and gold. For these two galaxies we also plot
the observed photometry as points offset in the wavelength direction so they
are easier to compare to the Milky Way values. The higher $\chi^{2}$ objects
are plotted with low opacity in pale blue. We also provide more detailed plots
of two separate ranges of the SEDs. The optical through near-IR regime is
depicted in 9(c), while the near- to mid-IR SED is depicted in 9(d). Portions
of the spectra that are based on observations are plotted with solid lines,
while modelled portions are plotted with dotted lines.

We also provide images for the four Brown et al. (2014) atlas galaxies that
fall below the $68\%$ limit in Fig. 10. At the top of each postage stamp image
we list the galaxy NGC number, as well as the $\chi^{2}$ value for comparing
photometry for that galaxy to our Milky Way SED. The tiles are presented in
order from smallest to largest $\chi^{2}$ value. At the bottom left of each
tile is a letter marking (A) or (B) which refers to the source for each given
image. (A) images are from ESA/Hubble (2021) and (B) images are from SDSS
(2021). The borders that surround each image match the colour coding in 9(b).

It is clear from this comparison that our GPR method produces results whose
spectral shape is comparable to observed galaxy SEDs. The small total number
of SEDs within the Brown et al. (2014) atlas that are consistent with the
Milky Way is no surprise; after cutting on inclination we are reduced to only
70 objects, the majority of which are early-type galaxies, leaving a limited
number of examples to cover the full range of star-forming galaxies. We
emphasise that there remains a need for more careful investigation and
modelling of the Milky Way SED we have obtained; this will be the topic of a
follow-up paper.

However, we will briefly comment on the galaxies from the Brown et al. (2014)
atlas whose SEDs are most consistent with the Milky Way. NGC 4138 is an Hubble
type SA(r)0 galaxy which contains an AGN and a star-forming ring. Using
estimates for mass from Jore et al. (1996) ($2.92\times
10^{10}\;\rm{M}_{\odot}$) and Kassin (2004) ($6.23\times
10^{10}\;\rm{M}_{\odot}$) and the star formation rate estimates from Wiegert &
English (2014) ($0.14\;\rm{M}_{\odot}\rm{yr}^{-1}$) and Brown et al. (2017)
($0.2\;\rm{M}_{\odot}\rm{yr}^{-1}$) yield a log sSFR of $\sim-11.3$ to
$-11.5$, significantly smaller than the Milky Way value ($-10.52$). NGC 4138
has a ${}^{0}(g-r)$ colour of $0.73\pm 0.05$ Brown et al. (2014), matching the
Milky Way value.

In contrast, NGC 3351/M95 is a galaxy of Hubble type SB(r)b, versus SBb or SBc
for the Galaxy; it contains a pseudo-bulge (Sandage & Bedke, 1994; Fisher et
al., 2009; Brown et al., 2014), much as the Milky Way is conjectured to
possess (see Bland-Hawthorn & Gerhard, 2016, and references therein).
According to measurements provided by Leroy et al. (2008); George et al.
(2019) NGC 3351 has a log specific star formation rate of $-10.43$ per year
(SFR $=0.940\;\rm{M}_{\odot}\rm{yr}^{-1}$), which is comparable to the log
sSFR of the Milky Way of $-10.52$. According to the measurements complied by
Brown et al. (2014), NGC 3351 has a ${}^{0}(g-r)$ colour of $0.74\pm 0.05$,
very similar to that of the Milky Way’s ($0.73\pm 0.05$, see Table 2). It
matches the Milky Way SED equally well as NGC 4138 at most wavelengths, with
the exception of the longest-wavelength W3 and W4 bands. However, we note that
the photometry for NGC 3351 from Dale et al. (2017) does not show the strong
red slope in the mid-IR seen in the SED atlas spectrum. If the Dale et al.
(2017) measurements were used in the mid-IR, the agreement with the Milky Way
SED would be significantly better.

Two other galaxies in the atlas have SEDs for which the $\chi^{2}$ value when
compared to the Milky Way SED are below the 68% significance level. NGC 5055
is a galaxy of Hubble type SAbc (Brown et al., 2014) with log sSFR of $-10.53$
as tabulated in Kennicutt et al. (2011). NGC 3265 is of Hubble type SA(rs)0
pec (Ann et al., 2015) with log sSFR of $-9.12$ (Kennicutt et al., 2011). The
set of objects whose SEDs are consistent with the Milky Way’s span a diverse
range of visual morphologies.

####  4.2.3 Impact of Physical Parameters on the Estimated SED of the Milky
Way

Figure 11: Isolating the contributions of each physical parameter to the SED.
Each panels shows the effect of varying one parameter while fixing the other
five parameters to the fiducial Milky Way values (see Section 2.3). We vary
most parameters over a range of $\pm 2\sigma$ from the Milky Way’s fiducial
values. In the case of axis ratio, we explore a wider range of values to
convey better the impact of inclination on the observed SED. The upper-left
panel depicts the different SEDs that would be inferred for different galaxy
axis ratios, with the expected results that an inclined Milky Way would look
much fainter than face-on. In the remainder of the panels we depict the log of
the predicted SED divided by the SED evaluated for fiducial Milky Way
parameters in order to make differences more clear; the plotted quantity is
thus $\Delta(\log{(\nu L_{\nu})})=\log{(\frac{\nu L_{\nu}}{{\nu
L_{\nu}^{\rm{MW}}}})}$. The predicted SED for the Milky Way value is plotted
with a dashed black line and open circles at the passbands. On the colour bars
the Milky Way’s value is marked by a horizontal white line. In $b/a$ the GP
captures the effects of dust reddening at higher inclinations. If the Milky
Way were to have a higher SFR, the SED would be brighter in both the UV and
IR. If the Milky Way’s disk were more extended we would observe an increased
UV brightness and decreased mid-IR brightness. We caution the reader that
predictions for disk scale lengths more than $1\sigma$ below the Milky Way
value may not be reliable. Relatively few galaxies of the mass of Milky Way
have sizes smaller than the Milky Way, causing the GPR to be poorly trained
for small $R_{d}$ values (Licquia et al., 2016). Changes in the Milky Way’s
$B/T$ on the SED would be minimal. As we decrease the intensity of the Milky
Way’s bar it seems to have a minimal effect, with a slightly increased UV
brightness. In general these follow our expectations of galaxy evolution which
we discuss further in Section 4.2.3.

Thus far we present results based upon a set of fiducial values (with
uncertainties) for the Milky Way’s stellar mass, star formation rate, axis
ratio, disk scale length, bulge-to-total mass ratio, and bar presence
(presented in Section 2.3). Because of the generality of the Gaussian process
regression fit, we can vary each of these parameters one at a time and test
its impact on the inferred SED.

In our previous analyses we randomly drew values from the distributions of
Milky Way properties and made predictions based on each of those draws, which
were then combined, as described in Section 3. This allowed us to incorporate
the uncertainties in the Milky Way’s physical properties into our results. In
this subsection, however, we will neglect these uncertainties in order to
isolate the effect of changing the central values for each parameter. We have
performed a similar analysis to the one presented below by sampling from Milky
Way uncertainties as a cross-check; the results are very similar so we do not
present them here for brevity.

In this analysis we keep the values for the five Milky Way parameters not
being studied at their fiducial mean. We then choose a discrete set of values
for the sixth parameter at which to evaluate the SED via GPR. In the case of
star formation rate, bulge-to-total mass ratio, and bar vote fraction we
select 8 values that lie evenly spaced between $\pm 2\sigma$ of the fiducial
mean value for the Milky Way (inclusive), in addition to evaluating at the
nominal value. For axis ratio we step through a wide range of possible galaxy
axis ratios instead of focusing around $0.9$ in order to capture the full
effects of inclination on the Galaxy SED. We have excluded mass from this
exercise, as changing mass would radically alter the normalisation of the
predicted SED.

For each of the values that we step through for the given test parameter, the
GPR is evaluated as before. We apply a $12\sigma$ cutoff on the training
sample for all parameters except the one being varied. By focusing on one
parameter at a time we can explore the impact each has on the predicted SED
for the Milky Way and can assess whether the GP is able to capture the
expected correlations between galaxy properties.

The results of this analysis are presented in Fig. 11. In the colour bars for
each panel the lighter shades correspond to smaller values for the parameter
being varied and darker shades correspond to larger values. In all panels the
fiducial value for the Milky Way is marked by a horizontal white line on the
colour bar.

The upper left panel shows the GPR-predicted SED for each axis ratio value
considered, with axes similar to Fig. 9. The SED evaluated with $b/a=0.9$, the
fiducial axis ratio value used for the Milky Way, is marked by open points. In
this panel one can see the effects of inclination reddening first hand, an
effect that has been seen repeatedly in analyses in spiral galaxies (e.g.,
Shao et al., 2007; Unterborn & Ryden, 2008; Maller et al., 2009; Masters et
al., 2010b). The cross-section for dust extinction and scattering generally
increases with decreasing wavelength, causing reddening effects to be the
strongest at the shortest wavelengths. Thus we expect the SEDs of galaxies to
appear redder the more inclined they are (Xiao et al., 2012). The increased
attenuation at higher inclinations for galaxies in the training sample causes
the GPR to predict a redder SED as the inclination increases (lower $b/a$). We
find a decreased brightness in the UV/optical and an increased brightness in
the IR as disks are viewed more edge-on, an effect also observed by SED
modellers (see e.g., Noll et al., 2009).

In the remainder of the panels we plot SEDs divided by the SED evaluated at
the fiducial Milky Way values, as effects are more subtle and would be
difficult to discern in an un-normalised plot. We include a plot of this type
based on varying axis ratio in the lower left panel, though we use a larger
y-range than the other normalised plots due to the large dynamic range
spanned.

An SED which matches the prediction for the fiducial MW values exactly would
fall along the horizontal line at $\log{(\nu L_{\nu})}-\log{(\nu
L_{\nu}^{\rm{MW}})}=0$. The photometric predictions for the Milky Way nominal
parameters hence correspond to the open circles along this line.

Star formation rate has the clearest effect on the SED, as seen in the top
middle panel. The amount of UV flux is a sensitive indicator of star formation
as it is dominated by hot, massive, short-lived stars. These stars contribute
to the flux at optical and near-IR wavelengths, but are subdominant there; but
in the mid-IR the SED responds strongly to star formation due to light from
hot stars that is reprocessed by dust. The GPR predictions reflect all of
these phenomena.

The upper right panel depicts the effect of disk scale length on the predicted
SED. We note that the predictions become unstable at shorter scale lengths due
to the small number of Milky Way-mass galaxies with radius smaller than our
Galaxy in the training samples (cf. Licquia et al. (2016)), so results for
disk scale lengths more than $1\sigma$ below the Milky Way value may not be
robust. Outside of that regime, it is clear that Milky Way-like galaxies with
shorter scale lengths exhibit significantly less flux in the UV than those
with longer scale lengths when $\rm{M}_{*}$, SFR, etc. are all held fixed,
along with smaller effects at optical-IR wavelengths.

A smaller disk, with other properties held fixed, could imply that the gas
within the disk is denser and dust columns are correspondingly greater. This
would in turn cause the SED to look fainter in the UV relative to the IR
compared to if the disk were more extended. We have tested the effect of
varying $b/a$ and $R_{d}$ simultaneously and find that we can compensate for a
change in one with a change in the other almost perfectly. Given the long-
standing scenario that inclination effects on the SEDs of disc galaxies are
driven by dust (Maller et al., 2009, e.g.), it is reasonable to hypothesize
that the effects of varying $R_{d}$ must relate to varying dust impact as
well.

The bottom middle panel shows the effect of varying only the bulge-to-total
ratio. Overall, the impact on the SED is small. We see that if the Milky Way
were to have a more massive bulge, it would appear slightly fainter in the UV.
This could reflect the fact that bulges have older stellar populations than
spiral disks; the effect may be subtle here as to first order the effect is
captured by variation in specific star formation rate. Bulges are also
susceptible to dust reddening (Tuffs et al., 2004; Driver et al., 2007), so a
larger bulge may also suffer from greater reddening effects at short
wavelengths.

The lower right panel shows the variation in SEDs as we change the bar vote
fraction. As this quantity increases we see a decrease in UV brightness and an
increase in mid-IR brightness. In general we speculate that most effects
caused by a bar may have been captured in the variation with other measured
properties (particularly star formation rate). The observed trend with
$p_{\rm{bar}}$ could cohere with a narrative of bar-induced star formation
suppression (see e.g., Tubbs, 1982; Haywood et al., 2016), which would affect
the UV bands the most and lead to a galaxy looking more red once the bar
instability has caused the consumption of cold gas in the disk (Masters et
al., 2010b). However, this would not explain why galaxies with stronger bars
are brighter in the mid-IR while simultaneously being fainter in the UV.
Alternatively, one could speculate that processes associated with bars could
modulate the properties of interstellar dust (e.g., by affecting the ability
of gas to cool and form molecular clouds); having more dust (at fixed
$\rm{M}_{*}$, SFR, and inclination) should cause a reduction of flux in the UV
and an increase in flux in the mid-IR, as observed here.

####  4.2.4 Exploration of Other Sources of Physical Parameter Measurements

As described in Section 2.2.1 and Section 2.2.2, there exist multiple options
for the values used for the stellar mass, star formation rate, and bulge-to-
total ratio for SDSS galaxies. In this subsection we discuss the impact on our
results of using different measurements of galaxy parameters or different
methods of defining our training samples.

##### Stellar Mass and Star Formation Rate

The GSWLC-M2 catalogue (Salim et al., 2016, 2018) includes estimates of
stellar masses and star formation rates computed based on the photometry