# Biogenic sulfur gases as biosignatures on temperate sub-Neptune waterworlds

Shang-Min Tsai Department of Earth and Planetary Sciences, University of
California, Riverside, CA, USA  Hamish Innes Atmospheric, Oceanic and
Planetary Physics, Department of Physics, University of Oxford, UK  Freie
Universität Berlin, Institute of Geological Sciences, Malteserstrasse 74-100,
12249 Berlin, Germany  Nicholas F. Wogan Space Science Division, NASA Ames
Research Center, Moffett Field, CA 94035  Virtual Planetary Laboratory,
University of Washington, Seattle, WA 98195  Edward W. Schwieterman Department
of Earth and Planetary Sciences, University of California, Riverside, CA, USA
Virtual Planetary Laboratory, University of Washington, Seattle, WA 98195
Blue Marble Space Institute of Science, Seattle, WA, USA

###### Abstract

Theoretical predictions and observational data indicate a class of sub-Neptune
exoplanets may have water-rich interiors covered by hydrogen-dominated
atmospheres. Provided suitable climate conditions, such planets could host
surface liquid oceans. Motivated by recent JWST observations of K2-18 b, we
self-consistently model the photochemistry and potential detectability of
biogenic sulfur gases in the atmospheres of temperate sub-Neptune waterworlds
for the first time. On Earth today, organic sulfur compounds produced by
marine biota are rapidly destroyed by photochemical processes before they can
accumulate to significant levels. Domagal-Goldman et al. (2011) suggest that
detectable biogenic sulfur signatures could emerge in Archean-like atmospheres
with higher biological production or low UV flux. In this study, we explore
biogenic sulfur across a wide range of biological fluxes and stellar UV
environments. Critically, the main photochemical sinks are absent on the
nightside of tidally locked planets. To address this, we further perform
experiments with a 3D GCM and a 2D photochemical model (VULCAN 2D (Tsai et
al., 2024)) to simulate the global distribution of biogenic gases to
investigate their terminator concentrations as seen via transmission
spectroscopy. Our models indicate that biogenic sulfur gases can rise to
potentially detectable levels on hydrogen-rich waterworlds, but only for
enhanced global biosulfur flux ($\gtrsim$20 times modern Earth’s flux). We
find that it is challenging to identify DMS at 3.4 $\mu m$ where it strongly
overlaps with CH4, whereas it is more plausible to detect DMS and companion
byproducts, ethylene (C2H4) and ethane (C2H6), in the mid-infrared between 9
and 13 $\mu m$.

††software: Exo-FMS (Lee et al., 2020, 2021), VULCAN (Tsai et al., 2017,
2021), Numpy (van der Walt et al., 2011)

##  1 Introduction

The Kepler mission revealed a large fraction of discovered planets with radii
between Earth and Neptune (Howard et al., 2012; Thompson et al., 2018; Bean et
al., 2021). Observed masses and radii indicate a population of sub-Neptunes
likely containing water-rich interiors, while evolution models also point to
the existence of these waterworlds on the population scale (Zeng et al., 2019;
Venturini et al., 2020; Luque & Pallé, 2022; Rogers et al., 2023). In
particular, several such planets orbiting M dwarf stars receive comparable
instellation111insolation from their host stars to Earth, placing them in
their planetary systems’ habitable zones (e.g., Table 1 in Madhusudhan et al.
(2021)). These temperate sub-Neptunes have sparked great interest in
atmospheric characterization and the assessment of habitability. Given
suitable climate conditions (Pierrehumbert, 2023; Innes et al., 2023), these
water-rich worlds could host surface water oceans below hydrogen-rich
atmospheres – recently referred to as “Hycean” (Hydrogen ocean) worlds
(Madhusudhan et al., 2021; Nixon & Madhusudhan, 2021). Promisingly, their
hydrogen-rich envelopes are more favorable for atmospheric characterization in
transmission observations due to large scale heights, compared to the higher
molecular weight N2- or CO2-dominated atmospheres.

Since the water detection in the atmosphere of the temperate sub-Neptune K2-18
b (Benneke et al., 2019; Tsiaras et al., 2019), follow-up analyses and
modeling efforts have sought to unveil the nature of K2-18 b and similar
temperate sub-Neptunes. Madhusudhan et al. (2020) explored the bulk
composition under different interior structures. Yu et al. (2021); Tsai et al.
(2021b) proposed using atmospheric compositional evolution to infer the
presence of a shallow surface, whereas Hu et al. (2021) investigated the
carbon inventories in a Hycean planet scenario. Recent JWST transit
observations revealed a strong detection of carbon dioxide (CO2) and methane
(CH4) on K2-18 b (Madhusudhan et al., 2023). Based on this basis and upper
limits of carbon monoxide (CO) and ammonia (NH3), Madhusudhan et al. (2023)
argued that the atmospheric characteristics of K2-18 b are more consistent
with a Hycean world scenario. However, several inconsistencies warrant
consideration: (1) It is challenging to accumulate $\%$ level of CH4 without
biological input or thermochemical recycling from the interior (Yu et al.,
2021; Tsai et al., 2021b). (2) On Hycean planets, CO can still be efficiently
produced to around 0.1 – 1 $\%$ level through photolysis of CO2, while the
quenched CO abundance is controlled by the interior temperature on Neptune-
like planet with a thick ($\gtrsim$ 100 bar) H2–atmosphere (Wogan et al.,
2024). (3) The cross-correlation method may be required to identify CO in JWST
NIRSpec data (Esparza-Borges et al., 2023), while current analysis based on a
single transit observation might not be sufficient to conclusively rule out
NH3 before the upcoming MIRI and revisit observations (JWST Cycle 1, GO-2722
and GO-2372).

Intriguingly, Madhusudhan et al. (2023) also reported a tentative detection of
dimethylsulfide (DMS), which is predominately produced by marine microbes on
Earth and regarded as a biosignature gas, especially for anoxic biospheres
(see Domagal-Goldman et al. (2011) and Schwieterman et al. (2018) for a
review). On modern Earth, DMS is the main biological source of sulfur, about
1.8–3.5 $\times$109 molecules cm-2s-1 (Seinfeld & Pandis, 2016; Cala et al.,
2023). DMS and methanethiol (CH3SH) are released from the degradation of
dimethylsulfoniopropionate (DMSP), an organosulfur compound produced by marine
phytoplankton. Microbial methylation and detoxification processes can also
actively produce DMS from hydrogen sulfide (H2S) (Li et al., 2023), allowing
hypothetical organisms to make use of the H2S reservoir in a H2-rich
environment on Hycean worlds. In addition, decayed organic substances in the
surface ocean can produce carbon disulfide (CS2) and carbonyl sulfide (OCS)
through photochemical processes. These biosulfur gases are not directly
associated with energy generation but are products or byproducts of
physiological responses, referred to as Type III Biosignatures in Seager et
al. (2013). On present-day Earth, these organic sulfur gases are readily
destroyed through photolysis or oxidation by OH radicals (Kettle et al., 2001;
Domagal-Goldman et al., 2011) and thus cannot accumulate to significant
concentrations exceeding the ppm level required for remote detection.

While great uncertainties persist regarding how biology could operate in an
H2-rich atmosphere, useful analogies can be drawn from anoxic environments on
Earth. As early microbial life emerged, the Archean atmosphere was more
reducing than present-day, with abundant CO2 and CH4 (e.g., Catling & Zahnle,
2020). As indicated by genomic analysis, organic sulfur cycling could date
back to the Proterozoic (Mateos et al., 2023), while sulfur-based metabolic
pathways likely emerged in Earth’s earliest biosphere in the Archean and may
have been common for biological evolution (Pilcher, 2003; House et al., 2003).
On Archean Earth-like worlds, most biogenic sulfur gases are rapidly destroyed
in the atmosphere similar to the present-day, except with the major sink being
the reaction with atomic O produced by photolysis (Domagal-Goldman et al.,
2011). Curiously, the potential sulfur-based biosignatures and chemical
pathways within a H2-dominated atmosphere have not yet been fully explored
with self-consistent photochemistry. Moreover, many close-in sub-Neptunes
around M stars might be tidally locked. On the permanent nightside, such
photochemical sinks of these organic sulfur gases do not exist. It is
therefore essential to quantify whether the biogenic gases can build up on the
terminators where transit observations probe.

In this Letter, we explore the global distribution and detectability of
methane- and sulfur-based biosignature gases on Hycean planets through
photochemical and climate modeling. While there is currently no robust
evidence that K2-18 b is a Hycean planet, we adopt its properties as a
paradigmatic example for Hycean exoplanets overall.

Figure 1: The average volume mixing ratios (VMRs) as a function of sulfur
biological flux (S${}_{\textrm{org}}$). We adopt the stellar spectrum of GJ
436 as an analogous star to K2-18 for our nominal Hycean K2-18 b model (top).
Additionally, we scaled the solar flux to match an equivalent flux (bottom).
Open circles are the average VMRs across the 1 to 10-4 bar pressure range,
whereas error bars denote the full span of VMRs within this pressure range.
The dashed line represents the limit of DMS given by surface deposition with
0.01 cm-1 dry deposition velocity. Combined with the deposition-free DMS
abundances in our model, the shaded blue region illustrates the upper and
lower bounds of DMS abundance. The black error bar marks the reported DMS VMR
from Madhusudhan et al. (2023) (M23). The orange band highlights the regime
where hydrocarbon haze precursors, C4H2 and C6H6, exceed $1\%$. Similarly, the
yellow band highlights the regime where elemental sulfur (S8) becomes
saturated and condenses.

##  2 Methods

Here, we model the atmospheric composition of a Hycean planet under a range of
biological surface emissions. Using the planetary parameters of K2-18 b as a
fiducial example, we vary the surface flux of biosulfur gases and stellar UV
flux to quantify the impacts on the global composition of Hycean planet
atmospheres.

###  2.1 1D climate and photochemical models

We iterate the radiative transfer (HELIOS; Malik et al. (2019a, b)) and
photochemical (VULCAN; Tsai et al. (2017, 2021)) models to obtain a self-
consistent 1D temperature profile of Hycean K2-18 b. The radiative transfer
model HELIOS has been modified to include the moist adiabatic lapse rate in
the convective adjustment steps to simulate atmospheres with condensing
components222https://github.com/exoclime/HELIOS/tree/development. For our
nominal Hycean planet, we assume a 1-bar H2–dominated atmosphere with elements
given by 100 $\times$ solar metallicity, following (Yu et al., 2021; Tsai et
al., 2021b). We fix the surface composition of CO2, N2, and H2O in the lower
boundary conditions. CO2 is set to 1$\%$ at the surface, based on the
retrieved CO2 abundance in Madhusudhan et al. (2023) and in accordance with
the carbon-inventory estimates in Kite & Ford (2018) and Hu et al. (2021).
Surface N2 is fixed to 0.006, corresponding to all nitrogen (from the 100
$\times$ solar metallicity) partitioning into N2. The water vapor is assumed
to be saturated above the ocean with a relative humidity of 50$\%$. The
spectrum of an analogous M2.5 star GJ 436 is adopted for both HELIOS and
VULCAN.

To achieve a marginally habitable surface temperature, we employ an ad hoc
surface albedo of 0.3, as previous studies have indicated that highly
reflective clouds are necessary to prevent a runaway greenhouse state on K2-18
b (Piette & Madhusudhan, 2020; Innes et al., 2023). Similarly, we do not
account for the possibility of convection inhibition, which would otherwise
raise the surface temperature beyond the supercritical point of water (Innes
et al., 2023; Pierrehumbert, 2023). This allowed us to model a habitable
atmospheric state supportive of liquid surface water. The temperature and
composition profiles are iterated between HELIOS and VULCAN until the
temperature change is about 1$\%$. The planetary parameters for our K2-18 b
Hycean model are listed in Table A1.

The self-consistent temperature and composition structures of our “lifeless”
K2-18 b Hycean planet (excluding biological surface fluxes) are illustrated in
Figure A1. We then fed the temperature profile into VULCAN to model the
compositional evolution in response to various biological surface fluxes. We
updated the S–N–C–H–O photochemical network by revising the reactions
involving H2CO (see Appendix B and Wogan et al. (2024)) and introduced
reactions related to dimethyl sulfide (DMS) and dimethyl disulfide
(DMDS)333https://github.com/exoclime/VULCAN/blob/master/thermo/SNCHO_DMS_photo_network_Tsai2024.txt.
We consider CH4 released from methanogenesis and a range of organic sulfur
gases biologically emitted from the surface, including CH3SH, CS2, DMS,
dimethyl disulfide (DMDS), and OCS, collectively termed “S${}_{\textrm{org}}$”
after Domagal-Goldman et al. (2011). During this stage, the temperature
profile from lifeless Hycean is held fixed and we neglect the thermal feedback
from the compositional variations due to biological input (mainly CH4 and
S${}_{\textrm{org}}$).

The rates of gases deposited to the ocean depend on the aerodynamics at the
atmosphere-ocean interface and the solubility of the gas (Hu et al., 2012;
Seinfeld & Pandis, 2016). Due to the uncertainty of parameterizing the dry
deposition velocities of several organic sulfur compounds, we assume zero dry
deposition for CH3SH, DMS, DMDS, and CS2, corresponding to the upper limit of
biosignature accumulation where biogenic sulfur compounds saturate the oceanic
mixed layer. The lower boundary conditions with modern Earth global average
S${}_{\textrm{org}}$ flux are described in Appendix A and Table A.

Figure 2: Left: The stellar spectra in the UV for GJ 436 (M2.5), a young (45
Myr) and active M star, Epsilon Eridani, and the Sun in the top panel. The
spectra of GJ 436 and Epsilon Eridani are from the MUSCLES survey (version
2.2) (France et al., 2016; Youngblood et al., 2016; Loyd et al., 2016), the
active M star spectrum is from HAZMAT (Peacock et al., 2020), and the solar
spectrum is from Gueymard (2018). The lower panel shows the UV photosphere
($\tau$ = 1) with different S${}_{\textrm{org}}$. The contributions from DMS,
DMDS, and CH3SH are further displayed for the 30$\times$ S${}_{\textrm{org}}$
model. Right: The upper panel shows the average VMRs of several sulfur species
as a function of attenuating total stellar UV flux for Earth’s
S${}_{\textrm{org}}$ flux. The bottom panel shows the span of abundance
averaged between 1 and 10-4 bar for 1$\times$ – 20$\times$
S${}_{\textrm{org}}$ flux around different stellar types

###  2.2 3D GCM 2D photochemical model

In the absence of photochemical destruction, biological sulfur compounds can
in principle build up on the night side of a tidally-locked Hycean planet.
Global circulation plays a central role in transporting species between the
dayside and nightside, where distinct chemical sources and sinks are at play
(Chen et al., 2018; Tsai et al., 2021b, 2023c). To account for the day-night
transport, we run VULCAN 2D (Tsai et al., 2021b, 2024) to model the equatorial
region, using temperature and wind structures from the 3D GCM ExoFMS (Lee et
al., 2021; Innes & Pierrehumbert, 2022) as input.

ExoFMS has previously been used to model the climate of a warmer K2-18 b
(surface temperature $>$ 500 K) without a liquid ocean and water condensation
(Innes & Pierrehumbert, 2022). In this study, we fit the shortwave and
longwave band opacities in ExoFMS to match our 1D HELIOS correlated-k results.
Therefore, the double-gray opacities in our Exo-FMS roughly represent the
opacities from the self-consistent 100$\times$ solar metallicity composition.
The input parameters are listed in Table A1 and the temperature and wind
structures can be found in Figure A2.

The temperatures and winds obtained from exo-FMS are averaged over the
equatorial region across 30∘ and divided into 32 longitude columns to set up
VULCAN 2D. In VULCAN 2D, the zonal winds from the GCM are directly implemented
for east-west advection, while vertical transport is parameterized by the eddy
diffusion coefficient ($K_{\textrm{zz}}$) derived from the mixing length
theory with $K_{\textrm{zz}}$ = 0.1 H $\times w_{\textrm{rms}}$, assuming
0.1$\times H$ (scale height) as the characteristic length (Smith, 1998;
Charnay et al., 2015). We use the same lower boundary conditions as the 1D
model uniformly across all longitudes.

##  3 Results

###  3.1 1D: S${}_{\textrm{org}}$ flux $\gtrsim$ 20$\times$ modern Earth value
is required to accumulate ppm-levels of DMS

Figure 1 illustrates the average mixing ratios of several carbon-bearing and
organic sulfur species in response to increases of biological production of
sulfur gases. With modern Earth’s S${}_{\textrm{org}}$ flux, surface DMS
reaches abundances of around 1–30 ppb, before being destroyed above 0.01 bar
in our Hycean K2-18 b model. The DMS concentrations on Hycean K2-18 b are
slightly higher than modern and Archean Earth’s values, which are around
0.01–1 ppb levels (Zhang et al., 2020; Domagal-Goldman et al., 2011). As the
S${}_{\textrm{org}}$ flux continues increasing, there is an abrupt transition
between 20–30 times S${}_{\textrm{org}}$ flux in our
S${}_{\textrm{org}}$-deposition-free setup. After this threshold, the
atmosphere becomes rich in organic sulfur species and poor in CO, while CH3SH,
DMS, and DMDS jump from sub-ppm level to above 0.1 $\%$. In the more
conservative scenario where the ocean is not saturated, DMS is expected to be
limited by the surface deposition and increased linearly with increasing
S${}_{\textrm{org}}$, as shown by the dashed line in Figure 1. Once
S${}_{\textrm{org}}$ exceeds 20–30$\times$ Earth’s value, H2S begins
accumulating in the upper atmosphere, as a result of photochemistry reducing
sulfur back to its thermochemically favored state. This occurs when CH3SH
builds up to sufficiently high abundances that its photodissociation becomes
significant, producing H2S via the following pathway

$\displaystyle\begin{aligned} \ce{CH3SH&->[h\nu]CH3+SH}\\\
\ce{CH3+SH&->CH4+S}\\\ \ce{S+H2&->[M]+H2S}\\\ \vglue 3.0pt\cr\end{aligned}$
(1)

In fact, the ultimate fate of S${}_{\textrm{org}}$ species in a Hycean
atmosphere is either photochemical oxidation into CH3SO or reduction back to
hydrocarbons and H2S, where CH3SO is mainly produced by DMDS reacting with
atomic O. The path to H2S becomes more dominant over CH3SO as
S${}_{\textrm{org}}$ flux increases and makes the atmosphere more reducing.
Although SO2 is also photochemically produced, it is highly soluble and
deposition into the ocean limits abundances to insignificant levels.

Notably, the atmosphere also becomes rich in hydrocarbons as
S${}_{\textrm{org}}$ increases, as pointed out in Domagal-Goldman et al.
(2011). CO follows an opposite trend to S${}_{\textrm{org}}$, showing an
abrupt decline after 20 times modern Earth’s S${}_{\textrm{org}}$ flux. In
this carbon-rich regime, the high abundances of haze precursors, such as C4H2
and C6H6 (Tsai et al., 2021), imply that the atmosphere would be covered by
photochemical hazes (Arney et al., 2018). A main pathway to produce organic
haze precursors in the upper atmosphere initiated by DMS photolysis is

$\displaystyle\begin{aligned} \ce{DMS(CH3SCH3)&->[h\nu]CH3+CH3S}\\\
\ce{CH3+CH3&->[M]C2H6}\\\ \ce{C2H6&->[M]C2H4+2H}\\\
\ce{C2H4&->[h\nu]C2H2+2H}\\\ \ce{C2H2&->[h\nu]C2H+H}\\\
\ce{C2H2+C2H&->C4H2+H}\\\ \ce{C4H2+H&->[M]C4H3}\\\ \ce{C2H2+C4H3&->[M]C6H5}\\\
\ce{C6H5+H2&->C6H6+H}\\\ \vglue 3.0pt\cr\end{aligned}$ (2)

For DMS at low concentrations below ppm-level, its primary sink is direct
photodissociation. After the transition to where DMS has a high concentration,
its photodissociation halts in the lower atmosphere due to shielding by DMDS,
CH3SH, and self-shielding. The reaction with atomic H becomes the main sink.
This shielding effect from S${}_{\textrm{org}}$ also reduces the photolysis of
CO2 causing the drop of CO when the flux $\gtrsim$ 20$\times$
S${}_{\textrm{org}}$. The overall impact of shielding across this transition
is evident in the top left panel of Figure 2, where DMS, DMDS, and CH3SH
significantly elevate the near-UV (NUV) photosphere and subsequently attenuate
the NUV flux reaching the surface.

For an equivalent solar flux, the same average abundances of hydrocarbon and
sulfur gases as a function of biological sulfur flux are shown in the bottom
of Figure 1. One major difference is that the stronger NUV irradiation
generates atomic H down to the surface, which significantly suppresses DMS and
other organic sulfur gases. DMS remains below 10 ppb level for
S${}_{\textrm{org}}$ flux up to 200 times modern Earth’s value. Moreover, the
concentration of DMS is confined near the surface, whereas for an M-dwarf host
star, DMS can extend vertically above 0.01 bar. Since S${}_{\textrm{org}}$
never reaches significant abundances under a Sun-like star, shielding of CO2
photolysis does not occur. Consequently, both CH4 and CO progressively
increase with S${}_{\textrm{org}}$ flux. Our model results with an active M
star show qualitatively similar trends to those with the equivalent solar flux
owing to its enhanced UV flux.

Figure 3: The equatorial volume-mixing-ratio distribution of DMS (left) and
DMDS (middle) in our Hycean K2-18 b model with modern Earth
S${}_{\textrm{org}}$ flux, computed by VULCAN 2D. The substellar point is
located at 180∘ longitude. The right panel depicts the vertical distribution
of several species across the equatorial region, where different shades
correspond to different longitudinal locations. The vertical profiles of DMS
excluding zonal transport are shown in red dotted (morning terminator, 73∘
westward sub-stellar point) and dashed-dotted (evening terminator, 73∘
eastward sub-stellar point) lines for comparison.

###  3.2 2D: DMS is homogenized by global circulation

The temperature and wind structures simulated by our Hycean K2-18 b GCM can be
found in Figure A2. The circulation exhibits an eastward zonal jet in the
troposphere and day-night thermally driven cells in the stratosphere. In the
upper atmosphere above 10-4 bar, the flow transitions to retrograde, likely
due to vertical momentum transport (Tsai et al., 2014). The mean tropospheric
zonal wind below 0.01 bar is about 20 m/s, consistent with the jet speed
estimated from wave balance (Hammond et al., 2020). This zonal wind speed
translates to a horizontal transport timescale of $\sim$ 10 days.

Figure 3 illustrates the equatorial distribution of DMS and DMDS for modern
Earth S${}_{\textrm{org}}$ production. DMS exhibits a rather uniform abundance
below 10-4 bar. As the zonal winds change to retrograde in the upper
atmosphere, DMS extends to higher altitudes around the evening terminator due
to the westward transport from the nightside. On the other hand, DMDS is
photochemically produced by CH3SH via

$\displaystyle\begin{aligned} \ce{CH3SH&->[h\nu]CH3S+CH3}\\\
\ce{CH3S+CH3S&->DMDS(CH3S2CH3)}\end{aligned}$ (3)

in the stratosphere between 10-2–10-4 bar and follows the thermally driven
cells. This relatively uniform distribution of DMS can be understood by
comparing relevant timescales. The photochemical lifetime of DMS on the
dayside ($\tau_{DMS}^{day}$) is determined by
$\tau_{\textrm{DMS}}^{\textrm{day}}=1/k_{\textrm{DMS}}$ , where
$k_{\textrm{DMS}}$ is its photolysis rate. $\tau_{DMS}^{day}$ is $\sim$
108–109 s in the troposphere below 10-2 bar. On the nightside, the replenish
timescale ($\tau_{DMS}^{night}$) is given by $\sim$ $H/K_{\textrm{zz}}^{2}$
$\sim$ 108 s. Given that $\tau_{DMS}^{day}$ and $\tau_{DMS}^{night}$ are
comparable to each other and considerably longer than the horizontal transport
timescale of $\sim$ 106 s, the S${}_{\textrm{org}}$ species on the nightside
is still constrained by the photochemical processes on the dayside through
global transport. We note that our 2D model already assumes an optimistic
scenario with a uniform S${}_{\textrm{org}}$ flux across the nightside. In
reality, the nightside S${}_{\textrm{org}}$ flux might be much lower due to
the lack of stellar energy input to power a (hypothetical) photosynthetic
biosphere. Overall, the DMS abundance predicted by our 2D photochemical model
is broadly consistent with that from 1D model, as indicated by the right panel
in Figure 3.

Figure 4: Synthetic transmission spectra of Hycean K2-18 b 1D models for
1$\times$ and 20$\times$ modern Earth biological sulfur flux. The JWST
NIRISS/SOSS and NIRSpec/G395H data from (Madhusudhan et al., 2023) are plotted
in the left panel. The simulated noise by Pandexo (Batalha et al., 2017) with
the JWST/MIRI LRS for 5 transits without random noise is plotted in the right
panel (gray error bars). The 20$\times$ S${}_{\textrm{org}}$ spectrum but with
DMS opacity removed is shown for comparison. The “no S${}_{\textrm{org}}$”
case includes a flux for CH4 equivalent to Modern Earth’s non-anthropogenic
emission (Domagal-Goldman et al., 2011), but omits fluxes for all biogenic
sulfur gases. The corresponding abundance profiles used to produce the spectra
can be found in Figure A3.

###  3.3 Transmission spectra

We simulated synthetic transmission spectra of our cloud-/haze-free Hycean
K2-18 b models using PLATON (Zhang et al., 2019, 2020), including opacity
sources of CH4, CO,CO2, C2H2, H2O,HCN, NH3, O2, NO, C2H4, C2H6, H2CO, N2, NO2,
H2S, COS, SH, SO2, (CIA) of H2–H2 and H2–He. We further include cross-sections
of CH3SH, DMS, and DMS from the Pacific Northwest National Laboratories (PNNL)
database (Johnson et al., 2004), which are limited to 1 bar pressure and room
temperature.

The transmission spectra for none, 1$\times$ and 20$\times$ modern Earth
S${}_{\textrm{org}}$ flux are compared in Figure 4. For the JWST NIRISS and
NIRSpec coverage below 5 $\mu m$, all of our Hycean K2-18 b models with
methanogenic flux produce CH4 absorption features between 1 and 4 $\mu m$ that
are broadly consistent with the JWST data. However, it is challenging to
distinguish the DMS absorption contribution within the NIRSpec coverage, even
with 20$\times$ S${}_{\textrm{org}}$ flux. The signal of DMS around 3.5 $\mu
m$ for 20$\times$ S${}_{\textrm{org}}$ flux does not exceed 10 ppm, smaller
than NIRSpec’s precision of 30–50 ppm at R $\sim$ 50 (Madhusudhan et al.,
2023). The most promising diagnostics window of DMS is in the mid-infrared, as
indicated in the right panel of Figure 4. For the 20$\times$
S${}_{\textrm{org}}$ case, DMS and the photochemical byproduct C2H4 contribute
to an absorption of about 50 ppm at 9–11 $\mu m$, in addition to the C2H6
feature at 7$\mu m$. The 7–$\mu m$ feature of C2H6 also overlaps with CS2,
making it difficult to differentiate the contribution of CS2. Our transmission
spectra reveal no discernible H2S spectral signatures. While H2S can reach
above ppm levels with outgassing and photochemical sources (e.g. pathway (1)),
its broad opacities that overlap with several species make it challenging to
detect.

Since ethane (C2H6) is efficiently produced from CH4 and has strong mid-
infrared absorption bands, C2H6 shows a prominent feature at 7 and 11–13 $\mu
m$ with or without the presence of S${}_{\textrm{org}}$ flux. By constraining
[CH4]/[C2H6], one can potentially estimate the relative contribution of
methanogenic CH4 flux relative to S${}_{\textrm{org}}$ flux (Domagal-Goldman
et al., 2011). Since our model likely overpredicts C6H6 in the gas phase due
to the lack of conversion to higher-order hydrocarbons and particles, we did
not include C6H6 opacities when generating the spectra in Figure 4. However,
it is interesting to keep in mind that the presence ($\gtrsim$ 50 ppm) of C6H6
would exhibit detectable features within NIRSpec wavelength range, in addition
to CH4 and CO2 (see Figure A4). For biogenic flux $\gtrsim$ 30$\times$ Earth’s
S${}_{\textrm{org}}$ flux, the atmospheric mean molecular weight exceeds 8 amu
in addition to the expected aerosols that would flatten transmission spectral
features, making it immediately distinguishable for $\lesssim$ 20$\times$
cases. We summarize the identifiable biosignature features with
1$\times$–20$\times$ S${}_{\textrm{org}}$ flux in Table 1.

Table 1: Detectability of potential biosignature gases with corresponding wavelength ranges on Hycean K2-18 b based on our models with 1$\times$–20$\times$ S${}_{\textrm{org}}$ flux.  ✓ denotes detectable ($\gtrsim$ 50 ppm) and  ✗ denotes nondetectable ($\lesssim$ 10 ppm), whereas ? means challenging, with shallow features ($\gtrsim$ 20 ppm but $\lesssim$ 50 ppm) and only present under enhanced 20$\times$ S${}_{\textrm{org}}$ flux. Species | 1–5 $\mu m$ | 5–15 $\mu m$
---|---|---
CH4 | ✓ | ✓
CO2 | ✓ | ✓
H2S | ✗ | ✗
CH3SH | ✗ | ✗
CS2 | ✗ | ?
DMS | ✗ | ?
DMDS | ✗ | ✗
C2H4 | ✗ | ?
C2H6 | ✗ | ✓

##  4 Discussion

Since K2-18 b requires reflective clouds and hazes to sustain water oceans, we
explore the potential UV shielding effects due to clouds and hazes in the
bottom panel of Figure 2. Organic sulfur species only start to build up when
the stellar UV is reduced to less than 20$\%$. For K2-18 b specifically,
Leconte et al. (2024) highlighted the observed methane features deeper than
$\sim$ 100 ppm (their Fig. 9) are not consistent with highly reflective
clouds. We argue that our photochemical results should remain robust to UV
attenuation by clouds and hazes, unless muted spectral features indicate
scattering by high-altitude aerosols can significantly reduce stellar UV flux.

The predicted requirement of 20$\times$ S${}_{\textrm{org}}$ flux to be
potentially detectable is equivalent to about 0.004 g m-2 based on the DMS lab
production rate (Seager et al., 2013), which is within the plausible range of
surface biomass density. Our predictions for the biogenic sulfur gases in
Hycean atmospheres are broadly consistent with previous studies of the Archean
Earth Domagal-Goldman et al. (2011) and Archean-like atmospheres on Trappist-1
planets (Meadows et al., 2023). While Meadows et al. (2023) concluded CH3SH to
be the most prominent S${}_{\textrm{org}}$ feature in N2-dominated
atmospheres, we find efficient CH3SH conversion into H2S in a H2 atmosphere.
Consequently, DMS appears to constitute the most detectable organic sulfur
molecule in the mid-infrared for Hycean worlds.

On modern Earth, the oxidation of DMS might participate in the sulfuric cloud
formation (Charlson et al., 1987) and affect climate feedback. It has been
suggested that this biologically driven feedbacks could help stabilize the
Earth’s climate (Lovelock, 1989), although the precise pathways and quantified
contributions from DMS to aerosols still remain elusive. Similarly, on a
Hycean world with highly elevated S${}_{\textrm{org}}$ flux
($\gtrsim$30$\times$ modern Earth’s value, see Figure 1.), hydrocarbon and
elemental sulfur hazes (S8) are expected to form. We speculate a similar
negative feedback loop can take place in a Hycean world, where photochemical
hazes derived from organic sulfur flux can increase albedo and prevent a
runaway state.

Lastly, we reiterate that constructing complete line lists for DMS, DMDS, and
CH3SH to correctly account for the pressure-broadening effects would improve
the interpretation of spectral data. Future work extending the study
encompassing the diversity of halomethanes, such as methyl chloride and methyl
bromide (e.g., Leung et al., 2022), will expand our understanding of the
biosignature gases on sub-Neptune waterworlds.

##  5 Conclusion

We find a biogenic sulfur flux of approximately 20 times higher than modern
Earth’s is required for DMS to reach detectable levels (i.e. above $\sim$ppm)
on a K2-18b-like Hycean world. At the terminators, our 2D photochemical model
shows only minor abundance enhancements in DMS compared to 1D model
predictions. On the other hand, JWST/MIRI could detect C2H6 around 7 $\mu m$
and the broad feature at 9–13 $\mu m$ contributed from DMS, C2H4, and C2H6
with 20$\times$ S${}_{\textrm{org}}$ flux, based on our fiducial Hycean model.
However, identifying DMS signatures within NIRSpec’s wavelength range is
demonstrated to be challenging, despite DMS appearing as the most promising
biogenic sulfur signature directly from the S${}_{\textrm{org}}$ source. The
moderate threshold for biological production suggests that the search for
biogenic sulfur gases as one class of potential biosignature is plausible for
Hycean worlds.

The authors thank Michaela Leung for the useful comments on an early draft of
this article. The authors also thank Stephen Klippenstein and Julie I. Moses
for the illuminating discussion on H2CO kinetics and for sharing the ab initio
calculations. S.-M.T. thanks Eva-Maria Ahrer, Guangwei Fu, Jake Taylor, and
Lili Alderson for helpful discussions on data resampling and feature
identification. S.-M.T. and E.W.S. acknowledge support from National
Aeronautics and Space Administration (NASA) through Exobiology Grant No.
80NSSC20K1437 and Interdisciplinary Consortia for Astrobiology Research (ICAR)
Grant Nos. 80NSSC23K1399, 80NSSC21K0905, and 80NSSC23K1398. H.I. is funded by
the European Union (ERC, DIVERSE, 101087755 and ERC, EXOCONDENSE, 740963).
Views and opinions expressed are however those of the author(s) only and do
not necessarily reflect those of the European Union or the European Research
Council Executive Agency. Neither the European Union nor the granting
authority can be held responsible for them. N.F.W. was supported by the NASA
Postdoctoral Program.

##  Appendix A Model input parameters

Table A1 lists the input parameters for Hycean K2-18 b with the HELIOS
radiative transfer model and the ExoFMS general circulation model. HELIOS
includes H2O, CH4, CO, CO2, NH3, HCN, C2H2, collision-induced absorption (CIA)
of H2–H2 and H2–He as opacity sources. We ignore the H2O self continuum
absorption and likely underestimate the tropospheric temperature. ExoFMS uses
double-gray opacities for its radiative transfer. The optical depth in the
infrared ($\tau_{\text{IR}}$) and shortwave ($\tau_{\text{SW}}$) bands is
given by:

$d\tau_{\text{i}}=\quantity(\kappa_{\text{i,d}}(1-q)+\kappa_{\text{i,w}}q)\quantity(f_{\text{i}}+2(1-f_{\text{i}}))\quantity(\frac{p}{p_{s}}))\frac{dp}{g},$
(A1)

where the $i$ subscript represents either the SW or IR band, $p$ is the
pressure and $p_{s}$ the surface pressure. The opacities $\kappa_{i,w}$ and
$\kappa_{i,d}$ represent the opacities due to the water vapour component and
dry component (all other gases) of the atmosphere respectively. The factor
$f_{i}$ accounts for pressure dependence of opacity, due to collision-induced
absorption and pressure broadening effects. To reproduce the inversion in the
HELIOS calculations, a second band was added in the SW part of the spectrum,
receiving a fraction $\alpha$ of the total incoming flux. The optical depth in
this band was given by:

$\tau_{SW,2}=\frac{\kappa_{\text{SW,2}}p}{g}.$ (A2)

The water mass concentration, $q$, used in this calculation is taken from the
1D HELIOS calculation and remains constant throughout the simulation (though
in reality this should be spatially inhomogeneous and linked to the local
saturation vapour pressure at a given temperature and pressure). The value of
$\kappa_{IR,w}$ was taken to match the runaway greenhouse OLR of $\approx 280$
Wm-2 (as in Innes et al., 2023). The other values were tuned to match the
HELIOS temperature-pressure profile. The parameters used are shown in Table
A2.

Table A1: Planetary parameters for Hycean K2-18 b. Symbol | Value | Unit | Description
---|---|---|---
Rp | 1.66 $\times$ 107 | m | Planetary radius
a | 0.1591 | AU | Orbital distance
Tirr | 1370 | K | Irradiation temperature
Tint | 0 | K | Internal temperature
g | 12.43 | m s-2 | Gravitational acceleration
cp | 2216 | J K-1 kg-1 | Specific heat capacity
$\Omega_{\rm p}$ | 2.2137 $\times$ 10-6 | rad s-1 | Planetary rotation rate
Table A2: Radiative parameters for the ExoFMS double-gray radiation scheme. Symbol | Value | Unit | Description
---|---|---|---
$\kappa_{IR,d}$ | $1.15\text{\times}{10}^{-3}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ | IR dry opacity
$\kappa_{\text{IR,w}}$ | $1.0\text{\times}{10}^{-2}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ | IR water opacity
$\kappa_{SW,d}$ | $2.0\text{\times}{10}^{-5}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ | SW dry opacity
$\kappa_{\text{SW,w}}$ | $4.0\text{\times}{10}^{-4}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ | SW water opacity
$\kappa_{\text{SW,2}}$ | $5.0\text{\times}{10}^{-2}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{kg}}^{-1}$ | SW second band opacity
$\alpha$ | $1.0\text{\times}{10}^{-2}$ | Dimensionless | Fraction of instellation in SW second band
$f_{\text{IR}}$ | 0.9 | Dimensionless | Controls pressure dependence of $\tau_{\text{IR}}$
$f_{\text{SW}}$ | 0.2 | Dimensionless | Controls pressure dependence of $\tau_{\text{SW}}$
Table A3: Lower boundary conditions for photochemical modeling of inhabited
Hycean K2-18 b.

Species | Surface Emission (molecules cm-2 s-1) | Vdep (cm s-1)
---|---|---
CH4${}^{\text{a}}$ | 7$\times$ 1010 | 0
CH3SH${}^{\text{a}}$ | 8.3$\times$ 108 | 0
DMS${}^{\text{a}}$ | 4.2$\times$ 109 | 0${}^{\text{d}}$
CS2${}^{\text{a}}$ | 1.4$\times$ 107 | 0
H2S${}^{\text{b}}$ | 2$\times$ 108 | 0.015
SO2${}^{\text{b}}$ | 9$\times$ 109 | 1
COS${}^{\text{b}}$ | 5.4$\times$ 107 | 0.003
H2O2${}^{\text{c}}$ | 0 | 1
CH3S${}^{\text{a}}$ | 0 | 0.01
HSO${}^{\text{a}}$ | 0 | 1
S${}^{\text{a}}$ | 0 | 1
SO${}^{\text{a}}$ | 0 | 0.0003
${}^{\text{a}}$ Domagal-Goldman et al. (2011)
${}^{\text{b}}$ Seinfeld & Pandis (2016) |  |
${}^{\text{c}}$ Hauglustaine et al. (1994) |  |
${}^{\text{d}}$ For the scenario that the surface layer of the ocean is saturated. The measured deposition velocity of DMS on land is 0.064–0.28 cm s-1 (Judeikis & G. Wren, 1977). |  |

Figure A1: The self-consistent temperature (left) and composition (right)
profiles generated by iterative runs between HELIOS and VULCAN. Sulfur
kinetics is not included in the VULCAN runs here for simplicity.

Figure A2: The top left panel shows the temperatures (color scale) and winds
(arrows) on the equatorial plane averaged across $\pm$30∘ latitudes from our
K2-18 b GCM (substellar point located at 180∘ longitude). The top right panel
displays the vertical temperature profiles around the equator corresponding to
different longitudinal locations, where the temperature profile from the 1D
radiative transfer calculation (HELIOS) is plotted in black for comparison.
The bottom two panels show the temperatures and winds at 0.1 bar and 0.1 mbar
level.

Figure A3: The abundance profiles of the Hycean K2-18 b 1D models with
1$\times$ (solid) and 20$\times$ (dashed) modern Earth values.

Figure A4: Synthetic transmission spectra of Hycean K2-18 b similar to Figure
4, but including C6H6 absorption for 20$\times$ S${}_{\textrm{org}}$ flux. The
spectra without C6H6 opacity and with 10$\%$ of C6H6, representing depletion
from haze formation, are shown for comparison.

##  Appendix B Updated H2CO reaction rates

Formaldehyde (H2CO) is one of the main intermediate species in CO2-CO-CH4
interconversion. Forming and breaking the double bond between C and O in H2CO
may facilitate an important step that controls the overall timescale (Tsai et
al., 2018). There are two three-body reactions involving H2CO that we adopt
the thermal dissociation rate coefficients due to the lack of low-pressure
rate limits for recombination:

$\ce{CH2OH->[M]H+H2CO}$ (B1)

and

$\ce{CH3O->[M]H+H2CO}$ (B2)

The recombination rates are then reversed by thermochemical data. In the
default chemical networks in VULCAN (e.g. NCHO_photo_network), the adopted
rate coefficients either underestimate the activation energy or are only valid
at high temperatures. As a result, these reversed rates tend to overpredict
the CO2–CH4 conversion efficiency around room temperatures. In this study, we
have adopted the low-pressure rate coefficients for B1 and B2 from Tsang &
Hampson (1986) and Tsang (1987), which are more consistent with the ab initio
calculations conducted by Klippenstein (2023). The original and updated rate
coefficients are summarized in Table A4. This update of H2CO kinetics
significantly lowers the CH4 abundance in our lifeless case (without
methanogenic flux; Figure A1)) from $\sim$10-5 to $\sim$10-8, which is
consistent with the calculations in Wogan et al. (2024).

Table A4:  Reactions |  Rate Coefficients${}^{\text{a}}$ |  Reference
---|---|---
Original
CH2OH -¿[M] H + H2CO |  $k_{0}=1.66\times 10^{-10}\exp(-12630/T)$
$k_{\infty}=3\times 10^{9}\exp(-14640/T)$ |  Cribb et al. (1992)
Tsuboi et al. (1981)
CH3O -¿[M] H + H2CO |  $k_{0}=9\times 10^{-11}\exp(-6790/T)$
$k_{\infty}=1.56\times 10^{15}T^{-0.39}\exp(-13300/T)$ |  Baulch et al. (1994)
Curran (2006)
Updated
CH2OH -¿[M] H + H2CO |  $k_{0}=7.48\times 10^{1}T^{-2.5}\exp(-17200/T)$
$k_{\infty}=4.53\times 10^{34}T^{-7.11}\exp(-22200/T)$ |  Tsang (1987)
Xu et al. (2015)
CH3O -¿[M] H + H2CO |  $k_{0}=6.51\times 10^{13}T^{-6.65}\exp(-16700/T)$
$k_{\infty}=3.17\times 10^{24}T^{-4.25}\exp(-13100/T)$ |  Tsang & Hampson (1986)
Xu et al. (2015)
${}^{\text{a}}$In low pressure rate coefficient, $k_{0}$ (cm6 molecules-2
s-1), and high pressure rate coefficient, $k_{\infty}$ (cm3 molecules-1 s-1).

## References

  * Arney et al. (2018) Arney, G., Domagal-Goldman, S. D., & Meadows, V. S. 2018, Astrobiology, 18, 311, doi: 10.1089/ast.2017.1666
  * Batalha et al. (2017) Batalha, N. E., Mandell, A., Pontoppidan, K., et al. 2017, PASP, 129, 064501, doi: 10.1088/1538-3873/aa65b0
  * Baulch et al. (1994) Baulch, D. L., Cobos, C. J., Cox, R. A., et al. 1994, Journal of Physical and Chemical Reference Data, 23, 847, doi: 10.1063/1.555953
  * Bean et al. (2021) Bean, J. L., Raymond, S. N., & Owen, J. E. 2021, Journal of Geophysical Research (Planets), 126, e06639, doi: 10.1029/2020JE006639
  * Benneke et al. (2019) Benneke, B., Wong, I., Piaulet, C., et al. 2019, ApJ, 887, L14, doi: 10.3847/2041-8213/ab59dc
  * Cala et al. (2023) Cala, B. A., Archer-Nicholls, S., Weber, J., et al. 2023, Atmospheric Chemistry & Physics, 23, 14735, doi: 10.5194/acp-23-14735-2023
  * Catling & Zahnle (2020) Catling, D. C., & Zahnle, K. J. 2020, Science Advances, 6, eaax1420, doi: 10.1126/sciadv.aax1420
  * Charlson et al. (1987) Charlson, R. J., Lovelock, J. E., Andreae, M. O., & Warren, S. G. 1987, Nature, 326, 655, doi: 10.1038/326655a0
  * Charnay et al. (2015) Charnay, B., Meadows, V., & Leconte, J. 2015, ApJ, 813, 15, doi: 10.1088/0004-637X/813/1/15
  * Chen et al. (2018) Chen, H., Wolf, E. T., Kopparapu, R., Domagal-Goldman, S., & Horton, D. E. 2018, ApJ, 868, L6, doi: 10.3847/2041-8213/aaedb2
  * Cribb et al. (1992) Cribb, P. H., Dove, J. E., & Yamazaki, S. 1992, Combustion and Flame, 88, 169, doi: https://doi.org/10.1016/0010-2180(92)90050-Y
  * Curran (2006) Curran, H. J. 2006, International Journal of Chemical Kinetics, 38, 250, doi: https://doi.org/10.1002/kin.20153
  * Domagal-Goldman et al. (2011) Domagal-Goldman, S. D., Meadows, V. S., Claire, M. W., & Kasting, J. F. 2011, Astrobiology, 11, 419, doi: 10.1089/ast.2010.0509
  * Esparza-Borges et al. (2023) Esparza-Borges, E., López-Morales, M., Adams Redai, J. I., et al. 2023, ApJ, 955, L19, doi: 10.3847/2041-8213/acf27b
  * France et al. (2016) France, K., Loyd, R. O. P., Youngblood, A., et al. 2016, ApJ, 820, 89, doi: 10.3847/0004-637X/820/2/89
  * Gueymard (2018) Gueymard, C. A. 2018, Solar Energy, 169, 434 , doi: https://doi.org/10.1016/j.solener.2018.04.067
  * Hammond et al. (2020) Hammond, M., Tsai, S.-M., & Pierrehumbert, R. T. 2020, ApJ, 901, 78, doi: 10.3847/1538-4357/abb08b
  * Hauglustaine et al. (1994) Hauglustaine, D. A., Granier, C., Brasseur, G. P., & MéGie, G. 1994, J. Geophys. Res., 99, 1173, doi: 10.1029/93JD02987
  * House et al. (2003) House, C. H., Runnegar, B., & Fitz-Gibbon, S. T. 2003, Geobiology, 1, 15, doi: 10.1046/j.1472-4669.2003.00004.x
  * Howard et al. (2012) Howard, A. W., Marcy, G. W., Bryson, S. T., et al. 2012, ApJS, 201, 15, doi: 10.1088/0067-0049/201/2/15
  * Hu et al. (2021) Hu, R., Damiano, M., Scheucher, M., et al. 2021, 921, L8, doi: 10.3847/2041-8213/ac1f92
  * Hu et al. (2012) Hu, R., Seager, S., & Bains, W. 2012, Astrophys. J., 761, doi: 10.1088/0004-637X/761/2/166
  * Innes & Pierrehumbert (2022) Innes, H., & Pierrehumbert, R. T. 2022, ApJ, 927, 38, doi: 10.3847/1538-4357/ac4887
  * Innes et al. (2023) Innes, H., Tsai, S.-M., & Pierrehumbert, R. T. 2023, ApJ, 953, 168, doi: 10.3847/1538-4357/ace346
  * Johnson et al. (2004) Johnson, T. J., Sams, R. L., & Sharpe, S. W. 2004, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5269, Chemical and Biological Point Sensors for Homeland Defense, ed. I. Sedlacek, Arthur J., R. Colton, & T. Vo-Dinh, 159–167, doi: 10.1117/12.515604
  * Judeikis & G. Wren (1977) Judeikis, H. S., & G. Wren, A. 1977, Atmospheric Environment (1967), 11, 1221, doi: https://doi.org/10.1016/0004-6981(77)90099-3
  * Kettle et al. (2001) Kettle, A. J., Rhee, T. S., von Hobe, M., et al. 2001, J. Geophys. Res., 106, 12,193, doi: 10.1029/2000JD900630
  * Kite & Ford (2018) Kite, E. S., & Ford, E. B. 2018, ApJ, 864, 75, doi: 10.3847/1538-4357/aad6e0
  * Klippenstein (2023) Klippenstein, S. J. 2023, Private Communication
  * Leconte et al. (2024) Leconte, J., Spiga, A., Clément, N., et al. 2024, arXiv e-prints, arXiv:2401.06608, doi: 10.48550/arXiv.2401.06608
  * Lee et al. (2021) Lee, E. K. H., Parmentier, V., Hammond, M., et al. 2021, MNRAS, 506, 2695.  https://arxiv.org/abs/2106.11664
  * Lee et al. (2020) Lee, G. K., Casewell, S. L., Chubb, K. L., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 4674, doi: 10.1093/mnras/staa1882
  * Leung et al. (2022) Leung, M., Schwieterman, E. W., Parenteau, M. N., & Fauchez, T. J. 2022, ApJ, 938, 6, doi: 10.3847/1538-4357/ac8799
  * Li et al. (2023) Li, C.-Y., Cao, H.-Y., Wang, Q., et al. 2023, The ISME Journal, 17, 1184, doi: 10.1038/s41396-023-01430-z
  * Lovelock (1989) Lovelock, J. E. 1989, Reviews of Geophysics, 27, 215, doi: 10.1029/RG027i002p00215
  * Loyd et al. (2016) Loyd, R. O. P., France, K., Youngblood, A., et al. 2016, ApJ, 824, 102, doi: 10.3847/0004-637X/824/2/102
  * Luque & Pallé (2022) Luque, R., & Pallé, E. 2022, Science, 377, 1211, doi: 10.1126/science.abl7164
  * Madhusudhan et al. (2020) Madhusudhan, N., Nixon, M. C., Welbanks, L., Piette, A. A. A., & Booth, R. A. 2020, ApJ, 891, L7, doi: 10.3847/2041-8213/ab7229
  * Madhusudhan et al. (2021) Madhusudhan, N., Piette, A. A. A., & Constantinou, S. 2021, ApJ, 918, 1, doi: 10.3847/1538-4357/abfd9c
  * Madhusudhan et al. (2023) Madhusudhan, N., Sarkar, S., Constantinou, S., et al. 2023, ApJ, 956, L13, doi: 10.3847/2041-8213/acf577
  * Malik et al. (2019b) Malik, M., Kempton, E. M. R., Koll, D. D. B., et al. 2019b, ApJ, 886, 142, doi: 10.3847/1538-4357/ab4a05
  * Malik et al. (2019a) Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019a, The Astronomical Journal, 157, 170, doi: 10.3847/1538-3881/ab1084
  * Mateos et al. (2023) Mateos, K., Chappell, G., Klos, A., et al. 2023, Science Advances, 9, eade4847, doi: 10.1126/sciadv.ade4847
  * Meadows et al. (2023) Meadows, V., Lincowski, A., & Lustig-Yaeger, J. 2023, in American Astronomical Society Meeting Abstracts, Vol. 55, American Astronomical Society Meeting Abstracts, 125.04
  * Nixon & Madhusudhan (2021) Nixon, M. C., & Madhusudhan, N. 2021, MNRAS, 505, 3414, doi: 10.1093/mnras/stab1500
  * Peacock et al. (2020) Peacock, S., Barman, T., Shkolnik, E. L., et al. 2020, ApJ, 895, 5, doi: 10.3847/1538-4357/ab893a
  * Pierrehumbert (2023) Pierrehumbert, R. T. 2023, ApJ, 944, 20, doi: 10.3847/1538-4357/acafdf
  * Piette & Madhusudhan (2020) Piette, A. A. A., & Madhusudhan, N. 2020, ApJ, 904, 154, doi: 10.3847/1538-4357/abbfb1
  * Pilcher (2003) Pilcher, C. B. 2003, Astrobiology, 3, 471, doi: 10.1089/153110703322610582
  * Rogers et al. (2023) Rogers, J. G., Schlichting, H. E., & Owen, J. E. 2023, ApJ, 947, L19, doi: 10.3847/2041-8213/acc86f
  * Schwieterman et al. (2018) Schwieterman, E. W., Kiang, N. Y., Parenteau, M. N., et al. 2018, Astrobiology, 18, 663, doi: 10.1089/ast.2017.1729
  * Seager et al. (2013) Seager, S., Bains, W., & Hu, R. 2013, ApJ, 777, 95, doi: 10.1088/0004-637X/777/2/95
  * Seinfeld & Pandis (2016) Seinfeld, J. H., & Pandis, S. N. 2016, Atmospheric chemistry and physics: from air pollution to climate change (John Wiley & Sons, Inc.)
  * Smith (1998) Smith, M. D. 1998, Icarus, 132, 176 , doi: https://doi.org/10.1006/icar.1997.5886
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38, doi: 10.3847/1538-4365/aab4f9
  * Tsai et al. (2014) Tsai, S.-M., Dobbs-Dixon, I., & Gu, P.-G. 2014, ApJ, 793, 141, doi: 10.1088/0004-637X/793/2/141
  * Tsai et al. (2021b) Tsai, S.-M., Innes, H., Lichtenberg, T., et al. 2021b, The Astrophysical Journal Letters, 922, L27, doi: 10.3847/2041-8213/ac399a
  * Tsai et al. (2018) Tsai, S.-M., Kitzmann, D., Lyons, J. R., et al. 2018, ApJ, 862, 31, doi: 10.3847/1538-4357/aac834
  * Tsai et al. (2017) Tsai, S.-M., Lyons, J. R., Grosheintz, L., et al. 2017, Astrophys. J. Suppl. Ser., 228, 1, doi: 10.3847/1538-4365/228/2/20
  * Tsai et al. (2021) Tsai, S.-M., Malik, M., Kitzmann, D., et al. 2021, The Astrophysical Journal, 923, 264, doi: 10.3847/1538-4357/ac29bc
  * Tsai et al. (2023c) Tsai, S.-M., Moses, J. I., Powell, D., & Lee, E. K. H. 2023c, ApJ, 959, L30, doi: 10.3847/2041-8213/ad1405
  * Tsai et al. (2024) Tsai, S.-M., Parmentier, V., Mendonça, J. M., et al. 2024, arXiv e-prints, arXiv:2310.17751, doi: 10.48550/arXiv.2310.17751
  * Tsang (1987) Tsang, W. 1987, Journal of Physical and Chemical Reference Data, 16, 471, doi: 10.1063/1.555802
  * Tsang & Hampson (1986) Tsang, W., & Hampson, R. F. 1986, Journal of Physical and Chemical Reference Data, 15, 1087, doi: 10.1063/1.555759
  * Tsiaras et al. (2019) Tsiaras, A., Waldmann, I. P., Tinetti, G., Tennyson, J., & Yurchenko, S. N. 2019, Nature Astronomy, 3, 1086, doi: 10.1038/s41550-019-0878-9
  * Tsuboi et al. (1981) Tsuboi, T., Katoh, M., Kikuchi, S., & Hashimoto, K. 1981, Japanese Journal of Applied Physics, 20, 985, doi: 10.1143/JJAP.20.985
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Venturini et al. (2020) Venturini, J., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, A&A, 643, L1, doi: 10.1051/0004-6361/202039141
  * Wogan et al. (2024) Wogan, N. F., Batalha, N. E., Zahnle, K. J., et al. 2024, ApJ, 963, L7, doi: 10.3847/2041-8213/ad2616
  * Xu et al. (2015) Xu, Z. F., Raghunath, P., & Lin, M. C. 2015, Journal of Physical Chemistry A, 119, 7404, doi: 10.1021/acs.jpca.5b00553
  * Youngblood et al. (2016) Youngblood, A., France, K., Loyd, R. O. P., et al. 2016, ApJ, 824, 101, doi: 10.3847/0004-637X/824/2/101
  * Yu et al. (2021) Yu, X., Moses, J. I., Fortney, J. J., & Zhang, X. 2021, ApJ, 914, 38, doi: 10.3847/1538-4357/abfdc7
  * Zeng et al. (2019) Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proceedings of the National Academy of Science, 116, 9723, doi: 10.1073/pnas.1812905116
  * Zhang et al. (2019) Zhang, M., Chachan, Y., Kempton, E. M. R., & Knutson, H. A. 2019, PASP, 131, 034501, doi: 10.1088/1538-3873/aaf5ad
  * Zhang et al. (2020) Zhang, M., Chachan, Y., Kempton, E. M. R., Knutson, H. A., & Chang, W. H. 2020, ApJ, 899, 27, doi: 10.3847/1538-4357/aba1e6
  * Zhang et al. (2020) Zhang, M., Park, K.-T., Yan, J., et al. 2020, Progress in Oceanography, 186, 102392, doi: https://doi.org/10.1016/j.pocean.2020.102392


# Resource Sharing in Energy Communities: A Cooperative Game Approach
††thanks: This work was supported in part by the National Science Foundation
under Award 2218110 and the Power Systems and Engineering Research Center
(PSERC) Research Project M-46.

Ahmed S. Alahmed and Lang Tong, ({asa278<EMAIL_ADDRESS>
0000-0002-4715-4379 0000-0003-3322-2681 School of Electrical and Computer
Engineering, Cornell University, Ithaca, USA

###### Abstract

We analyze the overall benefits of an energy community cooperative game under
which distributed energy resources (DER) are shared behind a regulated
distribution utility meter under a general net energy metering (NEM) tariff.
Two community DER scheduling algorithms are examined. The first is a community
with centrally controlled DER, whereas the second is decentralized letting its
members schedule their own DER locally. For both communities, we prove that
the cooperative game’s value function is superadditive, hence the grand
coalition achieves the highest welfare. We also prove the balancedness of the
cooperative game under the two DER scheduling algorithms, which means that
there is a welfare re-distribution scheme that de-incentivizes players from
leaving the grand coalition to form smaller ones. Lastly, we present five ex-
post and an ex-ante welfare re-distribution mechanisms and evaluate them in
simulation, in addition to investigating the performance of various community
sizes under the two DER scheduling algorithms.

###### Index Terms:

Cooperative games, DER aggregation, energy communities, net metering, payoff
allocations.

##  I Introduction

The proliferating penetration levels of behind-the-meter (BTM) DER under NEM
surfaced multiple issues including distribution system operator (DSO) revenue
shortfalls, giving rise to successor NEM tariffs that reward DSO customers at
lower rates for energy exports to the grid [1]. In this work, a cooperative
game-theoretic formulation is developed to capture the interaction of energy
prosumers in an energy community coalition and whether forming the coalition
can provide its members benefits higher than those achieved under the DSO’s
successor NEM tariffs. Another goal is to evaluate several ex-post allocation
mechanisms and compare them to the ex-ante Dynamic-NEM (D-NEM) mechanism in
[2] that uses pricing signals to induce community members’ to achieve the
community’s maximum welfare.

Cooperative game theory has been extensively used in literature to model the
interaction of players and welfare/cost re-distribution in energy communities,
as it allows for the analysis of coalition formation, and the allocation of
benefits in a fair and efficient manner [3, 4, 5, 6, 7, 8, 9]. Various ex-
post111By ex-post, we refer to the re-distribution of coalition’s realized
costs/welfare after its resources have been scheduled and exercised. solution
concepts have been proposed with a focus on the fairness and computational
complexity of cost allocation, and whether the allocation is stabilizing,
i.e., no player or a subset of players would find it better off to leave the
grand coalition and form smaller ones. A stabilizing nucleolus-based
allocation was proposed in [3, 5, 8]. Shapley and Nash bargaining-based
solutions were examined in [6, 8, 4, 7], and a proportional rule based cost
allocation was considered in [3, 4]. A mechanism that explicitly incorporates
whether players are causing or mitigating cost to the coalition in allocating
costs was proposed in [9].

Whereas some of the proposed ex-post allocation schemes were shown to belong
to the cooperative game’s core, including the allocation in [9], which also
satisfies the cost causation principle and has a low computation cost, the
scheduling of the community DER is decoupled from the welfare/cost re-
distribution problem. More specifically, the aforementioned well-established
allocation schemes (e.g., Shapley, nucleolus, proportional, etc.) distribute
the coalition value ex-post, which necessitates the existence of a central
planner who schedules the coalition DER to maximize sharing benefits. This is
unlike the D-NEM mechanism, proposed in [2], which employs an ex-ante pricing
rule that induces players to achieve the maximum social welfare and cost
causation principle, while stabilizing the coalition.

To this end, the contributions of this work involve proposing two local energy
coalitions among prosumers, where flexible demands and renewable distributed
generation (DG) are optimized through two different DER scheduling algorithms:

  1. 1.

Centralized scheduling: A model that takes full control of players’ DER to
maximize the community’s welfare.

  2. 2.

Decentralized scheduling: A model that adopts types of the DSO NEM tariff and
lets players schedule their own DER.

For both energy scheduling algorithms, we formulate a cooperative game and
show that the value function, characterized by the coalition welfare, is
superadditive and the game is balanced. For both scheduling algorithms, we
evaluate five different ex-post payoff allocation mechanisms to distribute the
coalition welfare and show their performance compared to the ex-ante D-NEM
mechanism in [2]. We show that, unlike D-NEM, the ex-post allocations do not
stabilize the centrally scheduled community, and require a central planner to
maximize the community’s social welfare.

##  II Problem Formulation

We consider a set of $H$ electricity customers indexed by $i\in{\cal
H}:=\\{1,2,\ldots,H\\}$, who may have flexible loads, BTM renewable DG, and
central renewable DG. The customers face a regulated DSO that implements a NEM
tariff [1]. This work analyzes the problem of energy coalitions in the form of
energy communities, under which the members collaboratively share their DER
behind the DSO’s revenue meter (Fig.1). As depicted in Fig.1, we analyze two
DER scheduling algorithms, where under the first a central controller
schedules all DER by controlling the local home energy management systems
(HEMS), whereas in the second, the central controller is absent, and each
prosumer is independently scheduled by its local HEMS in a decentralized
manner. Clearly, the privacy of members is preserved only under the latter
control scheme.

In the sequel, we mathematically model customers’ DER, payment and surplus
function, followed by formulating the scheduling algorithms under both
scheduling algorithms.

Figure 1: Energy community coalition framework. Player’s flexible consumption
and renewables are $d_{i},r_{i}\in\mathbb{R}_{+}$, respectively, and their net
consumption is $z_{i}\in\mathbb{R}$. The arrows point to the positive
direction of energy flow.

###  II-A Member Resources, Payment, and Surplus

####  II-A1 Flexible Consumption, Renewable DG, and Net Consumption

Every customer $i\in{\cal H}$ is assumed to have $K$ flexible loads indexed by
$k\in\mathcal{K}=\\{1,2,\ldots,K\\}$. The energy consumption bundle of the
devices is denoted by

$\bm{d}_{i}=\left(d_{i1},\ldots,d_{iK}\right)\in\mathcal{D}_{i}:=\\{\bm{d}_{i}:\underline{\bm{d}}_{i}\preceq\bm{d}_{i}\preceq\overline{\bm{d}}_{i}\\}\subseteq\mathbb{R}_{+}^{K},$

where $\underline{\bm{d}}_{i}$ and $\overline{\bm{d}}_{i}$ are the consumption
bundle’s lower and upper limits of customer $i$, respectively.

The random renewable DG of every customer $i\in{\cal H}$, is denoted by
$r_{i}\in\mathbb{R}_{+}$. Customers may own shares of the central generation,
which is already accounted for in $r_{i}$.

The net consumption $z_{i}\in\mathbb{R}$ of every prosumer $i\in{\cal H}$ is
defined as $z_{i}:=\bm{1}^{\top}\bm{d}_{i}-r_{i}$, where $z_{i}\geq 0$ and
$z_{i}<0$ represent a net-consuming and net-producing customers, respectively.
To ensure network integrity regardless of customers’ net-consumption, the DSO
imposes dynamic operating envelopes $\underline{z}_{i}\leq
z_{i}\leq\overline{z}_{i}$, where $\underline{z}_{i}\leq 0$ and
$\overline{z}_{i}\geq 0$ are the export and import envelopes, respectively
[10]. For every $i\in{\cal H}$, we assume
$\bar{z}_{i}\geq\mathbf{1}^{\top}\underline{\bm{d}}_{i}-r_{i}$, and
$\underline{z}_{i}\leq\mathbf{1}^{\top}\overline{\bm{d}}_{i}-r_{i}$ to avoid
infeasibility of the prosumer’s decision.

####  II-A2 Customer Payment and Surplus

Every $i\in{\cal H}$ customer is billed under the utility’s NEM X regime [11]
based on its net consumption and the NEM X parameter $\pi=(\pi^{+},\pi^{-})$
as

$P^{\pi}(z_{i})=\pi^{+}[z_{i}]^{+}-\pi^{-}[z_{i}]^{-},$ (1)

where $\pi^{+},\pi^{-}\in\mathbb{R}_{+}$ are the retail rate, and export rate,
respectively, and $[x]^{+}:=\max\\{0,x\\},[x]^{-}:=\max\\{0,-x\\}$ are the
positive and negative parts of $x\in\mathbb{R}$. We assume
$\pi^{-}\leq\pi^{+}$ as the majority of NEM successor policies consider to
better align the value of excess generation with DSO deferred costs, and avoid
risk-free price arbitrage [1, 2].

The surplus $S_{i}^{\pi}:\mathbb{R}^{K}\rightarrow\mathbb{R}$ of every
$i\in{\cal H}$ is
($S_{i}^{\pi}(\bm{d}_{i},z_{i}):=U_{i}(\bm{d}_{i})-P^{\pi}(z_{i})$), where
$U_{i}(\bm{d}_{i})$ is the utility function of consuming $\bm{d}_{i}$. For all
$i\in{\cal H}$, $U_{i}(\cdot)$ is assumed to be additive (i.e.,
$U_{i}(\bm{d}_{i})=\sum_{k\in\mathcal{K}}U_{ik}\left(d_{ik}\right)$), strictly
concave, non-decreasing, and continuously differentiable with
$U_{ik}(\bm{0})=0,\forall k\in\mathcal{K}$.

###  II-B Community Energy Management Models

Let ${\cal N}\subseteq{\cal H}$ represent an energy coalition of customers,
with cardinality $|{\cal N}|$, who pool their DER to improve their overall
benefits and jointly share the coalition welfare (Fig.1). Here, we present the
two energy community control structures.

####  II-B1 Community with Centralized Scheduling

Under this scheme, the operator of coalition ${\cal N}\subseteq{\cal H}$
maximizes community welfare by fully controlling customers’ resources. Thus,
the central controller solves a welfare maximization program

$\displaystyle\underset{\\{\bm{d}_{i}\\}_{i=1}^{|{\cal
N}|},\\{z_{i}\\}_{i=1}^{|{\cal N}|}}{\text{ maximize}}$
$\displaystyle~{}~{}W({\cal N}):=\sum_{i\in{\cal
N}}S_{i}^{\pi}(\bm{d}_{i},z_{i})$ (2a) subject to
$\displaystyle\sum_{i\in{\cal N}}P^{\pi}(z_{i})=P^{\pi}(\sum_{i\in{\cal
N}}z_{i})$ (2b) $\displaystyle\text{for all }~{}i=1,\ldots,N$ $\displaystyle
z_{i}=\bm{1}^{\top}\bm{d}_{i}-r_{i}$ (2c) $\displaystyle\underline{z}_{i}\leq
z_{i}\leq\overline{z}_{i}$ (2d)
$\displaystyle\underline{\bm{d}}_{i}\preceq\bm{d}_{i}\preceq\overline{\bm{d}}_{i},$
(2e)

where (2b) is the operator’s profit-neutrality constraint. Denote by
$\bm{d}_{i}^{{\cal N}^{\ast}}$ and $z_{i}^{{\cal N}^{\ast}}$ the optimal
consumption, and net-consumption, respectively, $\forall i\in{\cal N}$, that
solve (2). The optimal welfare of (2) is therefore $W^{\ast}({\cal
N}):=\sum_{i\in{\cal N}}S_{i}^{\pi}(\bm{d}^{{\cal N}^{\ast}}_{i},z^{{\cal
N}^{\ast}}_{i})$.

An ex-post payoff allocation is needed to distribute the coalition welfare, as
will be described in Sec.III.

####  II-B2 Community with Decentralized Scheduling

Under this scheme, the operator of coalition ${\cal N}\subseteq{\cal H}$
adopts the DSO’s NEM X and lets its members schedule their own resources as if
they were facing the DSO, i.e., every $i\in{\cal N}$ solves

$\displaystyle\underset{\bm{d}_{i},z_{i}}{\text{ maximize}}$
$\displaystyle~{}~{}W(i):=S_{i}^{\pi}(\bm{d}_{i},z_{i})$ (3a) subject to
$\displaystyle(\ref{eq:PpiNetCons})-(\ref{eq:Ppidi}).$ (3b)

For all $i\in{\cal N}$, let $\bm{d}_{i}^{\dagger}$ and $z_{i}^{\dagger}$
denote, respectively, the optimal consumption, and net-consumption that solve
(3) and result in the customer maximum welfare
$W^{\dagger}(i):=S_{i}^{\pi}(\bm{d}_{i}^{\dagger},z^{\dagger}_{i})$.
Evidently, the total welfare of the decentralized community is suboptimal.
Nonetheless, the ex-ante D-NEM in [2] achieves the maximum welfare, while
adopting a decentralized scheduling scheme. D-NEM is briefly described in
Sec.III-D.

##  III Energy Community Cooperative Game

Here, we introduce the cooperative energy community game and its properties
under the two DER scheduling algorithms.

###  III-A Coalitions and Value Functions

An energy community coalition is formed when any subset of prosumers ${\cal
N}$ decide to cooperate by sharing their DER, and jointly face the DSO. The
set of all possible coalitions is defined as the power set of ${\cal
H},2^{\cal H}:=\\{{\cal N}:{\cal N}\subseteq{\cal H}\\}$. ${\cal H}$ and
$\emptyset$ are the grand coalition and empty coalition, respectively.

An $|{\cal H}|$-player cooperative game is defined by the pair (${\cal
H},\nu$), where $\nu:2^{\cal H}\rightarrow\mathbb{R}$ is a characteristic
function that assigns value $\nu({\cal N})$ to each coalition ${\cal
N}\subseteq{\cal H}$. The value function $\nu({\cal N})$ is characterized by
the welfare of the energy community coalition ${\cal N}$

$\nu({\cal N}):=W({\cal N}),$ (4)

with $\nu({\cal N})=0$ if ${\cal N}=\emptyset$, and $\nu({\cal H})$ being the
value of the grand coalition. Note that the value of a singelton coalition
$\nu(i)=W(i)$ is simply the welfare of a standalone customer who autonomously
faces the DSO’s NEM X.

Next, we formulate the value function under both the centralized and
decentralized energy management schemes.

####  III-A1 Community with Centralized Scheduling

Here, the $|{\cal N}|$ players’ DER are scheduled by the community manager,
which results in the maximum social welfare in (2). From the definition in
(4), the cooperative game value function under this scheme is
$\nu^{\ast}({\cal N})=W^{\ast}({\cal N})$.

It turns out that the value function under the centralized scheme is
superadditive, which means that the grand coalition achieves the highest
welfare.

######  Lemma 1.

The cooperative game under centralized DER control $({\cal H},\nu^{\ast})$ is
superadditive, i.e., $\nu^{\ast}({\cal N})+\nu^{\ast}({\cal
T})\leq\nu^{\ast}({\cal N}\cup{\cal T}),~{}~{}\forall{\cal N},{\cal
T}\subseteq{\cal H},{\cal N}\cap{\cal T}=\emptyset$.

Proof: See Appendix V-A. $\blacksquare$

####  III-A2 Community with Decentralized Scheduling

Here, customers continue scheduling their DER, as if they were facing the DSO,
therefore, their coalition welfare is $\sum_{i\in{\cal
N}}S_{i}^{\pi}(\bm{d}^{\dagger}_{i},z^{\dagger}_{i})$. Because we enforce a
profit-neutral operator, we adjust this welfare by adding a term that ensures
the re-distribution of any excess profit:

$\displaystyle W^{\dagger}({\cal N}):=\sum_{i\in{\cal
N}}W^{\dagger}(i)+\big{(}\sum_{i\in{\cal
N}}P^{\pi}(z_{i}^{\dagger})-P^{\pi}(z_{{\cal N}}^{\dagger})\big{)},$ (5)

where $z_{{\cal N}}^{\dagger}:=\sum_{i\in{\cal N}}z_{i}^{\dagger}$. The second
term on the right-hand side of (5) is the difference between the payment the
community operator collects from its members, and the payment it submits to
the DSO. From (4), the cooperative game value function is $\nu^{\dagger}({\cal
N})=W^{\dagger}({\cal N})$.

Interestingly, and as formally stated in Lemma 2, although the players were
neither centrally scheduled, as in (2), nor incentivized by prices to achieve
the maximum welfare, as in [2], the value function of their cooperation is
superadditive.

######  Lemma 2.

The cooperative game under decentralized DER control $({\cal
H},\nu^{\dagger})$ is superadditive, i.e., $\nu^{\dagger}({\cal
N})+\nu^{\dagger}({\cal T})\leq\nu^{\dagger}({\cal N}\cup{\cal
T}),~{}~{}\forall{\cal N},{\cal T}\subseteq{\cal H},~{}{\cal N}\cap{\cal
T}=\emptyset$.

Proof: See Appendix V-B. $\blacksquare$

###  III-B Payoff Allocation, Imputation, and Core

We showed that both the community under both scheduling schemes has a
superadditive value function, which raises the normative question of how we
should assign the coalition value to each player of the cooperative game.

######  Definition 1 (Payoff allocation and imputation).

Let $\bm{\psi}\in\mathbb{R}^{H}$ denote the payoff allocation vector for
coalition ${\cal N}\subseteq{\cal H}$, with $\psi_{i}$ being the payoff to
prosumer $i\in{\cal H}$. No payoff is assigned to prosumers outside the
coalition, i.e., $\psi_{i}=0,~{}\forall i\notin{\cal N}$.

The allocation $\bm{\psi}$ is an imputation if: 1) the grand coalition welfare
is wholly distributed to the members (efficiency), i.e., $\sum_{i\in{\cal
H}}\psi_{i}=\nu({\cal H})$, and 2) individuals are better off joining the
coalition (individual rationality), i.e., $\psi_{i}\geq\nu(i),~{}\forall
i\in{\cal H}$. Therefore, the set of all imputations is defined as

${\cal I}:=\big{\\{}\bm{\psi}\in\mathbb{R}^{H}:\sum_{i\in{\cal
H}}\psi_{i}=\nu({\cal H}),~{}\psi_{i}\geq\nu(i),\forall i\in{\cal
H}\big{\\}}.$ (6)

It is not hard to show that for any superadditive game, ${\cal I}$ is
nonempty. Despite the favorable properties of a payoff allocation that is an
imputation $\bm{\psi}\in{\cal I}$, it does not guarantee group rationality,
which ensures that no subset of prosumers would be better-off if they jointly
withdraw from the grand coalition to form their own coalition. Therefore, we
further require the payoff allocation to be in the core, which is the
fundamental solution concept for cooperative games.

######  Definition 2 (Core and coalition stability).

The core ${\cal C}$ for the cooperative game $({\cal H},\nu)$ is defined as
the set of all payoffs that are imputations and ensures that no coalition can
have a value greater than the sum of its members’ payoffs

${\cal C}:=\big{\\{}\bm{\psi}\in{\cal I}:\sum_{i\in{\cal
N}}\psi_{i}\geq\nu({\cal N}),~{}\forall{\cal N}\in 2^{\cal H}\big{\\}}.$ (7)

The coalition is stable if $\bm{\psi}\in{\cal C}$.

###  III-C Balanced Games

Balanced games, which gives a necessary and sufficient condition for the
existence of a stabilizing allocation (nonempty core) in a cooperative game
are defined next.

######  Definition 3 (Balanced game).

A map $\alpha:2^{{\cal H}}\rightarrow[0,1]$ is called a balanced map if
$\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal N})\bm{1}_{{\cal N}}(i)=1$, where
$\bm{1}_{{\cal N}}(i)=1$ if $i\in{\cal N}$ and $\bm{1}_{{\cal N}}(i)=0$ if
$i\notin{\cal N}$. The cooperative game (${\cal H},\nu$) is balanced if for
any balanced map $\alpha$, $\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal
N})\nu({\cal N})\leq\nu({\cal H})$.

The following Theorems show that the cooperative games under both DER
scheduling algorithms are balanced.

######  Theorem 1 (Balancedness under central DER scheduling).

The energy community cooperative game (${\cal H},\nu^{\ast}$) is balanced.

Proof: See Appendix V-D. $\blacksquare$

######  Theorem 2 (Balancedness under decentral DER scheduling).

The energy community cooperative game (${\cal H},\nu^{\dagger}$) is balanced.

Proof: See Appendix V-E. $\blacksquare$
From the Bondareva-Shapley Theorem [12], the balancedness of the cooperative
game under the decentralized and centralized control schemes is a necessary
and sufficient condition for the non-emptiness of their core. The
superadditivity and balancedness of the cooperative game imply the existence
of a stable payoff allocation of the coalitional welfare that is satisfactory
for every member of the grand coalition.

###  III-D Payoff Allocation Mechanisms

We present five ex-post payoff allocations and compare them to D-NEM in [2,
10], which is an ex-ante allocation mechanism. For all of the ex-post
allocation schemes, we ensure that the utility of consumption $U_{i}(\cdot)$
is wholly assigned to each member $i$, as such value is not transferable.
Hence, only the coalition payment $P^{\pi}(z_{\cal N})$ is re-distributed.

####  III-D1 Equal Division Allocation

The equal division allocation equally splits the payment to the players [13].
Therefore, the allocation to each player under this scheme is

$\psi^{(1)}_{i}=U_{i}-P^{\pi}(z_{\cal N})/|{\cal N}|,\quad\forall i\in{\cal
N}.$

####  III-D2 Egalitarian Allocation

This allocation assigns to player $i$ all of his individual payment plus an
equal share of any possible cost savings achieved by forming the coalition
[13]

$\psi^{(2)}_{i}=U_{i}-\left(P^{\pi}(z_{i})+\frac{P^{\pi}(z_{\cal
N})-\sum_{j=1}^{N}P^{\pi}(z_{j})}{|{\cal N}|}\right),\forall i\in{\cal N}.$

####  III-D3 Proportional Allocation

The proportional allocation assigns payments to each member that is
proportionate to the cost they impose on the coalition

$\psi_{i}^{(3)}=U_{i}-\frac{W(i)}{W({\cal N})}P^{\pi}(z_{\cal N}),\quad\forall
i\in{\cal N},$

####  III-D4 Allocation in [9]

The ex-post allocation in [9] assigns payment to players based on their net
consumption and the aggregate net consumption of the community as

$\psi_{i}^{(4)}=U_{i}-\begin{cases}\pi^{+}z_{i}&\text{ if }z_{{\cal N}}\geq
0\\\ \pi^{-}z_{i}&\text{ if }z_{{\cal N}}<0\end{cases},\quad\forall i\in{\cal
N}.$

####  III-D5 Shapley Value Allocation

Shapley value distributes payoff among the players based on their marginal
contributions to the coalition [12] as, for every $i\in{\cal N}$,

$\displaystyle\psi_{i}^{(5)}=U_{i}-\sum_{{\cal N}\in 2^{\cal H},i\in{\cal
N}}\frac{(|{\cal H}|-|{\cal N}|)!(|{\cal N}|-1)!}{|{\cal H}|!}(P^{\pi}(z_{\cal
N})-P^{\pi}(z_{{\cal N}\setminus i})).$

####  III-D6 Dynamic NEM [2, 10]

Unlike the allocation mechanisms in Sec. III-D1–III-D5, which are ex-post and
reliant on a centralized energy management scheme for social welfare
maximization, the welfare-maximizing D-NEM employs an ex-ante and threshold-
based pricing mechanism $\pi^{\mbox{\tiny DNEM}}$ that sets the community
price based on the available renewables, as

$\pi^{\mbox{\tiny DNEM}}(\bm{r})=\begin{cases}\pi^{+}&\text{ if
}\sum_{i\in{\cal N}}r_{i}<f_{\cal N}(\pi^{+})\\\ \pi^{z}(\bm{r})&\text{ if
}\sum_{i\in{\cal N}}r_{i}\in[f_{\cal N}(\pi^{+}),f_{\cal N}(\pi^{-})]\\\
\pi^{-}&\text{ if }\sum_{i\in{\cal N}}r_{i}<f_{\cal N}(\pi^{-}).\end{cases}$

The computation of the thresholds $(f_{\cal N}(\pi^{+}),f_{\cal N}(\pi^{-}))$
and the price $\pi^{z}(\bm{r})\in[\pi^{-},\pi^{+}]$ are described in [2] and
[10] without and with operating envelopes, respectively. The work in [2, 10]
proved that the optimal response of members under D-NEM is equivalent to the
optimal response under the centralized control scheme, i.e.,
$(\bm{d}_{i}^{{\cal N}\ast},z_{i}^{{\cal N}\ast}),\forall i\in{\cal N}$.
Therefore, the member payoff (surplus) under D-NEM,

$\psi_{i}^{(6)}=U_{i}(\bm{d}_{i}^{{\cal N}\ast})-\pi^{\mbox{\tiny
DNEM}}(\bm{r})\cdot z_{i}^{{\cal N}\ast},$

is attained rather than distributed ex-post, as members’ react to the D-NEM
price $\pi^{\mbox{\tiny DNEM}}(\bm{r})$.

##  IV Numerical Results

We evaluated the performance of community scheduling algorithms and the ex-
post and ex-ante allocation mechanisms by constructing a hypothetical
residential energy community coalition of various sizes using a one-year
consumption and DER data of 20 households.222Using PecanStreet data for
households in Austin, TX in 2018. The community faces a NEM policy with an
hourly billing period, a ToU-based retail rate with $\pi^{+}_{\mbox{\tiny
ON}}=\$0.40$/kWh and $\pi^{+}_{\mbox{\tiny OFF}}=\$0.20$/kWh as on- and off-
peak prices, respectively, and an export rate $\pi^{-}$ based on the wholesale
market price.333Using the averaged 2018 real-time wholesale prices in Texas
(ERCOT). The DSO operating envelopes were assumed to be homogeneous with
$\overline{z}_{i}=-\underline{z}_{i}=6kW,\forall i\in\mathcal{N}$.

For every $i\in\mathcal{N}$, the household’s consumption preferences are
characterized by a quadratic concave and non-decreasing utility function of
the form, for every device $k\in\mathcal{K}$,

$U_{ik}(d_{ik})=\left\\{\begin{array}[]{ll}\alpha_{ik}d_{ik}-\frac{1}{2}\beta_{ik}d_{ik}^{2},&0\leq
d_{ik}\leq\frac{\alpha_{ik}}{\beta_{ik}}\\\
\frac{\alpha_{ik}^{2}}{2\beta_{ik}},&d_{ik}>\frac{\alpha_{ik}}{\beta_{ik}},\end{array}\right.$
(8)

where $\alpha_{ik},\beta_{ik}$ are parameters that are learned and calibrated
using historical retail prices,444We used Data.AustinTexas.gov historical
residential rates in Austin, TX. consumption,555We used pre-2018 PecanStreet
data for households in Austin, TX. and elasticities for each load type (see
appendix D in [1]).

Figure 2 shows the average daily welfare gain (left) and normalized difference
(right) of the centralized (2) and decentralized (3) scheduling algorithms
over the aggregate maximum welfare of standalone DSO customers, as the
coalition size is increased by randomly adding households from the 20
households data. For every coalition size, 1000 Monte Carlo runs are
implemented over the randomly selected households.

Two main observations are in order. First, as the coalition size increased,
under both scheduling algorithms, the welfare difference and gain increased,
which is due to the superadditivity of the value function under both schemes
(Lemmas 1–2). The gain, however, saturated as coalition size increased, which
was at 2.45% and 1.55% for centralized and decentralized scheduling,
respectively. Also, as expected, the welfare difference, hence gain, of the
centralized scheduling dominates decentralized scheduling. Second, we observe
that although DER was decentrally scheduled under D-NEM, the community
achieved welfare gains equivalent to centralized scheduling. Moreover, re-
distributing the welfare of larger coalitions in Fig.2 is more computationally
intensive under well-established ex-post allocations such as Shapley and the
nucleolus; an issue that is mitigated under D-NEM.

Figure 2: Average daily welfare normalized difference and gain (%).

Table I delves into the payoff allocation of the five presented schemes in
Sec.III-D and D-NEM under both centrally and decentrally scheduled coalitions,
and presents how often the individual rationality axiom was violated, i.e.,
number of hours whereby individual rationality fails $x$ over number of total
hours, hence $x/(24\times 365\times|{\cal N}|)$. Under decentralized
scheduling, only equal division and proportional allocations failed in
achieving individual rationality, which necessarily mean that they do not
belong to the core of the game.

Under centralized scheduling, however, all of the ex-post allocation schemes
failed in achieving individual rationality, meaning that they are not
stabilizing allocations. Only the ex-ante D-NEM achieved individual
rationality. We also observe that the frequency of violating individual
rationality was consistently higher under the centralized scheduling over the
decentralized one, which is because the community payment is subadditive only
under the decentralized scheduling. Lastly, Table I shows that the frequency
of violating individual rationality by the ex-post allocations increased when
$|{\cal N}|$ increased.

TABLE I: Percentage of time of violating individual rationality (%).

$|{\cal N}|$ | DER Schedule | Eq. division | Egalitarian | Proportional | [9] | Shapley | D-NEM
---|---|---|---|---|---|---|---
4 | Centralized | 55.3 | 5.0 | 27.8 | 2.5 | 5.5 | 0
Decentralized | 55.1 | 0 | 27.6 | 0 | 0 | 0
10 | Centralized | 61.1 | 5.1 | 30.1 | 3.0 | 5.6 | 0
Decentralized | 60.4 | 0 | 29.8 | 0 | 0 | 0

The daily welfare gain of each player in a coalition of size $|{\cal N}|=6$
under both DER management algorithms and under the ex-post and ex-ante
allocation mechanisms is shown in Fig.3. One can see that the egalitarian,
[9], Shapley, and D-NEM allocations presented similar patterns of welfare
distribution. D-NEM has only six boxplots because the gain is the same under
both DER controls. Under both DER control algorithms, individual rationality
was violated for all 6 customers under the equal division and proportional
allocations. All 6 customers were individually rational under decentrally
controlled coalitions and egalitarian, [9], and Shapley allocations. These
three allocations however fail to achieve individual rationality
(particularly, customers 4–6) under the centralized DER control case, which
maximizes the community welfare. The ex-ante D-NEM mechanism is the only
policy that achieved individual rationality when the community welfare was
maximized.

Figure 3: Daily welfare gain (%) of players ($|{\cal N}|=6$) under
decentralized† and centralized∗ DER controls. Top panels (left to right):
equal division, egalitarian, proportional. Bottom panels (left to right): [9],
Shapley, D-NEM.

##  V Conclusion

In this paper, a cooperative game theoretic approach is employed to analyze
resource sharing and welfare distribution in energy community coalitions. We
showed that when the community DER is centrally controlled to maximize the
community welfare, and when the DER of each player are locally controlled to
maximize each player’s surplus, the cooperation value function is
superadditive and the game is balanced. The payoff of the game is then
allocated using five different ex-post payoff allocation mechanisms, which are
compared to an ex-ante and welfare maximizing mechanism. In simulations, we
evaluated the welfare gains of the energy coalition under both scheduling
algorithms and the performance of the ex-post and ex-ante allocation
mechanisms.

## References

  * [1] A. S. Alahmed and L. Tong, “Integrating distributed energy resources: Optimal prosumer decisions and impacts of net metering tariffs,” _SIGENERGY Energy Inform. Rev._ , vol. 2, no. 2, p. 13–31, Aug. 2022. [Online]. Available: https://doi.org/10.1145/3555006.3555008
  * [2] ——, “Dynamic net metering for energy communities,” 2023. [Online]. Available: https://arxiv.org/abs/2306.13677
  * [3] Y. Yang, G. Hu, and C. J. Spanos, “Optimal sharing and fair cost allocation of community energy storage,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 5, pp. 4185–4194, Sep. 2021.
  * [4] P. Chakraborty, E. Baeyens, and P. P. Khargonekar, “Cost causation based allocations of costs for market integration of renewable energy,” _IEEE Transactions on Power Systems_ , vol. 33, no. 1, pp. 70–83, 2018.
  * [5] L. Han, T. Morstyn, and M. McCulloch, “Incentivizing prosumer coalitions with energy management using cooperative game theory,” _IEEE Transactions on Power Systems_ , vol. 34, no. 1, pp. 303–313, 2019.
  * [6] A. Fleischhacker, C. Corinaldesi, G. Lettner, H. Auer, and A. Botterud, “Stabilizing energy communities through energy pricing or PV expansion,” _IEEE Trans. on Smart Grid_ , vol. 13, no. 1, pp. 728–737, 2022.
  * [7] S. Cui, Y.-W. Wang, Y. Shi, and J.-W. Xiao, “Community energy cooperation with the presence of cheating behaviors,” _IEEE Transactions on Smart Grid_ , vol. 12, no. 1, pp. 561–573, 2021.
  * [8] N. Vespermann, T. Hamacher, and J. Kazempour, “Access Economy for Storage in Energy Communities,” _IEEE Transactions on Power Systems_ , vol. 36, no. 3, pp. 2234–2250, 2021.
  * [9] P. Chakraborty, E. Baeyens, P. P. Khargonekar, K. Poolla, and P. Varaiya, “Analysis of solar energy aggregation under various billing mechanisms,” _IEEE Tran. on Smart Grid_ , vol. 10, no. 4, pp. 4175–4187, 2019.
  * [10] A. S. Alahmed, G. Cavraro, A. Bernstein, and L. Tong, “Operating-envelopes-aware decentralized welfare maximization for energy communities,” in _2023 59th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , 2023, pp. 1–8.
  * [11] A. S. Alahmed and L. Tong, “On net energy metering X: Optimal prosumer decisions, social welfare, and cross-subsidies,” _IEEE Transactions on Smart Grid_ , vol. 14, no. 02, pp. 1652–1663, 2023.
  * [12] L. S. Shapley, “On balanced sets and cores,” _Naval Research Logistics Quarterly_ , vol. 14, no. 4, pp. 453–460, 1967. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/nav.3800140404
  * [13] H. P. Young, _Cost Allocation: Methods, Principles, Applications_. Amsterdam: North Holland Publishing Co., March 1985. [Online]. Available: https://pure.iiasa.ac.at/id/eprint/2596/

## Appendix: Proofs

###  V-A Proof of Lemma 1

We need to show that for all ${\cal N},{\cal T}\subseteq{\cal H},~{}{\cal
N}\cap{\cal T}=\emptyset$,

$\displaystyle\nu^{\ast}({\cal N})+\nu^{\ast}({\cal T})$
$\displaystyle\leq\nu^{\ast}({\cal N}\cup{\cal T}).$ (9)

Note that

$\displaystyle\nu^{\ast}({\cal N})$ $\displaystyle=W^{\ast}(\bm{d}^{{\cal
N}^{\ast}},\bm{z}^{{\cal N}^{\ast}})=\sum_{i\in{\cal N}}U_{i}(\bm{d}^{{\cal
N}^{\ast}}_{i})-P^{\pi}(\sum_{i\in{\cal N}}z^{{\cal N}^{\ast}}_{i})$ (10)
$\displaystyle\nu^{\ast}({\cal T})$ $\displaystyle=W^{\ast}(\bm{d}^{{\cal
T}^{\ast}},\bm{z}^{{\cal T}^{\ast}})=\sum_{i\in{\cal T}}U_{i}(\bm{d}^{{\cal
T}^{\ast}}_{i})-P^{\pi}(\sum_{i\in{\cal T}}z^{{\cal T}^{\ast}}_{i}),$ (11)

where for any coalition ${\cal M}\in 2^{{\cal H}},\bm{d}^{{\cal
M}^{\ast}}:=\\{\bm{d}^{{\cal M}^{\ast}}_{i}\\},\bm{z}^{{\cal
M}^{\ast}}:=\\{z^{{\cal M}^{\ast}}_{i}\\},\forall i\in{\cal M}$, and the
second term in (10)-(11) above is from leveraging constraint (2b). If
coalitions ${\cal N}$ and ${\cal T}$ form one coalition ${\cal N}\cup{\cal T}$
and continue scheduling, as before, their new coalition’s value becomes

$\displaystyle\nu({{\cal N}}\cup{{\cal T}})=W(\bm{d}^{{\cal N}^{\ast}\cup{\cal
T}^{\ast}},\bm{z}^{{\cal N}^{\ast}\cup{\cal T}^{\ast}})$
$\displaystyle=\sum_{i\in{\cal N}\cup{\cal T}}U_{i}(\bm{d}^{{\cal
N}^{\ast}\cup{\cal T}^{\ast}}_{i})-P^{\pi}(\sum_{i\in{\cal N}\cup{\cal
T}}z^{{\cal N}^{\ast}\cup{\cal T}^{\ast}}_{i})$ $\displaystyle=\sum_{i\in{\cal
N}}U_{i}(\bm{d}^{{\cal N}^{\ast}}_{i})+\sum_{i\in{\cal T}}U_{i}(\bm{d}^{{\cal
T}^{\ast}}_{i})-P^{\pi}(\sum_{i\in{\cal N}}z^{{\cal
N}^{\ast}}_{i}+\sum_{i\in{\cal T}}z^{{\cal T}^{\ast}}_{i}).$

Note that because

$\displaystyle\nu({{\cal N}}\cup{{\cal T}})-\nu^{\ast}({\cal
N})-\nu^{\ast}({\cal T})$ $\displaystyle=P^{\pi}(\sum_{i\in{\cal N}}z^{{\cal
N}^{\ast}}_{i})+P^{\pi}(\sum_{i\in{\cal T}}z^{{\cal
T}^{\ast}}_{i})-P^{\pi}(\sum_{i\in{\cal N}}z^{{\cal
N}^{\ast}}_{i}+\sum_{i\in{\cal T}}z^{{\cal T}^{\ast}}_{i})$
$\displaystyle\stackrel{{\scriptstyle\text{(Lem.\ref{lem:AdditiveDecentral})}}}{{\geq}}0,$

we have

$\nu({{\cal N}}\cup{{\cal T}})\geq\nu^{\ast}({\cal N})+\nu^{\ast}({\cal T}).$
(12)

Now, by solving (2) for coalition ${\cal N}\cup{\cal T}$ we have

$\displaystyle\nu^{\ast}({\cal N}\cup{\cal T})$
$\displaystyle=W^{\ast}(\bm{d}^{{({\cal N}\cup{\cal
T})}^{\ast}},\bm{z}^{{({\cal N}\cup{\cal T})}^{\ast}})$ $\displaystyle\geq
W(\bm{d}^{{\cal N}^{\ast}\cup{\cal T}^{\ast}},\bm{z}^{{\cal N}^{\ast}\cup{\cal
T}^{\ast}})=\nu({{\cal N}}\cup{{\cal T}})$
$\displaystyle\stackrel{{\scriptstyle\text{(\ref{eq:lem1_3})}}}{{\geq}}\nu^{\ast}({\cal
N})+\nu^{\ast}({\cal T}),$

which completes the proof. $\blacksquare$

###  V-B Proof of Lemma 2

Given that the players’ face the same NEM X tariff design regardless of the
community size, their optimal consumption and net consumption schedules
($\bm{d}^{\dagger}_{i},z^{\dagger}_{i}$) stay the same. Therefore, for any
${\cal N},{\cal T}\subseteq{\cal H}$, we have

$\displaystyle\nu^{\dagger}({\cal N}\cup{\cal T})$
$\displaystyle=W^{\dagger}({\cal N}\cup{\cal T})=\sum_{i\in{\cal N}\cup{\cal
T}}U_{i}(\bm{d}^{\dagger}_{i})-P^{\pi}(z^{\dagger}_{{\cal
N}}+z^{\dagger}_{{\cal T}}).$ $\displaystyle\nu^{\dagger}({\cal N})$
$\displaystyle=W^{\dagger}({\cal N})=\sum_{i\in{\cal
N}}U_{i}(\bm{d}^{\dagger}_{i})-P^{\pi}(z^{\dagger}_{{\cal N}})$
$\displaystyle\nu^{\dagger}({\cal T})$ $\displaystyle=W^{\dagger}({\cal
T})=\sum_{i\in{\cal T}}U_{i}(\bm{d}^{\dagger}_{i})-P^{\pi}(z^{\dagger}_{{\cal
T}}),$

where for any coalition ${\cal M}\in 2^{{\cal H}}$, $z^{\dagger}_{{\cal
M}}:=\sum_{i\in{\cal M}}z_{i}^{\dagger}$. Given that the utility function
terms cancel out, the proof is complete if we show the following

$\displaystyle\nu^{\dagger}({\cal N})+\nu^{\dagger}({\cal T})$
$\displaystyle\leq\nu^{\dagger}({\cal N}\cup{\cal T}),$

which can reformulated to

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal T}})\leq
P^{\pi}(z^{\dagger}_{{\cal N}})+P^{\pi}(z^{\dagger}_{{\cal T}}).$ (13)

where $z^{\dagger}_{{\cal N}\cup{\cal T}}:=z^{\dagger}_{{\cal
N}}+z^{\dagger}_{{\cal T}}$. We prove (13) using the following four exhaustive
cases.

  * •

Case 1 ($z^{\dagger}_{\cal N}\geq 0,z^{\dagger}_{\cal T}\geq 0$): This leads
to $z^{\dagger}_{{\cal N}\cup{\cal T}}\geq 0$. The LHS and RHS in (13) are,
respectively:

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal T}})$
$\displaystyle=\pi^{+}(z^{\dagger}_{\cal N}+z^{\dagger}_{\cal T})$
$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}})+P^{\pi}(z^{\dagger}_{{\cal
T}})$ $\displaystyle=\pi^{+}z^{\dagger}_{\cal N}+\pi^{+}z^{\dagger}_{\cal
T}=\pi^{+}(z^{\dagger}_{\cal N}+z^{\dagger}_{\cal T}),$

hence (13) holds.

  * •

Case 2 ($z^{\dagger}_{\cal N}\leq 0,z^{\dagger}_{\cal T}\leq 0$): This leads
to $z^{\dagger}_{{\cal N}\cup{\cal T}}\leq 0$. The LHS and RHS in (13) are,
respectively:

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal T}})$
$\displaystyle=\pi^{-}(z^{\dagger}_{\cal N}+z^{\dagger}_{\cal T})$
$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}})+P^{\pi}(z^{\dagger}_{{\cal
T}})$ $\displaystyle=\pi^{-}z^{\dagger}_{\cal N}+\pi^{-}z^{\dagger}_{\cal
T}=\pi^{-}(z^{\dagger}_{\cal N}+z^{\dagger}_{\cal T}),$

hence (13) holds.

  * •

Case 3 ($z^{\dagger}_{\cal N}\geq 0,z^{\dagger}_{\cal T}\leq 0$): We have two
sub-cases here

    * –

Case 3.1 ($z^{\dagger}_{{\cal N}\cup{\cal T}}\geq 0$): The LHS and RHS in (13)
are, respectively:

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal T}})$
$\displaystyle=\pi^{+}(z^{\dagger}_{\cal N}+z^{\dagger}_{\cal T})$
$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}})+P^{\pi}(z^{\dagger}_{{\cal
T}})$ $\displaystyle=\pi^{+}z^{\dagger}_{\cal N}+\pi^{-}z^{\dagger}_{\cal T},$

By computing the difference

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal
T}})-P^{\pi}(z^{\dagger}_{{\cal N}})-P^{\pi}(z^{\dagger}_{{\cal
T}})=(\pi^{+}-\pi^{-})z^{\dagger}_{\cal T}\leq 0,$

we verify that (13) holds, as $\pi^{-}\leq\pi^{+}$.

    * –

Case 3.2 ($z^{\dagger}_{{\cal N}\cup{\cal T}}\leq 0$): Following the same
steps in Case 3.1, we get

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal
T}})-P^{\pi}(z^{\dagger}_{{\cal N}})-P^{\pi}(z^{\dagger}_{{\cal
T}})=(\pi^{-}-\pi^{+})z^{\dagger}_{\cal N}\leq 0,$

which verifies that (13) holds, as $\pi^{-}\leq\pi^{+}$.

  * •

Case 4 ($z^{\dagger}_{\cal N}\leq 0,z^{\dagger}_{\cal T}\geq 0$): We have two
sub-cases here

    * –

Case 4.1 ($z^{\dagger}_{{\cal N}\cup{\cal T}}\geq 0$): The LHS and RHS in (13)
are, respectively:

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal T}})$
$\displaystyle=\pi^{+}(z^{\dagger}_{\cal N}+z^{\dagger}_{\cal T})$
$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}})+P^{\pi}(z^{\dagger}_{{\cal
T}})$ $\displaystyle=\pi^{-}(z^{\dagger}_{\cal N})+\pi^{+}(z^{\dagger}_{\cal
T}),$

By computing the difference

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal
T}})-P^{\pi}(z^{\dagger}_{{\cal N}})-P^{\pi}(z^{\dagger}_{{\cal
T}})=(\pi^{+}-\pi^{-})z^{\dagger}_{\cal N}\leq 0,$

we verify that (13) holds, as $\pi^{-}\leq\pi^{+}$.

    * –

Case 4.2 ($z^{\dagger}_{{\cal N}\cup{\cal T}}\leq 0$): Following the same
steps in Case 4.1, we get

$\displaystyle P^{\pi}(z^{\dagger}_{{\cal N}\cup{\cal
T}})-P^{\pi}(z^{\dagger}_{{\cal N}})-P^{\pi}(z^{\dagger}_{{\cal
T}})=(\pi^{-}-\pi^{+})z^{\dagger}_{\cal T}\leq 0,$

which verifies that (13) holds, as $\pi^{-}\leq\pi^{+}$.

The four cases prove (13), which says that the payment function under the
decentralized scheduling ($\bm{d}^{\dagger}_{i},z^{\dagger}_{i}$) is
subadditive, hence $\nu^{\dagger}({\cal N})+\nu^{\dagger}({\cal
T})\leq\nu^{\dagger}({\cal N}\cup{\cal T}).$ $\blacksquare$

######  Lemma 3.

For any coalition ${\cal N}\in 2^{{\cal H}}$, the value function under
centralized DER scheduling $\nu^{\ast}({\cal N})=W^{\ast}(\bm{d}^{{\cal
N}^{\ast}},\bm{z}^{{\cal N}^{\ast}})$ obeys

$W^{\ast}(\beta\bm{d}^{{\cal N}^{\ast}},\beta\bm{z}^{{\cal
N}^{\ast}})\geq\beta W^{\ast}(\bm{d}^{{\cal N}^{\ast}},\bm{z}^{{\cal
N}^{\ast}}),$ (14)

where $\beta\in[0,1]$.

###  V-C Lemma 3 and Proof of Lemma 3

Recall the welfare maximization problem under centralized DER scheduling

$\displaystyle W^{\ast}(\bm{d}^{{\cal N}^{\ast}},\bm{z}^{{\cal
N}^{\ast}}):=\underset{\\{\bm{d}_{i}\\}_{i=1}^{|{\cal
N}|},\\{z_{i}\\}_{i=1}^{|{\cal N}|}}{\text{ maximize}}$
$\displaystyle~{}~{}W({\cal N}):=\sum_{i\in{\cal
N}}S_{i}^{\pi}(\bm{d}_{i},z_{i})$ subject to $\displaystyle\sum_{i\in{\cal
N}}P^{\pi}(z_{i})=P^{\pi}(\sum_{i\in{\cal N}}z_{i})$ $\displaystyle\text{for
all }~{}i=1,\ldots,N$ $\displaystyle z_{i}=\bm{1}^{\top}\bm{d}_{i}-r_{i}$
$\displaystyle\underline{z}_{i}\leq z_{i}\leq\overline{z}_{i}$
$\displaystyle\underline{\bm{d}}_{i}\preceq\bm{d}_{i}\preceq\overline{\bm{d}}_{i},$

which can be more compactly written as

$\displaystyle W^{\ast}(\bm{d}^{{\cal N}^{\ast}},\bm{z}^{{\cal
N}^{\ast}}):=\underset{(\bm{d},\bm{z})\in\Omega^{{\cal N}}}{\text{
maximize}}\sum_{i\in{\cal N}}U_{i}(\bm{d}_{i})-P^{\pi}(\sum_{i\in{\cal
N}}z_{i}),$

where we used the constraint (2b), the notation $\bm{d}^{{\cal
M}^{\ast}}:=\\{\bm{d}^{{\cal M}^{\ast}}_{i}\\},\bm{z}^{{\cal
M}^{\ast}}:=\\{z^{{\cal M}^{\ast}}_{i}\\},\forall i\in{\cal M}\in 2^{{\cal
H}}$, and the feasible set notation

$\displaystyle\Omega^{{\cal N}}:=\\{$
$\displaystyle(\bm{d}_{1},\ldots,\bm{d}_{N}),(z_{1},\ldots,z_{N}):z_{i}=\bm{1}^{\top}\bm{d}_{i}-r_{i},$
$\displaystyle\underline{z}_{i}\leq
z_{i}\leq\overline{z}_{i},\underline{\bm{d}}_{i}\preceq\bm{d}_{i}\preceq\overline{\bm{d}}_{i}\text{
for all }i=1,\ldots,N\\}.$

Given that the constraints are linear/affine, we can scale the decision
variables by $\beta$ without changing the feasible set $\Omega^{{\cal N}}$.

Now, we solve

$\displaystyle W^{\ast}(\beta\bm{d}^{{\cal N}^{\ast}},\beta\bm{z}$ ${}^{{\cal
N}^{\ast}})=\underset{(\bm{d},\bm{z})\in\Omega^{{\cal
N}}}{\text{maximize}}\bigg{(}\sum_{i\in{\cal
N}}U_{i}(\beta\bm{d}_{i})-P^{\pi}(\sum_{i\in{\cal N}}\beta z_{i})\bigg{)}$
$\displaystyle\stackrel{{\scriptstyle\text{(A)}}}{{=}}\underset{(\bm{d},\bm{z})\in\Omega^{{\cal
N}}}{\text{maximize}}\bigg{(}\sum_{i\in{\cal N}}U_{i}(\beta\bm{d}_{i})-\beta
P^{\pi}(\sum_{i\in{\cal N}}z_{i})\bigg{)}$
$\displaystyle\stackrel{{\scriptstyle\text{(B)}}}{{\geq}}\underset{(\bm{d},\bm{z})\in\Omega^{{\cal
N}}}{\text{maximize}}\bigg{(}\beta\sum_{i\in{\cal N}}U_{i}(\bm{d}_{i})-\beta
P^{\pi}(\sum_{i\in{\cal N}}z_{i})\bigg{)}$
$\displaystyle=\beta\cdot\underset{(\bm{d},\bm{z})\in\Omega^{{\cal
N}}}{\text{maximize}}\bigg{(}\sum_{i\in{\cal
N}}U_{i}(\bm{d}_{i})-P^{\pi}(\sum_{i\in{\cal N}}z_{i})\bigg{)}$
$\displaystyle=\beta W^{\ast}(\bm{d}^{{\cal N}^{\ast}},\bm{z}^{{\cal
N}^{\ast}}),$

where step (A) is because, for any $x\in\mathbb{R}$ and $\beta\in[0,1]$,

$\displaystyle P^{\pi}(\beta x)$ $\displaystyle=\pi^{+}[\beta
x]^{+}-\pi^{-}[\beta x]^{-}$
$\displaystyle=\beta\pi^{+}[x]^{+}-\beta\pi^{-}[x]^{-}=\beta P^{\pi}(x),$

and step (B) is due to the concavity, additivity and $U_{i}(\bm{0})=0$
properties of $U_{i}(\cdot)$, which yield, for all $i\in{\cal N},k\in{\cal
K}$,

$\displaystyle U_{ik}(\beta d_{ik})$ $\displaystyle=U_{ik}(0+\beta d_{ik})$
$\displaystyle\geq(1-\beta)U_{ik}(0)+\beta U_{ik}(d_{ik})=\beta
U_{ik}(d_{ik}),$

which completes the proof. $\blacksquare$

###  V-D Proof of Theorem 1

We prove Theorem 1 by leveraging Lemma 1 and Lemma 3. The balancedness of the
cooperative game under centralized DER scheduling (${\cal H},\nu^{\ast}$) is
proved, if we show

$\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal N})\nu^{\ast}({\cal
N})\leq\nu^{\ast}({\cal H}).$ (16)

By reformulating the LHS of (16), we get

$\displaystyle\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal N})\nu^{\ast}({\cal
N})=\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal N})W^{\ast}(\bm{d}^{{\cal
N}^{\ast}},\bm{z}^{{\cal N}^{\ast}})$
$\displaystyle\stackrel{{\scriptstyle\text{(Lem.\ref{lem:homogenity})}}}{{\leq}}\sum_{{\cal
N}\in 2^{\cal H}}W^{\ast}(\alpha({\cal N})\bm{d}^{{\cal
N}^{\ast}},\alpha({\cal N})\bm{z}^{{\cal N}^{\ast}})$
$\displaystyle\stackrel{{\scriptstyle\text{(Lem.\ref{lem:AdditiveCentral})}}}{{\leq}}W^{\ast}(\sum_{{\cal
N}\in 2^{\cal H}}\alpha({\cal N})\bm{d}^{{\cal N}^{\ast}},\sum_{{\cal N}\in
2^{\cal H}}\alpha({\cal N})\bm{z}^{{\cal N}^{\ast}})$
$\displaystyle=W^{\ast}(\sum_{{\cal N}\in 2^{{\cal H}}}\alpha({\cal
N})\bm{1}_{{\cal N}}(i)\bm{d}^{{\cal H}^{\ast}},\sum_{{\cal N}\in 2^{{\cal
H}}}\alpha({\cal N})\bm{1}_{{\cal N}}(i)\bm{z}^{{\cal H}^{\ast}})$
$\displaystyle\stackrel{{\scriptstyle\text{(Def.\ref{def:BalancedGames})}}}{{=}}W^{\ast}(\bm{d}^{{\cal
H}^{\ast}},\bm{z}^{{\cal H}^{\ast}})=\nu^{\ast}({\cal H}),$

which completes the proof. $\blacksquare$

###  V-E Proof of Theorem 2

We prove the balancedness of the game (${\cal H},\nu^{\dagger}$) using the
subadditivity of the payment function under decentralized scheduling
($\bm{d}^{\dagger}_{i},z^{\dagger}_{i}$) shown in Lemma 2. From Def.3 we need
to show that

$\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal N})\nu^{\dagger}({\cal
N})\leq\nu^{\dagger}({\cal H}).$ (17)

By reformulating the LHS of (17), we get

$\displaystyle\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal N})\nu^{\dagger}({\cal
N})$ $\displaystyle=\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal
N})W^{\dagger}({\cal
N})\stackrel{{\scriptstyle\text{(\ref{eq:Welfare1Mod})}}}{{=}}\sum_{{\cal
N}\in 2^{\cal H}}\alpha({\cal N})\bigg{(}\sum_{i\in{\cal
N}}U_{i}(\bm{d}_{i}^{\dagger})-P^{\pi}(z_{{\cal N}}^{\dagger})\bigg{)}$
$\displaystyle=\sum_{i\in{\cal H}}\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal
N})\bm{1}_{{\cal N}}(i)U_{i}(\bm{d}_{i}^{\dagger})-\sum_{{\cal N}\in 2^{\cal
H}}\alpha({\cal N})P^{\pi}(z_{{\cal N}}^{\dagger})$
$\displaystyle\stackrel{{\scriptstyle\text{(Def.\ref{def:BalancedGames})}}}{{=}}\sum_{i\in{\cal
H}}U_{i}(\bm{d}_{i}^{\dagger})-\sum_{{\cal N}\in 2^{\cal H}}\alpha({\cal
N})P^{\pi}(z_{{\cal N}}^{\dagger})$
$\displaystyle\stackrel{{\scriptstyle\text{(A)}}}{{=}}\sum_{i\in{\cal
H}}U_{i}(\bm{d}_{i}^{\dagger})-\sum_{{\cal N}\in 2^{\cal
H}}P^{\pi}(\alpha({\cal N})z_{{\cal N}}^{\dagger})$
$\displaystyle\stackrel{{\scriptstyle\text{(Lem.\ref{lem:AdditiveDecentral})}}}{{\leq}}\sum_{i\in{\cal
H}}U_{i}(\bm{d}_{i}^{\dagger})-P^{\pi}(\sum_{{\cal N}\in 2^{\cal
H}}\alpha({\cal N})z_{{\cal N}}^{\dagger})$ $\displaystyle=\sum_{i\in{\cal
H}}U_{i}(\bm{d}_{i}^{\dagger})-P^{\pi}(\sum_{i\in{\cal H}}\sum_{{\cal N}\in
2^{\cal H}}\alpha({\cal N})\bm{1}_{{\cal N}}(i)z_{i}^{\dagger})$
$\displaystyle\stackrel{{\scriptstyle\text{(Def.\ref{def:BalancedGames})}}}{{=}}\sum_{i\in{\cal
H}}U_{i}(\bm{d}_{i}^{\dagger})-P^{\pi}(z_{{\cal
H}}^{\dagger})\stackrel{{\scriptstyle\text{(\ref{eq:Welfare1Mod})}}}{{=}}W^{\dagger}({\cal
H})=\nu^{\dagger}({\cal H}),$

where in step (A) we used the homogeneity of $P^{\pi}(\cdot)$ under
decentralized DER scheduling, i.e., for ${\cal N}\in 2^{{\cal H}}$ and scalar
$\beta\geq 0$,

$P^{\pi}(\beta z_{{\cal N}}^{\dagger})=\pi^{+}[\beta z_{\cal N}]^{+}-[\beta
z_{\cal N}]^{-}=\beta P^{\pi}(z_{{\cal N}}^{\dagger}),$

which completes the proof. $\blacksquare$