# The extremal generalised Randić index for a given degree range

John Haslegrave Lancaster University, UK, email<EMAIL_ADDRESS>
###### Abstract

O and Shi proved that the Randić index of any graph $G$ with minimum degree at
least $\delta$ and maximum degree at most $\Delta$ is at least
$\frac{\sqrt{\delta\Delta}}{\delta+\Delta}|G|$, with equality if and only if
the graph is $(\delta,\Delta)$-biregular. In this note we give a short proof
via a more general statement. We apply the latter to classify the graphs in
any given degree range which minimise (or maximise) the generalised Randić
index for any exponent, and describe the transitions between different types
of behaviour precisely.

The Randić index, introduced in 1975 [7], is one of the most-studied
topological graph indices: see [5] and references therein. It is defined as

$R(G)=\sum_{uv\in E(G)}(d(u)d(v))^{-1/2}.$ (1)

Among connected $n$-vertex graphs, or more generally those with all degrees
positive, the Randić index is minimised by the star, with value $\sqrt{n-1}$
[2]. However, for bounded-degree graphs it is $\Theta(n)$, and O and Shi [6],
improving on an earlier bound due to Aouchiche, Hansen and Zheng [1], gave the
best-possible lower bound for graphs with all degrees in a specified range.
The maximum value of the Randić index in a given range is less interesting,
being attained by any regular graph independent of degree.

Here we give a short proof of a more general statement that immediately
implies the result of O and Shi, but applies to any degree-based topological
index (which we define below). We then apply our result to the generalised
Randić index $R_{\alpha}(G)$, introduced by Bollobás and Erdős [2], where the
exponent $-1/2$ in (1) is replaced by a fixed real number $\alpha$. We give a
complete classification of the minimal value of the generalised Randić index
for graphs in any given degree range and for any value of $\alpha$, which is
divided into three regimes. We remark that similar, but more complicated,
methods were used to analyse two specific degree-based indices by Deng,
Balachandran and Elumalai [4], building on work of Dalfó [3].

For integers $r,s\geq 1$, we say a graph is $(r,s)$-biregular if every vertex
has degree either $r$ or $s$, and every edge has one vertex of each degree.
(Note that $(r,r)$-biregular is synonymous with $r$-regular.) We define a
degree-based (topological) index to be a graph invariant of the form

$F(G)=\sum_{uv\in E(G)}f(d(u),d(v))$ (2)

for some function $f$ satisfying $f(x,y)\equiv f(y,x)$. The following result
gives sharp bounds on the extremal values for graphs with degrees in a
specified range (noting that maximising $F(G)$ is equivalent to minimising
$-F(G)$).

######  Theorem 1.

Fix integers $\Delta\geq\delta\geq 1$ and let $F$ satisfy (2). Let
$a,b\in\\{\delta,\ldots,\Delta\\}$ with $a\leq b$ be chosen to minimise
$\frac{ab\,f(a,b)}{a+b}$. Then any graph $G$ with minimum degree at least
$\delta$ and maximum degree at most $\Delta$ satisfies
$F(G)\geq\frac{abf(a,b)}{a+b}|G|$. Equality occurs if $G$ is
$(a,b)$-biregular, and this is the only case when equality occurs provided the
choice of $a,b$ is unique.

###### Proof.

Consider the edge weighting $w:E(G)\to\mathbb{R}^{+}$ given by
$w(uv)=1/d(u)+1/d(v)$. Note that

$\sum_{e\in E(G)}w(e)=\sum_{u\in V(G)}\sum_{v\in N(u)}1/d(u)=\sum_{u\in
V(G)}1=|G|.$

For any edge $uv$, where without loss of generality $d(u)\leq d(v)$, we have

$\frac{f(d(u),d(v))}{w(uv)}=\frac{d(u)d(v)f(d(u),d(v))}{d(u)+d(v)}\geq\frac{ab\,f(a,b)}{a+b},$
(3)

with equality if, and provided $a,b$ are uniquely determined only if, $d(u)=a$
and $d(v)=b$. Thus

$F(G)\geq\sum_{e\in
E(G)}\frac{ab\,f(a,b)}{a+b}w(e)=\frac{ab\,f(a,b)}{a+b}|G|,$

with equality if and only if equality holds in (3) for every edge. ∎

For $f(a,b)=(ab)^{-1/2}$, corresponding to the Randić index, we have

$\frac{ab\,f(a,b)}{a+b}=(x+1/x)^{-1},$

where $x=\sqrt{b/a}$. Since $x+1/x$ is increasing for $x\geq 1$, this is
minimised (uniquely) when $(a,b)=(\delta,\Delta)$ and maximised when $a=b$.
Thus we immediately obtain the following corollary.

######  Corollary 2.

If $G$ has minimum degree at least $\delta$ and maximum degree at most
$\Delta$ then

$\frac{\sqrt{\delta\Delta}}{\delta+\Delta}|G|\leq R(G)\leq\frac{1}{2}|G|,$

with the lower bound being achieved if and only if $G$ is
$(\delta,\Delta)$-biregular and the upper bound being achieved if and only if
every component of $G$ is regular.

We conclude by applying Theorem 1 to the generalised Randić index
$R_{\alpha}(G)$.

Figure 1: The three regions having different minimising graphs for the
generalised Randić index.

######  Theorem 3.

Fix integers $\Delta>\delta>0$ and write $c=\Delta/\delta$. If $G$ has minimum
degree at least $\delta$ and maximum degree at most $\Delta$, then:

  1. 1.

for $\alpha\leq\log_{c}\Bigl{(}\frac{1+c}{2c}\Bigr{)}$, we have
$R_{\alpha}(G)\geq\frac{\Delta^{1+2\alpha}}{2}|G|$, with equality if $G$ is
$\Delta$-regular;

  2. 2.

for
$\log_{c}\Bigl{(}\frac{1+c}{2c}\Bigr{)}\leq\alpha\leq\log_{c}\Bigl{(}\frac{2}{1+c}\Bigr{)}$,
we have $R_{\alpha}(G)\geq\frac{(\delta\Delta)^{1+\alpha}}{\delta+\Delta}|G|$,
with equality if $G$ is $(\delta,\Delta)$-biregular;

  3. 3.

for $\log_{c}\Bigl{(}\frac{2}{1+c}\Bigr{)}\leq\alpha$ we have
$R_{\alpha}(G)\geq\frac{\delta^{1+2\alpha}}{2}|G|$, with equality if $G$ is
$\delta$-regular.

These are the only graphs where equality occurs unless $\alpha$ lies on the
boundary between two regimes, in which case both equality cases apply.

Additionally, for $\alpha<-1/2$ we have
$R_{\alpha}(G)\leq\frac{\delta^{1+2\alpha}}{2}|G|$, with equality if and only
if $G$ is $\delta$-regular, whereas for $\alpha>-1/2$ we have
$R_{\alpha}(G)\leq\frac{\Delta^{1+2\alpha}}{2}|G|$, with equality if and only
if $G$ is $\Delta$-regular.

###### Proof.

Consider optimising the function $g(x,y)=\frac{x+y}{(xy)^{1+\alpha}}$ for
$x,y\in[\delta,\Delta]$, and write $c=\delta/\Delta$. Note that

$\frac{\partial g(x,y)}{\partial x}=\frac{-\alpha
x-(1+\alpha)y}{x^{2+\alpha}y^{1+\alpha}}\quad\text{and}\quad\frac{\partial^{2}g(x,y)}{\partial
x^{2}}=\frac{\alpha(1+\alpha)x+(1+\alpha)(2+\alpha)y}{x^{3+\alpha}y^{1+\alpha}},$

and similarly for $y$.

For $\alpha\leq-1$, $g$ is increasing in both coordinates and so is maximised
for $x=y=\Delta$. Similarly, $g$ is maximised for $x=y=\delta$ if $\alpha\geq
0$.

We now consider the case $-1\leq\alpha\leq-1/2$. Here $g$ is an increasing
function of the larger coordinate, and a convex function of the smaller. It
follows that the maximum occurs only when the larger coordinate is $\Delta$
and the smaller is either $\delta$ or $\Delta$. Comparing $g(\delta,\Delta)$
with $g(\Delta,\Delta)$, we see that the former is larger when $1+c\geq
2c^{1+\alpha}$, which is equivalent to
$\alpha\geq\log_{c}\Bigl{(}\frac{1+c}{2c}\Bigr{)}$.

The case $-1/2\leq\alpha\leq 0$ is similar, with inequalities reversed and
$\delta,\Delta$ exchanged. Thus we compare $g(\delta,\Delta)$ with
$g(\delta,\delta)$, with the former being larger if $1+c\geq 2c^{-\alpha}$,
which is equivalent to $\alpha\leq\log_{c}\Bigl{(}\frac{2}{1+c}\Bigr{)}$.

By Theorem 1, $R_{\alpha}(G)\geq M^{-1}|G|$, where $M$ is the maximum value of
$g(x,y)$, and this, together with consideration of equality cases, gives the
claimed lower bounds.

For the upper bounds, we instead minimise $g$. Here, for any $\alpha<-1/2$ we
have $g$ increasing in $x$ for $x\geq\frac{-1-\alpha}{\alpha}y$, and so from
any $(x,y)\neq(\delta,\delta)$ the value of $g(x,y)$ can be strictly decreased
by reducing one coordinate. Thus the minimum uniquely occurs at
$g(\delta,\delta)$, and similarly for $\alpha>-1/2$ it uniquely occurs at
$g(\Delta,\Delta)$. ∎

## Acknowledgments

I am grateful to Saieed Akbari for bringing the result of [6] to my attention.

## References

  * [1] M. Aouchiche, P. Hansen and M. Zheng,  Variable neighborhood search for extremal graphs. XIX. Further conjectures and results about the Randić index.  MATCH Commun. Math. Comput. Chem. 58 (2007), 83–102.
  * [2] B. Bollobás and P. Erdős,  Graphs of extremal weights.  Ars Combin. 50 (1998), 225–233.
  * [3] C. Dalfó,  On the Randić index of graphs.  Discrete Math. 342:10 (2019), 2792–2796.
  * [4] H. Deng, S. Balachandran and S. Elumalai,  Some tight bounds for the harmonic index and the variation of the Randić index of graphs.  Discrete Math. 342:7 (2019), 2060–2065.
  * [5] X. Li and Y. Shi,  A survey on the Randić index.  MATCH Commun. Math. Comput. Chem. 59 (2008), 127–156.
  * [6] S. O and Y. Shi,  Sharp bounds for the Randić index of graphs with given minimum and maximum degree. Discrete Appl. Math. 247 (2018), 111–115.
  * [7] M. Randić,  Characterization of molecular branching.  J. Amer. Chem. Soc. 97 (1975), 6609–6615.


# Exotic states induced by co-evolving connection weights and phases in
complex networks

S. Thamizharasan1, V. K. Chandrasekar2, M. Senthilvelan1, Rico Berner3,4,
Eckehard Schöll3 and D. V. Senthilkumar5 1 Department of Nonlinear Dynamics,
School of Physics, Bharathidasan University, Tiruchirappalli - 620 024, Tamil
Nadu, India,
2Centre for Nonlinear Science & Engineering, Department of Physics, School of
Electrical & Electronics Engineering, SASTRA Deemed University, Thanjavur -613
401, Tamil Nadu, India.
3 Institut für Theoretische Physik, Technische Universität Berlin,
Hardenbergstrasse 36, 10623 Berlin, Germany
4 Institut für Physik, Humboldt-Universität zu Berlin, Newtonstraße 15, 12489
Berlin, Germany
5School of Physics, Indian Institute of Science Education and Research,
Thiruvananthapuram-695 551, Kerala, India.

###### Abstract

We consider an adaptive network, whose connection weights co-evolve in
congruence with the dynamical states of the local nodes that are under the
influence of an external stimulus. The adaptive dynamical system mimics the
adaptive synaptic connections common in neuronal networks. The adaptive
network under external forcing displays exotic dynamical states such as
itinerant chimeras whose population density of coherent and incoherent domains
co-evolves with the synaptic connection, bump states and bump frequency
cluster states, which do not exist in adaptive networks without forcing. In
addition the adaptive network also exhibits partial synchronization patterns
such as phase and frequency clusters, forced entrained, and incoherent states.
We introduce two measures for the strength of incoherence based on the
standard deviation of the temporally averaged (mean) frequency and on the mean
frequency in order to classify the emergent dynamical states as well as their
transitions. We provide a two-parameter phase diagram showing the wealth of
dynamical states. We additionally deduce the stability condition for the
frequency-entrained state. We use the paradigmatic Kuramoto model of phase
oscillators which is a simple generic model that has been widely employed in
unraveling a plethora of cooperative phenomena in natural and man-made
systems.

###### pacs:

##  I Introduction

The paradigm of complex networks hosts a plethora of complex collective
dynamical states observed in a wide variety of natural and technological
systems sbvl2006 ; aaad2008 . Investigations of the dynamics on and of complex
networks have been an active area of research for more than a decade tgbb2008
; tntv2012 ; trb2021 . In particular, adaptive networks are sophisticated
complex networks, a basis for smart systems that self-adapt, in which the
connection weights co-evolve along with the dynamical states of the network.
For instance, neurophysiological experiments have revealed that spike timing
differences between pre and post synaptic neurons determine the evolution of
synaptic connections hmjl1997 . Other examples include the reaction rate
evolving as a function of the state variable in chemical systems sjsk2001 ,
state-dependent plasticity in epidemics, biological and social systems
tgcjdd2006 ; fseo2008 ; lhck2018 . Such co-evolution of (synaptic) network
connections is thought to provide a basis for higher order brain functions.
Recently, adaptive networks have been intensively investigated and shown to
exhibit exotic intriguing dynamical states including chimera states tata2009 ;
cbp2011 ; rgaa2011 ; lt2014 ; xzsb2015 ; llovp2016 ; dvk2016 ; scpi2017 ;
dvk2017 ; vajaa2018 ; mhm2019 ; vrrb2019 ; rbes2019 ; rbjf2019 ; rbjs2020 .

In particular, phenomena such as self-organization of hierarchical
multilayered structures resulting in multifrequency clusters or chimera states
are peculiar to the adaptive networks dvk2017 ; rbes2019 . Specifically, a
slow adaptation mechanism is identified as pivotal for the emergence of stable
multi-clusters. A special type of metastable cluster states called ”double
antipodal” clusters have been shown to play an important role in achieving the
asymptotic dynamics of the adaptive networks rbjf2019 . Very recently, partial
synchronization patterns like phase clusters and more complex stable clusters
such as splay states, antipodal states and double antipodal states have been
found to emerge due to the delicate balance between adaption and multiplexing
complex networks, which are otherwise unstable in single layer networks
rbjs2020 . The results provide exclusive evidence that network adaptation
provides a mechanism giving rise to a variety of novel dynamical scenarios
that are, however, far from our current understanding. Therefore, there is a
need for much in-depth investigations.

External stimuli are an inevitable factor to be accounted for along with the
adaptive networks as there exist a plethora of phenomena evoked by the
external force. For instance, synaptic connections themselves evolve in
response to the external stimulus in brain networks gpcn1994 . Other examples
include circadian rhythms mcm1982 ; jcd2003 ; rgf2005 , resonances izk2002 ;
izk2008 ; aaz2000 ; mr2008 ; eiv2003 ; asp1997 ; czjk2001 , event-related
desynchronization and synchronization gpcn1994 . Hence, it is very natural to
extend the study on adaptive networks to systems with external stimulus. Phase
diagrams of the forced Kuramoto model have been studied in great detail
ssyk1986 ; jaal2005 ; far2016 . In this work, we consider an adaptive network
of phase oscillators, where the coupling weights co-evolve along with phases,
mimicking the spike-timing-dependent synaptic connections between neurons.
Under the influence of an external forcing, we will show the existence of
several exotic self-organizing dynamical states such as conventional chimera
states characterized by static coherent and incoherent groups, an only
recently discovered type of adaptive chimera termed as itinerant chimera of
two distinct types, characterized by co-evolving coherent and incoherent
groups along with the adaptive coupling, bump states characterized by a
localized inactive coherent and coexisting incoherent group, and a bump
frequency cluster among other intriguing collective dynamical patterns. We use
the Kuramoto order parameter, spatio-temporal plots, instantaneous phase
profiles and mean frequency profiles of the oscillators to illustrate the
various self-organizing dynamical states. We also introduce the strength of
incoherence $S$, based on the standard deviation of the mean frequencies, and
$\hat{S}$, based on the mean phase frequencies along with the Kuramoto order
parameters to classify the dynamical states and their transitions in two-
parameter phase diagrams. We also deduce the stability condition for the
frequency-entrained state.

Figure 1: First and second rows depict the order parameters $R_{(1,2)}(t)$ and
space-time plot in the range $t\in(26400,26500)$, while the third and fourth
rows depict a snapshot of instantaneous phases taken at $t=26475$, and
temporally averaged mean frequencies, respectively. The mean is taken over
$t\in(24000,26500)$. (a)-(d) Two-cluster state for $\alpha=0.40$ and $f=0.0$,
(e)-(h) chimera state for $\alpha=0.40$ and $f=0.78$, and (i)-(l) incoherent
state for $\alpha=1.54$ and $f=0.20$. Parameters: $\lambda=1$,
$\varepsilon=0.005$.

The structure of the paper is as follows. In section II, we introduce our
model and discuss its constituents in detail. We unravel the different
dynamical states exhibited by the globally coupled identical phase oscillators
under the influence of the external forcing in Sec. III. We introduce the
strength of incoherence to quantify and classify the different dynamical
states in Sec. IV. Dynamical transitions in one- and two-parameter space are
discussed in Sec. V. Finally, we provide a summary and conclusion in Sec. VI.

##  II Model

We consider a network of globally coupled identical phase oscillators under
the influence of an external forcing represented by

$\displaystyle\dot{\psi_{i}}$ $\displaystyle=$
$\displaystyle\sigma-\frac{1}{N}\sum_{j=1}^{N}\Gamma_{ij}(\psi_{i}-\psi_{j})+fG(\psi_{i},\Omega),$
(1)

where, $i=1,2,3,\ldots,N$, $\psi_{i}\in[0,2\pi)$ is the phase of the $i^{th}$
oscillator, $\sigma$ is the identical natural frequency of the oscillators,
$f$ and $\Omega$ are the strength and frequency of the external drive,
respectively. When a coupled oscillator model with external forcing term is
reduced to a phase oscillator model, the forcing term will turn out to depend
on the phases of the oscillators. The latter form corresponds to the active
rotator model that can be regarded as a paradigmatic model for neurons with
type-1 excitability similar to Morris-Lecar models bljgo2004 . The coupling
function is given by $\Gamma_{ij}(\psi)=k_{ij}\sin(\psi+\alpha)$, where
$k_{ij}$ is the coupling weight from $j^{th}$ to $i^{th}$ oscillator and
$\alpha$ quantifies the phase-lag induced by a small coupling delay mmav2018 .
It is to be noted that $\alpha$ maybe also seen as the manifestation of the
nonisochronocity parameter of the Stuart-Landau oscillator in the reduced
phase model jkbook . The adaptation rule for $k_{ij}$ is given by

$\displaystyle\dot{k}_{ij}=\varepsilon\Lambda_{ij}(\psi_{i}-\psi_{j}),$ (2)

where we restrict the coupling values to the region $|k_{ij}(t)|\leq 1$ and
the time scale $\varepsilon^{-1}$ is considered to be very much longer than
that of the phase oscillators. We choose $2\pi$-periodic functions for the
evolution of $k_{ij}$ and $G(\psi,\Omega)$. For simplicity, $\Lambda(\psi)$
and $G(\psi,\Omega)$ take the form $\Lambda(\psi)=\sin(\psi)$ and
$G(\psi,\Omega)=\sin(\psi-\Omega t)$ with the lowest-order Fourier mode.

Note that the evolution of the coupling strength has been chosen as

$\displaystyle\dot{k}_{ij}=\varepsilon\sin(\psi_{i}-\psi_{j}+\beta)$ (3)

to account for the co-evolution of the coupling weights and the dynamical
states of the nodes of the network tata2009 . Specifically, neurophysiological
experiments indicate that the change in the strength of the synaptic coupling
between neurons depends on the relative timing of the pre- and postsynaptic
spikes gqbmmp1998 . Hence, it is natural that the dynamics of the coupling
weights depend on the relative timing of the oscillators. In particular,
models of phase oscillators have been successfully used to explain a variety
of dynamical mechanisms induced by spike-timing dependent plasticity llovp2016
; vrrb2019 .

Depending on the value of $\beta$, the dynamics of the coupling weights mimics
Hebbian-like function, spike-timing dependent plasticity-like function and
anti-Hebbian-like function tata2009 . In particular, when $\beta=0$,
$\dot{k}_{ij}=\varepsilon\sin(\psi_{i}-\psi_{j})$, this situation is
essentially the same as that in the case of spike-timing dependent plasticity.
We consider the case of $\beta=0$ and investigate the effect of external
stimuli in terms of an external forcing on the emerging dynamical states.

A two-cluster state, a coherent state with a fixed phase relation and an
incoherent state with frustration have been reported with the above evolution
equation for the coupling weights in the Kuramoto model in the appropriate
ranges of $\beta$ tata2009 . When the connection weights are fixed and simply
provide heterogeneity in the interaction strength between the Kuramoto
oscillators without any specific evolution rule, it was shown that
synchronization can be completely inhibited when the weights are strongly
anti-correlated, otherwise the synchronization transition observed in the
standard Kuramoto model has been generalized to this case of heterogeneous
coupling ghpdhz2008 .

Introducing $\psi_{i}=\theta_{i}+\Omega t$, the globally coupled identical
phase oscillators with the external forcing can be represented in the rotating
frame as

$\displaystyle\dot{\theta_{i}}$ $\displaystyle=$
$\displaystyle\lambda-\frac{1}{N}\sum_{j=1}^{N}k_{ij}\sin(\theta_{i}-\theta_{j}+\alpha)+f\sin(\theta_{i}),$
(4)

with $\lambda=\sigma-\Omega$. The latter form of the phase oscillator model is
also known from studies on active rotators ifsy2020 . Throughout the paper, we
consider $N=100$ oscillators, and initial conditions for which $\theta_{i}$
and $k_{ij}$ are uniformly distributed in the interval $[0,2\pi)$ and
$(-1,1),\forall j$, respectively, $\lambda=1$ and $\varepsilon=0.005$.

##  III Dynamical states

The Kuramoto order parameter

$\displaystyle
R_{l}=\left|\frac{1}{N}\sum_{j=1}^{N}e^{il\theta_{j}}\right|,~{}l=1,2$ (5)

characterizes the nature of the dynamical states. $R_{1}=1$ holds for the
completely synchronized state, while $R_{2}=1$ characterizes two cluster
(antipodal) states, the latter including the completely synchronized state.
The order parameters, spatio-temporal evolution, instantaneous phases and mean
frequencies of all oscillators are depicted in the first to fourth row,
respectively, of Figs. 1, 2 and 5. The unity value of $R_{2}$ in Fig. 1(a) for
$\alpha=0.4$ and $f=0.0$ elucidates that the network exhibits a cluster state,
here a two-cluster state. However, a small but non-vanishing value of $R_{1}$
indicates that the two clusters are not of equal size, which is also
corroborated by the space-time plot and the snapshot of instantaneous phases
in Figs. 1(b) and 1(c), respectively. Note that their mean frequency profile
$\omega_{j}=\langle\dot{\theta}_{j}\rangle$, where $\langle\cdot\rangle$
represents a long time average, depicted in Fig. 1(d), shows that the
oscillators are frequency entrained. The instantaneous phases are obtained at
$t=26475$ in all the figures of the manuscript unless otherwise specified.

Figure 2: Same plots as in Fig. 1 for different parameters $\alpha$ and $f$.
(a)-(d) itinerant chimera of first kind (IC1) for $\alpha=1.54$ and $f=0.78$,
(e)-(h) bump state $\alpha=0.10$, and $f=0.90$ and (i)-(l) itinerant chimera
of second kind (IC2) for $\alpha=1.54$ and $f=0.99$. Figure 3: Phase snapshots
for an itinerant chimera of the first kind (IC1) at different instants of
time. (a) $t=24004$, (b) $t=24007$, (c) $t=24013$ and (d) $t=24034$. The
parameter values are the same as in Fig. 2 for IC1.

The network of globally coupled phase oscillators also exhibits a chimera
state upon increasing the forcing amplitude to $f=0.78$ for the same $\alpha$
(see the middle column of Fig. 1). The order parameters $R_{(1,2)}$ in Fig.
1(e) vary within a small range, which is attributed to the coexistence of
coherent and incoherent domains, characterizing the chimera state. Distributed
frequencies among the oscillators constituting the incoherent domain result in
the time-dependent variations of the order parameters. The coherent and
incoherent domains of phase oscillators are well visible in the space-time
plot and the snapshot of the instantaneous phases depicted in Figs. 1(f) and
1(g), respectively. The mean frequency profile (see Fig. 1(h)) showing
distinct groups of phase oscillators with distributed and entrained
frequencies is yet another evidence for the existence of chimera in the
network of phase oscillators. In addition, we have further analyzed the
emergence of chimera states with respect to the number of coupled oscillators
$N$. Remarkably, we observe chimera states even for a small number of
oscillators, see Appendix C.

Random evolution of the instantaneous phase of the oscillators is observed for
$\alpha=1.54$ and $f=0.20$. Fluctuating values of $R_{(1,2)}$ near zero in
Fig. 1(i) confirm the incoherent nature of the phase evolution. Further, the
space-time plot and the snapshot of the instantaneous phases in Figs. 1(j) and
1(k), respectively, corroborate the random nature of the phase oscillators.
Randomness is also observed in the mean frequency profile (see Fig. 1(l))
corroborating the chaotic distribution of their mean frequencies.

The chimera states reported so far in the literature are mostly characterized
by coexisting coherent and incoherent domains with constant population
density. However, turbulent and rotating chimeras are characterized by
irregular and periodic temporal evolution of the coherent domains fpkswh2016 ;
oeoek2019 . In this work, we show the emergence of a recently discovered and
thus rather less understood type of chimera characterized by coexisting
coherent and incoherent groups, whose population densities co-evolve with the
connection weights and the phase of the oscillators.

In general, a chimera state refers to a co-existing group of coherent and
incoherent dynamical states arising out of an ensemble of identical systems.
Various types of chimera states including frequency chimera, amplitude
chimera, spiral chimera, traveling chimera, breathing chimera, etc., have been
reported in the literature fpsj2021 ; oeoek2019 . The coherent and incoherent
groups characterizing the chimera state can be qualitatively identified from
the snapshot of the instantaneous phases of the oscillators, snapshot of the
mean frequency and even from their spatiotemporal patterns. Qualitative
measures such as local order parameter, strength of incoherence have also been
in use to identify and characterize the various chimera states. More
discussion on the specific types of chimera states, their physical relevance,
experimental realizations and their characterizations can be found in the
recent reviews on chimera states fpsj2021 ; smbkb2019 ; es2020 .

Figure 4: Evolution of the coherent and incoherent groups of the itinerant
chimera. Red (light gray) corresponds to the coherent and blue (dark gray) to
the incoherent groups. The parameter values are the same as in Fig. 2. (a)
itinerant chimera of the first kind and (b) itinerant chimera of the second
kind.

The small variations of the order parameters $R_{1}$ and $R_{2}$ in Fig. 2(a)
for $f=0.78$ indicate the coexistence of coherent and incoherent domains
characterizing a chimera state. Phases evolving in concurrence with the
adaptive coupling are evident from the space-time plot in Fig. 2(b). Evolution
of connection weights is included in the appendix. Randomly distributed phases
without any phase correlation may be observed from the snapshot of the
instantaneous phases in Fig. 2(c). Nevertheless, there exists a set of
oscillators with nearly identical phases at any given point of time
constituting the coherent group while the rest of the oscillators with a drift
in their phases constitute the incoherent group which can be best observed for
a suitable sorting of the indices $i=1,\ldots,N$ as illustrated in Fig. 2(b)
and 2(c). The sorting is chosen according to increasing mean frequencies (see
Fig. 2(d).) It is to be noted that the oscillators spontaneously switch
between the coherent and incoherent domains as the oscillators self-organize
in accordance with the adaptive coupling. Such a self-organizing pattern is
reported as “traveling of the oscillators from one domain to another or as
traveling of the chimera core across the network”’ and is referred to as
“itinerant chimera” dvkvvk2019 . Since the dynamical nature of the self-
organizing patterns in the first column of Fig. 2 exactly resembles that of
the itinerant chimera, we call it an itinerant chimera of first kind. The
existence of coherent and incoherent domains of the itinerant chimera of first
kind is also evident from the mean frequency profile depicted in Fig. 2(d).

Figure 5: Same characteristic plots as in Fig. 1 for different parameters
$\alpha$ and $f$. (a)-(d) frequency cluster state for $\alpha=0.40$ and
$f=0.83$, (e)-(h) frequency cluster state with bump frequency cluster state
$\alpha=0.60$ and $f=0.95$ and (i)-(l) forced entrainment state for
$\alpha=0.40$ and $f=1.40$.

Phase snapshots of the itinerant chimera of first kind at four different
instants are illustrated in Figs. 3(a)-3(d). The oscillators constituting the
coherent domain are depicted as a filled triangle, while the oscillators
constituting the incoherent domain are indicated by filled circles. Further,
the temporal evolution of the oscillators constituting coherent (marked in
red/light gray) and incoherent (marked in blue/dark gray) groups of the
itinerant chimera of the first kind is depicted in Fig. 4(a), which resembles
the characteristics of the itinerant chimera reported by Kasatkin et al
dvkvvk2019 .

We observe a bump state upon increasing the forcing to $f=0.9$, where a
majority of the oscillators are entrained to the external forcing for
$\alpha=0.1$, while the rest evolve independently displaying an active state.
The coexistence of a coherent quiescent state and an incoherent desynchronized
oscillating state is known as a bump state mocl2007 ; crloeo2020 ; hsda2020 ;
ifoeo2021 . The corresponding order parameters are depicted in Fig. 2(e),
which oscillate near unity as a majority of the oscillators are locked to the
external forcing. The spatio-temporal evolution of this dynamical state with
coherent domain constituted by the oscillators entrained with the external
forcing and incoherent domain constituted by independently evolving phase
oscillators is shown in Fig. 2(f). The snapshot of the instantaneous phases
(see Fig. 2(g)) and the mean frequency profile (see Fig. 2(h)) also
corroborate the coexistence of two distinct domains. Note that the phase
oscillators with near zero mean frequency and the oscillators with a finite
value of average frequencies together constitute the bump state. Let
$\theta_{i}=\theta^{s}$, $i=1,2,\ldots,N_{s}$ be the oscillators locked to the
external force and $\theta_{i}=\theta_{i}^{d}$, $i=N_{s}+1,N_{s}+2,\ldots,N$
be the oscillators evolving independently, then the instantaneous angular
velocity corresponding to the coherent quiescent state and the incoherent
desynchronized oscillating state can be written as

$\displaystyle\dot{\theta^{s}}$
$\displaystyle=\lambda-\frac{N_{s}}{N}k_{1}\sin\alpha-\frac{1}{N}\sum_{j=N_{s}+1}^{N}k_{sj}\sin(\theta^{s}-\theta_{j}^{d}+\alpha)+f\sin\theta^{s},$
(6a) $\displaystyle\dot{\theta_{i}^{d}}$
$\displaystyle=\lambda-\frac{N_{s}}{N}k_{2}\sin(\theta_{i}^{d}-\theta^{s}+\alpha)-\frac{1}{N}\sum_{j=N_{s}+1}^{N}k_{ij}\sin(\theta_{i}^{d}-\theta_{j}^{d}+\alpha)+f\sin\theta_{i}^{d},$
(6b)

where
$k_{1}=\frac{1}{N_{s}}\sum_{j=1}^{N_{s}}k_{ij},~{}\forall~{}i=\\{1,2,\ldots,N_{s}\\}$
and
$k_{2}=\frac{1}{N_{s}}\sum_{j=1}^{N_{s}}k_{ij},~{}\forall~{}i=\\{N_{s}+1,N_{s}+2,\ldots,N\\}$.
Indeed, we have confirmed from the numerical analysis that the entrained
oscillators as well as those evolving independently satisfy the above
equations.

The itinerant chimera of the second kind is depicted in the third column of
Fig. 2 for $\alpha=1.54$ and $f=0.99$. The essential difference between
itinerant chimera of the first and second kind is that the latter is
characterized by a high degree of phase-synchronization at certain intervals
of time as quantified by the epochs of the order parameters $R_{1,2}$ near
unity in Fig. 2(i). We can speculate that such a chimera state with a large
scale phase entrainment may be related to neuropathological states with
characteristic large-scale synchronization (epochs of abnormal/massive
synchronization) such as in epileptic seizures pjuws2006 ; mfrb2020 ; es2020 .
Small range variations of the order parameter in Fig. 2(i) indicate the
chimera state. Relatively high degree of phase coherence coexisting with a low
degree of phase coherence, compared to that of the itinerant chimera of the
first kind, is evident from the space-time plot in Fig. 2(j) and from the
snapshot of the instantaneous phases in Fig. 2(k). The indices of the globally
coupled oscillators are reordered suitably to visualize the coherent and
incoherent domains clearly according to increasing mean frequency. It is to be
noted that the range of the mean frequency distribution is rather narrow in
the range $(0,0.2)$ for the itinerant chimera of the second kind (see Fig.
2(l)), whereas that of the itinerant chimera of the first kind is distributed
in a much broader range $(0,1.0)$ corroborating their degree of phase
coherence. The plots of phase snapshots of the itinerant chimera of the second
kind at four different instants of time are illustrated in Figs. 6(a)-6(d).
The oscillators constituting the coherent group are depicted by a filled
triangle, while the oscillators constituting the incoherent group are
indicated by filled circles. It can be seen from these figures that the number
of oscillators constituting coherent and incoherent groups vary as they self-
organize in concurrence with the adaptive coupling. Further, the temporal
evolution of the oscillators constituting coherent (marked in red/light gray)
and incoherent (marked in blue/dark gray) groups of the itinerant chimera of
the second kind is depicted in Fig. 4(b).

Figure 6: Phase snapshots of the itinerant chimera of the second kind at
different instants of time. (a) $t=24040$, (b) $t=24055$, (c) $t=24280$ and
(d) $t=24317$. The parameter values are the same as in Fig. 2 for the
itinerant chimera of the second kind.

The network of globally coupled phase oscillators with external forcing
exhibits a frequency cluster state for $\alpha=0.40$ and $f=0.83$. The order
parameters $R_{1}$ and $R_{2}$ corresponding to the frequency cluster vary
between the values $0.25$ and $1$ (see Fig. 5(a)). The oscillators split into
two clusters of different frequency of oscillations (see Fig. 5(b))
constituting the frequency cluster state, which is also evident from the
snapshot of the instantaneous phases in Fig. 5(c) and from the mean frequency
profile depicted in Fig. 5(d). Similar partially synchronized patterns like
phase clusters and more complex stable clusters are shown to emerge due to the
delicate balance between adaptation and multiplexing in complex networks
rbjs2020 ; dvk2016 .

In Fig. 5(e-h), we show another frequency cluster for $\alpha=0.60$ and
$f=0.95$. Here, both order parameters $R_{1}$ and $R_{2}$ are uniformly
oscillating between $0.5$ and $1$ as shown in Fig. 5(e). The spatio-temporal
plot (see Fig. 5(f)) elucidates the coexistence of a cluster with fully
pronounced spike oscillations and a cluster with small amplitude sub-threshold
oscillations induced by the other cluster, which we call the bump frequency
cluster state. The instantaneous phases (Fig. 5(g)) and the mean frequency
profile (Fig. 5(h)) illustrate that phases of one group of oscillators are
nearly locked with the frequency of the external forcing, while those of the
second group undergo a complete cycle indicating the bump frequency cluster
state.

A forced entrained state is observed for $\alpha=0.40$ and $f=1.40$, where all
the oscillators are entrained by the external forcing. Both order parameters
acquire unity values (see Fig. 5(i)) corroborating the forced entrained state.
The oscillator phases acquire the same phase as that of the external forcing
as depicted in the space-time plot in Fig. 5(j) and in the snapshot of the
instantaneous phases in Fig. 5(k). The mean frequency profile in Fig. 5(l)
also confirms the frequency entrained state.

The coupling weights reach the frozen states asymptotically for two cluster
and forced entrained states (see Appendix A), while for all other states the
coupling weights evolve with time. The time variation of the coupling weights
are larger in the chimera, incoherent state, itinerant chimera of first and
second kind but we found very negligible variations in frequency cluster, bump
state and bump frequency cluster state. Further one to one correlation of the
dynamical states with its corresponding evolution of the coupling weight
matrices require more in depth statistical analysis of the latter. In the
following, we will introduce the strength of incoherence rgvkc2014  to
quantify the different dynamical states and to classify them in one- and two-
parameter space plots.

Figure 7: Time-averaged order parameters ($\langle R_{1}\rangle$ and $\langle
R_{2}\rangle$) and the strength of incoherence ($S$ and $\hat{S}$) as a
function of the external forcing strength $f$ elucidating the dynamical
transitions for the parameter values (a) $\alpha=0.30$ and (b) $\alpha=1.54$.
Abbreviations: TC - two cluster state, CHI - chimera state, ICS - incoherent
state, IC1 - itinerant chimera states of first and IC2 - second kind, BS -
bump state, FC - frequency cluster state, BFC - bump frequency cluster state,
FE - forced entertainment state.

##  IV Strength of incoherence

In order to define a measure of incoherence, we divide the $N$ oscillators
into $M$ (even) bins each of size $n=N/M$. We estimate the local standard
deviation ($\sigma_{m}$) of the frequencies at each bin as

$\displaystyle\sigma_{m}$
$\displaystyle=\sqrt{\frac{1}{n}\sum_{j=n(m-1)+1}^{mn}\left[\omega_{j}-\bar{\omega}_{m}\right]^{2}},$
(7) $\displaystyle\qquad\qquad\qquad\qquad\qquad~{}m=1,2,\dots,M.$

with $\bar{\omega}_{m}=\frac{1}{n}\sum_{j=n(m-1)+1}^{mn}\omega_{j}$. Note that
$\omega_{j}$ is the time averaged frequency of the $j$th oscillator. Now, the
strength of incoherence $S$ is calculated using the formula rgvkc2014

$\displaystyle S=1-\frac{\sum_{m=1}^{M}s_{m}}{M},\qquad
s_{m}=\Theta(\delta-\sigma_{m}),$ (8)

where $\Theta$ is the Heaviside step function and $\delta$ is a predefined
threshold. The local standard deviation $\sigma_{m}$ becomes zero if the
network of oscillators is frequency-entrained as in the two-cluster state (see
Fig. 1(a)-1(c)). Since $\sigma_{m}<\delta$, $s_{m}$ becomes unity for all $m$.
As a consequence, the strength of incoherence for the phase clusters with
frequency entrainment will be $S=0$. In contrast, for the evolution of phases
with random frequencies such as the incoherent state (see Fig. 1(i)-1(l)),
$\sigma_{m}$ takes a finite value larger than the predefined threshold
$\delta$ and hence $s_{m}=0$ for all $m$ resulting in the strength of
incoherence $S=1$ (refer Appendix B for more details). On the other hand, for
the case of conventional chimeras, frequencies of the phase oscillators that
constitute the coherent domain are entrained, whereas those of the incoherent
domain are completely uncorrelated, and hence $S$ takes a value between zero
and unity depending on the number of oscillators constituting the coherent
domain. Since the phases evolve with random frequencies for the itinerant
chimeras, $S$ for these states also takes the value unity similar to the
incoherent state. However, $S$ for the mean frequency profile of the itinerant
chimera of the first kind (see Fig. 2(d)) and of the second kind (see Fig.
2(l)), where the oscillator indices are reordered appropriately, takes a value
between zero and unity in analogy with the chimera state. In the bump state,
the coherent domain is entrained with the frequency of the external forcing
while the oscillators of the incoherent domain evolve with random frequencies
and hence the strength of incoherence for the bump state is given by $0<S<1$.
For the frequency clusters, there is a discontinuous jump in the frequency and
$s_{m}=0$ only in the bin where such discontinuity exists, while $s_{m}=1$ in
the other bins. Hence the strength of incoherence takes the value $1/M$ close
to zero due to frequency entrainment in the $M-1$ bins note . For a suitable
choice of $M$ the value of $S$ might be exactly zero. In any case, $S$ depends
on the choice of $M$), whereas in the case of the bump and chimera states $S$
takes a large value between $0$ and $1$ as there is an appreciable number of
bins with a large standard deviation. The natural frequencies of the phase
oscillators are completely entrained to that of the stimulus in a forced
entrained state (see Fig. 5(g)-5(i)), and hence $s_{m}=1$ in all the $M$ bins
and as a consequence this state is characterized by $S=0$.

States that are entrained to the external stimulus such as forced entrained
and bump frequency cluster states are characterized by near zero mean of the
time-averaged frequencies and can be more clearly distinguished using the mean
of the time-averaged frequencies in each bin. In particular, the conventional
chimera state, where the strength of the incoherence takes a value $0<S<1$,
can be clearly distinguished using the mean of the time-averaged frequencies
in each bin instead of their standard deviation in Eq. (8). For this purpose,
we define a modified strength of incoherence using the mean of the time-
averaged frequencies in each bin as

$\displaystyle\hat{S}=1-\frac{\sum_{m=1}^{M}\hat{s}_{m}}{M},\qquad\hat{s}_{m}=\Theta(\delta-\bar{\omega}_{m}),$
(9)

where $\bar{\omega}_{m}$ in each bin for the two-cluster state has a finite
value and hence the corresponding $\hat{s}_{m}=0$. As a consequence
$\hat{S}=1$ for the two-cluster state, whereas $S=0$ for this state. For the
same reason, $\hat{S}=1$ for the chimera, incoherent, and itinerant chimera of
the first and second kind, and forced entrained states. Note that for the
forced entrainment states $S=\hat{S}=0$ because of the complete entrainment
with the stimulus. Bump and the bump frequency cluster states are
characterized by intermediate values between $0$ and $1$ for $S$ and $\hat{S}$
as a fraction of the phase oscillators are entrained with the stimulus while
the rest have a finite time-averaged frequency.

In the next section, we classify the dynamical states using the time-averaged
order parameters $\langle R_{1}\rangle$, $\langle R_{2}\rangle$, along with
the strengths of incoherence $S$ and $\hat{S}$ in one- and two-parameter
phases.

##  V Map of regimes

Transitions in parameter space among the different dynamical states are
depicted in Figs. 7(a) and 7(b) as a function of the strength of the external
stimulus for two different $\alpha=0.3$ and $\alpha=1.54$, respectively. The
degree of the strength of incoherence and that of the time-averaged order
parameters have been used to clearly distinguish the dynamical states and
their transitions. Two-cluster state (TC) exists in the range of
$f\in(0,0.38)$ (see Fig. 7(a)) and there is a transition from two-cluster to
itinerant chimera states of the first kind (IC1) which prevails in the range
$f\in[0.38,0.73)$. Upon increasing $f$ further, there is a transition from IC1
to forced entrained (FE) state via chimera (CHI), frequency cluster (FC), bump
frequency cluster (BFC), and bump states (BS). Chimera states exist in the
range $f\in[0.73,0.84)$, frequency cluster states in the range
$f\in[0.84,0.89)$, bump frequency cluster states in the range
$f\in[0.89,0.98)$ and $f\in[1.24,1.30)$, bump states in the range
$f\in[0.98,1.24)$ and forced entrainment states in the range $f\in[1.30,1.5)$.
The transition from incoherent to forced entrainment states via itinerant
chimera states of the first (IC1) and second kind (IC2) is depicted in Fig.
7(b) for $\alpha=1.54$.

Figure 8: Two-parameter phase diagram in the ($\alpha$, $f$) parameter space
depicting various collective dynamical states of the adaptive network with
external periodic forcing obtained from one initial condition. The dynamical
states in the parameter space indicated by TC, ICS, IC1, CHI, FC, BS, BFC,
IC2, FE corresponds to two-cluster, incoherent state, itinerant chimera of the
first kind, chimera, frequency cluster, bump state, bump frequency cluster,
itinerant chimera of the second kind and forced entrained state, respectively.
The dashed line corresponds to the condition $f\geq\lambda$ with $\lambda=1$
for the forced entrained state.

Complementing the analysis in Fig. 7, the dynamical scenarios in the
$(\alpha,f)$ parameter space are shown in Fig. 8. Here the initial conditions
are distributed such that for all $i,j=1\dots,N$ the coupling weights $k_{ij}$
are uniformly distributed between $-1$ to $1$, while the phases are uniformly
distributed between $[0,2\pi)$ in the entire explored range of $\alpha$ and
$f$. For low values of $f$, that is approximately $f\in(0,0.38)$, there is a
transition from TC state to incoherent state (ICS) as a function of $\alpha$.
An itinerant chimera of first kind is observed in the entire explored range of
$\alpha$ in a narrow range of $f\in(0.38,0.56)$. Moreover, the itinerant
chimera of first kind re-emerges in the range of $f\in(0.56,0.96)$ via the
chimera state (CHI) for larger values of $\alpha$. The spread of the bump
state (BS) increases as a function of $\alpha$ and $f$ until $\alpha=0.42$,
which then decreases until $f=1.20$. In a narrow range of $f\in(0.75,0.79)$,
there is a transition from the bump state to an itinerant chimera of first
kind via a chimera state for increasing $\alpha$. In the range of
$\alpha\in(0.11,0.81)$ there is a very fine region of frequency cluster states
separating the chimera state and bump frequency cluster state. Then there is a
transition from a bump state to an itinerant chimera of first kind via bump
frequency and chimera state, which is then followed by a transition from a
bump state to an itinerant chimera of second kind via a bump frequency cluster
and chimera state as a function of $\alpha$ in the appropriate ranges of $f$.
There is also a transition from a bump state to forced entrainment state as a
function of $\alpha$ for larger values of $f$. The forced entrainment state
prevails in the entire parameter space for $f>1.4$. The dashed line
corresponds to the stability condition of the forced entrained state above
which this state is stable, and it can be deduced as shown in the following.

Figure 9: Same as Fig. 8 but for different initial conditions.

In the forced entrained state, all the oscillators are entrained to the same
phase $\theta_{i}=\theta^{*}$ given by

$\displaystyle\theta^{*}=\sin^{-1}\left(\frac{\sum_{j}k_{ij}\sin(\alpha)-N\lambda}{Nf}\right),$
(10)

which possesses two solutions only if the corresponding coupling weights
fulfill the condition that $\frac{1}{N}\sum^{N}_{j=1}k^{*}_{ij}=\eta$ is
independent of $i$ and further the condition $|\eta\sin(\alpha)/N-\lambda|\leq
f$ is met. With $\sin^{-1}:[-1,1]\to[-\pi/2,\pi/2]$ the two solutions are
given by $\theta^{*}_{1}=\theta^{*}$ and $\theta^{*}_{2}=\pi-\theta^{*}$. Note
that for the forced entrained state, we have $\dot{k}_{ij}=0$ for all $i,j$.
The stability of the fixed point can be deduced using linear stability
analysis. The diagonal components ($DF_{ii}$) of the Jacobian matrix $J$ are

$\displaystyle DF_{ii}$ $\displaystyle=$
$\displaystyle\left(\frac{k_{ii}}{N}-\eta\right)\cos\alpha+f\cos(\theta_{i}^{*})$

The off-diagonal components ($DF_{ij}$) of the Jacobian matrix $J$ are

$\displaystyle DF_{ij}$ $\displaystyle=$
$\displaystyle\frac{k^{*}_{ij}}{N}\cos\alpha.$

The characteristic equation of the Jacobian matrix $J$ is given by

$\displaystyle(\mu+\eta\cos\alpha-f\cos(\theta^{*}))^{N}=0,$

which leads to the N-degenerate eigenvalues

$\displaystyle\mu_{k}=-\eta\cos\alpha\pm
f\sqrt{1-\bigg{(}\frac{\eta\sin\alpha-\lambda}{f}\bigg{)}^{2}},$

for $k=0,\dots,N-1$ and “+” and “-” corresponding to the solutions
$\theta^{*}_{1}$ and $\theta^{*}_{2}$, respectively. The forced locked state
emerge via a fold bifurcation at

$\displaystyle\sqrt{\bigg{(}\eta\cos\alpha\bigg{)}^{2}+\bigg{(}\eta\sin\alpha-\lambda\bigg{)}^{2}}=\pm
f.$ (11)

Note that for a suitable choice of initial conditions one may observe
$\frac{1}{N}\sum_{j}k^{*}_{ij}=\eta=0$. Hence, the condition for the existence
of a stable forced entrained state simplifies to $|\lambda|\leq f$. To show
the validity of the derived condition, the dynamical scenarios in the
$(\alpha,f)$ parameter space for uniformly distributed initial conditions that
lead to $\eta=0$ are explored and the results are depicted in Fig. 9. The
dynamical transitions are almost similar to those in Fig. 8 except for an
increase or decrease in the spread of the observed dynamical states. It is to
be noted that the stability condition given by $f\geq|\lambda|$ with
$\lambda=1$ (indicated by the dashed line), above which the forced entrained
state is stable, exactly matches with the simulation results for the chosen
set of initial conditions, in accordance with the condition derived from the
stability analysis.

The multistability between the observed dynamical states is depicted in the
two-parameter phase diagram in the ($\alpha$, $f$) parameter space in Fig. 10
for 10 different initial conditions. The dynamical states are the same as in
Fig. 8. Distinct dynamical states are represented as shades and patterns as
shown in Fig. 10. It is evident from the figure that the adaptive network with
external periodic forcing exhibits rich multistability.

Figure 10: Two-parameter phase diagram in the ($\alpha$, $f$) parameter space
depicting various multistable states of the adaptive network with external
periodic forcing. The dynamical states are the same as in Fig. 8, illustrating
a wealth of multistability among the distinct dynamical states for $10$
different initial conditions.

##  VI Conclusion

Using the order parameters $R_{1}$ and $R_{2}$ along with the space-time
plots, snapshots of instantaneous phases, and mean frequencies of the
collective dynamical states of an adaptive network under the influence of an
external periodic forcing, we have unraveled exotic types of chimera state,
namely, itinerant chimeras including bump states and bump frequency cluster
states. The population density of the oscillators constituting the coherent
and incoherent groups of the itinerant chimera co-evolves along with the
connection weights, which in turn co-evolve locally, and contribute to the
diversity of the frequency distribution of the resulting collective dynamical
state. In addition, the adaptive network is also found to admit
synchronization patterns such as phase and frequency clusters, forced
entrained state including the incoherent nature of the phase dynamics. We have
also defined two versions of the strength of incoherence to quantify and
corroborate the distinct dynamical states and their transitions. An analytical
stability condition for the forced entrained state is also found to agree well
with the simulation results for appropriate initial conditions of the coupling
weight matrix. In view of the recent profound relevance of chimera states in
complex neuronal patterns rmjw2010 ; lmmays1977 ; mtjwb2016 ; pjuws2006 ;
smbkb2019 ; es2020 , we speculate that the observed chimera states may provide
clues about the underlying dynamical mechanism of several neuronal and
neuropathological states because of their dynamic frequency distribution
depending on the co-evolving connection weights and phases. Further, as a
major motivation for the employed model hails from neuroscience, it is
important to investigate the emergence of the observed dynamical patterns in
relevant biological models and to explore their role in the self-organization
of adaptive dynamical patterns of neuronal assemblies.

###### Acknowledgements.

The work of V.K.C. is supported by the SERB-DST-MATRICS Grant No.
MTR/2018/000676 and DST-CRG Project under Grant No. CRG/2020/004353. The work
of V. K. C. is also supported by DST, New Delhi for computational facilities
under the DST-FIST programme (SR/FST/PS-1/2020/135)to the Department of
Physics. RB and ES acknowledge support by DFG (Project No. 440145547). The
work of MS forms a part of a research project sponsored by Council of
Scientific and Industrial Research (CSIR) under the Grant No. 03(1397)/17/EMR-
II. DVS is supported by the DST-SERB-CRG Project under Grant No.
CRG/2021/000816.

## Appendix A: Evolution of coupling weights

Coupling weights of distinct dynamical states at different time instance are
depicted in Figs. A1 \- A3. Evolution of coupling weights for two cluster,
chimera and incoherent states are depicted in Fig. A1. Evolution of coupling
weights for itinerant chimera of the first kind, bump state and itinerant
chimera of the second kind are shown in Fig. A2, while that of frequency
cluster, bump frequency cluster state and forced entrained state are
illustrated in Figs A3. Note that except for the entrained states such as two-
cluster state, frequency cluster and frequency entrainment state (where the
structure of the coupling matrix does not change qualitatively), the inter-
cluster coupling weights of all other observed states evolve in time.

Figure A1: Coupling weight matrix at different times. First row represents the
two-cluster state for $\alpha=0.40$ and $f=0.00$, second row represents the
chimera state $\alpha=0.40$ and $f=0.78$ and third row represents the
incoherent state for $\alpha=1.54$ and $f=0.20$. Figure A2: Coupling weight
matrix at different times. First row represents the itinerant chimera of the
first kind for $\alpha=1.54$ and $f=0.78$, second row represents the bump
state $\alpha=0.10$ and $f=0.90$ and third row represents the itinerant
chimera of the second kind for $\alpha=1.54$ and $f=0.99$. Figure A3: Coupling
weight matrix at different times. First row represents the frequency cluster
state for $\alpha=0.40$ and $f=0.83$, second row represents bump frequency
cluster state $\alpha=0.60$ and $f=0.95$ and third row represents the forced
entrainment state for $\alpha=0.40$ and $f=1.40$.

## Appendix B: Estimation of the strength of incoherence

Spatiotemporal plots of the chimera state for three different values of the
strength of the external forcing $f=0.73,0.77$ and $0.81$ are depicted in the
first row of Fig. A4. Other parameter are fixed as $\alpha=0.4$ and
$\varepsilon=0.005$. Time averaged frequencies of the oscillators, averaged
over the time interval $t\in(24000,26500)$, are depicted in the second row of
Fig. A4. Both the spatiotemporal plot and the time averaged frequencies
clearly illustrate the co-existing coherent and incoherent groups of the
chimera states. Now, $N=100$ oscillators, indexed in the same order as in the
spatiotemporal and the time averaged frequencies plots, are divided into
$m=20$ bins with $n=5$ oscillators in each bin. The standard deviation
$\sigma_{m}$ for each bin is calculated using (4).
$s_{m}=\Theta(\delta-\sigma_{m})$ in each bin depicted in the third row of
Fig. A4. Since, $s_{m}=1$ for the coherent groups and $s_{m}=0$ for incoherent
groups and consequently, the strength of incoherence take values $S=0$ and $1$
for coherent and incoherent groups, respectively, which is estimated using
(8). The strength of incoherence for the chimera states in Fig. A4 are also
included in the bottom row.

Figure A4: Chimera states for different $f$. First row depicts the time
evolution of phases in the range $t\in(26400,26500)$. Second and third row
depicts the time averaged frequencies in the range $t\in(24000,26500)$ and the
strength of incoherence, respectively. (a)-(c) $f=0.73$, (d)-(f) $f=0.77$ and
(g)-(i) $f=0.81$. Other parameter values are fixed as $\alpha=0.4$ and
$\varepsilon=0.005$.

## Appendix C: Chimeras for smaller $N$

Chimera states for $N=6,12,30$ and $50$ are depicted in Fig. A5(a)- A5(d),
respectively, to illustrate the emergence of chimera states even for smaller
$N$. Other parameter values are fixed as $\alpha=0.4,f=0.73$ and
$\varepsilon=0.005$.

Figure A5: Chimera states for different $N$. (a) $N=6$, (b) $N=12$, (c) $N=30$
and (d) $N=50$. Other parameter values are fixed as $\alpha=0.4,f=0.73$ and
$\varepsilon=0.005$.

## References

  * (1) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D. -U. Hwang, Phys. Rep. 424, 175-308 (2006).
  * (2) A. Arenas, A. Diaz-Guilera, J. Kurths, Y. MAreno, and C. Zhou, Phys. Rep. 469, 93-153 (2008).
  * (3) T. Gross, and B. Blasius, J. R. Soc. Interface 5, 259 (2008).
  * (4) Tamás Nepusz and Tamás Vicsek, Nature Physics 8, 568 (2012).
  * (5) R. Berner, Patterns of Synchrony in Complex Networks of Adaptively Coupled Oscillators, Springer Theses, Springer, Cham (2021).
  * (6) H. Markram, J. Lübke, and B. Sakmann, Science 275 213 (1997); L. F. Abbott, and S. Nelson, Nat. Neurosci. 3 1178 (2000); N. Caporale, and Y. Dan, Ann. Rev. Neurosci. 31 25 (2008).
  * (7) S. Jain, and S. Krishna, Proc. Natl. Acad. Sci. 98, 543 (2001).
  * (8) T. Gross, C. J. D. D’ Lima, and B. Blasius, Phys. Rev. Lett. 96, 208701 (2006); T. Gross, and B. Blasius, J. R. Soc. Interface 5 259 (2008).
  * (9) F. Sorrentino, and E. Ott, Phys. Rev. Lett. 100, 114101 (2008); F. Vazquez, V. M. Eguiluz, and M. S. Miguel, Phys. Rev. Lett. 100 108702 (2008).
  * (10) L. Horstmeyer, C. Kuehn, and S. Thurner, Phys. Rev. E 98, 042313 (2018).
  * (11) T. Aoki, and T. Aoyagi, Phys. Rev. Lett. 102, 034101 (2009); Phys. Rev. E 84, 066109 (2011).
  * (12) C. B. Picallo and H. Riecke, Phys. Rev. E 83, 036206 (2011).
  * (13) R. Gutiérrez, A. Amann, S. Assenza, J. Gómez-Gardeñes, V. Latora, and S. Boccaletti, Phys. Rev. Lett. 107, 234103 (2011).
  * (14) L. Timms and L.Q. English,, Phys. Rev. E 89, 032906 (2014).
  * (15) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Phys. Rev. Lett. 114, 038701 (2015).
  * (16) L. Lücken, O. V. Popovych, P. Tass, and S. Yanchuk, Phys. Rev. E 93, 032210 (2016).
  * (17) S. Chakravartula, P. Indic, B. Sundaram, and T. Killingback, PLoS ONE 12, e0178975 (2017).
  * (18) D. V. Kasatkin, S. Yanchuk, E. Schöll, and V. I. Nekorkin, Phys. Rev. E 96, 062211 (2017).
  * (19) D. V. Kasatkin, and V. I. Nekorkin, Quant. Electron. 58, 877 (2016); Chaos 28, 093115 (2018); Eur. Phys. J. Spec. Top. 227, 1051 (2018).
  * (20) V. Avalos-Gaytán, J.A.Almendral, I. Leyva, F. Battiston, V.Nicosia, V. Latora, and S. Boccaletti, Phys. Rev. E 97, 042301 (2018).
  * (21) M. H. Matheny et al, Science. 363, 1057 (2019).
  * (22) V. Röhr, R. Berner, E. L. Lameu, O. V. Popovych, and S. Yanchuk, PLoS ONE 14, e0225094 (2019).
  * (23) R. Berner, E. Schöll, and S. Yanchuk, SIAM J. Appl. Dyn. Syst 18, 2227 (2019).
  * (24) R. Berner, J. Fialkowski, D. V. Kasatkin, V. I. Nekorkin, S. Yanchuk, and E. Schöll, Chaos 29, 103134 (2019).
  * (25) R. Berner, J. Sawicki, and E. Schöll, Phys. Rev. Lett. 124, 088301 (2020).
  * (26) G. Pfurtscheller and C. Neuper, Neurosci. Lett. Suppl. 174, 93 (1994); C. M. Krause, H. Lang, M.. Laine, and B. Pörn, Electroencephalogr. Clin. Neurophysiol. 98, 319 (1996); L. Leocani, C. Toro, P. Manganotti, P. Zhuang, and M. Hallet, Electroencephalogr. Clin. Neurophysiol. 104, 199 (1997).
  * (27) M.C. Moore-Ede, F.M. Sulzman, and C.A. Fuller, The Clocks That Time Us (Harvard University Press, Cambridge, 1982).
  * (28) J.C. Dunlap, J.J. Loros, and P.J. Decoursey, Chronobiology: Biological Timekeeping (Sinauer Associates, 2003).
  * (29) R.G. Foster and L. Kreitzman, Rhythms of Life: The Biological Clocks that Control the Daily Lives of Every Living Thing (Yale University Press, New Haven, 2005).
  * (30) I. Z. Kiss, Y. Zhai, J. L. Hudson, Science 296, 1676 (2002).
  * (31) I. Z. Kiss, Y. Zhai, J. L. Hudson, Physical Review E 77, 046204 (2008).
  * (32) A. A. Zaikin, J. Kurths, and L. Schimansky-Geier, Phys. Rev. Lett. 85, 227 (2000).
  * (33) M. Rusconi, A. Zaikin, N. Marwan, and J. Kurths, Phys. Rev. Lett. 100, 128101 (2008).
  * (34) E. I. Volkov, E. Ullner, A. A. Zaikin and J. Kurths, Phys. Rev. E 68, 026214 (2003).
  * (35) A. S. Pikovsky and J. Kurths, Phys. Rev. Lett. 78, 775 (1997).
  * (36) C. Zhou, J. Kurths, and B. Hu, Phys. Rev. Lett. 87, 098101 (2001).
  * (37) S. Shinomoto, and Y. Kuramoto, Progress of Theoretical Physics 75, 1105 (1986).
  * (38) J. A. Acebron, L. L. Bonilla, C. J. Perez Vincente, F. Ritort, and R. Spigler Rev. Mod. Physics 77, 1-49 (2005).
  * (39) F. A. Rodrigues, Thomas K. DM. Peron, Peng Ji, and J. Kurths, Phys. Rep. 610, 1-98 (2016).
  * (40) B. Lindner, J. Garc a-Ojalvo, A. B. Neiman, L. and Schimansky-Geier, Phys. Rep. 392, 321 (2004).
  * (41) M. Madadi Asl, A. Valizadeh, and P. A. Tass, , Front. Physiol. 9, 1849 (2018).
  * (42) G. Q. Bi and M. M. Poo, J. Neurosci. 18, 10464 (1998).
  * (43) G. H. Paissan and D. H. Zanettee, Physica D 237 818 (2008).
  * (44) A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences (Cambridge Univer- sity Press, Cambridge, UK, 2001).
  * (45) I. Franović, S. Yanchuk, S. Eydam, I. Bacic, and M. Wolfrum, Chaos 30, 8, 083109 (2020).
  * (46) F. P. Kemeth, S. W. Haugland, L. Schmidt, I. G. Kevrekidis, and K. Krischer, Chaos 26, 094815 (2016).
  * (47) O. E. Omel’chenko, and E. Knobloch, New J. Phys. 21, 093034 (2019).
  * (48) F. Parastesh, S. Jafari, H. Azarnoush, Z. Shahriari, Z. Wang, S. Boccaletti and M. Perc, Phys. Rep. 898, 1-144 (2021).
  * (49) D. V. Kasatkin, V. V. Klinshov, and V. I. Nekorkin, Phys. Rev. E 99, 022203 (2019).
  * (50) M. Owen, C. Laing, and S. Coombes, New J. Phys. 9, 378 (2007).
  * (51) C. R. Laing, and O. E. Omelćhenko, Chaos 30, 043117 (2020).
  * (52) H. Schmidt, and D. Avitabile, Chaos 30, 033133 (2020).
  * (53) I. Franović, O. E. Omelćhenko, and M. Wolfrum, (arXiv:2104.03025v1, 2021).
  * (54) M. Gerster, R. Berner, J. Sawicki, A. Zakharova, A.S koch, J. Hlinka, K. Lehnertz, and E. Schöll, Chaos 30, 123130 (2020).
  * (55) P. J. Uhlhaas, and W. Singer, Neuron 52, 155 (2006).
  * (56) R. Gopal, V. K. Chandrasekar, A. Venkatesan, and M. Lakshmanan, Phys. Rev. E 89, 052914 (2014).
  * (57) Note that one needs to be very careful about the choice of $M$.
  * (58) R. Ma, J. Wang, and Z. Liu, Europhys. Lett. 91, 40060 (2010).
  * (59) L. M. Mukhametov, A. Y. Supin, and I. G. Polyakova, Brain Res. 134, 581 (1977).
  * (60) M. Tamaki, J. W. Bang, T. Watanabe, and Y. Sasaki, Curr. Biol. 26, 1190 (2016).
  * (61) S. Majhi, B. K. Bera, D. Ghosh, and M. Perc, Phys. Life Rev. 28, 100 (2019).
  * (62) E. Schöll, Chimeras in Physics and Biology: Synchronization and Desynchronization of Rhythms, Nova Acta Leopoldina 425, 67 (2020).


For CivilComments, the identity annotation is a value between $(0, 1)$ (it is the average rating among multiple raters), and we include an example into the underrepresented group only if the rating $>$ 0.99 (i.e. all raters agree about the identity) following [Koh et al., 2021]. However, we do note that this leads to a under coverage of the group membership, as many comments with plausible identity mentions are not included into the group identity labels.

Model. For tabular experiments, we use a 2-layer Dense ResNet model with 128 hidden units and pre-activation dropout rate = 0.1, using a random-feature Gaussian process with hidden dimension 256 as the output layer [Liu et al., 2022] (In the preliminary experiments, we tried larger models with update to 6-layers and 1024 hidden units, and did not observe significant improvement). For language experiments, we used BERT$_{\mbox{small}}$ mode initialized from the official pre-trained checkpoint released at BERT GitHub page[Turc et al., 2019][<https://github.com/google-research/bert>]. In each active learning round, we train the Dense ResNet model with Adam optimizer with learning rate 0.1, batch size 256 and maximum epoch 200; and train the BERT model with AdamW optimizer (learning rate 1e-5) for 6 epochs with batch size 16.

Experiment Protocol. Boxes with thick outlines (Underrepresentation Label, Introspective Training Method, Sampling Signal) indicates the experiment components where the methods differ.

Active Learning Protocol.

<Ref> visualizes the experiment protocol. As shown, in each stage, we first (optionally) trains a cross validated ensemble to estimate the under-representation labels, where we split the data into 10 cross-validation splits, and train ensemble members on 1 split and predict the rest of the 9 splits. We then use the ensemble's in-sample and out-of-sample predictions to compute the underrepresentation label $\hat{b}_i$ (<Ref>), and conduct introspective training (<ref>) to generate the final active sampling signals for 8 rounds to generate the final sampled data (red box). At the end of round 8, we estimate the underrepresentation label for the final sampled data, and send it to the final model for reweighted training to generate the full accuracy-fairness frontier. The sampling model is always a 10-member ensemble of neural Gaussian process (introduced in <Ref>), and the final model is always a single DNN with architecture identical to the sampling model (i.e., 2-layer Dense ResNet for census income and BERT$_{small}$ for toxicity detection).

For both tasks, we randomly sample as small subset as the initial labelled dataset (2,500 out of 32,561 total training examples for census income, and 50,000 out of total 405,130 examples for toxicity detection), and use the rest of the training set as the unlabelled set for active learning. For each sampling round, the AL model acquires 1,500 examples for census income, and 15,000 examples for the toxicity detection, so the total sample reaches roughly half of the total training set size after 8 rounds.

In the final model training, we use the standard re-weighting objective [Liu et al., 2021]:

$$\sum_{(x,y)\not\in \hat{\gB}} L_{ce}(y, f(\rvx)) + \lambda \sum_{(x,y)\in \hat{\gB}} L_{ce}(y, f(\rvx))$$

where $\hat{\gB}$ is the set of underrepresented examples identified by the underrepresentation label, i.e., $(x_i,y_i) \in \hat{\gB}$ if $\hat{b}_i > t$. We vary the thresholds $t$ and the up-weight coefficient $\lambda$ over a 2D grid ($t \in \{0.05, 0.1, 0.15, ..., 1.0\}$ and $log(\lambda) \in \{0., 0.5, 1, 1.5, \dots, 10.\}$) to get a collection of model accuracy-fairness performances (i.e., accuracy v.s. worst-group accuracy), and use them to identify the Pareto frontier defined by this combination of data and reweighting signal.

Active Learning Signals.

In this work, we consider four types of active sampling signals. Recall that the sampling model (neural Gaussian process ensemble) is a K-member ensemble that generates three predictive quantities: (1) label probability $\{p_k(y|\rvx)\}_{k=1}^{10}$, (2) underrepresentation probability $\{p_k(b|\rvx)\}_{k=1}^{10}$ and (3) predictive variance $\{v(\rvx)\}_{k=1}^{10}$ (recall that $\E_k$ and $Var_k$ are the empirical mean and variance).

* Margin: The gap between the highest class probability and the second highest class probability for the output label. In the binary prediction context, this is equivalent to $2 * |p(y|\rvx) - 0.5|$, i.e., the gap between the mean predicted probability and the null value of 0.5. We use the mean predictive label probability of the ensemble, which leads to:

$$Margin(\rvx)=2 * |\E_k(p_k(y|\rvx)) - 0.5|.$$

* Predicted Underrepresentation: The mean predictive underrepresentation probability of the ensemble, which leads to:

$$Underrep(\rvx)= \E_k(p_k(b|\rvx)).$$

* Diversity: i.e., Ensemble Diversity (introduced in <Ref>). The variance of label predictions:

$$Diversity(\rvx)= Var_k(p_k(y|\rvx)).$$

* Variance: i.e., Predictive Variance (introduced in <Ref>). The mean of predictive variances:

$$Variance(\rvx)= \E_k(v_k(\rvx)).$$

§ NOISE-BIAS-VARIANCE DECOMPOSITION IN RIDGE REGRESSION

Consider fitting a ridge regression model $f(\rvx_i)=\beta^\top \rvx_i$ to the Gaussian observation data $y_i=\theta^{\top} \rvx_i + \epsilon, \epsilon \sim N(0, \sigma^2)$ under an imbalanced experiment design, where we have $|\gG|$ treatment groups and $n_g$ observations in each group. . Here, $\rvx_i = [1_{g_i = 1}, ..., 1_{g_i = |\gG|}]$ is a $|\gG| \times 1$ one-hot indicator of the membership of $g_i$ for each group in $\gG$, and $\theta=[\theta_1, \dots, \theta_{|\gG|}]$ is the true effect for each group. Then, under ridge regression, the noise-bias-variance decomposition for group $g$ is $E_D(L(y, f_D))={\color{OliveGreen} \sigma^2} + {\color{Maroon} \frac{(\lambda\theta_g)^{ 2}}{(n_g + \lambda)^2}} + {\color{RoyalBlue} \frac{\sigma^2 n_g}{(n_g + \lambda)^2}}$, where the regularization parameter $\lambda$  modulates a trade-off between the  bias and  variance terms. In <Ref>, we also treat the case of group-specific noise $\epsilon_i \stackrel{indep}{\sim} N(0, \sigma_g^2)$ .

§.§ Error Decomposition in a General Setting

We first derive the decomposition in a general setting with data $\{y_i, \phi_i\}_{i=1}^n$, where $\phi_i$ is the $d \times 1$ (fixed) features that follows a distribution $P(\phi)$. We consider a well-specified scenario where the data generation mechanism as:

$$ y_i = \tilde{y}_i + \epsilon, \quad \mbox{where} \quad \tilde{y}_i = \theta^\top \phi_i, \;\; \epsilon \stackrel{iid}{\sim} N(0, \sigma^2),$$

and $\theta_{d \times 1}=[\theta_1, \dots, \theta_{|\gG|}]$ is the true coefficient. Under ridge regression, we fit a linear model $f(\rvx_i)=\phi_i^\top\beta$ to the data by minimizing the following squared loss objective:

$$||\rvy_{n \times 1} - \Phi_{n \times d}\beta_{d \times 1} ||_2^2 + \lambda ||\beta||_2^2,$$

which gives rise to the following solution:

\begin{align}
\hat{\beta} = (\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top \rvy.
\label{eq:ridge_coef}
\end{align}

notice $\hat{\beta}$ is a random variable that depends on the data $\Phi_{n \times d} = [\phi_1^\top, \dots, \phi_n^\top] \stackrel{iid}{\sim} P(\phi)$. Notice that under squared loss, the ensemble predictors $\bar{f} = argmin_{f} E_\Phi[(f - \hat{\beta}^\top \rvx_i)^2]$ is simply the mean of individual predictors, i.e, $\bar{f} = E_\Phi(\hat{\beta}^\top \rvx_i)=\bar{\beta}^\top \rvx_i$, where $\bar{\beta}=E_\Phi(\hat{\beta})$.

Consequently, given a new observation $\{y, \phi\}$, the noise-bias-variance decomposition of $\hat{\beta}$ under squared loss is:

\begin{align}
E[(y - \Phi\hat{\beta})^2]
E_y[(y - \tilde{y})^2] + (\tilde{y} -  \bar{\beta}^\top \phi_i)^2 +
E_{\Phi}[\bar{\beta}^\top \phi_i  - \hat{\beta}^\top \phi_i ]^2
\nonumber    \\
\underbrace{\sigma^2\vphantom{\phi_i^\top[\theta -  \bar{\beta}]}}_{Noise} +
\underbrace{\phi_i^\top[\theta -  \bar{\beta}][\theta -  \bar{\beta}]^\top\phi_i}_{Bias} +
\underbrace{\phi_i^\top Var(\hat{\beta}) \phi_i}_{variance}.
\label{eq:error_decomposition_ridge_full}
\end{align}

As shown, to obtain a closed-form expression of the decomposition, we need to first derive the expressions of $Bias(\hat{\beta})=[\theta - \bar{\beta}]$ and $Var(\hat{\beta})$. Under the expression of the ridge predictor (<ref>), we have:

\begin{align*}
&= [\theta-\bar{\beta}] \\
&= \theta - E[(\Phi^\top \Phi + \lambda I_d)^{-1}\Phi^\top\Phi\theta]\\
&= E[I - (\Phi^\top \Phi + \lambda I_d)^{-1}\Phi^\top\Phi]\theta\\
&= \lambda * E[(\Phi^\top \Phi + \lambda I_d)^{-1}]\theta;\\
&= E(Var(\hat{\beta}|\Phi)) + Var(E(\hat{\beta}|\Phi)),
\end{align*}

\begin{align*}
&= E[(\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top Var(\rvy) \Phi(\Phi^\top\Phi + \lambda I_d)^{-1}]
\\
&= \sigma^2 E[(\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top\Phi(\Phi^\top\Phi + \lambda I_d)^{-1}] \\
&= \sigma^2 * E[(\Phi^\top\Phi + \lambda I_d)^{-1} - \lambda (\Phi^\top\Phi + \lambda I_d)^{-2}]. \\
&= E[S \theta\theta^\top S^\top] - E[S]\theta\theta^\top E[S^\top]
\end{align*}

where $S=(\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top\Phi$.

As shown, the above expression depends on the random-matrix moments
$E[(\Phi^\top \Phi + \lambda I_d)^{-1}]$,
$E[(\Phi^\top\Phi + \lambda I_d)^{-2}]$,
$E[(\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top\Phi]$ and
$E[S \theta\theta^\top S^\top]$.

§.§ Error Decomposition under Orthogonal Design

The above moments are in general difficult to solve due to the involvement of matrix inverse and product within the expectation. However, a closed-form expression is possible under an orthogonal design where $\phi_i = [1_{g_i = 1}, \dots, 1_{g_i = |\gG|}]$ is the one-hot vector of treatment group memberships. Then, denote $Diag[z_g]$ the diagonal matrix with diagonal elements $z_g$ and $[z_{g g'}]_{g g'}$ the full matrix whose $(g, g')$ element is $z_{g g'}$, we have:

\begin{align*}
\Phi^\top\Phi &= diag[n_g] \\
E[(\Phi^\top \Phi + \lambda I_d)^{-1}] &= diag[\frac{1}{n_g + \lambda}], \\
E[(\Phi^\top \Phi + \lambda I_d)^{-2}] &= diag[\frac{1}{(n_g + \lambda)^2}], \\
E[(\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top\Phi] &= diag[\frac{n_g}{n_g + \lambda}],
\end{align*}

\begin{align*}
E[(\Phi^\top\Phi + \lambda I_d)^{-1}\Phi^\top\Phi \theta\theta^\top \Phi^\top\Phi(\Phi^\top\Phi + \lambda I_d)^{-1}] &=
[\frac{n_g n'_g}{(n_g + \lambda)(n'_g + \lambda)} \theta_g \theta_{g'}]_{g g'}.
\end{align*}

We are now ready to derive the full decomposition (<ref>), without loss of generality, we assume $\phi_i$ belongs to group $g$. Then:

\begin{align*}
\phi_i^\top Bias(\hat{\beta})
&= \lambda * \phi_i^\top E[(\Phi^\top \Phi + \lambda I_d)^{-1}]\theta
= \frac{\lambda}{n_g + \lambda}\theta_g;\\
\phi_i^\top Var(\hat{\beta}) \phi_i
&= \sigma^2 * \phi_i^\top E[(\Phi^\top\Phi + \lambda I_d)^{-1} - \lambda (\Phi^\top\Phi + \lambda I_d)^{-2}] \phi_i; \\
&= \frac{\sigma^2}{n_g + \lambda} - \frac{\lambda \sigma^2}{(n_g + \lambda)^2}
= \frac{\sigma^2 n_g}{(n_g + \lambda)^2}.
\end{align*}

Consequently, we have the noise-bias-variance decomposition in (<ref>) as:

||ϕ_i^⊤Bias(β̂)||_2^2       = (λθ_g)^2/(n_g + λ)^2;

ϕ_i^⊤Var(β̂) ϕ_i       = σ^2 n_g/(n_g + λ)^2.

§.§ Error Decomposition under Orthogonal Design and Heterogeneous Noise

We now consider the case where $y_i \sim N(\theta^\top\phi_i, \sigma_g^2)$ follows a normal distribution with group-specific noise. Using the same decomposition as in <ref>, we see that:

\begin{align}
E[(y - \Phi\hat{\beta})^2]
E_y[(y - \tilde{y})^2] + (\tilde{y} -  \bar{\beta}^\top \phi_i)^2 +
E_{\Phi}[\bar{\beta}^\top \phi_i  - \hat{\beta}^\top \phi_i ]^2
\nonumber    \\
\underbrace{\sigma_g^2\vphantom{\phi_i^\top[\theta -  \bar{\beta}]}}_{Noise} +
\underbrace{\phi_i^\top[\theta -  \bar{\beta}][\theta -  \bar{\beta}]^\top\phi_i}_{Bias} +
\underbrace{\phi_i^\top Var(\hat{\beta}) \phi_i}_{variance}.
\label{eq:error_decomposition_ridge_full}
\end{align}

As shown, the nature of the bias and variance decomposition in fact does not change, and the noise component is now the group-specific variance $\sigma_g^2$. Therefore, by following the same derivation as in <Ref>, we have:

||ϕ_i^⊤Bias(β̂)||_2^2       = (λθ_g)^2/(n_g + λ)^2;

ϕ_i^⊤Var(β̂) ϕ_i       = σ^2 n_g/(n_g + λ)^2.

§ PROOF OF <REF>

Through introspective training, there is a guarantee on a model's bias-awareness based on its hidden representation and uncertainty estimates. At convergence, a well-trained model $f = (f_y, f_b)$ should satisfy the property that $p(b=1 | x) = \sigma(f_b(\rvx))$.

(I) (Bias-aware Hidden Representation)  We denote the odds for $\rvx$ belonging to the underrepresented group $\gB$ as $o_b(\rvx) = p(\rvx|b=1)/p(\rvx|b=0)$. Using Bayes' theorem, we derive the following:

\begin{align}
\nonumber
p(b|\rvx) &= \sigma(\beta^T h(\rvx) + \beta_0)
\\\nonumber
log \, \frac{p(b=1|\rvx)}{p(b=0|\rvx)} &= \beta^T h(\rvx) + \beta_0
\\\nonumber
log \, \frac{p(\rvx|b=1)p(b=1)}{p(\rvx|b=0)p(b=0)} &= \beta^T h(\rvx) + \beta_0
\\\nonumber
\beta^T h(\rvx) + \beta_0 &= log \, P(\rvx|b=1) - log P(\rvx|b=0) + log\frac{p(b=1)}{p(b=0)}
\\
\beta^T h(\rvx) + \beta_0 &= log \, o_b(\rvx) + log\frac{p(b=1)}{p(b=0)}
\label{eq:prop1_part1}
\end{align}

Hence, the hidden representation is aware of the likelihood ratio of whether an example $\rvx$ belongs to the underrepresented group, and the last-layer bias $\beta_0$ corresponds to the marginal likelihood ratio of the prevalence of the underrepresented groups $p(b=1)/p(b=0)$ .

(II) (Bias-aware Embedding Distance)
Next, we examine the embedding distance between two examples $(\rvx_1, \rvx_2)$, i.e., $||h(\rvx_1) - h(\rvx_2)||_2$.

The Cauchy-Schwarz inequality states that for two vectors $\rvu$ and $\rvv$ of the Euclidean space, $|\langle \rvu, \rvv \rangle| \leq ||\rvu|| \, ||\rvv||$. Hence, the distance between two embeddings can be expressed as $ \beta^T [h(\rvx_1) - h(\rvx_2)] \leq ||\beta||_2 \, ||h(\rvx_1) - h(\rvx_2)||_2$. Using this property and <Ref>, we derive the following:

\begin{align}
\nonumber
\beta^T h(\rvx_1) - \beta^T h(\rvx_2) &= log \, o_b(\rvx_1) - log \, o_b(\rvx_2) \\
\nonumber
\beta^T [h(\rvx_1) - h(\rvx_2)] &= log \, o_b(\rvx_1) - log \, o_b(\rvx_2) \\
\nonumber
log \, o_b(\rvx_1) - log \, o_b(\rvx_2) & \leq ||\beta||_2 \, ||h(\rvx_1) - h(\rvx_2)||_2 \\
\nonumber
\frac{1}{||\beta||_2} [log \, o_b(\rvx_1) - log \, o_b(\rvx_2)] & \leq ||h(\rvx_1) - h(\rvx_2)||_2 \\
\frac{1}{||\beta||_2} log \frac{o_b(\rvx_1)}{o_b(\rvx_2)} & \leq ||h(\rvx_1) - h(\rvx_2)||_2
\end{align}

Since the above inequality is invariant to the relative position of $(\rvx_1, \rvx_2)$, we also have: $\frac{1}{||\beta||_2} log \frac{o_b(\rvx_2)}{o_b(\rvx_1)} \leq ||h(\rvx_1) - h(\rvx_2)||_2$, which implies:

$$||h(\rvx_1) - h(\rvx_2)||_2 \geq \frac{1}{||\beta||_2} * max(log \frac{o_b(\rvx_1)}{o_b(\rvx_2)}, log \frac{o_b(\rvx_2)}{o_b(\rvx_1)}).$$

As shown, the distance between the hidden embeddings $h(\rvx_1)$, $h(\rvx_2)$ is lower-bounded by the log-odds ratio that a given example is in the underrepresented group. With this guarantee on the model's learned embedding distance, we expect the hidden features to be more diverse than when trained on the main task alone, since it needs to sufficient features to distinguish the underrepresented-group examples from those of the majority in the hidden space.

§ PERFORMANCE GUARANTEE FOR LOSS-BASED TAIL-GROUP DETECTION

In this section, we derive a lower bound for the group detection performance based on the self-play estimator introduced in <Ref>. We derive the result in a general setting without further assumptions on the data distribution. The goal here is to provide a broadly applicable, and mathematically rigorous account of how various aspects of the data distribution and the model behavior impacts the performance of the subgroup detection procedure, rather than deriving the tightest possible performance guarantee for a specific  family of model class or data distribution. Our proof technique develops a novel Cantelli-type lower bound on the tail probabilities of the form $P(F(l)>q)$ (where $F$ is a CDF and $l$ a random variable), and also a novel tight upper bound on $Var(F(l))$, which may be of independent interest.

Let's first establish some notations. Recall $P(\rvx, y, g)$ is the data-generating distribution (<Ref>). For the majority group $g=0$ and minority-group $g=1$, we denote $P_0(.)=P(.|g=0)$ and $P_1(.)=P(.|g=1)$ the data distribution for the majority and minority group, and recall $\gamma_0=P(G=0)$ and $\gamma_1=1-\gamma_0$ the prevalence of these two groups. Furthermore:

Group-specific error and its distribution. Given a loss function $L(.,.)$ of Bergman divergence family and a fixed, unfair model $f_y$ that violates the sufficiency criteria [Arjovsky et al., 2019, Shui et al., 2022], denote $l_0=L(\tilde{y}_0, f_y(\rvx_0))$ and $l_1=L(\tilde{y}_1, f_y(\rvx_1))$ the cross-validation  generalization gaps for the majority and the minority examples $(\tilde{y}_0, x_0) \sim P_0$ and $(\tilde{y}_1, x_1) \sim P_1$, where $\tilde{y}_g$ represents the true label without label noise. Notice that $(l_0, l_1)$ are random variables due to the randomness in $(P_0, P_1)$. To this end, also denote $(\mu_0, \mu_1), (\sigma_0^2, \sigma_1^2)$ the means and variances of the group-specific losses $(l_0, l_1)$. Due to the $f_y$'s violation of the sufficiency criteria (i.e., $E_0(y|f_y(\rvx)=t) \neq E_1(y|f_y(\rvx)=t)$), we expect the model $f_y$'s cross-validation loss is systematically worse for the minority groups, i.e., there's a systematic difference in group-specific loss $\mu_1 - \mu_0 = d >0$. Finally, denote $F(l_i)=P(l < l_i)$ the CDF for the distribution of the loss $l=L(y, f_y(\rvx)), (y, \rvx) \sim P$, and $(F_0, F_1)$ the CDF for the distributions of $l_0, l_1$, respectively.

Rank-based Estimator. We identify the minority group examples using a rank-based estimator. Specifically, recall $F(l_i)=P(l < l_i)$ is the CDF of the population loss, then the rank-based estimator for subgroup detection is:

\begin{align}
\hat{I}(g_i=1) = I(F(l_i) > q).
\label{eq:rank_estimator}
\end{align}

That is, we include a training example $\rvx_i$ into the introspective training only if the population quantile of its generalization error is higher than $q$, which is a user-specific threshold controlling the precision and recall of the estimator's identification performance.

Then, the below theorem describes how the estimator performance $P(g_i=1|F(l_i)>q)$ is related to the user-specified threshold $q$, and the characteristics of data distribution (i.e., the group prevalence $\gamma_0, \gamma_1$) as well as the classifier performance (in terms of the group-specific loss distribution $F_0, F_1$).

[Tail-group detection performance of rank-based estimator]
For $(l_0, l_1)$ a pair of random variables for minority- and majority-group examples. Denote $d=E(l_1 - l_0)>0$ and $\sigma^2=Var(l_1-l_0)$ the mean and variance of the between-group generalization gap.
Then, given a user-specified threshold $q \in (0, \frac{d^2}{d^2 + \sigma^2})$, the performance of the rank-based estimator $P(g=1|F(l)>q)$ is bounded by:

\begin{align}
P(g=1|F(l)>q) \geq (1-\gamma_0)^2 + \gamma_0 * \frac{1-\gamma_0}{1-q} * \frac{z^2}{z^2 + 1}
\quad \mbox{where} \quad
z=\frac{E[F_0(l_1)] - q}{\sqrt{Var[F_0(l_1)]}}.
\label{eq:rank_estimator_lower_bound_app}
\end{align}

Here, $F_0(l_1)=P(l_0 < l_1)$ is the majority-group loss CDF $F_0$ evaluated at the minority-group loss $l_1=L(\tilde{y}_1, f_y(\rvx_1))$ where $(\tilde{y}_1, \rvx_1) \sim P_1$.

Proof is at <Ref>.  Notice here due to the randomness in $l_1$, $F_0(l_1)$ is a random variable that follows a continuous distribution and is bounded within $F_0 \in [0, 1]$. Therefore $F_0(l_1)$ has valid moments $E[F_0(l_1)]$ and $Var[F_0(l_1)]$. We see that the performance bound (<ref>) is intuitively sensible: a lower majority group prevalence $\gamma_0$, a higher rank threshold $q$, and a higher  likelihood for loss dominance $F_0(l_1)=P(l_0 < l_1)$ all contributes a stronger identification performance. Here, we see that $z=\frac{E[F_0(l_1)] - q}{\sqrt{Var[F_0(l_1)]}}$ is a distribution dependent quantity that governs the difficulty of identifying the minority group.

Clearly, a larger magnitude of $z$ leads to a stronger guarantee in detection performance in Theorem <ref>. For the interested readers, the below result provides a lower bound for the magnitude of $z$ in terms of the characteristics of model performance (i.e., the moments of the group-specific loss distributions (e.g., $d=\mu_1 - \mu_2$ and $(\sigma^2_1, \sigma^2_2)$), which can be used to obtain a more precise understanding of $P(g=1|F(l)>q)$ in (<Ref>) in practice:

[Lower bound on the standardized likelihood of loss dominance $z$.]
Consider the standardized likelihood of loss dominance $z=\frac{E[F_0(l_1)] - q}{\sqrt{Var[F_0(l_1)]}}$, where the likelihood of loss dominance $F_0(l_1)=P(l_0 < l_1)$ is a random variable in terms of $l_1=L(\tilde{y}_1, f_y(\rvx_1)), (\tilde{y}_1, \rvx_1)\sim P_1$. We have:

\begin{align}
E[F_0(l_1)] &\geq \frac{d^2}{\sigma^2 + d^2} \quad \mbox{where} \quad
d=\mu_1 - \mu_0 > 0, \;
\sigma^2 = Var(l_1 - l_2)\leq \sigma_0^2 + \sigma_1^2,
\nonumber
\\
Var[F_0(l_1)] &\leq \frac{F_0(\mu_0)^2}{4} * \frac{\sigma_1^2}{\sigma_1^2 + d^2},
\label{eq:variance_lower_bound_app}
\end{align}

where $F_0(\mu_0)=P(l_0 \leq \mu_0)$ is the probability of the majority-group loss $l_0$ smaller than its mean $\mu_0 =E(l_0)$. Therefore, $z$ can be bounded by:

\begin{align}
z=\frac{E[F_0(l_1)] - q}{\sqrt{Var[F_0(l_1)]}}
\geq
\frac{2}{F_0(\mu_0)}
* \sqrt{\frac{\sigma_1^2 + d^2}{\sigma_1^2}} * (\frac{d^2}{(\sigma_1^2 + \sigma_2^2) + d^2} - q).
\label{eq:standardized_likelihood_lower_bound_app}
\end{align}

The proof is in <Ref>. It relies on a novel upper bound of CDF variance (i.e., <Ref>), which is important for guaranteeing a high magnitude of $z$ and consequently a tighter bound for group detection performance in <Ref>.

The combination of Theorems <ref>-<ref> provides us an opportunity to quantitatively understand of the group detection performance $P(g=1|F(l)>q)$ in terms of the estimator configuration (i.e., the user-specified threshold $q$), data distribution (i.e., the majority-group prevalence $\gamma_0$), and the degree of unfairness of the vanilla model $f_y$ (in terms of the distributions of the group-specific losses). For example, consider a setting with majority-group prevalence $\gamma_0=0.85$, expected  between-group loss gap $d=1$, group-specific variances $\sigma_1=\sigma_2=0.15$ and  $P(l_0 < \mu_0)=0.5$, using a percentile threshold $q=0.9$, the rank-based estimator $I(F(l)>q)$ has group detection performance $P(g=1|F(l)>q)>0.9$.

It is also worth commenting that, following the previous discussion, we see that the success of group identification relies on valid estimation of model's generalization error (with respect to the true label). To this end, the K-fold cross-validated ensemble procedure we used in this work is known to produce unbiased and low-variance estimates for the expected generalization loss [Blum et al., 1999, Kumar et al., 2013].

Finally, we highlight that the above results are broadly applicable and derived under weak, nearly assumption-free conditions. A even tighter performance bounds can be obtained by making further assumptions on the family of data distributions or model class, which is outside the scope of the current work. Another interesting direction is to incorporate the finite-sample estimation error of K-fold cross validation into the analysis [Blum et al., 1999, Kumar et al., 2013, Bayle et al., 2020], although this necessitates a careful treatment of model's generalization behavior (i.e., loss stability) in the subgroup setting, which we will pursue in the future work.

§ PROOF OF <REF>

Let $r\br{\alpha}=\br{\alpha n}^{-p}$ and $\ragg\br{n} = \tau n^{-q} + \delta$. As discussed in section <ref>, we assume that the group-specific risk decays as

\[E[R\br{\hat{f}_{\alpha, n}|G=g}] = c_g r\br{\alpha_g} + \ragg\br{n}\]

where $\hat{f}_{\alpha, n}$ is the minimizer of the $\gamma$-weighted empirical risk on a dataset of size $n$ with allocation $\alpha$ and the expectation is with respect to the randomness of the dataset used to train $\hat{f}$.

Under the above assumption, the optimization problem for $\alpha$ looks like

\[\min_{\alpha \in \Delta^{|\gG|}} \omega \br{\sum_g \gamma_g \br{c_g r\br{\alpha_g} + \ragg\br{n}}} + (1-\omega) \max_g \br{c_g r\br{\alpha_g} + \ragg\br{n}}\]

The term $\ragg\br{n}$ is a constant that does not impact the optimal solution. Hence, we drop it and focus on the problem

\[\min_{\alpha \in \Delta^{|\gG|}} \omega \br{\sum_g \gamma_g \br{c_g r\br{\alpha_g}}} + (1-\omega) \max_g \br{c_g r\br{\alpha_g}}\]

We begin by noting that the objective is strictly convex and the domain of optimization is bounded. Hence, there exists a unique global optimum.

Turning the constraint $\sum_{g \in \gG} \alpha_g = 1$ into a Lagrangian, we obtain

\[\min_{\alpha \geq 0} \omega \br{\sum_g \gamma_g \br{c_g r\br{\alpha_g}}} + (1-\omega) \max_g \br{c_g r\br{\alpha_g}} + \lambda\br{\sum_{g \in \gG} \alpha_g - 1}\]

where $\lambda \in \R$ is a Lagrange multiplier.

At the optimum, we know that the Lagrangian should have $0$ within its subdifferential wrt $\alpha$  [Boyd & Vandenberghe, 2004]. The subdifferential is given by

\[\left\{\begin{pmatrix}\br{\omega \gamma_1 + (1-\omega)\mu_1} c_1 r^\prime\br{\alpha_1} + \lambda \\ \br{\omega\gamma_2 + (1-\omega)\mu_2} c_2 r^\prime\br{\alpha_2} + \lambda\\ \vdots \\ \br{\omega\gamma_{|\gG|} + (1-\omega)\mu_{|\gG|}} c_{|\gG|} r^\prime\br{\alpha_{|\gG|}} + \lambda \end{pmatrix}  \text{ where } \mu \in \Delta^{|\gG|} \text{ is such that } \mu_g > 0 \iff g \in \argmax_{g^\prime \in \gG} c_{g^\prime} r\br{\alpha_{g^\prime}}\right\}\]

where $r^\prime$ denotes the derivative of $r$. Since $0$ belongs to the subdifferential at the optimum, there must exist $\mu$ satisfying the constraints above such that

\[\alpha^\star_g = \rpinv\br{-\frac{\lambda}{c_g\br{\omega \gamma_g + (1-\omega) \mu_g}}} \quad \forall g \in \gG\]

where $\rpinv$ is the inverse of $r^\prime$. Since $r^\prime$ is a homogeneuous function of its argument, so is its inverse, and $\lambda$ can be eliminated to enforce the constraint $\sum_g \alpha_g=1$, so that the optimal solution $\alpha^\star_g$ is

\[\alpha^\star_g = \frac{\rpinv\br{-\br{c_g\br{\omega \gamma_g + (1-\omega) \mu_g}}^{-1}}}{\sum_{g^\prime \in \gG} \rpinv\br{-\br{c_{g^\prime}\br{\omega \gamma_{g^\prime} + (1-\omega) \mu_{g^\prime}}}}}\]

Let $s$ denote the denominator and $\theta_g = \omega + (1-\omega) \frac{\mu_g}{\gamma_g}$ so that

\[\alpha^\star_g = \frac{\rpinv\br{-\br{c_g\gamma_g \theta_g}^{-1}}}{s}\]

and $\theta$ must satisfy:

\[\sum_{g \in \gG} \theta_g \gamma_g = 1, \theta_g \geq \omega \quad \forall g \in \gG, \theta_g > \omega\iff g \in \argmax_{g^\prime \in \gG} c_{g^\prime} r\br{\alpha_{g^\prime}^\star}\]

Using the fact that $r\br{t}=\br{nt}^{-p}$, we have that

\[\argmax_{g^\prime} c_{g^\prime} r\br{\alpha_{g^\prime}^\star} = \argmin_{g^\prime} c_{g^\prime}^{-\frac{1}{p}} \alpha_{g^\prime}^\star = \argmin_{g^\prime} c_{g^\prime}^{\frac{1}{p+1}-\frac{1}{p}} \br{\gamma_g \theta_g}^{\frac{1}{p+1}} = \argmin_{g^\prime} c_{g^\prime}^{-\frac{1}{p}} \gamma_g \theta_g\]

If we sort groups in ascending order according to $c_g^{-\frac{1}{p}} \gamma_g$ to obtain $g_1, g_2, \ldots$, a value of $\theta$ satisfying the conditions above can be computed as follows:

Initialize $\theta_g = \omega  \quad \forall g \in \gG$

Set $k=1$

Until $\sum_{g} \theta_g \gamma_g = 1$ or $k=|\gG|$ repeat:

\[
\begin{array}{ll}
\text{a)} &  l = \min(k+1, |\gG|) \\
\text{b)} & \text{ Set } \theta_{g_j}=\theta_{g_j} t \text{ (for } j \leq k \text{) for the largest $t \geq 1$ such that (i) } \sum_{g} \theta_g \gamma_g = 1 \text{ or (ii) } c_{g_1}^{-\frac{1}{p}}\gamma_{g_1}\theta_{g_1}= \ldots = c_{g_{l}}^{-\frac{1}{p}}\gamma_{g_{l}} \theta_{g_{l}} \\
\text{c)} & k=k+1
\end{array}
\]

It is easy to see that this algorithm must terminate as it can go for at most $|\gG|$ rounds. Further, by construction, we have $\theta_{g} \geq \omega \forall g \in \gG$ since we start at these values and only scale up any of the $\theta_g$.

At the $j$-th iteration of the loop,  we have that

\[c_{g_1}^{-\frac{1}{p}} \gamma_{g_1}\theta_{g_1}=\ldots=c_{g_j}^{-\frac{1}{p}} \gamma_{g_j}\theta_{g_j}\]

In the $j$-th iteration, we scale all $\theta_{g_j}$ ($m=1,\ldots,j+1$) up by the same factor $t$ by the same amount so that we achieve the above for all groups upto $g_{j+1}$, or we hit the constraint $\sum_g \theta_g\gamma_g=1$. If the former happens, we begin the next iteration with the same invariant. If the latter happens, we have obtained a $\theta$ that satisfies $\sum_g \theta_g \gamma_g=1$ and $\theta_{g_i} > \omega $ for $i < k$ and $g_i \in \argmin_{g^\prime} c_{g^\prime}^{-\frac{1}{p}} \gamma_g \theta_g$ for $i < k$.

If we reach the iteration where $k=|\gG|$, we can simply scale up all the $\theta_g$ by the same amount until $\sum_g \theta_g\gamma_g=1$ is satisfied and we would have that all $c_g^{-\frac{1}{p}} \gamma_g\theta_g$ is equal for all $g \in \gG$ and all $\theta_g > \omega$.

§ PROOF OF THEOREMS <REF> AND <REF>