d}\bm{u}\Biggr{]}=\\\
&=\sum_{m=0}^{\infty}q^{*}_{M}(m)\left(\int_{\mathbb{G}}\exp\left\\{\sum_{j=1}^{d}-s_{j}f_{j}(x)\right\\}P_{0}({\rm
d}x)\prod_{j=1}^{d}H^{*}_{j}({\rm d}s_{j})\right)^{m}\\\
&=\sum_{m=0}^{\infty}q^{*}_{M}(m)\left(\prod_{j=1}^{d}\int_{0}^{\infty}\int_{\mathbb{X}}\frac{\exp\left\\{-s_{j}\left(f_{j}(x)+u_{j}\right)\right\\}}{\psi_{j}(u_{j})}H_{j}({\rm
d}s_{j})P_{0}({\rm d}x)\right)^{m},\end{split}$ (S45)

where the second equality follows by definition of $H^{*}_{j}$.

The Laplace transform derived in Equation (S45) must coincide with the second
term of Equation (S43). To show this, we focus on the expected value in the
numerator of the second term in Equation (S43), since the denominator will
follow by setting $f_{j}=0$. We use Lemma 1 for $\Phi^{!}_{\bm{s},\bm{y}}$
with $f(\bm{w},x)=\sum_{j=1}^{d}w_{j}(u_{j}+f_{j}(x))$ to get

$\begin{split}&\mathbb{E}\Biggl{[}\exp\left\\{-\sum_{j=1}^{d}\int_{\mathbb{G}}w_{j}\left(u_{j}+f_{j}(x)\right)\Phi^{!}_{\bm{s},\bm{y}}({\rm
d}\bm{w},{\rm d}x)\right\\}\Biggr{]}\\\
&\quad=\sum_{m=0}^{\infty}q_{M}(m+k)\frac{(m+k)!}{m!}\frac{1}{\mathbb{E}\left[M^{(K)}\right]}\left(\prod_{j=1}^{d}\int_{0}^{\infty}\int_{\mathbb{X}}e^{-w_{j}(u_{j}+f_{j}(x))}H_{j}({\rm
d}w_{j})P_{0}({\rm d}x)\right)^{m}.\end{split}$ (S46)

Therefore, the denominator equals to

$\begin{split}\mathbb{E}\left[\exp\left\\{-\sum_{j=1}^{d}\int_{\mathbb{G}}w_{j}u_{j}\Phi^{!}_{\bm{s},\bm{y}}({\rm
d}\bm{w},{\rm
d}x)\right\\}\right]=\sum_{m=0}^{\infty}q_{M}(m+k)\frac{(m+k)!}{m!}\frac{1}{\mathbb{E}[M^{(K)}]}\left(\prod_{j=1}^{d}\psi_{j}(u_{j})\right)^{m}.\end{split}$
(S47)

To conclude, we multiply and divide the numerator by
$\prod_{j=1}^{d}\left(\psi_{j}(u_{j})\right)^{m}$ and then take the ratio
between equations (S46) and (S47), simplifying common terms. The identity with
Equation (S45) follows after recognizing the definition of $q^{*}_{M}(m)$,
which has been stated above. ∎

###  S2.8 Full conditional distribution of $\mathbf{U}_{n}$

As a byproduct of Theorem 3, we also obtain the full conditional distribution
of vector $\bm{U}_{n}$.

######  Corollary 3.

The conditional distribution of $\bm{U}_{n}$ given a realization of
$(\theta_{1},\dots,\theta_{d})$ in $K$ distinct values
$(\theta^{**}_{1},\dots,\theta^{**}_{K})$ with counts
$(\bm{n}_{1},\dots,\bm{n}_{d})$ depends on the distinct values only through
their counts and it has density on $(\mathbb{R}^{+})^{d}$ given by

$f_{\bm{U}_{n}}(\bm{u}\mid\bm{n}_{1},\dots,\bm{n}_{d})\propto\Psi(K,\bm{u})\prod_{j=1}^{d}\left\\{\frac{u_{j}^{n_{j}-1}}{\Gamma(n_{j})}\prod_{k=1}^{K}\kappa_{j}(u_{j},n_{j,k})\right\\}$

###### Proof.

The statement follows from Equation (S42), conditioning on $\bm{\theta}$ and
computing the expected value as done in Equation (S47). ∎

###  S2.9 Proof of Equation (12)

We refer to the case when $q_{M}$ is chosen to be a $1-$shifted Poisson
distribution with parameter $\Lambda$, that is

$\mathbb{P}\left(M=m\mid\Lambda\right)\ =\ q_{M}(m)\ =\
\frac{e^{-\Lambda}\Lambda^{m-1}}{(m-1)!},\quad\text{ for }m=1,2,\dots$

In such a case, the $\Psi(K,\bm{u})$ function, defined in Equation (S15) for
$K\geq 0$ and $\bm{u}\in(\mathbb{R})^{d}$ equals

$\begin{split}\Psi(K,\bm{u})&=\sum_{m=0}^{\infty}\frac{(m+K)!}{m!}\frac{e^{-\Lambda}\Lambda^{m+K-1}}{(m+K-1)!}\prod_{j=1}^{d}\left(\psi_{j}(u_{j})\right)^{m}\\\
&=e^{-\Lambda}\Lambda^{K-1}\sum_{m=0}^{\infty}\frac{(m+K)\Lambda^{m}}{m!}\bm{\psi}(\bm{u})^{m},\end{split}$

where $\bm{\psi}(\bm{u})=\prod_{j=1}^{d}(\psi_{j}(u_{j})$. We note that, in
the case under consideration, $\psi_{j}(u_{j})=1/(u_{j}+1)^{\gamma_{j}}<1$,
hence, the series does converge. As a consequence, we can split the series in
two parts

$\begin{split}\Psi(K,\bm{u})&=e^{-\Lambda}\Lambda^{K-1}\left(\sum_{m=1}^{\infty}\frac{\Lambda^{m}}{(m-1)!}\bm{\psi}(\bm{u})^{m}+K\sum_{m=0}^{\infty}\frac{\Lambda^{m}}{(m)!}\bm{\psi}(\bm{u})^{m}\right),\end{split}$

where we point out that the first sum starts from $m=1$ because we simplified
$m/m!$, which is exactly equal to $0$ when $m=0$. Then, we change of index in
the first sum ($h=m-1$)

$\begin{split}\Psi(K,\bm{u})&=e^{-\Lambda}\Lambda^{K-1}\left(\Lambda\bm{\psi}(\bm{u})\sum_{h=0}^{\infty}\frac{1}{h!}\left(\Lambda\bm{\psi}(\bm{u})\right)^{h}+K\sum_{m=0}^{\infty}\frac{\Lambda^{m}}{(m)!}\prod_{j=1}^{d}\left(\psi_{j}(u_{j})\right)^{m}\right).\end{split}$

The result follows after recognizing the exponential series

$\Psi(K,\bm{u})=e^{-\Lambda}\Lambda^{K-1}\left(\Lambda\bm{\psi}(\bm{u})e^{\Lambda\bm{\psi}(\bm{u})}+Ke^{\Lambda\bm{\psi}(\bm{u})}\right)=\Lambda^{K-1}\left(K+\Lambda\bm{\psi}(\bm{u})\right)e^{-\Lambda\left(1-\bm{\psi}(\bm{u})\right)}.$

###  S2.10 Proof of Theorem 4

###### Proof.

Consider a sequence of $n$ observations divided in $d=2$ groups of length
$n_{1}$ and $n_{2}$, respectively. We want to compute the probability of
having $K_{(n_{1},n_{2})}=K$ global clusters under a $\operatorname{Vec-
FDP}(q_{M},\bm{\gamma},\Lambda)$ prior.
Firstly, we remind the following formula (Charalambides, 2002, Theorem 8.16)

$\lvert C(n,K;-\gamma,-\rho)\rvert\ =\
\frac{1}{K!}\sum\binom{n}{r_{1},\ldots,r_{K}}\prod_{j=1}^{K}(\gamma)_{r_{j}}$
(S48)

where $(x)_{n}$ denotes the Pochhammer symbol, i.e.,
$(x)_{n}=\frac{\Gamma(x+n)}{\Gamma(x)}$. In general, this is defined for any
$n\in\mathbb{R}^{+}$ but, for positive and integer values of $n$, it reduces
to the rising factorial of $x$ of order $n$. Moreover, in Equation (S48) the
sum is extended over all the vectors $(r_{1},\ldots,r_{K})$ of positive
integers such that $\sum_{j=1}^{K}r_{j}=n$.
The probability under study $\mathbb{P}(K_{(n_{1},n_{2})}=K)$ can be written
as follows in terms of the pEPPF function,

$\mathbb{P}(K_{n_{1},n_{2}}=K)=\sum_{(\star)}\frac{1}{K!}\prod_{j=1}^{2}\binom{n_{j}}{n_{j,1},\ldots,n_{j,K}}\cdot\Pi_{K}^{(n)}(\bm{n}_{1},\bm{n}_{2})$

where the sum $(\star)$ is extended over all the vectors
$(n_{1,1},\ldots,n_{1,K})$ and $(n_{2,1},\ldots,n_{2,K})$ of non-negative
integers satisfying the following constraints

$n_{j,k}\geq
0,\;k\in{1,\ldots,K},\;j\in\\{1,2\\},\;\sum_{k=1}^{K}n_{j,k}=n_{j}\;j\in\\{1,2\\}\text{
and }n_{1,k}+n_{2,k}\geq 1\text{ as }k=1,\ldots,K.$

By exploiting the expression of the pEPPF we get

$\begin{split}\mathbb{P}(K_{(n_{1},n_{2})}=K)&=\sum_{\bar{M}=0}^{\infty}\frac{(\bar{M}+K)!}{\bar{M}!}q_{M}(\bar{M}+K)\prod_{j=1}^{2}\frac{1}{(\gamma_{j}(\bar{M}+K))_{n_{j}}}\\\
&\qquad\times\sum_{(\star)}\frac{1}{K!}\prod_{j=1}^{2}\binom{n_{j}}{n_{j,1},\ldots,n_{j,K}}\cdot\prod_{j=1}^{2}\prod_{k=1}^{K}(\gamma_{j})_{n_{j,k}}.\end{split}$
(S49)

We now focus on the evaluation of the sums over the partitions in (S49). This
may be calculated by counting how many zeros we have in the vectors
$(n_{j,1},\ldots,n_{j,K})$. We denote by $r_{1}$ (resp. $r_{2}$) the number of
possible zeros in the first vector (resp. second vector). Thanks to the
validity of the condition $n_{1,k}+n_{2,k}\geq 1$, we know that if $n_{1,k}$
then $n_{2,k}\geq 1$ and vice versa. Moreover we observe that one has
$\binom{K}{r_{1}}$ ways to choose the zeros in the first vector and
$\binom{K-r_{1}}{r_{2}}$ possibilities to choose the zeros in the second one.
Thus, we get

$\displaystyle\sum_{(\star)}\frac{1}{K!}\prod_{j=1}^{2}\binom{n_{j}}{n_{j,1},\ldots,n_{j,K}}\cdot\prod_{j=1}^{2}\prod_{k=1}^{K}(\gamma_{j})_{n_{j,k}}$
$\displaystyle\qquad=\sum_{r_{1}=0}^{K}\sum_{r_{2}=0}^{K-r_{1}}\binom{K}{r_{1}}\binom{K-r_{1}}{r_{2}}\sum_{(\star_{1})}\sum_{(\star_{2})}\frac{1}{K!}\prod_{j=1}^{2}\binom{n_{j}}{n_{j,1},\ldots,n_{j,K-r_{j}}}\cdot\prod_{j=1}^{2}\prod_{k=1}^{K-r_{j}}(\gamma_{j})_{n_{j,k}}$

where

$(\star_{j})=\left\\{(n_{j,1},\ldots,n_{j,K-r_{j}}):\;n_{j,k}\geq
1,k=1,\ldots,K-r_{j},\,\sum_{k=1}^{K-r_{j}}n_{j,k}=n_{j}\right\\}$

as $j=1,2$. We now exploit (S48) to obtain

$\begin{split}&\sum_{(\star)}\frac{1}{K!}\prod_{j=1}^{2}\binom{n_{j}}{n_{j,1},\ldots,n_{j,K}}\cdot\prod_{j=1}^{2}\prod_{k=1}^{K}(\gamma_{j})_{n_{j,k}}\\\
&\qquad\times\sum_{r_{1}=0}^{K}\sum_{r_{2}=0}^{K-r_{1}}\prod_{j=1}^{2}\lvert
C(n_{j},K-r_{j};-\gamma_{j})\rvert\end{split}$ (S50)

If we now replace (S50) in (S49), the statement of the theorem follows. ∎

##  S3 Predictive distributions

In this section, we provide and prove the formula for the predictive
distribution for a new client entering in a restaurant $j$. This has already
been stated in Equation (16) for the special case of a $\operatorname{Vec-
FDP}$ prior. Instead, in this section we refer to the more general
$\operatorname{Vec-NIFPP}$. More precisely, consider a realization
$\left(\bm{\theta}_{1},\dots,\bm{\theta}_{d}\right)$ from
$\left(P_{1},\dots,P_{d}\right)\sim\operatorname{Vec-NIFPP}(q_{M},H,P_{0})$
with $K$ distinct values $\left(\theta^{**}_{1},\dots,\theta^{**}_{K}\right)$
and a partition $\rho=\left\\{C_{1},\dots,C_{K}\right\\}$ with counts
$\left(\bm{n}_{1},\dots,\bm{n}_{d}\right)$ satisfying the constraints in
Equation (7). Following the approach of James et al. (2009), Favaro and Teh
(2013) and Argiento and De Iorio (2022) we work conditionally to
$\bm{U}_{n}=\bm{u}$. Then, for each group, say $j$, we have

$\begin{split}\mathbb{P}(\theta&{}_{j,n_{j}+1}\in\cdot\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d},\bm{u})\\\
&\propto\sum_{k=1}^{K}\frac{\kappa_{j}(u_{j},n_{j,k}+1)}{\kappa_{j}(u_{j},n_{j,k})}\delta_{\theta^{**}_{k}}(\cdot)\
+\ \kappa_{j}(u_{j},1)\left(\prod_{h\neq
j}\kappa_{h}(u_{h},0)\right)\frac{\Psi(K+1,\bm{u})}{\Psi(K,\bm{u})}P_{0}(\cdot)\end{split}$
(S51)

Without loss of generality, consider a new client entering the first
restaurant, where $n_{1}$ customers are currently seated. This $(n_{1}+1)$-th
client can $(i)$ choose an already existing table, or $(ii)$ sit to a non-
existing table. Note that, as detailed in Section 3.2, we distinguish between
empty and non-existing tables. The latter refer to tables that do not appear
in any restaurant of the franchise, while we admit the possibility of having
tables with no clients in some restaurant as long as at least one customer in
the whole franchise is seated in the corresponding table, see Figure S2.

As far as the proof of Equation S51 is concerned, considering case $(i)$: we
compute the probability of sitting to table $k$, $k=1,\dots,K$, where dish
$\theta^{**}_{k}$ is served. This probability equals

$\begin{split}\mathbb{P}\bigl{(}\theta_{1,n_{1}+1}=\theta^{**}_{k}\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d}\bigr{)}\
=\ \Pi^{(n+1)}_{K}(\widetilde{\bm{n}}_{1},\dots,\bm{n}_{d}),\end{split}$

where $\widetilde{\bm{n}}_{1}=(n_{1,1},\dots,n_{1,k}+1,\dots,n_{1,K})$. If we
further condition to $\bm{U}_{n}=u$, we obtain

$\begin{split}\mathbb{P}\bigl{(}\theta_{1,n_{1}+1}&=\theta^{**}_{k}\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d},\bm{u}\bigr{)}\\\
&\propto\Psi(\bm{u},K)\left(\prod_{j=2}^{K}\frac{(u_{j})^{n_{j}-1}}{\Gamma(n_{j})}\prod_{l=1}^{K}\kappa_{j}(u_{j},n_{j,l})\right)\frac{(u_{1})^{n_{1}}}{\Gamma(n_{1}+1)}\kappa_{1}(u_{1},n_{1,k}+1)\prod_{l\neq
k}\kappa_{1}(u_{1},n_{1,l})\\\
&=\Psi(\bm{u},K)\left(\prod_{j=2}^{K}\frac{(u_{j})^{n_{j}-1}}{\Gamma(n_{j})}\prod_{l=1}^{K}\kappa_{j}(u_{j},n_{j,l})\right)\frac{(u_{1})^{n_{1}}}{\Gamma(n_{1}+1)}\frac{\kappa_{1}(u_{1},n_{1,k}+1)}{\kappa_{1}(u_{1},n_{1,k})}\prod_{l=1}^{K}\kappa_{1}(u_{1},n_{1,l}),\end{split}$

where the second equality follows since we just multiplied and divided by
$\kappa_{1}(u_{1},n_{1,k})$.

Now, considering case $(ii)$: we are interested in the probability of sitting
to a new, non-existing, table where a new dish $\theta^{**}_{K+1}$ drawn from
$P_{0}$ is served. This equals to

$\begin{split}\mathbb{P}\bigl{(}\theta_{1,n_{1}+1}\in\cdot\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d}\bigr{)}\
=\
\Pi^{(n+1)}_{K+1}(\bm{n}^{\dagger}_{1},\dots,\bm{n}^{\dagger}_{d}),\end{split}$

where $\bm{n}^{\dagger}_{1}=(n_{1,1},\dots,n_{1,K},1)$ and
$\bm{n}^{\dagger}_{j}=(n_{j,1},\dots,n_{j,K},0)$ for each $j=2,\dots,d$. As
before, we further condition to $\bm{U}_{n}=u$ and we obtain

$\begin{split}\mathbb{P}\left(\theta_{1,n_{1}+1}\in\cdot\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d},\bm{u}\right)&\propto\
\Psi(\bm{u},K+1)\left(\prod_{j=2}^{K}\frac{(u_{j})^{n_{j}-1}}{\Gamma(n_{j})}\kappa_{j}(u_{j},0)\prod_{l=1}^{K}\kappa_{j}(u_{j},n_{j,l})\right)\\\
&\qquad\times\frac{(u_{1})^{n_{1}}}{\Gamma(n_{1}+1)}\kappa_{1}(u_{1},1)\prod_{l=1}^{K}\kappa_{1}(u_{1},n_{1,l})P_{0}(\cdot).\end{split}$

We can collect the previous expression in the following compact form

$\begin{split}&\mathbb{P}\left(\theta_{1,n_{1}+1}\in\cdot\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d},\bm{u}\right)\propto\\\
&\Psi(\bm{u},K)\sum_{k=1}^{K}\left\\{\left(\prod_{j=2}^{K}\frac{(u_{j})^{n_{j}-1}}{\Gamma(n_{j})}\prod_{l=1}^{K}\kappa_{j}(u_{j},n_{j,l})\right)\frac{(u_{1})^{n_{1}}}{\Gamma(n_{1}+1)}\frac{\kappa_{1}(u_{1},n_{1,k}+1)}{\kappa_{1}(u_{1},n_{1,k})}\prod_{l=1}^{K}\kappa_{1}(u_{1},n_{1,l})\delta_{\theta^{**}}(\cdot)\right\\}\\\
&+\Psi(\bm{u},K+1)\left(\prod_{j=2}^{K}\frac{(u_{j})^{n_{j}-1}}{\Gamma(n_{j})}\kappa_{j}(u_{j},0)\prod_{l=1}^{K}\kappa_{j}(u_{j},n_{j,l})\right)\frac{(u_{1})^{n_{1}}}{\Gamma(n_{1}+1)}\kappa_{1}(u_{1},1)\prod_{l=1}^{K}\kappa_{1}(u_{1},n_{1,l})P_{0}(\cdot).\end{split}$

The latter is a mixture (with unnormalized weights) of Dirac point masses
centered in the distinct values $\theta^{**}_{k}$ and the prior $P_{0}$.
Finally, we cancel out the common term in the unnormalized weights and we get

$\begin{split}\mathbb{P}\bigl{(}\theta&{}_{1,n_{1}+1}\in\cdot\mid\bm{\theta}_{1},\dots,\bm{\theta}_{d},\bm{u}\bigr{)}\\\
&\propto\Psi(\bm{u},K)\sum_{k=1}^{K}\left\\{\frac{\kappa_{1}(u_{1},n_{1,k}+1)}{\kappa_{1}(u_{1},n_{1,k})}\delta_{\theta^{**}}(\cdot)\right\\}+\Psi(\bm{u},K+1)\left(\prod_{j=2}^{K}\kappa_{j}(u_{j},0)\right)\kappa_{1}(u_{1},1)P_{0}(\cdot).\end{split}$

Dividing by $\Psi(\bm{u},K)$, Equation Equation S51 is obtained.

Figure S2 provides a graphical representation of the the Chinese restaurant
franchise process for Vec-NIFPP based on a sample of six, five and five
observations in three groups, respectively. The first (left panel) and the
second client (middle panel) enter the first restaurant and sit to different
tables. Tables having the same color must be prepared in all other
restaurants. Finally, (right panel) one of the possible configurations when
all client are seated. Global clustering is obtained by merging tables having
the same color.

Figure S2: Chinese restaurant franchise process representation for
$\operatorname{Vec-NIFPP}$ based on a sample of $6$, $5$ and $5$ observations
in $d=3$ groups, respectively.

##  S4 MCMC Algorithm for HMFM

From this point on, we place our attention to the HMFM model introduced in
Section 4.

###  S4.1 Hyperparameters setup

As mentioned in Section 4, we rely on the sample equivalence principle
(Diaconis and Ylvisaker, 1979) to design a suitable reparametrization of the
prior

$\Lambda\mid\gamma\sim\operatorname{Gamma}\left(a_{\Lambda}+da_{\gamma},b_{\Lambda}+b_{\gamma}\sum_{j=1}^{d}\gamma_{j}\right),$

with an easier interpretation. Indeed, it is straightforward to show that

$\mathbb{E}[\Lambda\mid\gamma]=\left(\frac{1}{\mathbb{E}[\Lambda]}\frac{a_{\Lambda}/a_{\gamma}}{a_{\Lambda}/a_{\gamma}+d}+\frac{1}{\widehat{\Lambda}}\frac{d}{a_{\Lambda}/a_{\gamma}+d}\right)^{-1},$

where $\mathbb{E}[\Lambda]$ is the prior mean of $\Lambda$, i.e.,
$a_{\Lambda}/b_{\Lambda}$ and
$\widehat{\Lambda}=\frac{a_{\gamma}}{b_{\gamma}}\frac{d}{\sum_{j=1}^{d}\gamma_{j}}$
is the maximum likelihood estimator for $\Lambda$ based on a sample
$\gamma_{1},\dots,\gamma_{d}$. Hence, $a_{\Lambda}/a_{\gamma}$ represents the
size of a prior sample and it can be compared to the size $d$. Let
$\Lambda_{0}$ and $\text{V}_{\Lambda}$ be the prior expected value and
variance of $\Lambda$, respectively, and $\gamma_{0}$ be a common prior guess
for the $\gamma_{j}$’s, so that
$\frac{1}{d}\sum_{j=1}^{d}\gamma_{j}=\gamma_{0}$. Then, assuming
$\mathbb{E}[\Lambda]=\widehat{\Lambda}$, we obtain

$a_{\gamma}=\frac{1}{d}\frac{\Lambda^{2}_{0}}{\text{V}_{\Lambda}},\quad
b_{\gamma}=\frac{a_{\gamma}}{\gamma_{0}\Lambda_{0}},\quad
a_{\Lambda}=\frac{\Lambda_{0}^{2}}{\text{V}_{\Lambda}},\quad
b_{\Lambda}=\frac{\Lambda_{0}}{\text{V}_{\Lambda}}.$

We now provide an intuition for choosing $\Lambda_{0}$, $\text{V}_{\Lambda}$
and $\gamma_{0}$. We would like our mixture to be sparse, which, in practice,
means that the expected number of clusters is (much) smaller than the expected
number of mixture components. Under the HMFM model, this translates to small
values of $\gamma_{j}$’s. Looking at each season separately, we note that each
element of the Vec-FDP is a finite DP. Hence, it is easy to compute the prior
distribution for each local number of cluster $K_{j}$, using Theorem 4 with
$n_{2}=0$, and compare it with the law of $M$. Setting a small number of
clusters for each season implies being a priori noninformative about the
heterogeneity across different seasons.

###  S4.2 Conditional Algorithm

The conditional algorithm for the HMFM model, introduced in Section 4.1,
provides full Bayesian inference on both the mixing parameters
$(P_{1},\dots,P_{d})$ and the clustering structure $\rho$. We use the
following notation, $\mathscr{L}\left(X\mid\text{rest}\right)$, to indicate
the distribution of a random variable $X$ conditionally to everything but
itself. The algorithm is a blocked Gibbs sampler whose state space is
$(\bm{c}_{1},\dots,\bm{c}_{d},\bm{\tau},\bm{S},M,\bm{U},\Lambda,\bm{\gamma})$.
The algorithm is composed of two main steps.

#### Step 1:

Metropolis-within-Gibbs algorithm to update one parameter at a time.

  * 1.1

Update cluster allocations $\bm{c}$. This step is done independently for each
group $j$ and for each individual $i$ by sampling $c_{j,i}$ from a discrete
distribution such that

$\mathbb{P}(c_{j,i}=m\mid\text{rest})\propto
S_{j,m}f(y_{j,i}\mid\tau_{m}),\quad\text{for }m=1,\dots,M.$

After sampling the whole vector, compute the number of allocated components
$K$.

  * 1.2

Update mixture parameters $\bm{\tau}=\left(\tau_{1},\dots,\tau_{M}\right)$.
From Theorem 3, the full conditional distribution further factorizes in two
independent parts, one for the first $K$ allocated components and one for the
remaining $M^{*}$ non-allocated components. In particular, for $k=1,\dots,K$,
sample $\tau_{k}$ independently from

$\mathbb{P}\left(\tau_{k}\in
d\tau_{k}\mid\text{rest}\right)\propto\left\\{\prod_{(j,i):c_{j,i}=k}f(y_{j,i}\mid\tau_{k})\right\\}P_{0}(d\tau_{k})$

while for $m=K+1,\dots,M$, sample $\tau_{m}$ independently from the prior
$P_{0}$.

  * 1.3

Update the unnormalized mixture weights $\bm{S}$ differently for the allocated
and non-allocated part of the process. In particular, for $k=1,\dots,K$,
sample $S_{j,k}$ from $\operatorname{Gamma}(\gamma_{j}+n_{j,k},U_{j}+1)$,
independently for each $j=1,\dots,d$, while for $m=K+1,\dots,M$, sample
$S_{j,m}$ from $\operatorname{Gamma}(\gamma_{j},U_{j}+1)$.

  * 1.4

Update $\Lambda$ from the mixture of gamma densities

$\begin{split}&\frac{K\left(b^{*}+1-\bm{\psi}(\bm{u})\right)}{(a^{*}-1)(\bm{\psi}(\bm{u})+K(b^{*}+1)}\operatorname{Gamma}\left(\Lambda;a^{*}+K-1,b^{*}+1-\bm{\psi}(\bm{u})\right)\\\
&\qquad+\frac{(a^{*}+K-1)(\bm{\psi}(\bm{u})}{(a^{*}-1)\bm{\psi}(\bm{u})+K(b^{*}+1)}\operatorname{Gamma}\left(\Lambda;a^{*}+K,b^{*}+1-\bm{\psi}(\bm{u})\right)\end{split}$

where $a^{*}=a_{\Lambda}+da_{\gamma};\
b^{*}=b_{\Lambda}+b_{\gamma}\sum_{j=1}^{d}\gamma_{j}.$

#### Step 2:

To avoid conditioning to $\bm{U}$ when updating $\bm{\gamma}$ and $M$, we
write the full conditional distribution of $\Delta_{2}$ as
$\mathscr{L}\left(\bm{U}\mid\text{rest}\right)\mathscr{L}\left(M,\bm{\gamma}\mid\Delta_{1},\bm{y}\right)$.

  * 2.1

Update $U_{j}$ from $\operatorname{Gamma}(n_{j},T_{j})$, independently for
each $j=1,\dots,d$.

  * 2.2

Update $M=K+M^{*}$, where $K$ has been updated in Step 1.1 and $M^{*}$ is
sampled from

$\mathbb{P}\left(M^{*}=m^{*}\mid\bm{\gamma},\Delta_{2},\bm{y}\right)\propto\frac{m^{*}+K}{m^{*}!}\Lambda^{m^{*}}\prod_{j=1}^{d}\frac{\Gamma\left(\gamma_{j}(m^{*}+K)\right)}{\Gamma\left(\gamma_{j}(m^{*}+K)+n_{j}\right)}.$

This step is accomplished using an adaptive Metropolis-Hastings step.

  * 2.3

Update $\gamma_{j}$ from a density proportional to

$\begin{split}\frac{\Gamma\left(\gamma_{j}(M^{*}+K)\right)}{\Gamma\left(\gamma_{j}(M^{*}+K)+n_{j}\right)}\prod_{k=1}^{K}\left(\frac{\Gamma\left(n_{j,k}+\gamma_{j}\right)}{\Gamma\left(\gamma_{j}\right)}\right)\operatorname{Gamma}(\gamma_{j};a_{\gamma},\Lambda
b_{\gamma}),\end{split}$

independently for each $j=1,\dots,d$. This update is performed via adaptive
Metropolis-Hastings step.

###  S4.3 Marginal Algorithm

We now turn our attention to the marginal algorithm for the HMFM model. The
algorithm is based on the Chinese restaurant franchise characterization of the
Vector of Finite Dirichlet process (see Section 3.2 in the main manuscript).
The state space of the Markov Chain is identified by the parameters
$\left(\bm{c}_{1},\dots,\bm{c}_{d},\bm{U},\Lambda,\bm{\gamma}\right)$. The
algorithm is composed of four steps:

#### Step 1:

Update cluster allocations $\bm{c}$. This is done in $n=\sum_{j=1}^{d}n_{j}$
sub-steps. For each restaurant $j=1,\dots,d$ and for each customers
$i=1,\ldots,n_{j}$, we remove observation $i$ from the cluster and denote by
$\rho_{j}^{-i}=\left\\{C_{j,1}^{-i},\ldots,C_{j,K^{-i}}^{-i}\right\\}$ the
resulting partition in $K^{-i}$ clusters. Then, the $i$-th observation is
assigned to a new cluster $C_{j,K^{-i}+1}^{-i}$ with probability proportional
to

$\begin{split}\mathbb{P}\left(i\in
C_{j,K^{-i}+1}^{-i}\mid\bm{u},\rho_{j}^{-i},\bm{y}\right)\propto\bm{\psi}(\bm{u})\gamma_{j}\Lambda\frac{K^{-i}+1+\Lambda\bm{\psi}(\bm{u})}{K^{-i}+\Lambda\bm{\psi}(\bm{u})}\mathscr{M}\left(y_{j,i}\right)\end{split}$

while it is assigned to each of the existing clusters $C_{j,l}^{-i}$, for
$l=1,\ldots,K^{-i}$, with probability proportional to

$\begin{split}\mathbb{P}\left(i\in
C_{j,l}^{-i}\mid\bm{u},\rho_{j}^{-i},\bm{y}\right)\propto\left(n_{j,k}^{-i}+\gamma_{j}\right)\frac{\mathscr{M}\left(\bm{y}_{C_{l}^{-i}\cup
i}\right)}{\mathscr{M}\left(\bm{y}_{C_{l}^{-i}}\right)}\end{split}$

where $\bm{\psi}(\bm{u})=\prod_{j=1}^{d}\psi_{j}\left(u_{j}\right)$ and, for
each global cluster $C_{l}\in\rho$, $\bm{y}_{C_{l}}$ is the vector of the
$y_{h,g}$, $h\in\\{1,\dots,d\\}$ and $g\in\\{1,\dots,n_{h}\\}$, such that
${h,g}\in C_{l}$ and

$\mathscr{M}\left(\bm{y}_{C_{l}}\right)=\int_{\mathbb{X}}\prod_{(h,g)\in
C_{l}}f\left(y_{h,g}\mid\theta\right)P_{0}(d\theta)$

is the marginal distribution of the data within the cluster $C_{l}$ with
sampling model $f\left(y_{h,g}\mid\theta\right)$ and prior $P_{0}(d\theta)$.

#### Step 2:

Update the whole vector $\bm{U}$ by sampling from a density proportional to

$\prod_{j=1}^{d}\left\\{U_{j}^{n_{j}-1}\frac{1}{(1+U_{j})^{n_{j}+K\gamma_{j}}}\right\\}\left(K+\Lambda\prod_{j=1}^{d}\frac{1}{(1+U_{j})^{\gamma_{j}}}\right)\exp\left\\{\Lambda\prod_{j=1}^{d}\frac{1}{(1+U_{j})^{\gamma_{j}}}\right\\}.$

#### Step 3:

Update the whole vector $\bm{\gamma}$ by sampling from a density proportional
to

$\begin{split}\prod_{j=1}^{d}&\left\\{\operatorname{Gamma}\left(\gamma_{j};a_{\gamma},\Lambda
b_{\gamma}\right)\frac{1}{(1+u_{j})^{K\gamma_{j}}}\prod_{m=1}^{K}\frac{\Gamma(\gamma_{j}+n_{j,m})}{\Gamma(\gamma_{j})}\right\\}\\\
&\times\left(K+\Lambda\prod_{j=1}^{d}\frac{1}{(1+U_{j})^{\gamma_{j}}}\right)\exp\left\\{\Lambda\prod_{j=1}^{d}\frac{1}{(1+U_{j})^{\gamma_{j}}}\right\\}\end{split}$

#### Step 4:

Update $\Lambda$ from the mixture of gamma densities

$\begin{split}&\frac{K\left(b^{*}+1-\bm{\psi}(\bm{u})\right)}{(a^{*}-1)(\bm{\psi}(\bm{u})+K(b^{*}+1)}\operatorname{Gamma}\left(\Lambda;a^{*}+K-1,b^{*}+1-\bm{\psi}(\bm{u})\right)\\\
&\qquad+\frac{(a^{*}+K-1)\bm{\psi}(\bm{u})}{(a^{*}-1)(\bm{\psi}(\bm{u})+K(b^{*}+1)}\operatorname{Gamma}\left(\Lambda;a^{*}+K,b^{*}+1-\bm{\psi}(\bm{u})\right)\end{split}$

where $a^{*}=a_{\Lambda}+da_{\gamma};\
b^{*}=b_{\Lambda}+b_{\gamma}\sum_{j=1}^{d}\gamma_{j}$.

Note that both Step 2 and 3 require to sample from $d$-dimensional
distributions which are not available in closed form and which can not be
factorized in the product of $d$ independent components. Joint updates must be
performed. We suggest to use adaptive Metropolis-Hastings techniques or
gradient based methods. In particular, we use here the Metropolis Adjusted
Langevin Algorithm (MALA) introduced by Roberts and Rosenthal (1998).

##  S5 Simulation study

In this section we present a simulation study to show our model’s ability to
perform clustering as well as density estimation. In particular, we
investigate the performance of the HMFM considering both the conditional
(HMFM-cond) and the marginal (HMFM-marg) samplers. Thus, we compare the HMFM
model with the HDP (Teh et al., 2006) mixture model. In the first scenario, we
also compare the HMFM against the MFM model that is assumed independently for
each group and fitted using the algorithm proposed by Argiento and De Iorio
(2022). The latter is implemented in the R package AntMAN (Ong et al., 2021).

The simulation study consists of three scenarios. The first experiment
considers $d=2$ groups which share a component. Thus, the example shows the
benefits of the hierarchical analysis against the independent group-specific
analysis. In the second scenario, again $d=2$ groups are considered, but
without any components shared across the groups. This example shows how the
HDP borrows too much information from the other group with respect to the
HMFM, leading to misleading results. Finally, the third experiment emphasizes
such differences in a more realistic situation, where $d=15$ groups are
considered; here, the HMFM outperforms the HDP in terms of both clustering and
density estimation. For all experiments, $50$ independent datasets have been
generated.

For a fair comparison, we fit all competing models by setting the same base
measure $P_{0}$. In particular, we follow the strategy explained in Section 4
and set a $\operatorname{Normal-InvGamma}$ prior with
$\mu_{0}=\text{mean}(\mathbf{y})$, $k_{0}=1/(\text{range}(\mathbf{y}))^{2}$,
$\nu_{0}=4$ and $\sigma^{2}_{0}=1/2$. Regarding the hyperparameters, for the
HDP we set default hyperprior $\alpha\sim\operatorname{Gamma}(1,1)$ and
$\gamma\sim\operatorname{Gamma}(1,0.1)$ (Teh et al., 2006). Also, we use
default hyperparameters for the independent analysis as suggested in the
AntMAN package, i.e., $\gamma\sim\operatorname{Gamma}(1,1)$ and
$\Lambda\sim\operatorname{Gamma}(1,1)$. Finally, for the HMFM, we follow the
strategy described in Section S4.1 to set the values of $\Lambda_{0}$,
$\text{V}_{\Lambda}$ and $\gamma_{0}$; the values are detailed for each
experiment below.

To assess the performance of recovering the true clustering, we compute the
Co-clustering error (CCE; Bassetti et al. 2020) and the Adjusted Rand Index
(ARI; Hubert and Arabie 1985). The CCE is based on the posterior probabilities
$\hat{\pi}_{lk}\in[0,1]$ of two data points $l$ and $k$ to belong to the same
cluster, i.e., the proportion of times that two observations have been
assigned to the same mixing component over the MCMC iterations. Specifically,
the error is defined as the average $L_{1}$ distance between the true pairwise
co-clustering matrix with elements $\pi_{lk}\in\\{0,1\\}$ and the estimated
co-clustering probability (similarity), namely

$\text{CCE}=\frac{1}{n}\sum_{l=1}^{n}\sum_{k=1}^{n}\ \lvert\
\pi_{lk}-\hat{\pi}_{lk}\ \rvert,$

where $n=\sum_{j=1}^{d}n_{j}$. Being a distance, lower values of CCE indicate
better fit, attaining zero in the absence of co-clustering error.

On the other hand, the ARI is computed using the estimated partition which is
obtained by minimizing the variation of information function (Wade and
Ghahramani, 2018; Dahl et al., 2022). The index ranges between zero (the true
and the estimated clustering do not agree on any pair of points) and one (the
two clustering are the same).

Regarding density estimation, we compute a Predictive Score (PS; Bassetti et
al. 2020) for each group $j=1,\dots,d$, that is the $L_{1}$ distance between
the group-specific density $f(y_{j,n_{j}+1})$ and the corresponding predictive
density $\hat{f}(y_{j,n_{j}+1}\mid\mathbf{Y})$, namely

$\text{PS}_{j}=\int_{-\infty}^{\infty}\ \lvert
f(y_{j,n_{j}+1})-\hat{f}(y_{j,n_{j}+1}\mid\mathbf{Y})\rvert\ dy_{j,n_{j}+1}.$

The PS is a distance so lower values indicate better fit, with zero being the
lowest possible value.

###  S5.1 Experiment 1

This experiments considers $d=2$ groups, both having two local clusters
($K_{1}=2$, $K_{2}=2$) one of which is shared; hence, the global number of
clusters is $K=3$. The mixing probabilities are set so that the shared
component has a lower value in the second group. Namely

$\begin{split}y_{1,i}\ &\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
0.5\text{N}(-3,0.1)+0.5\text{N}(0,0.5),\ \text{for }i=1,\dots,300\\\ y_{2,i}\
&\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
0.2\text{N}(0,0.5)+0.8\text{N}(1.75,1.5),\ \text{for }i=1,\dots,300\\\
\end{split}$ (S52)

Hence, the component specific parameters are $\mu=(-3,0,1.75)$ and
$\sigma=(\sqrt{0.1},\sqrt{0.5},\sqrt{1.5})$. We generate $n_{1}=300$ and
$n_{2}=300$ observations for each group. Note that the second group is defined
so that the two components strongly overlap. As a result, the shared component
is completely masked in this group. Nevertheless, the masked component can be
spotted by exploiting the sharing of information with the other group. Figure
S3 shows the empirical distribution of a dataset generated from model Equation
(S52) and the underlying densities.

Figure S3: Empirical distributions of a dataset simulated under Experiment 1.
Dots represent the observations while lines represent the underlying
densities. Colors relate to the mixing components.

For each simulated dataset, we fit the HDP, the indepedent group-specific MFM
and the proposed HMFM by setting, following the guidelines in Section S4.1,
$\Lambda_{0}=5$, $\text{V}_{\Lambda}=5$, $\gamma_{0}=0.5$. A clear advantage
of joint modelling is the possibility of deriving model-based clustering also
across different groups, which is not possible when we run independent
analyses. However, in the sake of a fair comparison with the MFM, clustering
performances are evaluated within each group, not globally. Table S1 presents
the results of model comparison based on the mean and standard deviation (in
brackets) of the ARI over the simulated datasets. The results show that all
methods, with negligible differences, are able to perfectly recover the true
clustering in the first group, where the two components are well separated.
Instead, the advantage of the sharing of information allowed by hierarchical
modeling is evident in the second group, where components overlap. Indeed, the
MFM fails to identify the presence of two different clusters in the second
group, gathering all observations together and obtaining an ARI value close to
zero. On the other hand, HDP and HMFM are able to borrow information from the
first group to recognize the presence of two clusters in the second group,
leading to higher values of ARI, with HMFM (both marginal and conditional
samplers) outperforming the HDP. To clarify, we compute the percentage of
simulated datasets for which each method is able to obtain a partition of the
second group consisting only of a single cluster. The results are reported in
the third column of Table S1 and show how the MFM consistently miss to
identify at least two clusters for the second group, showing the limitations
of the independent analyses against the hierarchical approach. Finally, we
point out that the HMFM is able to outperform the HDP, with slightly better
performances for the marginal sampler over the conditional one.

Table S1: The first two columns report the mean and the standard deviation (in brackets) of the ARI for the two groups over the $50$ simulated datasets under Experiment 1. The last column shows the percentage of times that each method gathers all observations of the second group in a single cluster, hence, lower values are better. | Group 1 | Group 2 | % 1 cluster
---|---|---|---
MFM | 0.992 (0.011) | 0.002 (0.037) | 98%
HDP | 0.992 (0.011) | 0.122 (0.183) | 68%
HMFM - Marginal | 0.990 (0.012) | 0.216 (0.191) | 42%
HMFM - Conditional | 0.991 (0.012) | 0.177 (0.190) | 52%

To conclude the clustering comparison, Figure S4 provides a picture of the
distribution of the CCE within each group over the simulated datasets. Not
surprisingly, differences are negligible in the first group while the MFM
commits an higher error in the second one, with also higher variability among
different datasets.

Figure S4:  Co-Clustering Error in the first (left) and second (right) group.
Boxplot are obtained over $50$ datasets simulated under Experiment 1.

Similar conclusions can be drawn in terms of density estimation. Figure S5
displays the boxplots for the group-specific PS. The figure shows that
differences are negligible in the first group, for which all methods provide a
good estimate of the density. However, the density for the second group
estimated using the MFM is worse than the one obtained by the HMFM and the
HDP.

Figure S5: Predictive Score in the first (left) and second (right) group.
Boxplot are obtained over $50$ datasets simulated under Experiment 1.

###  S5.2 Experiment 2

This experiment considers $d=2$ groups coming from a single component, not
shared across the groups, namely $K_{1}=1$, $K_{2}=1$ so that $K=2$. The
component specific parameters are $\mu=(0,1)$ and $\sigma=(1,1)$, hence the
sampling model is $y_{1,i}\ \stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
\text{N}(0,1)$ and $y_{2,i}\ \stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
\text{N}(1,1)$, each for $i=1,\dots,n_{j}$.

Figure S6: Empirical distributions of a dataset simulated under Experiment 2.
Dots represent the observations while lines represent the underlying
densities. Colors relate to the mixing components.

We repeat the experiment for an increasing number of observations,
$n=50,100,200$, which we equally divide in the two groups. See Figure S6 for
an example of the empirical distribution of a dataset with n=100 generated
from the model and the underlying densities. This experiment is designed to
assess whether the borrowing of information may lead to misleading results in
situations where groups do not share any features. Hence, for this scenario,
we confine the comparison between the HDP and the HMFM. The latter, is fitted
by setting $\Lambda_{0}=10$, $\text{V}_{\Lambda}=2$, $\gamma_{0}=0.01$.

Table S2 presents the results of model comparison based on the mean and
standard deviation (in brackets) of the ARI over the simulated datasets.
Although the scenario is simple and we would expect that both methods will
recover the true partition, the table shows that the HDP struggles to find the
underlying global clustering, especially for a small number of observations.
Rather, the HMFM is able to recover the true partition for all the sample
sizes.

Table S2:  Mean and standard deviation (in brackets) of the ARI for the two groups over the 50 simulated datasets under Experiment 2.  | n = 50 | n = 100 | n = 200
---|---|---|---
HDP | 0.515 (0.484) | 0.877 (0.310) | 0.957 (0.172)
HMFM - Marginal | 0.980 (0.141) | 0.999 (0.005) | 1 (0)
HMFM - Conditional | 0.978 (0.141) | 0.994 (0.033) | 1 (0)

Same conclusions can be drawn from the CCE, see Figure S7, which takes into
account all the posterior pairwise probabilities of observations to be
clustered together. Clearly, the HMFM outperforms the HDP in terms of
clustering estimation. Also, the proposed model provides better estimates of
the densities over the HDP, see Figure S8, although differences are here less
evident.

This simple example showcases that the HDP’s propensity to excessively share
information across groups can compromise the accurate recovery of the true
underlying partition. To understand why, we resort to the restaurant franchise
metaphor discussed in Section 3.2. According to the HDP, when a new customer
arrives, the probability of consuming a dish not yet served in that specific
restaurant, but available in other franchise restaurants, hinges on the
cumulative number of clients eating that dish across all restaurants. In
contrast, the HMFM uses information only about whether a dish is being
consumed or not in other restaurants. This experiment unveils a potential
issue in the HDP’s clustering mechanism; the concentration of all clients
within the first restaurant (group) eating the same dish increases the
probability of that dish being offered in the second restaurant (group), even
though no components are shared between them. Rather, this confusion is
circumvented by the HMFM’s clustering mechanism. We notice that this
phenomenon is directly tied to the choice of prior, and its impact diminishes
as the sample size increases, as shown in Table S2. However, even though the
HDP reaches the correct clustering for $n=200$, Figure S7 shows that the HMFM
model still outperforms the competing approach.

Figure S7:  Co-Clustering Error for different sample sizes. Boxplots are
obtained over $50$ datasets simulated under Experiment 2.

Figure S8:  Predictive score for the first group (top panels) and the second
group (bottom panels) for different sample sizes. Boxplot are obtained over
$50$ datasets simulated under Experiment 2.

###  S5.3 Experiment 3

The final experiment considers $d=15$ groups, each consisting of $n_{j}=30$
observations. Data have been generated so that for the first $12$ groups, we
consider the following component specific parameters $\mu_{1}=(-3,0,1)$,
$\sigma_{1}=(\sqrt{0.5},\sqrt{0.5},\sqrt{0.5})$. To increase the variability
and differences between the groups, for each $j=1,\dots,12$ we first sample
the number of local components $K_{j}\in\\{2,3\\}$ and then, conditionally to
$K_{j}$, we choose uniformly the components to be selected. Regarding the
weights, for each of the $12$ groups, we assign most of the mass to the second
component, if included, namely, $\pi_{j,2}=0.5$ if $K_{j}=3$ or
$\pi_{j,2}=2/3$ if $K_{j}=2$. The remaining components receive the remaining
mass split in equal parts. Finally, the last three groups all share two
components whose parameters are $\mu_{2}=(-1.5,1.5)$,
$\sigma_{2}=(\sqrt{0.5},\sqrt{0.5})$ with equal mixing weights. Then, the
global number of components is $K=5$. Figure S9 shows an example of data
simulated under this scenario.

Figure S9: Empirical distributions of a dataset simulated under Experiment 3.
Dots represent the observations while lines represent the underlying
densities. Colors relate to the mixing components.

For each simulated dataset, we fit the HDP and the HMFM; in particular, the
latter, is fitted by setting the hyperparameters $\Lambda_{0}=15$,
$\text{V}_{\Lambda}=3$, $\gamma_{0}=0.05$. Figure S10 shows the boxplot of the
ARI (left) and CCE (right) of the global clustering evaluated over the
simulated datasets. Both metrics show that the HFMF significantly outperforms
the HDP. The same conclusions can be drawn also inspecting the ARI and the CCE
withing each group (not shown here for brevity).

Figure S10:  ARI (left) and CCE (right) of the global clustering evaluated
over the datasets simulated under the Experiment 3.

Similar to the previous experiments, the proposed mixture model outperforms
the HDP also in terms on group-specific density estimation, even though the
differences are less evident than for the clustering.

Figure S11:  Predictive Score in all the groups. Boxplot are obtained over
$50$ datasets simulated under Experiment 3.

###  S5.4 Comparison of computational times

The goal of this section it to asses the scalability of the HDP and the HMFM
when the number of observation increases. In particular, we consider
$n=(24,50,100,200,400,800,1600)$ data, equally divided in $d=2$ groups. Data
are generated from a simple mixture model with well separated components,
enabling all methods to quickly converge to the true clustering. We run each
method for $10,000$ iterations, burning out half of them. The total execution
time is then divided for the number of iterations, so that the mean time for a
single iteration is computed and displayed in Figure S12.

Figure S12:  Execution time in log-scale for the HMFM, HDP and calculation of
Stirling numbers. The dashed line and the dotted line represent, respectively,
a linear and a quadratic increasing trend with the number of observations $n$.

The HDP sampling strategy is based on the direct assignment scheme suggested
by Teh et al. (2006). Differently from the traditional Chinese restaurant
franchise process where observations are firstly assigned to the tables and
then the menu (mixture components) are assigned to the tables, the direct
assignment sampler assigns data directly to the mixture components. The
computational bottleneck of this algorithm is the computation of the unsigned
Stirling numbers $s(n_{j,k},m)$ for each $j=1,\dots,d$, $k=1,\dots,K$ and for
all positive integers $m\leq n_{j,k}$, where $n_{j,k}$ is the number of
observations in group $j$ assigned to cluster $k$. The computational effort to
compute $s(n_{j,k},m)$ is proportional to $(n_{j,k})^{2}$, as shown in Figure
S12. As expected, for large $n$, the slope of the computational time needed to
calculate the Striling numbers equals the one of the HDP. In particular, such
slope is steeper than the one of HMFM (both marginal and conditional sampler)
which grows with a linear trend. This empirically verifies that our proposed
model is faster and scales much better when increasing the number of
observations.

##  S6 Additional results on shot put data analysis

Here we detail the complete specification of the HMFM model for the shot put
data:

$\begin{split}\bm{y}_{j,1},\dots,\bm{y}_{j,n_{j}}\mid\bm{S}_{j},\theta,M&\
\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
\sum_{m=1}^{M}\frac{S_{j,m}}{T_{j}}d\text{N}_{N_{j,i}}(\bm{y}_{j,i}\mid\mu_{m}\mathbf{1}_{N_{j,i}}+X_{j,i}\bm{\beta}_{j},\sigma^{2}_{m}\mathbf{I}_{N_{j,i}}),\quad\text{for
}j=1,\dots,d\\\ \bm{\beta}_{1},\dots,\bm{\beta}_{d}&\
\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
\text{N}_{d}\left(\bm{\beta}_{0},\Sigma_{0}\right)\\\
\tau_{1},\dots,\tau_{M}\mid M&\ \stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
\operatorname{Normal-InvGamma}\left(\mu_{0},k_{0},\nu_{0},\sigma_{0}\right)\\\
S_{j,1},\dots,S_{j,M}\mid M,\gamma_{j}&\
\stackrel{{\scriptstyle\mathrm{iid}}}{{\sim}}\
\operatorname{Gamma}(\gamma_{j},1),\quad\text{for }j=1,\dots,d\\\
M\mid\Lambda&\
\mathrel{\overset{\makebox[0.0pt]{\mbox{\tiny\leavevmode\nobreak\ }}}{\sim}}\
\operatorname{Pois}_{1}(\Lambda)\\\ \gamma_{1},\dots,\gamma_{d}\mid\Lambda&\
\sim\ \operatorname{Gamma}(a_{\gamma},\Lambda b_{\gamma})\\\ \Lambda&\ \sim\
\operatorname{Gamma}(a_{\Lambda},b_{\Lambda}),\end{split}$

where $d\text{N}_{N_{j,i}}(\bm{y}\mid\bm{\mu},\Sigma)$ is the density of a
multivariate normal distribution evaluated in $\bm{y}$ with mean $\bm{\mu}$
and covariance matrix $\Sigma$.

Figure S13: Shot put measurements collected throughout an athlete’s career for
four randomly selected athletes. Vertical dotted lines delimit seasons.

Table S3 presents the local clustering sizes, for each season and for each
cluster. Finally, Table S4 shows some summary statistics for the global
clusters and compare them with the overall dataset.

Table S3: Cluster cardinalities for each cluster (rows) and for each season (columns).  | S1 | S2 | S3 | S4 | S5 | S6 | S7 | S8 | S9 | S10 | S11 | S12 | S13 | S14 | S15
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Cluster 1 | 0 | 4 | 5 | 9 | 7 | 12 | 6 | 9 | 11 | 7 | 7 | 5 | 6 | 5 | 3
Cluster 2 | 14 | 10 | 14 | 13 | 19 | 16 | 21 | 9 | 11 | 14 | 4 | 9 | 3 | 6 | 4
Cluster 3 | 11 | 14 | 10 | 31 | 34 | 31 | 34 | 31 | 24 | 17 | 17 | 0 | 8 | 2 | 0
Cluster 4 | 0 | 14 | 20 | 45 | 54 | 51 | 38 | 35 | 46 | 24 | 10 | 20 | 0 | 0 | 9
Cluster 5 | 33 | 29 | 41 | 61 | 70 | 66 | 49 | 39 | 20 | 11 | 12 | 0 | 17 | 6 | 0
Cluster 6 | 0 | 0 | 14 | 11 | 7 | 11 | 19 | 0 | 0 | 0 | 4 | 1 | 0 | 0 | 0
Cluster 7 | 9 | 39 | 70 | 108 | 79 | 42 | 36 | 52 | 29 | 21 | 21 | 29 | 0 | 0 | 7
Cluster 8 | 29 | 85 | 109 | 79 | 72 | 59 | 48 | 23 | 30 | 29 | 4 | 0 | 7 | 12 | 0
Cluster 9 | 75 | 148 | 90 | 42 | 38 | 46 | 34 | 36 | 22 | 13 | 18 | 9 | 16 | 0 | 2
Cluster 10 | 200 | 43 | 30 | 4 | 3 | 7 | 9 | 21 | 18 | 22 | 12 | 16 | 0 | 0 | 6
Cluster 11 | 31 | 17 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 12 | 0 | 0 | 11 | 0
Cluster 12 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 2
Cluster 13 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Cluster 14 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 0
Cluster 15 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 3 | 0 | 0 | 0
Table S4: Cluster statistics for each global cluster. The first row represents the case where all data are pooled in a single cluster. Each subsequent row corresponds to a specific global cluster. The columns display the sample mean of various cluster statistics, differentiated by M (male) or F (female) athletes in the cluster. Cluster | SizeM | SizeF | MeanM | MeanF | VarM | VarF | AgeM | AgeF
---|---|---|---|---|---|---|---|---
Pooled | 1925 | 1683 | 1.12 | -1.29 | 1.51 | 2.16 | 24.47 | 22.77
1 | 42 | 54 | 3.47 | 1.50 | 0.28 | 0.48 | 27.70 | 26.08
2 | 77 | 90 | 2.85 | 0.72 | 0.28 | 0.38 | 27.52 | 24.91
3 | 147 | 117 | 2.26 | 0.17 | 0.28 | 0.42 | 25.79 | 25.05
4 | 234 | 132 | 1.70 | -0.39 | 0.29 | 0.42 | 25.80 | 23.54
5 | 260 | 194 | 1.18 | -1.05 | 0.21 | 0.33 | 24.41 | 22.34
6 | 42 | 25 | 1.14 | -1.02 | 0.86 | 0.98 | 24.24 | 21.69
7 | 312 | 228 | 0.72 | -1.56 | 0.28 | 0.39 | 23.84 | 22.17
8 | 316 | 270 | 0.24 | -2.08 | 0.26 | 0.34 | 22.71 | 21.85
9 | 277 | 309 | -0.37 | -2.68 | 0.28 | 0.34 | 22.18 | 21.38
10 | 178 | 212 | -1.07 | -3.30 | 0.23 | 0.33 | 22.05 | 21.03
11 | 35 | 41 | -1.81 | -3.56 | 0.29 | 1.20 | 21.40 | 25.78
12 | 0 | 6 | - | 2.63 | - | 0.32 | - | 28.05
13 | 2 | 1 | 1.08 | 0.00 | 13.05 | 1.83 | 26.96 | 28.00
14 | 1 | 1 | -0.46 | -4.19 | 2.19 | 0.47 | 30.80 | 31.00
15 | 2 | 3 | -4.76 | -4.20 | 18.10 | 2.60 | 36.33 | 25.42

## References

  * Argiento et al. (2020) Argiento, R., Cremaschi, A., and Vannucci, M. (2020), “Hierarchical normalized completely random measures to cluster grouped data,” Journal of the American Statistical Association, 115, 318–333.
  * Argiento and De Iorio (2022) Argiento, R. and De Iorio, M. (2022), “Is infinity that far? A Bayesian nonparametric perspective of finite mixture models,” The Annals of Statistics, 50, 2641–2663.
  * Baccelli et al. (2020) Baccelli, F., Błaszczyszyn, B., and Karray, M. (2020), Random measures, point processes, and stochastic geometry, Inria.
  * Bassetti et al. (2020) Bassetti, F., Casarin, R., and Rossini, L. (2020), “Hierarchical species sampling models,” Bayesian Analysis, 15, 809–838.
  * Beraha et al. (2023) Beraha, M., Argiento, R., Camerlenghi, F., and Guglielmi, A. (2023), “Normalized random meaures with interacting atoms for Bayesian nonparametric mixtures,” arXiv 2302.09034.
  * Beraha et al. (2021) Beraha, M., Argiento, R., Møller, J., and Guglielmi, A. (2021), “MCMC cmputations for Bayesian mixture models using repulsive point processes,” Journal of Computational and Graphical Statistics, 31, 1–37.
  * Camerlenghi et al. (2019) Camerlenghi, F., Lijoi, A., Orbanz, P., and Prünster, I. (2019), “Distribution theory for hierarchical processes,” The Annals of Statistics, 47, 67–92.
  * Charalambides (2002) Charalambides, C. A. (2002), Enumerative combinatorics, CRC Press.
  * Dahl et al. (2022) Dahl, D. B., Johnson, D. J., and Müller, P. (2022), “Search algorithms and loss functions for Bayesian clustering,” Journal of Computational and Graphical Statistics, 31, 1189–1201.
  * Daley and Vere-Jones (2008) Daley, D. J. and Vere-Jones, D. (2008), An introduction to the theory of point processes, Probability and its Applications, Springer, New York, second edition.
  * D’Angelo et al. (2023) D’Angelo, L., Canale, A., Yu, Z., and Guindani, M. (2023), “Bayesian nonparametric analysis for the detection of spikes in noisy calcium imaging data,” Biometrics, 79, 1370–1382.
  * De Blasi et al. (2015) De Blasi, P., Favaro, S., Lijoi, A., Mena, R. H., Prünster, I., and Ruggiero, M. (2015), “Are Gibbs-type priors the most natural generalization of the Dirichlet process?” IEEE Transactions on Pattern Analysis and Machine Intelligence, 37, 212–229.
  * Denti et al. (2023) Denti, F., Camerlenghi, F., Guindani, M., and Mira, A. (2023), “A common atoms model for the Bayesian nonparametric analysis of nested data,” Journal of the American Statistical Assocciation, 118, 405–416.
  * Diaconis and Ylvisaker (1979) Diaconis, P. and Ylvisaker, D. (1979), “Conjugate priors for exponential families,” The Annals of Statistics, 7, 269–281.
  * Dolmeta et al. (2023) Dolmeta, P., Argiento, R., and Montagna, S. (2023), “Bayesian GARCH modeling of functional sports data,” Statistical Methods & Applications, 32, 401–423.
  * Fang and Lai (1997) Fang, H. and Lai, T.-Y. (1997), “Co-kurtosis and capital asset pricing,” Financial Review, 32, 293–307.
  * Favaro and Teh (2013) Favaro, S. and Teh, Y. W. (2013), “MCMC for normalized random measure mixture models,” Statistical Science, 28, 335–359.
  * Friend and Westerfield (1980) Friend, I. and Westerfield, R. (1980), “Co-skewness and capital asset pricing,” The Journal of Finance, 35, 897–913.
  * Frühwirth-Schnatter and Malsiner-Walli (2019) Frühwirth-Schnatter, S. and Malsiner-Walli, G. (2019), “From here to infinity: sparse finite versus Dirichlet process mixtures in model-based clustering,” Advances in data analysis and classification, 13, 33–64.
  * Ghilotti et al. (2023) Ghilotti, L., Beraha, M., and Guglielmi, A. (2023), “Bayesian clustering of high-dimensional data via latent repulsive mixtures,” arXiv 2303.02438.
  * Griffin and Leisen (2017) Griffin, J. E. and Leisen, F. (2017), “Compound random measures and their use in Bayesian non-parametrics,” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 79, 525–545.
  * Hoff (2009) Hoff, P. D. (2009), A first course in Bayesian statistical methods, Springer.
  * Hubert and Arabie (1985) Hubert, L. and Arabie, P. (1985), “Comparing partitions,” Journal of Classification, 2, 193–218.
  * James et al. (2009) James, L. F., Lijoi, A., and Prünster, I. (2009), “Posterior analysis for normalized random measures with independent increments,” Scandinavian Journal of Statistics, 36, 76–97.
  * Kallenberg (2005) Kallenberg, O. (2005), Probabilistic symmetries and invariance principles, Probability and its Applications, Springer, New York.
  * Kallenberg (2017) — (2017), Random measures, theory and applications, volume 77 of Probability Theory and Stochastic Modelling, Springer, Cham.
  * Kallenberg (2021) — (2021), Foundations of modern probability, volume 99 of Probability Theory and Stochastic Modelling, Springer Cham.
  * Lijoi et al. (2007a) Lijoi, A., Mena, R. H., and Prünster, I. (2007a), “Controlling the reinforcement in Bayesian non-parametric mixture models,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69, 715–740.
  * Lijoi et al. (2007b) Lijoi, A., Mena, R. H., and Prünster, I. (2007b), “Bayesian nonparametric estimation of the probability of discovering new species,” Biometrika, 94, 769–786.
  * Lijoi et al. (2014) Lijoi, A., Nipoti, B., and Prünster, I. (2014), “Bayesian inference with dependent normalized completely random measures,” Bernoulli, 20, 1260–1291.
  * Malsiner-Walli et al. (2017) Malsiner-Walli, G., Frühwirth-Schnatter, S., and Grün, B. (2017), “Identifying mixtures of mixtures using Bayesian estimation,” Journal of Computational and Graphical Statistics, 26, 285–295.
  * Malsiner-Walli et al. (2016) Malsiner-Walli, G., Frühwirth-Schnatter, S., and Grün, B. (2016), “Model-based clustering based on sparse finite Gaussian mixtures,” Statistics and computing, 26, 303–324.
  * Miller and Harrison (2018) Miller, J. W. and Harrison, M. T. (2018), “Mixture models with a prior on the number of components,” Journal of the American Statistical Association, 113, 340–356.
  * Nobile (2004) Nobile, A. (2004), “On the posterior distribution of the number of components in a finite mixture,” The Annals of Statistics, 32, 2044–2073.
  * Ong et al. (2021) Ong, P., Argiento, R., Bodin, B., and De Iorio, M. (2021), AntMAN: Anthology of Mixture Analysis Tools, URL https://CRAN.R-project.org/package=AntMAN. R package version 1.1.0.
  * Page et al. (2013) Page, G., Barney, B., and McGuire, A. (2013), “Effect of position, usage rate, and per game minutes played on NBA player production curves,” Journal of Quantitative Analysis in Sports, 9, 337–345.
  * Petralia et al. (2012) Petralia, F., Rao, V., and Dunson, D. (2012), “Repulsive Mixtures,” in Pereira, F., Burges, C., Bottou, L., and Weinberger, K. (editors), Advances in Neural Information Processing Systems, volume 25, Curran Associates, Inc.
  * Regazzini et al. (2003) Regazzini, E., Lijoi, A., and Prünster, I. (2003), “Distributional results for means of normalized random measures with independent increments,” The Annals of Statistics, 31, 560–585.
  * Richardson and Green (1997) Richardson, S. and Green, P. J. (1997), “On Bayesian analysis of mixtures with an unknown number of components (with discussion),” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59, 731–792.
  * Roberts and Rosenthal (1998) Roberts, G. O. and Rosenthal, J. S. (1998), “Optimal Scaling of Discrete Approximations to Langevin Diffusions,” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 60, 255–268.
  * Rodríguez et al. (2008) Rodríguez, A., Dunson, D. B., and Gelfand, A. E. (2008), “The nested Dirichlet process,” Journal of the American Statistical Association, 103, 1131–1144.
  * Rousseau and Mengersen (2011) Rousseau, J. and Mengersen, K. (2011), “Asymptotic behaviour of the posterior distribution in overfitted mixture models,” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 73, 689–710.
  * Teh et al. (2006) Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M. (2006), “Hierarchical Dirichlet processes,” Journal of the American Statistical Association, 101, 1566–1581.
  * Van Havre et al. (2015) Van Havre, Z., White, N., Rousseau, J., and Mengersen, K. (2015), “Overfitting Bayesian mixture models with an unknown number of components,” PloS one, 10, e0131739.
  * Wade and Ghahramani (2018) Wade, S. and Ghahramani, Z. (2018), “Bayesian cluster analysis: Point estimation and credible balls (with discussion),” Bayesian Analysis, 13, 559–626.


# Observation of a pseudogap in the optical conductivity of underdoped
Ba1-xKxFe2As2

Y. M. Dai  LPEM, ESPCI-ParisTech, CNRS, UPMC, 10 rue Vauquelin, F-75231 Paris
Cedex 5, France  Beijing National Laboratory for Condensed Matter Physics,
National Laboratory for Superconductivity, Institute of Physics, Chinese
Academy of Sciences, P.O. Box 603, Beijing 100190, China  B. Xu  B. Shen
Beijing National Laboratory for Condensed Matter Physics, National Laboratory
for Superconductivity, Institute of Physics, Chinese Academy of Sciences, P.O.
Box 603, Beijing 100190, China  H. H. Wen  Beijing National Laboratory for
Condensed Matter Physics, National Laboratory for Superconductivity, Institute
of Physics, Chinese Academy of Sciences, P.O. Box 603, Beijing 100190, China
National Laboratory of Solid State Microstructures and Department of Physics,
Nanjing University, Nanjing 210093, China  J. P. Hu  Beijing National
Laboratory for Condensed Matter Physics, National Laboratory for
Superconductivity, Institute of Physics, Chinese Academy of Sciences, P.O. Box
603, Beijing 100190, China  Department of Physics, Purdue University, West
Lafayette, Indiana 47907, USA  X. G. Qiu<EMAIL_ADDRESS>Beijing
National Laboratory for Condensed Matter Physics, National Laboratory for
Superconductivity, Institute of Physics, Chinese Academy of Sciences, P.O. Box
603, Beijing 100190, China  R. P. S. M. Lobo<EMAIL_ADDRESS>LPEM, ESPCI-
ParisTech, CNRS, UPMC, 10 rue Vauquelin, F-75231 Paris Cedex 5, France

(August 27, 2024)

###### Abstract

We report the observation of a pseudogap in the _ab_ -plane optical
conductivity of underdoped Ba1-xKxFe2As2 ($x=0.2$ and 0.12) single crystals.
Both samples show prominent gaps opened by a spin density wave (SDW) order and
superconductivity at the transition temperatures $T_{\it SDW}$ and $T_{c}$,
respectively. In addition, we observe an evident pseudogap below
$T^{\ast}\sim$ 75 K, a temperature much lower than $T_{\it SDW}$ but much
higher than $T_{c}$. A spectral weight analysis shows that the pseudogap is
closely connected to the superconducting gap, indicating the possibility of
its being a precursor of superconductivity. The doping dependence of the gaps
is also supportive of such a scenario.

###### pacs:

74.25.Gz, 74.70.Xa, 78.30.-j

Among all the families of iron-pnictide superconductors discovered to date,
Kamihara et al. (2008); Rotter et al. (2008a); Sefat et al. (2008); Li et al.
(2009); Torikachvili et al. (2008); Tapp et al. (2008); Hsu et al. (2008);
Sales et al. (2009) the BaFe2As2 (Ba122) family is one of the most studied.
The parent BaFe2As2 composition is a poor Pauli-paramagnetic metal with a
structural and magnetic phase transition at 140 K. Rotter et al. (2008b)
Superconductivity arises with the suppression of magnetism which can be
achieved by applying pressure Torikachvili et al. (2008) or chemical
substitution. Rotter et al. (2008a); Sefat et al. (2008); Li et al. (2009) The
substitution of Ba with K atoms yields hole-doping Rotter et al. (2008a) with
a maximum $T_{c}\approx 39$ K and the substitution of Fe atoms by Co or Ni
atoms results in electron-doping Sefat et al. (2008); Li et al. (2009) with a
maximum $T_{c}\approx 25$ K. Extensive studies have been carried out in the
parent BaFe2As2, Hu et al. (2008); Akrap et al. (2009) electron-doped
BaFe${}_{2-x}A_{x}$As2 ($A$ = Co, Ni), Lobo et al. (2010); Tu et al. (2010);
Teague et al. (2011); Terashima et al. (2009) as well as optimally hole-doped
Ba0.6K0.4Fe2As2 Li et al. (2008); Shan et al. (2011); Ding et al. (2008)
compounds. However, the hole-underdoped regime of the phase diagram is
relatively unexplored.

This hole-underdoped region is arguably the most important regime because of
the following two reasons. First, the superconducting mechanism is deeply tied
with magnetism. The interplay between magnetism and superconductivity is
manifest in this regime. In a considerably large portion of the underdoped
regime, the SDW phase and superconductivity coexist. Park et al. (2009); Goko
et al. (2009); Aczel et al. (2008); Massee et al. (2009); Chia et al. (2010)
Second, in cuprates, the most exciting, yet puzzling, physics takes place in
the hole-underdoped regime. This regime thus is pivotal to the comparison
between iron-pnictides and cuprates.

Xu _et al._ have performed the surface sensitive angle-resolved photoemission
(ARPES) measurements on underdoped Ba1-xKxFe2As2. Xu et al. (2011) Their data
showed a distinct pseudogap coexisting with the superconducting gap and
suggested that both the pseudogap and superconductivity are driven by
antiferromagnetic fluctuations. However, one key issue in understanding the
origin of the pseudogap and, in particular, its relation to superconductivity
is the question of whether it shares electronic states with the
superconducting condensation. Yu et al. (2008) Infrared spectroscopy probes
the charge dynamics of bulk materials and the spectral weight analysis is a
powerful tool to address this issue.

We present broadband infrared spectroscopy measurements on two underdoped
Ba1-xKxFe2As2 ($x=0.2$ and 0.12) single crystals. In both samples, the opening
of the SDW gap and the superconducting gap was clearly observed on the optical
conductivity. In addition, another small gap opens below $T^{\ast}\sim$ 75 K,
closely resembling the famous pseudogap in the hole-underdoped cuprates. We
find that the SDW gap depletes the spectral weight available for the
superconducting condensate, which suggests that the SDW order competes with
superconductivity. However, both doping and temperature dependence of the
spectral weight inside the pseudogap indicate that it shares the same
electronic origin with the superconducting gap.

High quality Ba1-xKxFe2As2 single crystals were grown by the self-flux method
using FeAs as the flux. Luo et al. (2008)

Figure 1:  (color online) Left panel: Temperature dependence of the
resistivity of Ba1-xKxFe2As2 ($x=0.2$) single crystal (red solid line). A
steep superconducting transition can be seen at $T_{c}=19$ K. The blue solid
squares are values from the zero frequency extrapolation of the optical
conductivity. The inset shows the derivative of the resistivity $d\rho/dT$ as
a function of temperature. The sharp peak at 104 K in $d\rho/dT$ is associated
with the SDW transition. The right panel depicts the same curves for $x=0.12$
sample with $T_{c}=11$ K, and $T_{\it SDW}=121$ K.

The left panel of Fig. 1 shows the temperature dependence of the DC
resistivity [$\rho(T)$] for the Ba1-xKxFe2As2 ($x=0.2$) sample. The $\rho(T)$
curve is characterized by a steep superconducting transition at $T_{c}$ = 19
K. The inset shows the derivative of the resistivity $d\rho/dT$ as a function
of temperature. The SDW transition manifests itself as a sharp peak in
$d\rho/dT$ at $T_{\it SDW}=104$ K, which corresponds to a small kink on the
$\rho(T)$ curve. The right panel displays the same curves for the $x=0.12$
sample, which has $T_{c}=11$ K, and $T_{\it SDW}=121$ K.

The _ab_ -plane reflectivity [$R(\omega)$] was measured at a near-normal angle
of incidence on Bruker IFS113v and IFS66v/s spectrometers. An _in situ_ gold
overfilling technique Homes et al. (1993) was used to obtain the absolute
reflectivity of the samples. Data from 20 to 12000$\leavevmode\nobreak\
\textrm{cm}^{-1}$ were collected at 18 different temperatures from 5 to 300 K
on freshly cleaved surfaces. In order to use Kramers-Kronig analysis, we
extended the data to the visible and UV range (10000 to
55000$\leavevmode\nobreak\ \textrm{cm}^{-1}$) at room temperature with an
AvaSpec-2048 $\times$ 14 model fiber optic spectrometer.

Figure 2 shows the infrared reflectivity at selected temperatures for both
samples up to 1200$\leavevmode\nobreak\ \textrm{cm}^{-1}$. The inset in each
panel displays the reflectivity for the full measured range at 300 K. For the
$x=0.2$ sample, shown in the top panel, the reflectivity exhibits a metallic
response and approaches unity at zero frequency. Below $T_{\it SDW}=104$ K, a
substantial suppression of $R(\omega)$ at about 650$\leavevmode\nobreak\
\textrm{cm}^{-1}$ sets in and intensifies with the decreasing temperature.

Figure 2:  (color online) Reflectivity of Ba1-xKxFe2As2 single crystals below
1200$\leavevmode\nobreak\ \textrm{cm}^{-1}$ at various temperatures. Top
panel: $x=0.2$; Bottom panel: $x=0.12$. Inset: Reflectivity of full measured
range at 300 K.

Simultaneously, the low frequency reflectivity continues increasing towards
unity. This is a signature of a partial SDW gap on the Fermi surface. Below 75
K, defined as $T^{\ast}$ here, another suppression of $R(\omega)$ appears in a
lower energy scale ($\sim 150\leavevmode\nobreak\ \textrm{cm}^{-1}$) signaling
the opening of a second partial gap (pseudogap) with a smaller value. Upon
crossing the superconducting transition, which occurs at $T_{c}$ = 19 K, the
reflectivity below $\sim 150\leavevmode\nobreak\ \textrm{cm}^{-1}$ increases
indicating the opening of a superconducting gap. Similar features are observed
on $R(\omega)$ for the $x=0.12$ sample as shown in the bottom panel of Fig. 2.

The real part of the optical conductivity $\sigma_{1}(\omega)$ was determined
by Kramers-Kronig analysis of the measured reflectivity.

Figure 3:  (color online) Top panel: Optical conductivity of Ba1-xKxFe2As2
($x=0.2$) at selected temperatures. The inset shows the enlarged view of the
optical conductivity at low frequencies. The bottom panel displays the same
spectra for the $x=0.12$ sample.

Figure 3 shows $\sigma_{1}(\omega)$ at different temperatures for the two
samples. The zero frequency extrapolations of $\sigma_{1}(\omega)$ represent
the inverse dc resistivity of the sample, shown as blue solid squares in Fig.
1, which are in good agreement with the transport measurement. The top panel
of Fig. 3 shows $\sigma_{1}(\omega)$ for Ba1-xKxFe2As2 ($x=0.2$) below
1700$\leavevmode\nobreak\ \textrm{cm}^{-1}$. At 150 K and 125 K, hence above
$T_{\it SDW}$, a Drude-like metallic response dominates the low frequency
optical conductivity. Below $T_{\it SDW}$, $\sigma_{1}(\omega)$ below about
650 $\leavevmode\nobreak\ \textrm{cm}^{-1}$ is severely suppressed. Meanwhile,
it increases in a higher energy scale from 650 $\leavevmode\nobreak\
\textrm{cm}^{-1}$ to 1700 $\leavevmode\nobreak\ \textrm{cm}^{-1}$. The optical
conductivity for the normal state and that for the SDW state just below
$T_{\it SDW}$ show an intersection point at about 650 $\leavevmode\nobreak\
\textrm{cm}^{-1}$. As the temperature decreases, both the low energy spectral
suppression and the high energy bulge become stronger; and the intersection
point moves to a higher energy scale. This spectral evolution manifests the
behavior of the SDW gap in this material: transfer of low frequency spectral
weight to high frequencies. If we take the intersection points as an
estimative of the gap values, we can see that the gap increases with
decreasing temperature. Below $T^{\ast}\sim 75$ K, a second suppression in the
optical conductivity below roughly 110$\leavevmode\nobreak\ \textrm{cm}^{-1}$
with a bulge extending from about 110$\leavevmode\nobreak\ \textrm{cm}^{-1}$
to 250$\leavevmode\nobreak\ \textrm{cm}^{-1}$ sets in and develops with the
temperature decrease, implying the opening of the pseudogap on the Fermi
surface. The inset of Fig. 3 shows the enlarged view of the low temperature
optical conductivity at low frequencies, where the pseudogap is seen more
clearly. This pseudogap is unlikely due to the SDW transition as (i) it opens
at 75 K, well below $T_{\it SDW}$ and (ii) it redistributes spectral weight at
a different, smaller energy scale.

The superconducting transition at $T_{c}$ = 19 K implies the opening of a
superconducting gap. As shown in the inset of Fig. 3, this leads to the
reduction of the optical conductivity at low frequencies between 20 K and 5 K.
The spectral weight lost in the transition is recovered by the
$\delta(\omega)$ function at zero frequency representing the infinite DC
conductivity in the superconducting state. This $\delta(\omega)$ function is
not visible in the $\sigma_{1}(\omega)$ spectra, because only finite
frequencies are experimentally measured. Nevertheless, its weight can be
calculated from the imaginary part of the optical conductivity. Zimmers et al.
(2004); Dordevic et al. (2002) Note that, the spectral depletion in
$\sigma_{1}(\omega)$ due to the superconducting condensate extends up to
180$\leavevmode\nobreak\ \textrm{cm}^{-1}$. This is the same energy scale of
the pseudogap, hinting that the superconducting gap and the pseudogap share
the same electronic states, and may have the same origin.

In the $x=0.12$ sample, very similar features are observed, as shown in the
bottom panel of Fig. 3, but remarkable differences exist: (i) the SDW gap
opens at a higher temperature ($T_{\it SDW}=121$ K) and the gap value shifts
to a higher energy scale ($\sim 750\leavevmode\nobreak\ \textrm{cm}^{-1}$);
(ii) The low frequency spectral suppression due to the SDW gap is stronger,
indicating that a larger part of the Fermi surface is removed in the SDW
state; (iii) In contrast to the SDW gap, both the pseudogap and the
superconducting gap features are weaker in the more underdoped sample. The
evolution of the three gaps (SDW, pseudogap and superconducting) with doping
also suggests that the pseudogap and the superconducting gap may have a common
origin while the SDW is a competitive order to superconductivity.

In order to investigate the origin of the pseudogap and the relationship among
all gaps, we analyzed the data utilizing a restricted spectral weight, defined
as:

$SW_{\omega_{a}}^{\omega_{b}}=\int_{\omega_{a}}^{\omega_{b}}\sigma_{1}(\omega)d\omega,$
(1)

where $\omega_{a}$ and $\omega_{b}$ are lower and upper cut-off frequencies,
respectively. By choosing appropriate values for $\omega_{a}$ and
$\omega_{b}$, one can study the relations among different phase transitions.
Yu et al. (2008) When replacing $\omega_{a}$ by 0 and $\omega_{b}$ by
$\infty$, we fall back to the standard $f$-sum rule and the spectral weight is
conserved.

Figure 4:  (color online) Temperature dependence of the spectral weight,
SW${}_{\omega_{a}}^{\omega_{b}}$ =
$\int_{\omega_{a}}^{\omega_{b}}\sigma_{1}(\omega)d\omega$, between different
lower and upper cutoff frequencies for the $x=0.2$ sample. The vertical dashed
lines denote $T_{c}$, $T^{\ast}$ and $T_{\it SDW}$.

Figure 4 shows the temperature dependence of the $x=0.2$ sample spectral
weight, normalized by its value at 300 K, at different cut-off frequencies.
The vertical dashed lines denote $T_{c}$, $T^{\ast}$ and $T_{\it SDW}$. The
blue solid circles in the top panel of Fig. 4 are the normalized spectral
weight with cut-off frequencies $\omega_{a}=0$ and
$\omega_{b}=12000\leavevmode\nobreak\ \textrm{cm}^{-1}$ as a function of
temperature. Here, the weight of the zero frequency $\delta$-function is
included below $T_{c}$. Moreover, since the optical conductivity is measured
only down to 20 $\leavevmode\nobreak\ \textrm{cm}^{-1}$, we estimate the
spectral weight below that energy by fitting the low frequency normal state
optical conductivity to a Drude model. The upper cut-off frequency
($\omega_{b}=12000\leavevmode\nobreak\ \textrm{cm}^{-1}$) is high enough to
cover the whole spectrum responsible for the phase transitions in this
material. Hence the blue solid circles form a flat line at about unity,
indicating that the spectral weight is conserved.

The red solid circles in the top panel show the temperature dependence of the
normalized spectral weight with cut-off frequencies $\omega_{a}=0^{+}$ and
$\omega_{b}=650\leavevmode\nobreak\ \textrm{cm}^{-1}$. Here $0^{+}$ means that
the superfluid weight is not included. Above $T_{\it SDW}$, the continuous
increase of the normalized SW${}_{0^{+}}^{650}$ with decreasing $T$ is related
to the narrowing of the Drude band. This is the typical optical response of a
metallic material. A strong spectral weight suppression occurs at $T_{\it
SDW}$, which is the consequence of the SDW gap opening. At $T_{c}$, another
sharp drop of the spectral weight breaks in, indicating the superconducting
gap opening.

The temperature dependence of the normalized SW${}_{650}^{1700}$, shown as
green solid circles in the bottom panel of Fig. 4, provides clues about the
relation between the superconducting and the SDW gaps. Above $T_{\it SDW}$,
the material shows a metallic response which can be described by a Drude peak
centered at zero frequency. With the temperature decrease, the DC conductivity
increases and the scattering rate reduces. The continuous narrowing of the
Drude band induces a transfer of spectral weight from the mid-infrared to the
far infrared, resulting in the continuous decrease of the spectral weight
observed in the 650–1700$\leavevmode\nobreak\ \textrm{cm}^{-1}$ range. Below
$T_{\it SDW}$, the opposite behavior dominates the optical conductivity. The
SDW gap depletes the spectral weight below 650$\leavevmode\nobreak\
\textrm{cm}^{-1}$ and transfers it to the 650–1700$\leavevmode\nobreak\
\textrm{cm}^{-1}$ range, leading to the continuous increase of
SW${}_{650}^{1700}$ with decreasing $T$. This behavior continues into the
superconducting state and does not show any feature at $T_{c}$. These
observations indicate that the SDW and superconducting gaps are separate and
even act as competitive orders in this material.

If a partial gap is due to a precursor order of superconductivity, for example
preformed pairs without phase coherence, once the long range superconductivity
is established, a significant part of the spectral weight transferred to high
frequencies by the partial gap should be transferred back to low energies and
join the superconducting condensate. Ioffe and Millis (1999, 2000) Whereas, a
partial gap due to a competitive order to superconductivity depletes the low-
energy spectral weight and holds it in a high energy scale without
transferring it back to the superfluid weight below $T_{c}$. Yu et al. (2008)
From the normalized SW${}_{650}^{1700}$ _vs_ $T$ curve (green solid circles)
we note that no loss of spectral weight is observed at $T_{c}$. This means
that the spectral weight transferred to high frequencies by the SDW gap
remains in the high frequency scale and does not contribute to the
superconducting condensate. Therefore, the SDW acts as a competitive order to
superconductivity in this material.

Along these lines, the origin of the pseudogap and its relationship to
superconductivity can be revealed by a close inspection of the temperature
dependence of the normalized SW${}_{110}^{250}$, shown as pink solid circles
in the bottom panel. Above $T^{*}$, this curve shows the same feature as the
normalized SW${}_{0^{+}}^{650}$ _vs_ $T$ curve, _i.e._ , continuous increase
upon cooling down followed by a suppression at $T_{\it SDW}$ due to the SDW
gap opening. At $T^{*}$, the spectral weight in the
110–250$\leavevmode\nobreak\ \textrm{cm}^{-1}$ range reaches a minimum and
starts to increase with decreasing temperature. This is due to the opening of
the pseudogap. The pseudogap, opening at $T^{*}$, depletes the spectral weight
below 110$\leavevmode\nobreak\ \textrm{cm}^{-1}$ and retrieves it in the
110–250$\leavevmode\nobreak\ \textrm{cm}^{-1}$ frequency range, leading to the
increase of SW${}_{110}^{250}$ below $T^{*}$. An interesting phenomenom
happens to the pseudogap when the material undergoes the superconducting
transition. In contrast to the case of the SDW gap, a significant loss of
spectral weight in the 110–250$\leavevmode\nobreak\ \textrm{cm}^{-1}$
frequency range is observed below $T_{c}$. This observation indicates that the
spectral weight transferred to the 110–250$\leavevmode\nobreak\
\textrm{cm}^{-1}$ range by the pseudogap joins the superconducting condensate
when superconductivity is established. Hence, the pseudogap is likely a
precursor order with respect to superconductivity.

In summary, we measured the optical conductivity of two underdoped
Ba1-xKxFe2As2 ($x=0.2$ and 0.12) single crystals. In both samples, besides the
SDW gap and superconducting gap, the optical conductivity reveals another
small partial gap (pseudogap) opening below $T^{\ast}\sim$ 75 K an
intermediate temperature between $T_{SDW}$ and $T_{c}$. A spectral weight
analysis shows that the SDW gap diminishes the low energy spectral weight
available for the superconducting condensate while the pseudogap shares the
same electronic states with the superconducting gap. These observations,
together with the doping dependence of these gaps, suggest the SDW as a
competitive order and the pseudogap as a precursor to superconductivity.

We thank Li Yu, Lei Shan, Cong Ren and Zhiguo Chen for helpful discussion.
Work in Paris was supported by the ANR under Grant No. BLAN07-1-183876
GAPSUPRA. Work in Beijing was supported by the National Science Foundation of
China (No. 91121004) and the Ministry of Science and Technology of China (973
Projects No. 2011CBA00107, No. 2012CB821400 and No. 2009CB929102). We
acknowledge the financial support from the Science and Technology Service of
the French Embassy in China.

## References

  * Kamihara et al. (2008) Y. Kamihara, T. Watanabe, M. Hirano, and H. Hosono, J.Am.Chem.Soc. 130, 3296 (2008).
  * Rotter et al. (2008a) M. Rotter, M. Tegel, and D. Johrendt, Phys. Rev. Lett. 101, 107006 (2008a).
  * Sefat et al. (2008) A. S. Sefat, R. Jin, M. A. McGuire, B. C. Sales, D. J. Singh, and D. Mandrus, Phys. Rev. Lett. 101, 117004 (2008).
  * Li et al. (2009) L. J. Li, Y. K. Luo, Q. B. Wang, H. Chen, Z. Ren, Q. Tao, Y. K. Li, X. Lin, M. He, Z. W. Zhu, et al., New Journal of Physics 11, 025008 (2009).
  * Torikachvili et al. (2008) M. S. Torikachvili, S. L. Bud’ko, N. Ni, and P. C. Canfield, Phys. Rev. Lett. 101, 057006 (2008).
  * Tapp et al. (2008) J. H. Tapp, Z. Tang, B. Lv, K. Sasmal, B. Lorenz, P. C. W. Chu, and A. M. Guloy, Phys. Rev. B 78, 060505 (2008).
  * Hsu et al. (2008) F.-C. Hsu, J.-Y. Luo, K.-W. Yeh, T.-K. Chen, T.-W. Huang, P. M. Wu, Y.-C. Lee, Y.-L. Huang, Y.-Y. Chu, D.-C. Yan, et al., Proc. Natl. Acad. Sci. 105, 14262 (2008).
  * Sales et al. (2009) B. C. Sales, A. S. Sefat, M. A. McGuire, R. Y. Jin, D. Mandrus, and Y. Mozharivskyj, Phys. Rev. B 79, 094521 (2009).
  * Rotter et al. (2008b) M. Rotter, M. Tegel, D. Johrendt, I. Schellenberg, W. Hermes, and R. Pöttgen, Phys. Rev. B 78, 020503 (2008b).
  * Hu et al. (2008) W. Z. Hu, J. Dong, G. Li, Z. Li, P. Zheng, G. F. Chen, J. L. Luo, and N. L. Wang, Phys. Rev. Lett. 101, 257005 (2008).
  * Akrap et al. (2009) A. Akrap, J. J. Tu, L. J. Li, G. H. Cao, Z. A. Xu, and C. C. Homes, Phys. Rev. B 80, 180502 (2009).
  * Lobo et al. (2010) R. P. S. M. Lobo, Y. M. Dai, U. Nagel, T. Rõ om, J. P. Carbotte, T. Timusk, A. Forget, and D. Colson, Phys. Rev. B 82, 100506 (2010).
  * Tu et al. (2010) J. J. Tu, J. Li, W. Liu, A. Punnoose, Y. Gong, Y. H. Ren, L. J. Li, G. H. Cao, Z. A. Xu, and C. C. Homes, Phys. Rev. B 82, 174509 (2010).
  * Teague et al. (2011) M. L. Teague, G. K. Drayna, G. P. Lockhart, P. Cheng, B. Shen, H.-H. Wen, and N.-C. Yeh, Phys. Rev. Lett. 106, 087004 (2011).
  * Terashima et al. (2009) K. Terashima, Y. Sekiba, J. H. Bowen, K. Nakayama, T. Kawahara, T. Sato, P. Richard, Y. M. Xu, L. J. Li, G. H. Cao, et al., Proc. Natl. Acad. Sci. 106, 7330 (2009).
  * Li et al. (2008) G. Li, W. Z. Hu, J. Dong, Z. Li, P. Zheng, G. F. Chen, J. L. Luo, and N. L. Wang, Phys. Rev. Lett. 101, 107004 (2008).
  * Shan et al. (2011) L. Shan, Y.-L. Wang, J. Gong, B. Shen, Y. Huang, H. Yang, C. Ren, and H.-H. Wen, Phys. Rev. B 83, 060510 (2011).
  * Ding et al. (2008) H. Ding, P. Richard, K. Nakayama, K. Sugawara, T. Arakane, Y. Sekiba, A. Takayama, S. Souma, T. Sato, T. Takahashi, et al., EPL (Europhysics Letters) 83, 47001 (2008).
  * Park et al. (2009) J. T. Park, D. S. Inosov, C. Niedermayer, G. L. Sun, D. Haug, N. B. Christensen, R. Dinnebier, A. V. Boris, A. J. Drew, L. Schulz, et al., Phys. Rev. Lett. 102, 117006 (2009).
  * Goko et al. (2009) T. Goko, A. A. Aczel, E. Baggio-Saitovitch, S. L. Bud’ko, P. C. Canfield, J. P. Carlo, G. F. Chen, P. Dai, A. C. Hamann, W. Z. Hu, et al., Phys. Rev. B 80, 024508 (2009).
  * Aczel et al. (2008) A. A. Aczel, E. Baggio-Saitovitch, S. L. Budko, P. C. Canfield, J. P. Carlo, G. F. Chen, P. Dai, T. Goko, W. Z. Hu, G. M. Luke, et al., Phys. Rev. B 78, 214503 (2008).
  * Massee et al. (2009) F. Massee, Y. Huang, R. Huisman, S. de Jong, J. B. Goedkoop, and M. S. Golden, Phys. Rev. B 79, 220517 (2009).
  * Chia et al. (2010) E. E. M. Chia, D. Talbayev, J.-X. Zhu, H. Q. Yuan, T. Park, J. D. Thompson, C. Panagopoulos, G. F. Chen, J. L. Luo, N. L. Wang, et al., Phys. Rev. Lett. 104, 027003 (2010).
  * Xu et al. (2011) Y.-M. Xu, P. Richard, K. Nakayama, T. Kawahara, Y. Sekiba, T. Qian, M. Neupane, S. Souma, T. Sato, T. Takahashi, et al., Nat Commun 2, 392 (2011).
  * Yu et al. (2008) L. Yu, D. Munzar, A. V. Boris, P. Yordanov, J. Chaloupka, T. Wolf, C. T. Lin, B. Keimer, and C. Bernhard, Phys. Rev. Lett. 100, 177004 (2008).
  * Luo et al. (2008) H. Luo, Z. Wang, H. Yang, P. Cheng, X. Zhu, and H.-H. Wen, Superconductor Science and Technology 21, 125014 (2008).
  * Homes et al. (1993) C. C. Homes, M. Reedyk, D. A. Cradles, and T. Timusk, Applied Optics 32, 2976 (1993).
  * Zimmers et al. (2004) A. Zimmers, R. P. S. M. Lobo, N. Bontemps, C. C. Homes, M. C. Barr, Y. Dagan, and R. L. Greene, Phys. Rev. B 70, 132502 (2004).
  * Dordevic et al. (2002) S. V. Dordevic, E. J. Singley, D. N. Basov, S. Komiya, Y. Ando, E. Bucher, C. C. Homes, and M. Strongin, Phys. Rev. B 65, 134511 (2002).
  * Ioffe and Millis (1999) L. B. Ioffe and A. J. Millis, Science 285, 1241 (1999).
  * Ioffe and Millis (2000) L. B. Ioffe and A. J. Millis, Phys. Rev. B 61, 9077 (2000).