# The Progenitor of the Peculiar Galaxy NGC 3077 111Based on data collected at
Subaru Telescope, which is operated by the National Astronomical Observatory
of Japan.

Sakurako Okamoto Subaru Telescope, National Astronomical Observatory of Japan,
650 North A’ohoku Place, Hilo, HI 96720, U.S.A.  National Astronomical
Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo, 181-8588, JAPAN  The
Graduate University for Advanced Studies, Osawa 2-21-1, Mitaka, Tokyo
181-8588, Japan  Nobuo Arimoto  National Astronomical Observatory of Japan,
Osawa 2-21-1, Mitaka, Tokyo, 181-8588, JAPAN  The Graduate University for
Advanced Studies, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan  Astronomy
Program, Department of Physics and Astronomy, Seoul National University, 599
Gwanak-ro, Gwanak-gu, Seoul, 151-742, Korea  Annette M.N. Ferguson Institute
for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill,
Edinburgh, EH9 3HJ U.K.  Mike J. Irwin Institute of Astronomy, University of
Cambridge, Madingley Road, Cambridge CB3 0HA, U.K.  Rokas Žemaitis Institute
for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill,
Edinburgh, EH9 3HJ U.K.

###### Abstract

We present a study of the structural properties and metallicity distribution
of the nearby peculiar galaxy NGC 3077. Using data from our survey of the M81
Group with the Hyper Suprime-Cam on the Subaru Telescope, we construct deep
color-magnitude diagrams that are used to probe the old red giant branch
population of NGC 3077. We map these stars out to and beyond the nominal tidal
radius, which allows us to derive the structural properties and stellar
content of the peripheral regions. We show that NGC 3077 has an extended
stellar halo and pronounced “S-shaped” tidal tails that diverge from the
radial profile of the inner region. The average metallicity of the old
population in NGC 3077 is estimated from individual RGBs to be
$\rm{[M/H]}=-0.98\pm 0.26$, which decreases with the distance from the galaxy
center as $\rm{[M/H]}=-0.17$ dex $\rm{R_{h}}^{-1}$. The metallicity of the
S-shaped structure is similar to that of the regions lying at $r\sim
4\times\rm{R_{h}}(\sim 30$ kpc), indicating that the stellar constituents of
the tidal tails have come from the outer envelope of NGC 3077. These results
suggest that this peculiar galaxy was probably a rather normal dwarf
elliptical galaxy before the tidal interaction with M81 and M82. We also
examine the evidence in our dataset for the six recently-reported ultra-faint
dwarf candidates around NGC 3077. We recover a spatial overdensity of sources
coinciding with only one of these.

galaxies: groups: individual (M81) — galaxies: individual (NGC 3077) —
galaxies: stellar content — galaxies: structure

##  1 Introduction

NGC 3077 is a member galaxy of the M81 Group and is located at a distance of
3.8 Mpc (Sakai & Madore, 2001). This peculiar galaxy has been classified as an
Irr II (Sandage, 1961), Seyfert (Wisniewski & Kleinmann, 1968), I0pec (de
Vaucouleurs et al., 1976), and S0/Sa (Benacchio & Galletta, 1981). Sandage et
al. (1975) commented that NGC 3077 is one of the few galaxies which cannot be
accommodated by the Hubble scheme. The difficulty in classifying its
morphology comes from the influence of the interaction with M81 and M82, which
has distorted the original characteristics of NGC 3077. Price & Gullixson
(1989) concluded that NGC 3077 was a small elliptical galaxy before the
interaction, similar to the M31 dwarf satellite NGC 185.

Together with M81 and M82, NGC 3077 belongs to one of the most famous galaxy
triplets. Conspicuous streams of atomic hydrogen gas (HI) result from tidal
interactions in the past (Yun et al., 1994; Yun, 1999; de Blok et al., 2018)
and both NGC 3077 and M82 exhibit recent starburst activity in their central
areas. It is rather surprising that NGC 3077 has drawn much less attention
from researchers in the past compared to M82. Based on a study of the resolved
stellar populations of the M81 Group with Hyper Suprime-Cam (HSC) on the
Subaru Telescope, Okamoto et al. (2015) showed that the young stars closely
follow the HI distribution and can be found in a stellar stream between M81
and NGC 3077 and in numerous outlying stellar associations. In the central
starburst region of NGC 3077, young star clusters are also formed (Harris et
al., 2004). X-ray observations with CHANDRA revealed that the central area of
NGC 3077 is dominated by a hot intercluster medium surrounded by H$\alpha$
super-shells, and the super-babbles have not yet exploded in contrast to the
case of M82 (Ott et al., 2003). Molecular gas clouds, as traced by CO
(Heithausen & Walter, 2000; Meier et al., 2001; Walter et al., 2002, 2006),
and the dust (Walter et al., 2011; Heithausen, 2012) are concentrated on the
east side of NGC 3077, coincident with the densest HI (Walter et al., 2006).
Dust can also be seen in the main body of NGC 3077, but the dust mass there is
one order of magnitude smaller than that in the eastern area (Walter et al.,
2011). It is not clear whether the dense HI, CO, and the dust were tidally
stripped from NGC 3077 or not (Walter et al., 2011).

The HI tidal streams that link the three galaxies are believed to have formed
during the last gravitational interaction between the three galaxies. The
simulation by Yun (1999) showed that a close encounter between M81 and NGC
3077 likely happened 280 Myrs ago, and another close encounter between M81 and
M82 occurred 220 Myr ago. Yun (1999) failed to reproduce the stream between
NGC 3077 and M82 as the tidally disrupted outer disk of M81; instead, he was
required to introduce a gas disk into the model for NGC 3077. There is no
observational evidence, however, that suggests the HI stream between M82 and
NGC 3077 comes from the tidally stripped gas from NGC 3077. The central star-
forming region suggests that the gas tidally stripped from M81 was accreted
onto NGC 3077 and induced star bursts (Thronson et al., 1991). From nuclear
emission line studies, Heckman (1980) pointed out that the close companions of
M81 (NGC 3077, M82, and NGC 2976) have abnormally high abundances of heavy
elements for their absolute magnitude, suggesting that these galaxies might
have borrowed metal-rich gas from M81 during close encounters. As Price &
Gullixson (1989) and Thronson et al. (1991) pointed out, if NGC 3077 was
originally a dwarf spheroidal or elliptical galaxy, this would be incompatible
with it having a disk of HI mass $\sim 10^{9}\rm{M}_{\odot}$ as Yun (1999)
assumed in his simulation.

From studying the color-magnitude diagrams (CMDs) and the luminosity
functions, Okamoto et al. (2019) demonstrate that the star formation in
outlying young stellar systems around M81 and NGC 3077 has continued until
very recently, roughly 30 Myr ago. They also found that there are no over-
densities of old red giant branch (RGB) stars associated with the star-forming
region “Garland” and other young stellar systems, which led them to conclude
that these systems are genuine young stellar systems whose star formation was
triggered by the latest tidal interactions between M81, NGC 3077, and M82.
Although RGB stars at the location of Garland were identified in the Hubble
Space Telescope (HST) imagery by Sakai & Madore (2001), it is very likely that
these RGB stars are halo stars of NGC 3077. Recently Žemaitis et al. (2022)
discovered the giant stellar stream associated with F8D1, an M81 satellite
galaxy at roughly 2 degrees in projection from M81, which imply the tidal
interactions of the M81 Group are more complicated than previously thought.

The dynamical history of the M81-NGC 3077 pair, and the M81 Group in general,
remains unclear. Numerical simulations of the M81-NGC 3077 interaction by
Thomasson & Donner (1993) showed that spiral arms and HI streams connecting
M81 and NGC 3077 are the results of the interaction between the two galaxies.
Oehm et al. (2017) introduced dynamical friction in their M81 Group
simulation, which suggests that three inner galaxies in the M81 Group are
likely to merge within the next 1-2 Gyr. Yun (1999) also found in their
simulation that three galaxies finally merged.

Okamoto et al. (2015) discovered a huge extended S-shape structure in the
stellar component of NGC 3077, which extends towards the north (the direction
to M82) and the south. Recently Bell et al. (2022) identified five ultra faint
dwarf (UFD) candidates of $\rm{M_{V}}\sim-6$ around NGC 3077, three of which
seem to be embedded in this extended halo. Smercina et al. (2020) also showed
that there is a significant amount of metal-rich stars currently unbound from
NGC 3077, and M82 due to the recent interaction, which may increase halo
metallicity and mass of M81 in the near future.

In this paper, we focus on the spatial structure and metallicity of the old
resolved stellar population of NGC 3077. Using data from our Subaru/HSC survey
of the M81 Group, we map these properties across the main body of the galaxy
as well as in its tidal extensions. Section 2 describes the data and CMDs. The
structure and stellar populations are examined in Section 3 and Section 4,
respectively. In Section 5, we discuss the metallicity, the origin of NGC
3077, the survival timescale of S-shaped structure, and examine the UFD
candidates of Bell et al. (2022). Finally, Section 6 gives a summary and the
conclusions of our findings.

##  2 Data and Color-Magnitude Diagrams

Figure 1: De-reddened CMDs of stellar objects in NGC 3077 and in the field.
(a) Stars within the tidal radius ($\rm{R_{t}}$) of NGC 3077. Blue and magenta
polygons indicate the selection criteria for RGB and AGB stars, respectively.
PARSEC isochrones of age 10 Gyr and $\rm{[M/H]}=-2.2,-1.75,-1.3,-0.8,-0.5$ are
overlaid as blue solid lines. The error bars show the mean photometric error
of sources within $-1.5<(g-i)_{0}<3.0$. (b) Same as CMD (a) but with density
contours overlaid. Contours show number densities binned by 0.025 mag in the
x-axis and 0.04 mag in the y-axis, and the levels are drawn every 50 from 50
to 300 stars per bin.  (c) Stellar objects in the field, selected as regions
which lie outside of $4\times\rm{R}_{25}$ of M81 and M82, and outside of
$\rm{R_{t}}$ of the dwarf galaxies within the survey footprint. The field
counts are scaled so as to represent the same area as that contained within
$\rm{R_{t}}$ of NGC 3077. The contour levels are the same with CMD (c). (d)
The field subtracted Hess diagram of NGC 3077. The binned difference between
CMD (b) and CMD (c) is shown. PARSEC isochrones of age 10 Gyr are overlaid as
black solid lines.

The data we analyze here are drawn from a photometric survey of the M81 Group
that we have conducted with Hyper Suprime-Cam on the Subaru Telescope
(Miyazaki et al., 2012). Our survey covers 12 square degrees around M81 with
seven pointings and reaches roughly two magnitudes below the RGB tip at the
distance of M81. The observing strategy, details of data processing and
calibration have been previously described in Okamoto et al. (2015, 2019) and
Žemaitis et al. (2022). Throughout this study, we use the version of the de-
reddened point-source catalog that is full-described in Okamoto et al. (2019).
To briefly recap, the HSC pipeline is used for the bias correction, flat-
fielding, mosaicking, stacking, and calibration (Bosch et al., 2018; Ivezić et
al., 2019; Axelrod et al., 2010). For the dataset under study here, we then
used the IRAF implementation of DAOPHOT (Stetson, 1987) to obtain point spread
function fitting photometry of sources. Artificial star tests were performed
on some parts of the reduced images using the ADDSTAR task in DAOPHOT to add
stars every 0.5 mag from mag=20.0 to 28.0. The result is shown in Figure 2 of
Okamoto et al. (2019), which indicates the point source catalog is at least
50$\%$ complete to 26 mag in both passbands, except for the crowded region of
NGC 3077 center. We fitted a modified sigmoid function to the point-source
detection rate as functions of the magnitude and color, as given by the
artificial star tests. As seen by the hollow center in the left panel of
Figure 2, DAOPHOT failed to detect sources in the highly crowded central few
arc minutes of NGC 3077. To account for this, we only consider the
r$>5\arcmin$ area for the analysis in this paper. We used the DAOPHOT
parameters $\chi^{2}$ and sharpness to separate point sources from the
extended objects and noise-like detections. Specifically, we select objects as
stars whose $\chi^{2}$ and sharpness values lie within $3\times\sigma$ of the
mean values of the artificial stars of the same magnitude.

Figure 1 shows the CMDs of stellar objects in NGC 3077 and in a reference
field selected to be devoid of member stars of the M81 Group galaxies. Black
dots in panels (a) and (b) are de-reddened point sources that lie within the
tidal radius of NGC 3077 derived in Section 3. In panel (a), theoretical
PARSEC v1.2S isochrones (Bressan et al., 2012) of age 10 Gyr and metallicity
$\rm{[M/H]}=-2.2,-1.75,-1.3,-0.8,-0.5$, and of age 16, 32, 100 Myr with
$\rm{[M/H]}=-0.75$, are overlaid as blue solid lines. In doing this, we have
assumed the distance modulus of $(m-M)_{0}=27.93$ to NGC 3077 (Sakai & Madore,
2001). The error bars show the mean photometric error of stars in
$-1.5<(g-i)_{0}<3.0$. Panel (c) shows point sources that lie outside of
$4\times\rm{R}_{25}$ of M81 and M82 (de Vaucouleurs et al., 1991), and outside
of $\rm{R_{t}}$ of the dwarf galaxies that lie within the survey footprint
studied in Okamoto et al. (2019); these objects are expected to be
contaminants, although they may contain a trace amount of stars in the
extended halos of M81 and NGC 3077. These objects are sampled so as to reflect
the same area as that contained within $\rm{R_{t}}$ of NGC 3077. Panel (d) is
the Hess diagram showing the relative density of stars between panel (b) and
panel (c). The bin size is 0.1mag.

A well-populated RGB is visualized by the extended contour at $i_{0}>24$ and
at $(g-i)_{0}>1$ in Figure 1(b), which is not seen in Figure 1(c). The field-
subtracted Hess diagram in Figure 1(d) also shows this overdensity. Another
striking feature in Figure 1(b) is the young blue population vertically
distributed at $(g-i)_{0}<-0.4$. We select the bright RGB population using the
solid blue polygon shown in Figure 1(a). The boundaries of this have been
chosen to optimize the selection of genuine NGC 3077 stars while limiting the
number of foreground and background contaminants. We note that there is a
substantial asymptotic giant branch (AGB) population within the dotted magenta
area above the RGB selection box. As discussed in Okamoto et al. (2015), these
bright AGB stars are centrally concentrated in NGC 3077, which implies these
stars are intrinsic populations. Further investigation on the AGB population
of NGC 3077 will be done in the future with high-resolution NIR imagery from
Euclid.

##  3 The Outer Structure of NGC 3077

Figure 2: Left: The spatial distribution of RGB stars in and around NGC 3077.
The solid ellipses represent the King profile parameters $\rm{R_{c}}$ and
$\rm{R_{t}}$ as estimated using star counts over the whole area, while the
dashed lines show those parameters derived considering only stars in the outer
region. Right: A smoothed density map of RGB stars shown on a logarithmic
scale and reflecting densities starting at $0.4\sigma$ of the background
level. The kernel density is estimated with the bandwidth of $1\arcmin.0$.
Other M81 Group galaxies projected nearby are labeled with white letters. The
white solid and dashed lines represent $\rm{R_{t}}$ as the same as the left
panel. The gray solid and dashed lines show the break radii
$\rm{R_{B}}=14\arcmin$ and $\rm{R_{B,E}}=10.5\arcmin$ (see Section 5.3). White
and black star symbols indicate the locations of UFD candidates reported by
Bell et al. (2022) (see Section 5.4).

Using the RGB stars selected from the CMD, we proceed to analyze the structure
of NGC 3077 in the same manner as done in Okamoto et al. (2019). Briefly, this
involves deriving the centroid, position angle, and ellipticity using the
density-weighted first and second moments of the RGB star spatial distribution
and then using these parameters to construct the radial star count profile.

As seen in Figure 2, NGC 3077 has an extended and diffuse stellar halo and
pronounced “S-shaped” tidal tails. The halo exhibits a rather boxy shape,
which is not well-matched by the position angle and the ellipticity derived
from the fit to the full star count distribution (solid ellipses). We also
derive the average ellipticity and position angle of the star counts excluding
the inner r$=8.4\arcmin$ area of NGC 3077. The parameters derived in this way
provide a much better description of the outer structure of NGC 3077 as shown
by the dashed lines in Figure 2. The resulting structural properties are
listed in Table 1. Hereafter we use $\rm{R}_{{\rm{(out)}}}$ to denote radii
estimated using the ellipse parameters derived from only the stars in the
outer region.

Figure 3: Completeness-corrected radial number density profile of RGB stars in
and around NGC 3077. Upper panel: Azimuthally-averaged radial star count
profiles derived within elliptical annuli. The best-fitting King, exponential,
Plummer, and Sersic profiles are overlaid as magenta, orange, green, and blue
lines, respectively, although the exponential profile (orange line) is almost
overlapped with the Sersic profile (blue). The foreground/background
contamination level is shown as a dashed horizontal line. The profile derived
using the ellipse parameters derived from stars in the outer region only is
shown as a gray color. Lower panel: Radial profiles for the west and east
sides of NGC 3077 using structural parameters derived from stars in the outer
region are shown as blue and magenta colors, respectively.  Table 1: The
structural properties of NGC 3077

Galaxy | $\alpha$ | $\delta$ | P.A.(a) |  $\epsilon$(b) |  $\rm{R_{c}}$(c) |  $\rm{R_{t}}$(d) |  $\rm{R_{exp}}$(e) | b(f) | n(g) |  $\rm{R_{h}}$(h)
---|---|---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) | (deg) |  | (arcmin) | (arcmin) | (arcmin) | (arcmin) |  | (arcmin)
|  |  |  |  | (kpc) | (kpc) | (kpc) | (kpc) |  | (kpc)
NGC 3077 | $10:03:14.2$ | $+68:44:19.9$ | $46.3$ | $0.22$ | $4.6\pm 0.2$ | $28\pm 1$ | $3.6\pm 0.1$ | $6.7\pm 0.1$ | $1.0\pm 0.1$ | $6.0\pm 0.1$
|  |  |  |  | $5.1\pm 0.2$ | $32\pm 2$ | $4.0\pm 0.1$ | $7.5\pm 0.1$ |  | $6.7\pm 0.1$
NGC 3077 (outer) |  |  | $72.3$ | $0.52$ | $5.1\pm 0.2$ | $37.0\pm 1.2$ | $4.31\pm 0.07$ | $8.2\pm 0.2$ | $1.18\pm 0.06$ | $7.2\pm 0.1$
|  |  |  |  | $5.7\pm 0.2$ | $41\pm 1.3$ | $4.8\pm 0.1$ | $9.1\pm 0.2$ |  | $8.1\pm 0.1$

Note. — (a) Position angle from north to east. (b) Ellipticity
$\epsilon=1-\rm{b/a}$ where a is the scale length of the galaxy along the
major axis and b is that along the minor axis. (c) Core radius of King
profile. (d) Tidal radius of King profile. (e) Scalelength of the exponential
profile. (f) Plummer radius. (g) Sersic index. (h) Half light radius.

The completeness-corrected radial number density distribution of RGB stars is
fits with the standard King, exponential, Plummer, and Sersic models via
least-squares minimization. The best-fit profiles are overlaid as magenta,
orange, green, and blue lines, respectively, in the upper panel of Figure 3.
The foreground/background contamination level is shown as a dotted horizontal
line. Our imagery does not resolve stars in the very inner part of NGC 3077,
as evidenced by the hole in the left panel of Figure 2, hence we exclude the
inner $5\arcmin$ region from the fit. The rise in the profile beyond
$\rm{R_{t}}$ is due to contamination from the stellar components of M81 and
the other dwarf galaxies in the surrounding area. The star count profile is
generally well-described by all the models. That said, there is evidence for a
clear kink in the radial profile at around $\rm{r}\approx 14\arcmin$. The
radial profile derived using the structural parameters of the outer region is
shown in the gray color in the upper panel of Figure 3, where the kink is even
more obvious. This radius is indicated with the gray solid ellipse in the
right panel of Figure 2. Such features have been suggested to reflect the
break radius $\rm{R_{B}}$ at which a transient excess in the stellar density
profile appears due to tidal stripping of a dwarf galaxy (Johnston et al.,
2002; Peñarrubia et al., 2009). The stars in these features have gained energy
at the last pericentric passage and are moving outwards as they settle into
new, less bound orbits in the dwarf. We will use this $\rm{R_{B}}$ to estimate
the elapsed time since NGC 3077 last pericentric passage around M81 in Section
5.3.

Figure 4: Metallicity distribution of RGB stars within
$4\times\rm{R_{h}}_{\rm{(out)}}$ of NGC 3077 derived using the $(g-i)_{0}$
colour and assuming a 10 Gyr population. The black open histogram shows the
completeness-corrected MDF using all RGB stars. The shaded histogram is the
field-subtracted MDF. The solid line shows the fitted single Gaussian
distribution to the field-subtracted MDF. The mean metallicity error is shown
on the upper left side.

To further investigate the bump feature, we divide the area around NGC 3077
into eastern and western halves. The lower panel of Figure 3 shows the radial
stellar density of the eastern and the western sides. The radial star counts
of each are constructed by calculating the completeness-corrected average
number of RGB stars with $\Delta\rm{RA}>0$ ($\Delta\rm{RA}<0$) in elliptical
annuli, defined using the structural parameters derived from stars in the
outer region. The profile of the western side is highly contaminated by the
presence of M81, while the distribution in the eastern side shows a pronounced
bump at $\rm{R_{B,E}}\approx 10.5\arcmin$. The gray dashed ellipse in Figure 2
shows that this also corresponds to the location where the structure of the
dwarf starts to become increasingly elongated.

##  4 The Metallicity Distribution of the Old Population

Figure 5: The smoothed density map of RGB stars shown on a logarithmic scale
on which is overlaid dashed ellipses indicating $0.75\times\rm{R_{c}}$ and
$1.5\times\rm{R_{c}}$ radii and solid ellipses indicating
$2\times\rm{R_{h}}_{\rm{(out)}},3\times\rm{R_{h}}_{\rm{(out)}},4\times\rm{R_{h}}_{\rm{(out)}}$.
CMDs of stars located within the nine labeled regions are shown in Figure 6.
Figure 6: De-reddened CMDs of stellar objects within nine regions defined in
Figure 5. The CMD of the reference region is scaled to reflect the same areal
coverage as that of region 6 and region 8. The Parsec v1.2S isochrones of 10
Gyr with $\rm{[M/H]}=-2.2,-1.75,-1.3,-0.8,-0.5$ are overlaid as magenta solid
lines.

In this section, we examine the metallicity distribution of RGB stars in the
inner and outer regions of NGC 3077, which we have established exhibit
different structural properties in the previous section. The black open
histogram in Figure 4 is the metallicity distribution function (MDF) of all
RGB stars lying within $4\times\rm{R_{h}}_{\rm{(out)}}$ of the NGC 3077
center, corresponding to the outermost ellipse shown in Figure 5. We derive
the metallicity of individual RGB stars using the $(g-i)_{0}$ color and
assuming a 10 Gyr population in the same manner as those in Okamoto et al.
(2019).

The photometric completeness of each RGB star is corrected according to its
color and magnitude. Foreground and background contamination is estimated
using point sources drawn from the reference region shown in Figure 5, which
is scaled so as to represent the same area. The field contamination corrected
MDF is shown as the shaded histogram in Figure 4. The metallicity errors of
individual RGB stars are estimated by performing a Monte Carlo simulation with
N$=2000$ for each star, sampling from a Gaussian distribution with a width
equal to the photometric error. The estimated uncertainties increase from
$0.09$ to $0.48$ with the decreasing metallicity from $\rm{[M/H]}=-0.3$ to
$-2.5$ due to the narrow range of RGB color between metal-poor isochrones. The
mean metallicity corrected for this effect is shown as the error-weighted mean
value $\rm{[M/H]}_{\rm{W}}=-0.98\pm 0.24$. The error bar in Figure 4 indicates
the uncertainty at the mean metallicity. Our analysis uses all stars within
the RGB box but experiments showed that the results did not change if we
instead only used the brightest half of the RGB.

Figure 7: MDFs of RGB stars in regions 1 to region 8, as defined in Figure 5.
The shaded histograms are the field-subtracted MDFs. The solid lines show
single Gaussian fits to the distributions. The blue dashed lines in regions
2-8 are the scaled versions of the Gaussian fit to region 1 for reference
purposes. The region ID, mean metallicity, weighted mean metallicity, and
dispersion are listed in each panel.

To compare the MDFs in different areas in the galaxy, we divide NGC 3077 into
eight regions; (1) $0.75~{}\rm{R_{c}}<r<1.5~{}\rm{R_{c}}$, (2)
$1.5~{}\rm{R_{c}}<r<2~{}\rm{R_{h(out)}}$, (3)
$2~{}\rm{R_{h(out)}}<r<3~{}\rm{R_{h(out)}}$, (4)
$3~{}\rm{R_{h(out)}}<r<4~{}\rm{R_{h(out)}}$, (5) and (6) north part of the
S-shaped structure, (7) and (8) south part of the S-shaped as shown in Figure
5. The CMDs of stars within each region are shown in Figure 6. In regions 1 to
3, the dominant populations are seen to be RGB stars with $1.0<(g-i)_{0}<2.0$
and $i_{0}>24.0$. This RGB population diminishes in prominence moving from the
inner to outer regions, and unresolved background galaxies with $(g-i)_{0}\sim
0.4$ and $i_{0}\sim 25.8$ are the dominant feature in the outer “S-shaped”
structure CMDs (regions 5-8) and the reference CMD. Vertical features
corresponding to young populations lie at $(g-i)_{0}<-0.4$ and are visible in
regions 1 to 4. The only evidence of young stars in the “S-shaped” feature
comes from region 7, where the young stellar stream linking NGC 3077 and M81
was discovered (Okamoto et al., 2015). Figure 6 also shows that the luminous
intermediate-age AGB population, which is located above the RGB tip, is also
less prominent, moving from the center to the outer regions.

Figure 7 shows the MDF of the RGB populations in each region, derived in the
same manner as for Figure 4. The foreground and background counts in the
reference field are scaled so as to represent the same area as each region.
The blue dashed lines in the panels of regions 2-8 are the scaled versions of
the Gaussian fit to region 1 and are shown for reference. Moving from the
inner to the outermost parts of the main body (regions 1 to 4), the mean
metallicity and metallicity dispersion decrease with radius. Though it is
within the error bars, a mild radial metallicity gradient can be derived as
$\Delta\rm{[M/H]}=-0.02$ dex/kpc in physical scale or $-0.14$ dex
$\rm{R_{h(out)}}^{-1}$ as a function of the half-light radius using the
structural parameters derived for the outer region. In the “S-shaped” feature
(region 5-7), the contaminants dominate over the RGB population, but the
reference field-subtracted histograms still show a comparable or slightly
metal-poor population with that of region 4. This is consistent with this
structure being made up of stars that have been stripped from the outermost
part of NGC 3077 due to the tidal interaction with M81 and M82.

##  5 Discussion

###  5.1 Metallicity

In the previous section, we derive a mild metallicity gradient in the old RGB
population of NGC 3077 in the sense of decreasing metallicity with increasing
radius. Such gradients are common in Local Group dwarf galaxies (Taibi et al.,
2022). NGC 3077’s gradient, $-0.14$ dex $\rm{R_{h(out)}}^{-1}$ is comparable
with that of the dwarf elliptical NGC 185, $-0.2$ dex $\rm{R_{h}}^{-1}$ (Taibi
et al., 2022). The MDF is also similar to that of NGC 185 (Crnojević et al.,
2014). Although NGC 3077 has much larger radius of $\rm{R_{h(out)}}=8.1$ kpc
than NGC 185 (0.5 kpc) due to the recent gravitational interaction with M81
and M82, these results support the suggestion by Price & Gullixson (1989) that
NGC 3077 was once a normal dwarf elliptical, akin to NGC 185 before its
interaction.

The presence of metallicity gradients in dwarf galaxies has been examined in
high-resolution cosmological simulations of galaxy formation. Recently Mercado
et al. (2021) use the FIRE-2 simulation suite to investigate the presence of
metallicity gradients in dwarf galaxies of stellar mass $10^{5.5}$ to
$10^{8.8}~{}\rm{M}_{\odot}$. They show that inverse metallicity gradients in
the stellar component are common and are set by the competition between the
“puffing” process of old, metal-poor stars by feedback-driven potential
fluctuations in the early phases of formation and the late time accretion of
metal-rich gas that provides fuel for extended star formation. In this sense,
NGC 3077’s metallicity gradient may eventually be washed out due to new-born
metal-rich stars formed in the peripheral regions by gas displaced due to the
recent tidal interactions.

Figure 8 shows the stellar mass and metallicity of NGC 3077 compared to that
of other galaxies in the local universe. The M81 dwarf galaxies studied in
Okamoto et al. (2019) and Lianou et al. (2010), nearby dwarf galaxies from
McConnachie (2012) and early-type dwarf galaxies in the Virgo cluster from
Toloba et al. (2014) are plotted. We convert galaxies $\rm{M_{V}}$ in the
literature to the stellar mass, assuming the M/L of dE and dIrr in Woo et al.
(2008). The mass-metallicity relations of the Local Group galaxies (Kirby et
al., 2013), FIRE simulation suite (Ma et al., 2016), and SDSS galaxies
(Gallazzi et al., 2005) are also plotted in Figure 8. NGC 3077’s stellar mass
of $7.2\pm 1.2\times 10^{8}\rm{M}_{\odot}$ is calculated using RGB and AGB
stars within $4\times\rm{R_{h(out)}}$ (the outermost ellipse in Figure 5) to
avoid contamination of M81 and nearby satellites. The $i$-band and $g$-band
luminosities of individual RGB and AGB stars located between
$0.75\times\rm{R_{c(out)}}$ to $4\times\rm{R_{h(out)}}$ are summed up with the
photometric completeness correction. The total luminosity is corrected to
account for the flux of inner $0.75\times\rm{R_{c(out)}}$ area using the
extrapolated Sersic radial profile and for the flux of stars fainter than the
RGB selection box using the 10 Gyr and $\rm{[M/H]}=-1.3$ luminosity function,
then it is corrected the background/foreground contamination using the
reference field. The total magnitude $\rm{M_{V}}$ is calculated from
transforming $\rm{M_{g}}$ and $\rm{M_{i}}$ (Jordi et al., 2006), and then
converted to the stellar mass assuming $\rm{M_{V}}/L=1.6$ (Woo et al., 2008).
Considering just the old stellar components of NGC 3077, as we have done here,
it can be seen that NGC 3077 lies very close to the mass-metallicity relation
defined by dwarf galaxies and is very close to dwarf ellipticals in the Virgo
cluster.

Figure 8: The stellar mass-metallicity relation of nearby galaxies. NGC 3077
is shown with a black diamond. The other galaxies shown include the M81 dwarf
galaxies studied in Okamoto et al. (2019) and Lianou et al. (2010), nearby
dwarf galaxies from McConnachie (2012) and early-type dwarf galaxies in the
Virgo cluster from Toloba et al. (2014). Where necessary, we have assumed the
M/L of Woo et al. (2008) to derive stellar masses.

###  5.2 Origin of NGC 3077

Various pieces of evidence support the idea that NGC 3077 was an NGC 185-like
dwarf elliptical galaxy before the interaction with M81 and M82. This includes
the fact that its radial stellar density profile is well-fit by a King and
other models which are traditionally used to describe dE/dSph. It is also
located on the mass-metallicity relation defined by nearby galaxies and
inhabits a region of the diagram that is populated by early-type dwarfs.
Finally, it has a stellar MDF that is very similar to those of low-mass
elliptical galaxies (Crnojević et al., 2014). Taken together, it seems
unlikely that NGC 3077 was originally a late-type dwarf with an extended gas
disk as Yun (1999) assumed in his numerical simulation. However, without this
assumption, Yun (1999) could not reproduce the global distribution of HI in
the M81 Group. This implies that there may yet be missing ingredients in the
modeling of the interaction history of the M81 Group. The recent discovery of
a giant tidal stream from the dwarf galaxy F8D1 (Žemaitis et al., 2022)
suggests that it also had a recent passage close to M81 and may also have been
stripped of its gas at this time. Further modeling of how interactions can
explain the extended HI gas in the M81 Group is clearly warranted.

Figure 9: The oxygen abundance of M81 and NGC 3077 on the HI map taken from de
Blok et al. (2018). The metallicity of 21 regions around M81 from Patterson et
al. (2012) and at the center of NGC 3077 from Calzetti et al. (2004) are
color-coded according to the metallicity.

Another argument against the gas around NGC 3077 being of intrinsic origin
comes from consideration of the gas-phase metallicity. The gas oxygen
abundance of NGC 3077 is much higher than expected from its luminosity
(Heckman, 1980). Figure 9 shows the oxygen abundances that have been measured
for M81 and NGC 3077 superposed on the HI map of de Blok et al. (2018).
Calzetti et al. (2004) derived $12+\log(\rm{O/H})=8.9\pm 0.1$ for NGC 3077 by
applying the Kewley & Dopita (2002) calibration to measurements of the strong
emission lines. By using the same Kewley & Dopita (2002) calibration,
Patterson et al. (2012) derived oxygen abundances in the gas disk of M81 and
found a radial gradient going from $12+\log(\rm{O/H})\sim 9.0$ at the galaxy
center to $\sim 8.6$ at 30 kpc. This strongly suggests that the gas currently
in the inner regions of NGC 3077 was once part of the outer disk of M81 and
accreted onto NGC 3077 as a result of the tidal interaction. In this scenario,
the HI stream that links M81 and NGC 3077 represents the extension of a tidal
arm that resulted from NGC 3077 falling in on a prograde orbit within the
plane of M81’s disk (Thomasson & Donner, 1993). NGC 3077 sits at the edge of
this extended gaseous arm. Gas might have flowed inwards to trigger the
central starburst (Lelli et al., 2014). The concentration of HI, CO, and the
dust in the eastern side of NGC 3077, collectively referred to as the
“Garland”, would not originate from NGC 3077 itself. Instead, it is more
likely to have formed during the last close encounter, much like the other
young systems such as Arp’s Loop and Holmberg IX.

###  5.3 Survival Timescale of the S-shaped Structure

Okamoto et al. (2015) discovered that NGC 3077 has a very extended stellar
halo, from which a faint “S-shaped” structure emanates to the north and south.
Here we try to estimate when this feature was formed. Using N-body
simulations, Peñarrubia et al. (2009) investigated how the structure of a
dwarf satellite galaxy evolves during a tidal encounter with a massive host.
Interestingly, the extended shape of the NGC 3077’s halo looks almost
identical to their findings (top-middle panel in their Figure 1) for a dwarf
spheroidal system on a highly eccentric orbit. This striking resemblance
demonstrates that, at least qualitatively, the extended halo of NGC 3077 is
consistent with being formed from the tidal interaction between NGC 3077, and
M81. Furthermore, it also suggests the possibility of timing the last
pericentric passage by comparing the observed stellar density profile with the
simulated profiles at various snapshots in time. Peñarrubia et al. (2009)
showed that there is a transient excess in the radial stellar density profile
when a dwarf spheroidal galaxy undergoes tidal stripping. They proposed that
the radius at which this transient excess occurs, the “break radius”, can be
used as a clock to measure the elapsed time since the last pericentric
passage. In the case of NGC 3077, we have found in Section 3 that the break
occurs at $\log(\rm{R_{B}}/\rm{R_{c}})\sim 0.5$, where $\rm{R_{c}}$ is a core
radius and $\rm{R_{B}}$ is the break radius. According to the simulation of
Peñarrubia et al. (2009), there exists a tight relation between the last
pericentric passage, $\rm{t-t_{p}}$, and the “break radius” $\rm{R_{B}}$, that
is given by $\rm{R_{B}}=\rm{C}{\sigma_{0}}(\rm{t-t_{p}})$, where
$\rm{C}=0.55\pm 0.03$ and $\sigma_{0}$ is the stellar velocity dispersion at
the galaxy center. With our measured $\log(\rm{R_{B}}/\rm{R_{c}})\sim 0.5$, we
then obtain $\rm{t-t_{p}}=6.3~{}\rm{t_{cr}}$, where
$\rm{t_{cr}}=\rm{R_{c}}/\sigma_{0}$ is the crossing time of NGC 3077.
Therefore, the extended halo and its characteristic S-shape were formed
roughly 6 crossing times ago. This value is also very close to the one of
$5.7\times\rm{t_{cr}}$, which describes the top-middle panel of Figure 1 in
Peñarrubia et al. (2009).

To go further and estimate the actual elapsed time since the last pericentric
passage, we need to know the central velocity dispersion of NGC 3077.
Unfortunately, we could not find a measurement of this quantity in the
literature. When $\rm{R_{c}}$ and $\sigma_{0}$ are given in kpc and km/sec,
respectively, the crossing time is given as follows:

$\displaystyle\rm{t_{cr}}\sim\rm{R_{c}}[\rm{kpc}]/(\sigma_{0}[\rm{km/s}]\times
10^{9}[\rm{yr}])$ (1)

From the density profile, we obtain $\rm{R_{c}}=5.1\pm 0.2$ kpc (see Table 1).
If we assume the velocity dispersion of NGC 3077 is similar to that of NGC
185, $\sigma_{0}=24\pm 1$km/s (Geha et al., 2010), we obtain
$\rm{t_{cr}}=1.7\times 10^{8}$yr. Thus, the last close encounter with M81 that
triggered the extension of NGC 3077 is estimated as $6.3\rm{t_{cr}}\sim 1$Gyr
ago. While this is an approximative estimation, it might be reasonable since
Oehm et al. (2017) calculated close encounters between M81, M82, and NGC 3077
happened in recent few Gyrs. Mancillas et al. (2019) also showed tidal low
surface brightness features survive by 0.7-4 Gyrs after major and
intermediate-mass merger events.

###  5.4 Examining the Bell et al. (2022) Ultra-Faint Dwarf Candidates

Figure 10: The spatial distributions of RGB stars around the six UFD
candidates reported by Bell et al. (2022). The solid lines represent the half-
light radii $\rm{R_{h}}$ as reported by those authors. The panel number
corresponds to the digit in Figure 2.

Recently, Bell et al. (2022) conducted a search for ultra-faint dwarf (UFD)
candidates in the M81 Group, finding six of them. Remarkably, all of these
candidates lie in close proximity to NGC 3077 which is the smallest system of
the interacting triplet. Three of these sources, J0959+6837, J1000+6841, and
J1004+6835, are located well within NGC 3077’s halo structure (see black star
symbols in the right panel of Figure 2), while J1001+6907 and J1002+6903 lie
around NGC3077’s tidal radius (see white stars in Figure 2). For these five
systems, Bell et al. (2022) base their search on an independent dataset that
they have obtained with HSC.

We revisit the evidence for these sources using the dataset we have presented
here. Figure 10 shows the spatial distribution of RGB stars around the UFD
candidates, including the sixth one, J0954+6821, defined as a definite dwarf
satellite located to the south of M81 by Bell et al. (2022). This system lies
within our HSC M81 Group coverage but outwith the immediate area of NGC 3077
analyzed in this paper. The panel number (2) to (6) corresponds to the digit
in Figure 2. While there are a few RGB stars located within the half-light
radii in all cases except panel (5), we could not see any obvious over-
densities in panels (2) to (6) associated with these systems. For J0954+6821
(panel (1)), we confirm that there are some point sources of $i_{0}>24$ that
coincide with a region of unresolved diffuse light in our pixel frame. From
the ground-based imagery only, it is rather difficult to determine if these
are RGB stars at the M81 distance of 3.6 Mpc or if they are more distant red
supergiant (RSG) stars. If these sources are intrinsically bright RSGs
($\rm{M_{V}}\sim-6$), the distance of J0954+6821 is estimated as D $>10$ Mpc,
well behind the M81 Group. For further confirmation, deep high-resolution
images from HST or JWST are required.

##  6 Summary

We have presented the structure and metallicity distribution function of NGC
3077. Based on the M81 galaxy group survey with Hyper Suprime-Cam on the
Subaru Telescope, we have constructed the deep CMD that reached the old RGB
population of NGC 3077 and covered well beyond the tidal radius, which allows
us to derive the structural properties and stellar content of the peripheral
regions. NGC 3077 has the extended stellar halo and the “S-shaped” tidal
feature, which is not well-matched by the position angle and the ellipticity
derived from the fit to the full star count distribution. The average
metallicity of individual RGBs decreases with the distance from the NGC 3077
center. The metallicity at the S-shape area shows a comparable or slightly
metal-poor population with that of $r\sim 4\times\rm{R_{h(out)}}$ area,
implying this structure being made up of stars that have been stripped from
the outermost part of NGC 3077 due to the tidal interaction with M81 and M82.
The old stellar component of NGC 3077 is on the mass-metallicity relation of
nearby low-mass ellipticals and shows the mild metallicity gradient
$\rm{[M/H]}=-0.14$ dex $\rm{R_{h(out)}}^{-1}$, which is comparable with that
of NGC 185. These results, as well as the oxygen abundance of NGC 3077, and
M81, suggest that this peculiar galaxy had been a normal dwarf elliptical
galaxy before the interaction. The close encounters with M81 and M82 might
supply the gas to NGC 3077 center and induce the current central starburst.
The ongoing strong tidal effect on NGC 3077 is also stripping the stellar
constituent from the outer envelope, which is now seen as the “S-shaped”
structure.

We are grateful to the entire staff at Subaru Telescope and the HSC team. This
research is based on data collected at the Subaru Telescope, which is operated
by the National Astronomical Observatory of Japan. We are honored and grateful
for the opportunity of observing the Universe from Maunakea, which has the
cultural, historical, and natural significance in Hawaii. This paper makes use
of software developed for Vera C. Rubin Observatory. We thank the Rubin
Observatory for making their code available as free software at
http://pipelines.lsst.io/. The Pan-STARRS1 Surveys (PS1) have been made
possible through contributions of the Institute for Astronomy, the University
of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its
participating institutes, the Max Planck Institute for Astronomy, Heidelberg
and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns
Hopkins University, Durham University, the University of Edinburgh, Queen’s
University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las
Cumbres Observatory Global Telescope Network Incorporated, the National
Central University of Taiwan, the Space Telescope Science Institute, the
National Aeronautics and Space Administration under Grant No. NNX08AR22G
issued through the Planetary Science Division of the NASA Science Mission
Directorate, the National Science Foundation under Grant No. AST-1238877, the
University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos
National Laboratory. S.O. acknowledges support in part from JSPS Grant-in-Aid
for Scientific Research (18H05875, 20K04031, 20H05855). N.A. thanks the Brain
Pool program for financial support, which is funded by the Ministry of Science
and ICT through the National Research Foundation of Korea (2018H1D3A2000902).
For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## References

  * Axelrod et al. (2010) Axelrod, T., Kantor, J., Lupton, R. H., & Pierfederici, F. 2010, in Proc. SPIE, Vol. 7740, Software and Cyberinfrastructure for Astronomy, 774015, doi: 10.1117/12.857297
  * Bell et al. (2022) Bell, E. F., Smercina, A., Price, P. A., et al. 2022, ApJ, 937, L3, doi: 10.3847/2041-8213/ac8e5e
  * Benacchio & Galletta (1981) Benacchio, L., & Galletta, G. 1981, ApJ, 243, L65, doi: 10.1086/183444
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5, doi: 10.1093/pasj/psx080
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Calzetti et al. (2004) Calzetti, D., Harris, J., Gallagher, III, J. S., et al. 2004, AJ, 127, 1405, doi: 10.1086/382095
  * Crnojević et al. (2014) Crnojević, D., Ferguson, A. M. N., Irwin, M. J., et al. 2014, MNRAS, 445, 3862, doi: 10.1093/mnras/stu2003
  * de Blok et al. (2018) de Blok, W. J. G., Walter, F., Ferguson, A. M. N., et al. 2018, ApJ, 865, 26, doi: 10.3847/1538-4357/aad557
  * de Vaucouleurs et al. (1976) de Vaucouleurs, G., de Vaucouleurs, A., & Corwin, Jr., H. G. 1976, Second reference catalogue of bright galaxies. Containing information on 4,364 galaxies with references to papers published between 1964 and 1975.
  * de Vaucouleurs et al. (1991) de Vaucouleurs, G., de Vaucouleurs, A., Corwin, Jr., H. G., et al. 1991, Third Reference Catalogue of Bright Galaxies. Volume I: Explanations and references. Volume II: Data for galaxies between 0h and 12h. Volume III: Data for galaxies between 12h and 24h.
  * Gallazzi et al. (2005) Gallazzi, A., Charlot, S., Brinchmann, J., White, S. D. M., & Tremonti, C. A. 2005, MNRAS, 362, 41, doi: 10.1111/j.1365-2966.2005.09321.x
  * Geha et al. (2010) Geha, M., van der Marel, R. P., Guhathakurta, P., et al. 2010, ApJ, 711, 361, doi: 10.1088/0004-637X/711/1/361
  * Harris et al. (2004) Harris, J., Calzetti, D., Gallagher, III, J. S., Smith, D. A., & Conselice, C. J. 2004, ApJ, 603, 503, doi: 10.1086/381669
  * Heckman (1980) Heckman, T. M. 1980, A&A, 87, 142
  * Heithausen (2012) Heithausen, A. 2012, A&A, 543, A21, doi: 10.1051/0004-6361/201117861
  * Heithausen & Walter (2000) Heithausen, A., & Walter, F. 2000, A&A, 361, 500
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Johnston et al. (2002) Johnston, K. V., Choi, P. I., & Guhathakurta, P. 2002, AJ, 124, 127, doi: 10.1086/341040
  * Jordi et al. (2006) Jordi, K., Grebel, E. K., & Ammon, K. 2006, A&A, 460, 339, doi: 10.1051/0004-6361:20066082
  * Kewley & Dopita (2002) Kewley, L. J., & Dopita, M. A. 2002, ApJS, 142, 35, doi: 10.1086/341326
  * Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102, doi: 10.1088/0004-637X/779/2/102
  * Lelli et al. (2014) Lelli, F., Verheijen, M., & Fraternali, F. 2014, MNRAS, 445, 1694, doi: 10.1093/mnras/stu1804
  * Lianou et al. (2010) Lianou, S., Grebel, E. K., & Koch, A. 2010, A&A, 521, A43, doi: 10.1051/0004-6361/200913364
  * Ma et al. (2016) Ma, X., Hopkins, P. F., Faucher-Giguère, C.-A., et al. 2016, MNRAS, 456, 2140, doi: 10.1093/mnras/stv2659
  * Mancillas et al. (2019) Mancillas, B., Duc, P.-A., Combes, F., et al. 2019, A&A, 632, A122, doi: 10.1051/0004-6361/201936320
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * Meier et al. (2001) Meier, D. S., Turner, J. L., & Beck, S. C. 2001, AJ, 122, 1770, doi: 10.1086/323136
  * Mercado et al. (2021) Mercado, F. J., Bullock, J. S., Boylan-Kolchin, M., et al. 2021, MNRAS, 501, 5121, doi: 10.1093/mnras/staa3958
  * Miyazaki et al. (2012) Miyazaki, S., Komiyama, Y., Nakaya, H., et al. 2012, in Proc. SPIE, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, 84460Z, doi: 10.1117/12.926844
  * Oehm et al. (2017) Oehm, W., Thies, I., & Kroupa, P. 2017, MNRAS, 467, 273, doi: 10.1093/mnras/stw3381
  * Okamoto et al. (2015) Okamoto, S., Arimoto, N., Ferguson, A. M. N., et al. 2015, ApJ, 809, L1, doi: 10.1088/2041-8205/809/1/L1
  * Okamoto et al. (2019) Okamoto, S., Arimoto, N., Ferguson, A. M. N., et al. 2019, ApJ, 884, 128, doi: 10.3847/1538-4357/ab44a7
  * Ott et al. (2003) Ott, J., Martin, C. L., & Walter, F. 2003, ApJ, 594, 776, doi: 10.1086/377040
  * Patterson et al. (2012) Patterson, M. T., Walterbos, R. A. M., Kennicutt, R. C., Chiappini, C., & Thilker, D. A. 2012, MNRAS, 422, 401, doi: 10.1111/j.1365-2966.2012.20616.x
  * Peñarrubia et al. (2009) Peñarrubia, J., Navarro, J. F., McConnachie, A. W., & Martin, N. F. 2009, ApJ, 698, 222, doi: 10.1088/0004-637X/698/1/222
  * Price & Gullixson (1989) Price, J. S., & Gullixson, C. A. 1989, ApJ, 337, 658, doi: 10.1086/167137
  * Sakai & Madore (2001) Sakai, S., & Madore, B. F. 2001, ApJ, 555, 280, doi: 10.1086/321459
  * Sandage (1961) Sandage, A. 1961, The Hubble Atlas of Galaxies
  * Sandage et al. (1975) Sandage, A., Sandage, M., & Kristian, J. 1975, Galaxies and the Universe
  * Smercina et al. (2020) Smercina, A., Bell, E. F., Price, P. A., et al. 2020, ApJ, 905, 60, doi: 10.3847/1538-4357/abc485
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977
  * Taibi et al. (2022) Taibi, S., Battaglia, G., Leaman, R., et al. 2022, A&A, 665, A92, doi: 10.1051/0004-6361/202243508
  * Thomasson & Donner (1993) Thomasson, M., & Donner, K. J. 1993, A&A, 272, 153
  * Thronson et al. (1991) Thronson, Jr., H. A., Wilton, C., & Ksir, A. 1991, MNRAS, 252, 543, doi: 10.1093/mnras/252.4.543
  * Toloba et al. (2014) Toloba, E., Guhathakurta, P., Peletier, R. F., et al. 2014, ApJS, 215, 17, doi: 10.1088/0067-0049/215/2/17
  * Žemaitis et al. (2022) Žemaitis, R., Ferguson, A. M. N., Okamoto, S., et al. 2022, arXiv e-prints, arXiv:2209.09713.  https://arxiv.org/abs/2209.09713
  * Walter et al. (2006) Walter, F., Martin, C. L., & Ott, J. 2006, AJ, 132, 2289, doi: 10.1086/508273
  * Walter et al. (2002) Walter, F., Weiss, A., Martin, C., & Scoville, N. 2002, AJ, 123, 225, doi: 10.1086/324633
  * Walter et al. (2011) Walter, F., Sandstrom, K., Aniano, G., et al. 2011, ApJ, 726, L11, doi: 10.1088/2041-8205/726/1/L11
  * Wisniewski & Kleinmann (1968) Wisniewski, W. Z., & Kleinmann, D. E. 1968, AJ, 73, 866, doi: 10.1086/110721
  * Woo et al. (2008) Woo, J., Courteau, S., & Dekel, A. 2008, MNRAS, 390, 1453, doi: 10.1111/j.1365-2966.2008.13770.x
  * Yun (1999) Yun, M. S. 1999, in IAU Symposium, Vol. 186, Galaxy Interactions at Low and High Redshift, ed. J. E. Barnes & D. B. Sanders, 81
  * Yun et al. (1994) Yun, M. S., Ho, P. T. P., & Lo, K. Y. 1994, Nature, 372, 530, doi: 10.1038/372530a0


1533 LABEL:LastPageJun. 14, 2018Jul. 24, 2019

A preliminary version of this paper appeared in ICALP 2017.

# Subspace-Invariant AC${}^{\text{0}}$ Formulas

Benjamin Rossman  University of Toronto<EMAIL_ADDRESS>
###### Abstract.

We consider the action of a linear subspace $U$ of $\\{0,1\\}^{n}$ on the set
of $\mathrm{AC}^{0}$ formulas with inputs labeled by literals in the set
$\\{X_{1},\overline{X}_{1},\dots,X_{n},\overline{X}_{n}\\}$, where an element
$u\in U$ acts on formulas by transposing the $i$th pair of literals for all
$i\in[n]$ such that $u_{i}=1$. A formula is $U$-invariant if it is fixed by
this action. For example, there is a well-known recursive construction of
depth $d+1$ formulas of size $O(n{\cdot}2^{dn^{1/d}})$ computing the
$n$-variable parity function; these formulas are easily seen to be
$P$-invariant where $P$ is the subspace of even-weight elements of
$\\{0,1\\}^{n}$. In this paper we establish a nearly matching
$2^{d(n^{1/d}-1)}$ lower bound on the $P$-invariant depth $d+1$ formula size
of parity. Quantitatively this improves the best known
$\Omega(2^{\frac{1}{84}d(n^{1/d}-1)})$ lower bound for unrestricted depth
$d+1$ formulas [rossman2015average], while avoiding the use of the switching
lemma. More generally, for any linear subspaces $U\subset V$, we show that if
a Boolean function is $U$-invariant and non-constant over $V$, then its
$U$-invariant depth $d+1$ formula size is at least $2^{d(m^{1/d}-1)}$ where
$m$ is the minimum Hamming weight of a vector in $U^{\bot}\setminus V^{\bot}$.

###### Key words and phrases:

circuit complexity, formula size, AC0, parity, symmetric computation

Supported by NSERC and Sloan Research Fellowship.

##  1\. Introduction

There are two natural group actions on the set of literals
$\\{X_{1},\overline{X}_{1},\dots,X_{n},\overline{X}_{n}\\}$: the symmetric
group $S_{n}$ acts by permuting indices, while $Z_{2}^{n}$ acts by toggling
negations. These group actions extend to the set of $n$-variable Boolean
functions, as well as the set of $n$-variable Boolean circuits. Here we
consider bounded-depth circuits with unbounded fan-in AND and OR gates and
inputs labeled by literals, also known as $\mathrm{AC}^{0}$ circuits. If $G$
is subgroup of $S_{n}$ or $Z_{2}^{n}$ (or more generally of the group
$Z_{2}^{n}\rtimes S_{n}$ that they generate), we say that a function or
circuit is $G$-invariant if it is fixed under the action of $G$ on the set of
$n$-variable functions or circuits. Note that every $G$-invariant circuit
computes a $G$-invariant function, and conversely every $G$-invariant function
is computable by a $G$-invariant circuit.

We define the $G$-invariant circuit size of a $G$-invariant function $f$ as
the minimum number of gates in a $G$-invariant circuit that computes $f$. This
may be compared to the unrestricted circuit size of $f$, noting that $f$ can
be computed (possibly more efficiently) by circuits that are not
$G$-invariant. Several questions arise. What gap, if any, exists between the
$G$-invariant vs. unrestricted circuit size of $G$-invariant functions? Are
lower bounds on $G$-invariant circuit size easier to obtain, and do they
suggest new strategies for proving lower bounds for unrestricted circuits? Is
there a nice characterization of functions computable by polynomial-size
$G$-invariant circuits? The same questions may be asked with respect to
$G$-invariant versions of other complexity measures, such as formula
(leaf)size, as well as bounded-depth versions of both circuit and formula
size, noting that the action of $G$ on circuits preserves both depth and fan-
out.

The answer to these questions appears to be very different for subgroups of
$S_{n}$ and subgroups of $Z_{2}^{n}$. This is illustrated by considering the
$n$-variable parity function, which maps each element of $\\{0,1\\}^{n}$ to
its Hamming weight modulo $2$. This function is both $S_{n}$-invariant (it is
a so-called “symmetric function”) and $P$-invariant where $P\subset Z_{2}^{n}$
is the index-2 subgroup of even-weight elements in $Z_{2}^{n}$. The smallest
known circuits and formulas for $\textsc{parity}_{n}$ have size $O(n)$ and
leafsize $O(n^{2})$, respectively. These circuits and formulas turn out to be
$P$-invariant, as do the smallest known bounded-depth circuits and formulas
(which we describe in §2.3). In contrast, the $S_{n}$-invariant circuit size
of $\textsc{parity}_{n}$ is known to be exponential [anderson2016symmetric].

###  1.1. Invariance under subgroups of $S_{n}$

$G$-invariant circuit complexity for subgroups $G$ of the symmetric group
$S_{n}$ has been previously studied from the standpoint of Descriptive
Complexity, an area of research concerned with the characterization of
complexity classes in terms of definability in different logics
[immerman2012descriptive]. Here one considers Boolean functions that encode
isomorphism-invariant properties of relational structures. Properties of
$m$-vertex simple graphs, for instance, are identified with $G$-invariant
functions $\\{0,1\\}^{n}\to\\{0,1\\}$ of $n=\binom{m}{2}$ variables, each
corresponding to a potential edge, where $G$ is the group $S_{m}$ acting on
the set of potential edges. More generally, if $\sigma$ is a finite relational
signature $\sigma$, one considers the action of $S_{m}$ on
$n=\sum_{R\in\sigma}m^{\mathrm{arity}(R)}$ variables encoding the possible
$\sigma$-structures with universe $[m]$.

Denenberg et al [denenberg1986definability] showed that $S_{m}$-invariant
circuits of polynomial size and constant depth (subject to a certain
uniformity condition) capture precisely the first-order definable properties
of finite $\sigma$-structures. Otto [otto1996logic] introduced a certain limit
object of finite circuits (imposing uniformity in a different way) and showed
a correspondence between the logic $L^{\omega}_{\infty\omega}$ (infinitary
logic with a bounded number of variables) and $S_{m}$-invariant circuits of
polynomial size and arbitrary depth. Otto also gave characterizations of
fixed-point logic and partial-fixed-point logic in terms of $S_{m}$-invariant
Boolean networks. Recently, Anderson and Dawar [anderson2016symmetric] showed
a correspondence between fixed-point logic and polynomial-size
$S_{m}$-invariant circuits, as well as between fixed-point logic with counting
and polynomial-size $S_{m}$-invariant circuits in the basis that includes
majority gates.

Choiceless Polynomial Time [BGS99, BGS02, dawar2015symmetric,
rossman2010choiceless] provides a different example of a $G$-invariant model
of computation, where $G\subseteq S_{n}$ is the automorphism group of the
input structure. Invariance under subgroups of $S_{n}$ has been explored in
other settings as well, see for instance [ajtai1994symmetric,
riis2000generating, rudich1988optimal].

###  1.2. Invariance under subgroups of $Z_{2}^{n}$

This paper initiates a study of invariant complexity with respect to subgroups
of $Z_{2}^{n}$. Since our methods are linear algebraic, we shall henceforth
identify $Z_{2}^{n}$ with the $\mathbb{F}_{2}$-vector space $\\{0,1\\}^{n}$
under coordinate-wise addition modulo $2$, denoted $\oplus$. We identify
subgroups of $Z_{2}^{n}$ with linear subspaces $U$ of $\\{0,1\\}^{n}$. A
function $f:\\{0,1\\}^{n}\to\\{0,1\\}$ is $U$-invariant if $f(x)=f(x\oplus u)$
for all $x\in\\{0,1\\}^{n}$ and $u\in U$. Note that $U$-invariant functions
are in one-to-one correspondence with functions from the quotient space
$\\{0,1\\}^{n}/U$ to $\\{0,1\\}$.

Our focus is on bounded-depth circuits and formulas. Returning to the example
of the $P$-invariant function $\textsc{parity}_{n}$ (where $P$ is the even-
weight subgroup of $\\{0,1\\}^{n}$), there is a well-known recursive
construction of depth $d+1$ circuits for $\textsc{parity}_{n}$, which we
describe in §2.3. Roughly speaking, one combines a depth 2 circuit for
$\textsc{parity}_{n^{1/d}}$ with depth $d$ circuits for
$\textsc{parity}_{n^{(d-1)/d}}$ on disjoint blocks of variables. This produces
a depth $d+1$ circuit of size $O(n^{1/d}{\cdot}2^{n^{1/d}})$, which converts
to a depth $d+1$ formula of leafsize $O(n{\cdot}2^{dn^{1/d}})$. Up to constant
factors, these circuit and formulas are the smallest known computing
$\textsc{parity}_{n}$ and they are easily seen to be $P$-invariant, as we
explain in §2.3.

The main result of this paper gives a nearly matching lower bound of
$2^{d(n^{1/d}-1)}$ on the $P$-invariant depth $d+1$ formula size of
$\textsc{parity}{}_{n}$. This implies a $2^{n^{1/d}-1}$ lower bound on the
$P$-invariant depth $d+1$ circuit size, via the basic fact that every
($U$-invariant) depth $d+1$ circuit of size $s$ is equivalent to a
($U$-invariant) depth $d+1$ formula of size at most $s^{d}$. Quantitatively,
the lower bounds are stronger than the best known
$\Omega(2^{\frac{1}{10}n^{1/d}})$ and $\Omega(2^{\frac{1}{84}d(n^{1/d}-1)})$
lower bounds for unrestricted depth $d+1$ circuits [Hastad86] and formulas
[rossman2015average], respectively. Of course, $P$-invariance is a severe
restriction for circuits and formulas, so it is no surprise that the lower
bounds we obtain is stronger and significantly easier to prove. The linear-
algebraic technique in this paper is entirely different from the “switching
lemma” approach of [Hastad86, rossman2015average].

The general form of our lower bound is the following:

######  Theorem 1.

Let $U\subset V$ be linear subspaces of $\\{0,1\\}^{n}$, and suppose $F$ is a
$U$-invariant depth $d+1$ formula which is non-constant over $V$. Then $F$ has
size at least $2^{d(m^{1/d}-1)}$ where $m=\min\\{|x|:x\in U^{\bot}\setminus
V^{\bot}\\}$, that is, the minimum Hamming weight of a vector $x$ which is
orthogonal to $U$ and non-orthogonal to $V$.

Here size refers to the number of depth 1 subformulas, as opposed to leafsize.
Note that the bound in Theorem 1 does not depend on the dimension $n$ of the
ambient space. Also note that aforementioned $2^{d(n^{1/d}-1)}$ lower bound
for $\textsc{parity}{}_{n}$ follows from the case $U=P$ and $V=\\{0,1\\}^{n}$.
(Here $m=n$ is witnessed by the all-1 vector, which is an element of
$P^{\bot}\setminus(\\{0,1\\}^{n})^{\bot}$.)

We remark that, since $\lim_{d\to\infty}d(m^{1/d}-1)=\ln(m)$, Theorem 1
implies an $m^{\ln(2)}$ lower bound on the size of unbounded-depth formulas
which are $U$-invariant and non-constant over $V$. Theorem 1 also implies a
$2^{m^{1/d}-1}$ lower bound for depth $d+1$ circuits; however, we get no
nontrivial lower bound for unbounded-depth circuits, since
$\lim_{d\to\infty}m^{1/d}-1=0$.

##  2\. Preliminaries

Let $n$ range over positive integers. $[n]$ is the set $\\{1,\dots,n\\}$.
$\ln(n)$ is the natural logarithm and $\log(n)$ is the base-2 logarithm.

The Hamming weight of a vector $x\in\\{0,1\\}^{n}$, denoted $|x|$, is the
cardinality of the set $\\{i\in[n]:x_{i}=1\\}$. For vectors
$x,y\in\\{0,1\\}^{n}$, let $x\oplus y$ denote the coordinate-wise sum modulo
$2$ and let $\langle x,y\rangle$ denote the inner product modulo $2$.

Let $\mathcal{L}$ denote the lattice of linear subspaces of $\\{0,1\\}^{n}$.
For $U,V\in\mathcal{L}$, let $U+V$ denote the subspace spanned by $U$ and $V$.
Let $V^{\bot}$ denote the orthogonal complement
$V^{\bot}=\\{x\in\\{0,1\\}^{n}:\langle x,v\rangle=0\text{ for all }v\in V\\}$.
We make use of the following facts about orthogonal complements over finite
fields:

$\displaystyle\dim(V)+\dim(V^{\bot})=n,\qquad U\subseteq V\Longleftrightarrow
V^{\bot}\subseteq U^{\bot},$ $\displaystyle
V=(V^{\bot})^{\bot},\qquad(U+V)^{\bot}=U^{\bot}\cap V^{\bot},\qquad(U\cap
V)^{\bot}=U^{\bot}+V^{\bot}.$

###  2.1. AC${}^{\text{0}}$ formulas

We write $\mathcal{F}$ for the set of $n$-variable $\mathrm{AC}^{0}$ formulas
(with unbounded fan-in AND and OR gates and leaves labeled by literals).
Formally, let $\mathcal{F}=\bigcup_{d\in\mathbb{N}}\mathcal{F}_{d}$ where
$\mathcal{F}_{d}$ is the set of depth $d$ formulas, defined inductively:

  * •

$\mathcal{F}_{0}$ is the set
$\\{X_{1},\overline{X}_{1},\dots,X_{n},\overline{X}_{n}\\}\cup\\{0,1\\}$,

  * •

$\mathcal{F}_{d+1}$ is the set of ordered pairs

$\displaystyle\\{(\mathsf{gate},\mathcal{G}):\mathsf{gate}\in\\{\mathrm{AND},\mathrm{OR}\\}\text{
and }\mathcal{G}\text{ is a nonempty subset of }\mathcal{F}_{d}\\}.$

Every formula $F\in\mathcal{F}$ computes a Boolean function
$\\{0,1\\}^{n}\to\\{0,1\\}$ in the usual way. For $x\in\\{0,1\\}^{n}$, we
write $F(x)$ for the value of $F$ on $x$. For a nonempty set
$S\subseteq\\{0,1\\}^{n}$ and $b\in\\{0,1\\}$, notation $F(S)\equiv b$ denotes
that $F(x)=b$ for all $x\in S$. We say that $F$ is non-constant on $S$ if
$F(S)\not\equiv 0$ and $F(S)\not\equiv 1$.

The depth of $F$ is the unique $d\in\mathbb{N}$ such that
$F\in\mathcal{F}_{d}$. Leafsize is the number of depth 0 subformulas, and size
is the number of depth 1 subformulas. Inductively,

$\displaystyle\mathrm{leafsize}(F)$ $\displaystyle=\begin{cases}\hskip
0.8pt1&\text{if }F\in\mathcal{F}_{0},\\\ \hskip
0.8pt\sum_{G\in\mathcal{G}}\mathrm{leafsize}(G)&\text{if
}F=(\mathsf{gate},\mathcal{G})\in\mathcal{F}\setminus\mathcal{F}_{0},\end{cases}$
$\displaystyle\mathrm{size}(F)$ $\displaystyle=\begin{cases}0&\text{if
}F\in\mathcal{F}_{0},\\\ 1&\text{if }F\in\mathcal{F}_{1},\\\
\sum_{G\in\mathcal{G}}\mathrm{size}(G)\phantom{\text{leaf}}&\text{if
}F=(\mathsf{gate},\mathcal{G})\in\mathcal{F}\setminus(\mathcal{F}_{0}\cup\mathcal{F}_{1}).\end{cases}$

Clearly $\mathrm{size}(F)\leq\mathrm{leafsize}(F)$. Note that size is within a
factor $2$ of the number of gates in $F$, which is how one usually measures
the size of circuits. Our lower bound naturally applies to size, while the
upper bound that we present in §2.3 is naturally presented in terms of
leafsize.

###  2.2. The action of $\\{0,1\\}^{n}$

We now formally define the action of $\\{0,1\\}^{n}$ (as the group
$Z_{2}^{n}$) on the set $\mathcal{F}$. For $u\in\\{0,1\\}^{n}$ and
$F\in\mathcal{F}$, let $F^{u}$ be the formula obtained from $F$ by exchanging
literals $X_{i}$ and $\overline{X}_{i}$ for every $i\in[n]$ with $u_{i}=1$.
Formally, this action is defined inductively by

$\displaystyle F^{u}=\begin{cases}F&\text{if $F\in\\{0,1\\}$,}\\\ X_{i}\text{
(resp.\ $\overline{X}_{i}$)}&\text{if $F=X_{i}$ (resp.\ $\overline{X}_{i}$)
and $u_{i}=0$,}\\\ \overline{X}_{i}\text{ (resp.\ $X_{i}$)}&\text{if $F=X_{i}$
(resp.\ $\overline{X}_{i}$) and $u_{i}=1$,}\\\
(\mathsf{gate},\\{G^{u}:G\in\mathcal{G}\\})&\text{if
$F=(\mathsf{gate},\mathcal{G})$.}\end{cases}$

Note that $F^{u}$ has the same depth and size as $F$ and computes the function
$F^{u}(x)=F(x\oplus u)$ for all $x\in\\{0,1\\}^{n}$.

Let $U$ be a linear subspace of $\\{0,1\\}^{n}$ (i.e., subgroup of
$Z_{2}^{n}$). We say that an $\mathrm{AC}^{0}$ formula $F$ is:

  * •

$U$-invariant if $F^{u}=F$ (i.e., these are syntactically identical formulas)
for every $u\in U$,

  * •

semantically $U$-invariant if $F$ computes a $U$-invariant function (i.e.,
$F(x)=F(x\oplus u)$ for every $u\in U$ and $x\in\\{0,1\\}^{n}$).

Note that every $U$-invariant formula is semantically $U$-invariant, but not
conversely. For example, the formula
$(\mathrm{AND},\\{0,X_{1},\dots,X_{n}\\})$ computes the identically zero
function and is therefore semantically $U$-invariant (for any $U$); however,
this formula is not $U$-invariant (for any nontrivial $U$).

###  2.3. Upper bound

We review the smallest known construction of bounded-depth formulas for
$\textsc{parity}_{n}$ (see [Hastad86]) and observe that these formulas are
$P$-invariant where $P$ is the even-weight subspace of $\\{0,1\\}^{n}$.

######  Proposition 2.

For all $d,n\geq 1$, $\textsc{parity}_{n}$ is computable by $P$-invariant
depth $d+1$ formulas with either AND or OR as output gate and leafsize at most
$n{\cdot}2^{dn^{1/d}}$. If $n^{1/d}$ is an integer, this bound improves to
$n{\cdot}2^{d(n^{1/d}-1)}$.

######  Proof 2.1.

Define $\beta(d,n)$ by the following recurrence:

$\displaystyle\beta(1,n)$ $\displaystyle=\begin{cases}1&\text{if }n=1,\\\
\infty&\text{if
}n>1,\end{cases}\qquad\beta(d+1,n)=\min_{\begin{subarray}{c}k,n_{1},\dots,n_{k}\geq
1\,:\\\
n_{1}+\dots+n_{k}=n\end{subarray}}2^{k-1}\sum_{i=1}^{k}\beta(d,n_{i}).$

We will construct depth $d+1$ formulas of leafsize $\beta(d+1,n)$. If
$n^{1/d}$ is an integer, we get the bound $\beta(d+1,n)\leq
n{\cdot}2^{d(n^{1/d}-1)}$ by setting $k=n^{1/d}$ and
$n_{1}=\dots=n_{k}=n^{(d-1)/d}$. For arbitrary $d,n\geq 1$, we get the bound
$\beta(d+1,n)\leq n{\cdot}2^{dn^{1/d}}$ by setting $k=\lceil n^{1/d}\rceil$
and $n_{1},\dots,n_{k}\in\\{\lfloor n/k\rfloor,\lceil n/k\rceil\\}$. In
particular, note that $\beta(2,n)=n2^{n-1}$.

In the base case $d=1$, we have the brute-force DNF (OR-of-ANDs) and CNF (AND-
of-ORs) formulas of leafsize $n2^{n-1}$ for parityn. These formulas are
clearly $P$-invariant. Otherwise (if $d\geq 2$), fix the optimal choice of
parameters $k,n_{1},\dots,n_{k}$ for $\beta(d+1,n)$. Partition $[n]$ into sets
$J_{1}\sqcup\dots\sqcup J_{k}$ of size $|J_{i}|=n_{i}$. Let
$\textsc{parity}{}_{J_{i}}$ be the parity function over variables
$\\{X_{j}:j\in J_{i}\\}$ and let $P_{J_{i}}$ be the subspace
$\\{u\in\\{0,1\\}^{n}:\bigoplus_{j\in J_{i}}u_{j}=0\\}$.

By the induction hypothesis, for each $i\in[k]$ there exists a
$P_{J_{i}}$-invariant formula $G_{i}$ computing $\textsc{parity}{}_{J_{i}}$
with depth $d$ and leafsize at most $\beta(d,n_{i})$ and output gate AND. Let
$H_{i}$ be the formula obtained from $G_{i}$ by transposing literals $X_{j}$
and $\overline{X}_{j}$ for any choice of $j\in J_{i}$; note that $H_{i}$
computes $1-\textsc{parity}{}_{J_{i}}$. Let $F$ be the brute-force DNF formula
for parityk over variables $Y_{1},\dots,Y_{k}$. We first form a depth $d+2$
formula $F^{\prime}$ by replacing each literal $Y_{i}$ (resp.
$\overline{Y}_{i}$) in $F$ with the formula $G_{i}$ (resp. $H_{i})$. The two
layers of gates in $F^{\prime}$ below the output consist entirely of AND
gates; these two layers may be combined into a single layer, producing a
formula $F^{\prime\prime}$ of depth $d+1$. Since each variable $Y_{i}$ occurs
in $2^{k-1}$ literals of $F$, the leafsize of $F^{\prime\prime}$ is
$2^{k-1}\sum_{i=1}^{k}\beta(d,n_{i})$ as required.

Finally, to see that $F^{\prime\prime}$ is $P$-invariant, consider an even-
weight vector $u\in\\{0,1\\}^{n}$. Note that $u$ projects to an even-weight
vector in $\\{0,1\\}^{k}$ whose $i$th coordinate is $\bigoplus_{j\in
J_{i}}u_{i}$. Then $u$ acts on $F^{\prime\prime}$ by transposing subformulas
$G_{i}$ and $H_{i}$ for all $i\in[k]$ such that $\bigoplus_{j\in
J_{i}}u_{i}=1$; therefore, $P$-invariance of $F^{\prime\prime}$ follows from
$P_{\\{Y_{1},\dots,Y_{k}\\}}$-invariance of $F$. If we take $F$ to be a CNF
instead of a DNF, the same construction produces $F^{\prime\prime}$ with OR
instead of AND as its output gate. ∎

######  Remark 3.

$\textsc{parity}{}_{n}$ is known to be computable by $P$-invariant formulas of
depth $\lceil\log n\rceil+1$ and leafsize $O(n^{2})$ [tarui2010smallest,
Yab54]. The $n{\cdot}2^{dn^{1/d}}$ upper bound of Proposition 2 is therefore
slack, as this equals $n^{3}$ when $d=\log n$, whereas
$n{\cdot}2^{d(n^{1/d}-1)}=n^{2}$. We suspect that the upper bound of
Proposition 2 can be improved that $O(n{\cdot}2^{d(n^{1/d}-1)})$ for all
$d\leq\log n$, perhaps by a more careful analysis of the recurrence for
$\beta(d+1,n)$. Let us add that $\Omega(n^{2})$ is a well-known lower bound
for any depth, without the assumption of $P$-invariance
[khrapchenko1971complexity].

##  3\. Linear-algebraic lemmas

Recall that $\mathcal{L}$ denotes the lattice of linear subspaces of
$\\{0,1\\}^{n}$. Let $U,V,S,T$ range over elements of $\mathcal{L}$. If $U$ is
a subspace of $V$, recall that a projection from $V$ to $U$ is a linear map
$\rho:V\to U$ such that $\rho(u)=u$ for every $u\in U$. We begin by showing
that if $U$ is a codimension-$k$ subspace of $V$ (i.e., $\dim(V)-\dim(U)=k$),
then there there exists a projection $\rho:V\to U$ with “Hamming-weight
stretch” $k+1$.

######  Lemma 4.

If $U$ is a codimension-$k$ subspace of $V$, then there exists a projection
$\rho$ from $V$ to $U$ such that $|\rho(v)|\leq(k+1)|v|$ for all $v\in V$.


Johnny Tiu, Shankar Ramharack, Patrick Hosein 11institutetext: The University
of the West Indies, St. Augustine, Trinidad and Tobago, 11email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
# On the Optimal Placement of Electric Vehicle Charging Stations

Johnny Tiu  Shankar Ramharack  Patrick Hosein

###### Abstract

Increasing the adoption of Electric Vehicles (EV) is an integral part of many
strategies to address climate change and air pollution. However, Electric
Vehicle adoption rates are inhibited by several factors which reduce the
confidence of potential buyers of electric vehicles, namely range anxiety and
limited charging infrastructure. The latter concern can be addressed by
carefully planning the placement of Electric Vehicle charging stations to
sustainably meet long-term demand. We outline a methodology to optimally
allocate Electric Vehicle Charging stations using a novel approach to demand
estimation. Using this approach the first two authors placed second in the
Academic category of an International Shell competition which was created to
obtain solutions to this problem. They were also provided with a grant to
further study this approach and to develop a deployment strategy.

###### keywords:

EV demand forecasting · EV charging infrastructure · Optimization

##  1 Introduction

Electric Vehicles (EVs) are distinguished by the use of electric motors as a
primary driver for propulsion and inherently produce lower greenhouse gas
(GHG) emissions compared to traditional Internal Combustion Engine (ICE)
Vehicles. Thus, they play a major role in many countries’ goals in achieving
net-zero GHG emissions [1], with the United States of America setting a goal
of 50% of all new vehicles sold in 2030 to be zero-emission [2] and the United
Kingdom designating that no diesel/gasoline-powered vehicles to be sold by
2030 [3].

In addition, the purchase cost differential which is one of the major barriers
to EV adoption is range anxiety [4], which is a fear that their EV will run
out of charge before reaching their destination and is underpinned by the
concern of available charging infrastructure. These each present their own
issues but are ultimately intertwined and can be addressed through proper
forecasting of electrical demand and optimal placement of charging stations.
Various modelling techniques have previously been used for short, medium and
long-term electric demand forecasting such as structural time series and state
space models, probabilistic models [5] , and Artificial Neural Networks (ANNs)
[6] and Convolutional Neural Networks (CNNs) [7]. We instead introduce a
mixed-integer programming model for placement optimization together with a
unique approach to demand estimation.

##  2 Problem Description

The Shell.ai Hackathon 2022 challenged teams to optimally place Electric
Vehicle charging stations so that the configuration remains robust to
demographic changes. They provided the required data, the constraints to be
met and the objective to be achieved. Participants then had to develop a
solution that would maximize the proposed objective function.

###  2.1 Dataset Description

For an undisclosed location, a synthetic dimensionless dataset based on real
world data was distributed to contestants. This consisted of:

Demand History:

Yearly demand data between 2010 and 2018 for each demand point in a $64\times
64$ grid.

Existing Infrastructure:

Distribution of Slow charging stations (SCS), Fast Charging Stations (FCS) and
the maximum capacity for 100 supply points as of 2018. These are plotted in
Figure 1 showing the spatial distribution of supply points and charging demand
intensity.

Figure 1: Charging Demand Heatmap and Supply point locations

###  2.2 Objective

There are two objectives. Firstly, the EV charging station demand must be
forecast over time and, using these predicted estimates, additional charging
stations must be placed to satisfy this demand while minimizing the cost of
the infrastructure. The total cost comprised of three components:

Customer Dissatisfaction:

This is a penalty based on how far customers have to travel to charge their
vehicle and so represents charging anxiety. Let $d_{ij}$ denote the Euclidean
distance between the $i^{th}$ demand point and the $j^{th}$ supply point and
let $x_{ij}$ denote the amount of demand of the $i^{th}$ demand point that is
satisfied by the $j^{th}$ supply point, then the customer dissatisfaction cost
is given by

$C_{cd}\equiv\alpha\sum_{i}\sum_{j}x_{ij}d_{ij}$ (1)

Demand Mismatch:

This is the cost associated with incorrectly forecasting the EV charging
demand. Let $D_{i}$ denote the true demand and $\hat{D}_{i}$ the predicted
demand of demand point $i$ then

$C_{dm}\equiv\beta\sum_{i}|\hat{D}_{i}-D_{i}|$ (2)

Infrastructure Cost:

This is the cost of operating, maintaining and developing the EV charging
infrastructure. Let $N^{scs}_{j}$ denote the number of slow charging stations
and $N^{fcs}_{j}$ the number of fast charging stations at the $j^{th}$ supply
point then

$C_{ic}\equiv\gamma\sum_{j}(N^{scs}_{j}+rN^{fcs}_{j})$ (3)

The scaling parameters were given as $\alpha=1$, $\beta=25$, $\gamma=600$ and
$r=1.5$ which is the ratio of costs associated with developing Fast and Slow
charging stations.

###  2.3 Problem Constraints

The constraints we provided as follows:

  1. 1.

All elements of the demand supply matrix $\mathbf{X}$ must be non-negative.

  2. 2.

The number of slow $N^{scs}_{j}$ and fast $N^{fcs}_{j}$ chargers must be non-
negative.

  3. 3.

The total number of chargers placed at a supply point $j$ must not exceed the
number of parking spots which we denote by $p_{j}$.

  4. 4.

Charging stations can only be added (not removed) over time at a supply point.

  5. 5.

The total demand satisfied at a supply point cannot exceed the available
supply.

  6. 6.

The predicted demand at each demand point must be exactly satisfied.

##  3 Mathematical Models

The two aspects of the problem, demand estimation and placement optimization,
are independent and so can be solved independently. We provide the models that
were used for each and outline the solution approach.

###  3.1 EV Demand Prediction Model

For each demand point we are given actual demand over a number of years and we
need to predict the demand for the next two years. For each demand point we
can use the data from 2010 to 2018 to predict the value for 2019. The issue
was that the demand for a particular point may vary widely from year to year
which resulted in a poor prediction. We therefore decided to used points
surrounding the demand point to get a better estimate of its demand. For a
given year (with known demand) we instead estimate the demand for point $k$ as
follows:

$\tilde{D}_{k}=\frac{\sum_{j}\frac{D_{j}}{(1+{d_{jk})^{\kappa}}}}{\sum_{j}\frac{1}{1+d_{jk})^{\kappa}}}$
(4)

Note that when $\kappa=0$ we use the average demand over all demand points
while if $\kappa=\infty$ then $\tilde{D}_{k}=D_{k}$. This approach was adapted
from an approach used in [8].

These modified demands $\tilde{D}_{j}$ are computed for each year (2010-2017)
for a given value of $\kappa$ and then a third order polynomial regression
algorithm is used to predict demand for 2018. We then compared the true and
predicted demands for 2018 to obtain the Mean Square Error (MSE) for this
value of $\kappa$. This is repeated for different values of $\kappa$. We found
that $\kappa=4.7$ provided the lowest MSE. Finally we used this value of
$\kappa$ for years 2010-2018 to predict the demands for 2019 and then for
2020. These demands were then used in the charging station placement problem
which we next describe.

###  3.2 Charging Station Placement Optimisation Model

The optimal placement of chargers can be formulated as a Mixed-Integer
Programming Problem. Note that, since the true demands for 2019 and 2020 are
not known then we cannot include this cost in our model (although the cost was
used by the organizers when evaluating the solution. The optimization problem
and constraints can therefore be stated as:

$\min_{\mathbf{X},\vec{N^{scs}},\vec{N^{fcs}}}Z=\alpha\sum_{i,j}(x_{ij}d_{ij})+\gamma\sum_{j}(N^{scs}_{j}+rN^{fcs}_{j})$
(5)

Subject To:

$x_{ij}\geq 0\quad\forall\;i,j$ (6) $N^{scs}_{j}\in\\{0,1,2,\dots\\},\quad
N^{fcs}_{j}\in\\{0,1,2,\dots\\}\quad\forall\;j$ (7)
$N^{scs}_{j}+N^{fcs}_{j}\leq p_{j}\;\;\;\forall j$ (8)
$\sum_{i}x_{ij}\leq\Gamma_{scs}N^{scs}_{j}+\Gamma_{fcs}N^{fscs}_{j}\quad\forall\;j$
(9) $N^{scs}_{j}\geq N^{scs}_{j}(2018)\;\;\;N^{fcs}_{j}\geq N^{fcs}_{j}(2018)$
(10) $\sum_{j}x_{ij}=\hat{D}_{i}\quad\forall\;i$ (11)

where $r=1.5,\alpha=1$ and $\gamma=600$. The capacity of a slow charging
station is $\Gamma_{scs}=200$ and that of a fast charging station is
$\Gamma_{fcs}=400$. Figure 2 contains a network flow formulation of the
problem.

$\sum_{i}\tilde{D}_{i}$$\tilde{D}_{i}$$\tilde{D}_{0}$$\tilde{D}_{4095}$$\frac{r}{\Gamma_{fcs}}$$\text{cost=}\frac{1}{\Gamma_{scs}}$$j$$s$$f$$\text{flow=}x_{i}{}_{j}$$x_{s}{}_{j}$$x_{f}{}_{j}$$\text{cost=}d_{i}{}_{j}$$i$$T$$x_{Ts}$$99$
Figure 2: Network Flow Model

This minimisation problem was solved using the Mixed-Integer Programming
Solver in Julia using the Branch and Cut option and using the projected
demands obtained from Section 3.1 and data on the existing infrastructure. A
lower bound to this problem can be obtained by relaxing the integer
constraints and solving the resulting Linear Programming problem exactly. The
Mixed-Integer Programming solution (which may not be optimal) was within 0.01%
of the lower bound which means it was near-optimal.

##  4 Results and Leader board ranking

The information required to evaluate the demand prediction was only available
to the organizers and hence they alone could compute the demand mismatch cost.
However we have shown that, given demands, the placement optimization solution
is already very close to optimal (within 0.01% within the optimal objective
function value).

In terms of the demand estimation aspect. The prediction produced a result
that performed second in the academic category (meaning that one other
University had a better solution). Our approach outperformed methods employing
polynomial regressions combined with MILP, ANN, CNN and ARIMA models.

##  5 Conclusion and Future Work

We outlined an approach for determining the optimal placement of Electrical
Vehicle charging stations. This consisted of two subproblems. We first
estimated demand for charging stations based on historical data. We introduced
a novel approach for estimating demand that allows one to improve the
robustness of the estimate by using nearby samples. We formulated the charging
station placement optimization problem as a Mixed Integer Programming problem
and demonstrated that the solution obtained was within 0.01% of the optimal
solution. We are in the process of improving these algorithms using a grant
received from Shell. In the future we will report on any new results.

## References

  * [1] Programme, U.N.E.: Paris declaration on electro-mobility and climate change; call to action (2015).  URL https://wedocs.unep.org/20.500.11822/16838
  * [2] Whitehouse.gov: Fact sheet: The biden-harris electric vehicle charging action plan.  URL https://www.whitehouse.gov/briefing-room/statements-releases/2021/12/13/fact-sheet-the-biden-harris-electric-vehicle-charging-action-plan/
  * [3] GOV.UK: Uk electric vehicle infrastructure strategy.  URL https://www.gov.uk/government/publications/uk-electric-vehicle-infrastructure-strategy
  * [4] Lee, H., Clark, A.: Charging the future: Challenges and opportunities for electric vehicle adoption.  HKS (2018)
  * [5] Amini, M.H., Kargarian, A., Karabasoglu, O.: Arima-based decoupled time series forecasting of electric vehicle charging demand for stochastic power system operation.  Electric Power Systems Research 140, 378–390 (2016)
  * [6] Jahangir, H., Tayarani, H., Ahmadian, A., Golkar, M.A., Miret, J., Tayarani, M., Gao, H.O.: Charging demand of plug-in electric vehicles: Forecasting travel behavior based on a novel rough artificial neural network approach.  Journal of cleaner production 229, 1029–1044 (2019)
  * [7] Zhang, X., Chan, K.W., Li, H., Wang, H., Qiu, J., Wang, G.: Deep-learning-based probabilistic forecasting of electric vehicle charging load with a novel queuing model.  IEEE transactions on cybernetics 51(6), 3157–3170 (2021)
  * [8] Hosein, P.: A data science approach to risk assessment for automobile insurance policies.  International Journal of Data Science and Analytics pp. 1–12 (2023)