$1.6$ in Zhou (2009): since the covariance function $\Sigma(\mathbf{z}_{i})$
is homogeneous and $\mathbf{z}_{i}$, $i=1,...,n$ are i.i.d. samples from
$f(\mathbf{z})$. Note that $x_{1}(z_{1}),...,x_{n}(z_{n})$ are i.i.d samples
form distribution $\int f(x\mid\Sigma(\mathbf{z}))f(\mathbf{z})$, which is
assumed to be sub-Gaussian by the Assumption T5. Thus we have
$P(\mathbf{X}\notin\mathcal{G}_{n,p})\leq\exp\\{-a_{1}n\\}$, for $n$ larger
than $a_{4}s_{j}^{*}\log p$, where $a_{4}$ and $a_{1}$ are positive constants,
by the similar argument with Lemma 3 in Atchadé (2019). Therefore, the
following restricted eigenvalue conditions hold:
$\underaccent{\tilde}{\kappa}(2s^{*},\mathbf{X}_{-j}^{T}\mathbf{X}_{-j}/n)$
and $\tilde{\kappa}(2s^{*},\mathbf{X}_{-j}^{T}\mathbf{X}_{-j}/n)$ are
constants.

For the second conclusion, first fix $i$, note that all eigenvalues of
$\Sigma$ are uniformly upper and lower bounded by constants. Then by Hanson-
Wright inequality (Rudelson & Vershynin, 2013), we have

$P(\|X_{i}\|_{2}\geq c_{1}\sqrt{n})\leq 2e^{-c_{2}n},$

for some constants $c_{1},c_{2}>0$. Then by the union bound and $n\geq
cs^{*}\log p$, by choosing large enough constant
$c_{1}^{\prime},c_{2}^{\prime}$, we have

$P(\|X_{i}\|_{2}\geq c_{1}^{\prime}\sqrt{n})\leq 2e^{-c_{2}^{\prime}n+\log
p}\leq 2e^{-c_{3}n}.$

∎

######  Lemma C.4.

Under assumptions of Theorem 3, we have

$\begin{split}\operatorname{\mathbf{E}}_{-j}\Bigg{[}P\bigg{(}\int(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))d\theta_{j}(\mathbf{z})\leq&-\alpha\int\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))d\theta_{j}(\mathbf{z})\\\
&+\mathrm{D_{KL}}(\hat{q}_{\theta_{j}}\|p_{\theta_{j}})+\log\left(\frac{1}{\zeta}\right)\bigg{)}\Bigg{]}\geq
1-\zeta.\end{split}$ (27)

###### Proof.

From (11), we have

$\displaystyle\operatorname{\mathbf{E}}_{-j}\operatorname{\mathbf{E}}_{j}\exp\left\\{\alpha\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}\right\\}=\exp\left(-(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*}\mid\mathbf{X}_{-j})\right).$

Thus, for any $\zeta\in(0,1)$, we have

$\displaystyle\operatorname{\mathbf{E}}_{-j}\operatorname{\mathbf{E}}_{j}\left[\exp\left\\{\alpha\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}+(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})-\log(1/\zeta)\right\\}\right]\leq\zeta.$

Integrating both sides of this inequality with respect to the prior
distribution $p_{\theta_{j}}$ and interchanging the integrals using Fubini’s
theorem, we have

$\displaystyle\operatorname{\mathbf{E}}_{-j}\operatorname{\mathbf{E}}_{j}\int\exp\left\\{\alpha\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}+(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})-\log(1/\zeta)\right\\}p_{\theta_{j}}(\theta_{j}(\mathbf{z}))d\theta_{j}(\mathbf{z})\leq\zeta.$

Next we use the following result from Yang et al. (2020).If $\mu$ is a
probability measure and $h$ is a measurable function such that $e^{h}\in
L_{1}(\mu)$, then

$\displaystyle\log\int e^{h}d\mu=\underset{\rho\ll\mu}{\sup}\left[\int hd\rho-
D_{KL}(\rho\|\mu)\right].$ (28)

We set
$h=\alpha\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}+(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})-\log(1/\zeta)$
and $\rho=\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))$ in the above result
where $\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))$ is the variational
estimate of the fractional posterior distribution.

Thus, we get

$\displaystyle\operatorname{\mathbf{E}}_{-j}\operatorname{\mathbf{E}}_{j}\exp\left[\int\left\\{\alpha\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}+(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})-\log(1/\zeta)\right\\}\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))d\theta_{j}-D_{KL}(\hat{q}_{\theta_{j}}\|p_{\theta_{j}})\right]\leq\zeta.$

Let

$\tilde{A}(\mathbf{X}_{-j})=\int\left\\{\alpha\log\frac{p(x_{j}\mid\theta_{j}(\mathbf{z}),\mathbf{X}_{-j})}{p(x_{j}\mid\theta_{j}^{*},\mathbf{X}_{-j})}+(1-\alpha)d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})-\log\frac{1}{\zeta}\right\\}\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))d\theta_{j}(\mathbf{z})-D_{KL}(\hat{q}_{\theta_{j}}\|p_{\theta_{j}}),$
(29)

and

$A(\mathbf{X}_{-j})=\left\\{x_{j}:\tilde{A}(\mathbf{X}_{-j})\leq 0\right\\}.$
(30)

Now we apply Markov’s inequality to get a probability statement. Thus, the
required statement follows from the following:

$\displaystyle\operatorname{\mathbf{E}}_{-j}P(\exp\left(\tilde{A}(\mathbf{X}_{-j})\right)>1\mid\mathbf{X}_{-j})$
$\displaystyle\leq$
$\displaystyle\operatorname{\mathbf{E}}_{-j}\operatorname{\mathbf{E}}_{j}(\exp\left(A(X_{-j})\right)\mid\mathbf{X}_{-j})=\zeta$
$\displaystyle\operatorname{\mathbf{E}}_{-j}P(A(\mathbf{X}_{-j})\mid\mathbf{X}_{-j})$
$\displaystyle\geq$ $\displaystyle 1-\zeta.$

∎

######  Lemma C.5.

Under the assumptions of Theorem 3, let the precision matrix of the
$p$-variate data generating distribution be $s^{*}$-sparse and have
eigenvalues bounded away from $0$ and $\infty$. For any $\zeta\in(0,1)$ and
$n\geq a_{1}s^{*}\log p$, and any measure $q_{\theta_{j}}\in\Gamma$ such that
$q_{\theta_{j}}\ll p_{\theta_{j}}$, we have

$P\left(\int\frac{1}{n}d_{\alpha}(\theta_{j}(\mathbf{z}),\theta_{j}^{*})\hat{q}_{\theta_{j}}(\theta_{j}(\mathbf{z}))d\theta_{j}(\mathbf{z})\leq\frac{\alpha}{n(1-\alpha)}\Psi(q_{\theta_{j}})+\frac{1}{n(1-\alpha)}\log(1/\zeta)\right)\geq
1-\zeta-\frac{c_{2}}{p^{c_{1}+1}}-\exp\\{-a_{2}n\\}.$

for some positive constants $a_{1}$, $D$ and $a_{2}$ , $c_{1},c_{2}$.

###### Proof.

Similar with Lemma C.2, we need to provide upper bound for

$2\|\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}\tilde{\beta}_{j}(\mathbf{z})-\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}{\beta}_{j}(\mathbf{z})\|^{2}+\|\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}\tilde{\beta}_{j}(\mathbf{z})-\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}{\beta^{*}_{j}}(\mathbf{z}_{l})\|^{2}+2\langle\mathrm{W}_{l}\mathbf{X}_{-j}(\tilde{\beta}_{j}(\mathbf{z})-\beta_{j}(\mathbf{z})),\varepsilon_{j}\rangle$

By the similar argument with Lemma 3 in Atchadé (2019), the following
restricted eigenvalue conditions hold:
$\underaccent{\tilde}{\kappa}(2s^{*},\mathbf{X}_{-j}^{T}\mathbf{X}_{-j}/n)$
and $\tilde{\kappa}(2s^{*},\mathbf{X}_{-j}^{T}\mathbf{X}_{-j}/n)$ are
constants. Therefore, with probability at least $1-\exp(-a_{1}n)$, we have
$\mathbf{X}\notin\mathcal{G}_{n,p}$, this gives us

$\|\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}\tilde{\beta}_{j}(\mathbf{z})-\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}{\beta^{*}_{j}}(\mathbf{z}_{l})\|^{2}\leq
cn\|\mathrm{W}\|_{2}\|\tilde{\beta}_{j}(\mathbf{z})-{\beta^{*}_{j}}(\mathbf{z}_{l})\|_{2}^{2}\leq
cn\tau^{-1}n^{-1}s_{j}^{*}\leq cn\epsilon^{2},$

given that $\tau^{-1}\leq c\sqrt{\log n}$.

In addition, by definition of
$\tilde{q}_{{\theta}_{j}}(\theta_{j}(\mathbf{z}))$, similarly we have

$2\|\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}\tilde{\beta}_{j}(\mathbf{z})-\mathrm{W}_{l}^{1/2}\mathbf{X}_{-j}{\beta}_{j}(\mathbf{z})\|^{2}\leq
cn\tau^{-1}\|\tilde{\beta}_{j}(\mathbf{z})-{\beta^{*}_{j}}(\mathbf{z}_{l})\|_{2}^{2}\leq
cn\epsilon^{2}.$

For the last term
$2\langle\mathrm{W}_{l}\mathbf{X}_{-j}(\tilde{\beta}_{j}(\mathbf{z})-\beta_{j}(\mathbf{z})),\varepsilon_{j}\rangle$,
first we have

$2\langle\mathrm{W}_{l}\mathbf{X}_{-j}(\tilde{\beta}_{j}(\mathbf{z})-\beta_{j}(\mathbf{z})),\varepsilon_{j}\rangle\leq
2\|\tilde{\beta}_{j}(\mathbf{z})-\beta_{j}(\mathbf{z})\|_{1}\|\mathbf{X}^{T}_{-j}\mathrm{W}_{l}\varepsilon_{j}\|_{\infty}\leq
c\sqrt{s_{j}^{*}}\tau\epsilon\|\mathbf{X}^{T}_{-j}\mathrm{W}_{l}\varepsilon_{j}\|_{\infty}.$

Note that $\mathbf{X}^{T}_{-j}\mathrm{W}_{l}\varepsilon_{j}$ is a $p-1$
dimensional Gaussian vector, with and with probability great than
$1-\exp(-a_{2}n)$, the scale of each component of the Gaussian vector is
bounded by $\sigma^{*}\sqrt{n}/\tau$ multiplied constant, by maximal
inequality of Gaussian random vector, we have

$P(\|\mathbf{X}^{T}_{-j}\mathrm{W}_{l}\varepsilon_{j}\|_{\infty}\geq t)\leq
e^{-\frac{\tau^{2}t^{2}}{2\sigma^{*2}n}+\log p},$

and we can choose $t=\sqrt{(c_{0}+2)n\log p}/\tau$ for a constant $c_{0}>0$,
then the probability upper bound becomes $p^{-(c_{0}+1)}$. Therefore, we have

$2\langle\mathrm{W}_{l}\mathbf{X}_{-j}(\tilde{\beta}_{j}(\mathbf{z})-\beta_{j}(\mathbf{z})),\varepsilon_{j}\rangle\leq
c\sqrt{s_{j}^{*}}\tau\epsilon\sqrt{n\log p}/\tau\leq
cn\epsilon^{2}+s_{j}^{*}\log p\leq cn\epsilon^{2},$

where the in the second inequality we use $2ab\leq a^{2}+b^{2}$. ∎

###  D Steps of the algorithm

The following algorithm provides the steps for covariate-dependent graph
estimation of $\mathbf{X}\in\mathbb{R}^{n\times p}$.

We select the bandwidth hyperparameter $\tau\in\mathbb{R}^{n}$ using a 2-step
approach for density estimation discussed in Dasgupta et al. (2020); Abramson
(1982); Van Kerm (2003). Under this approach, bandwidths are initialized using
Silverman’s rule of thumb, and the density is subsequently refined by updating
the bandwidth values. We follow this methodology to estimate the density of
$\mathbf{z}$, and use the updated bandwidths from the second step for $\tau$.

We next fix $x_{j}$ as the response and consider the task of performing $n$
weighted spike-and-slab regressions with $\mathbf{X}_{-j}$ as predictors using
the weights calculated using the bandwidth $\tau$ and the covariates. Each of
these regressions requires the specification of three hyperparameters:
$\pi,\sigma^{2}$, and $\sigma^{2}_{\theta}$. To select the hyperparameters, we
use a hybrid of model averaging and grid search. We first generate candidate
grids of $\pi,\sigma^{2}$, and $\sigma^{2}_{\theta}$ values. We denote the
grid of $\pi$ candidates by $\Theta_{\pi}$, and the Cartesian product between
the grid of $\sigma^{2}$ and $\sigma^{2}_{\theta}$ candidates as $\Theta$.

Next, for each $\pi\in\Theta_{\pi}$, we fit a spike-and-slab regression
weighted with respect to individual $l$ for each
$(\sigma^{2},\sigma^{2}_{\theta})\in\Theta,l\in 1,...,n$. We make a global
selection of $\sigma^{2}$ and $\sigma^{2}_{\theta}$ for each of the
$\pi\in\Theta_{\pi}$ such that the sum of the ELBO across all $n$ weighted
regressions is maximized. This grid search produces $\lvert\Theta_{\pi}\rvert$
models per individual. For each of these models, we calculate a model
averaging weight by taking the softmax over the ELBOs and use these to average
over the variational approximations to the posterior quantities to construct
the final model. Finally, to obtain the graph estimate, we symmetrize the
posterior inclusion probabilities from the final model and threshold them at
$0.5$.

###  E Discrete Covariate Simulation Study

For the discrete covariate, we perform experiments in which we vary the data
dimensionality, the distribution of the covariate levels, and the strength of
the signal in the ground-truth precision structures. As with the continuous
covariate, we perform 50 trials per experiment. We compare the performance of
W-PL to mgm (Haslbeck & Waldorp, 2020), as well as to the method of Carbonetto
et al. (2012) applied in a pseudo-likelihood fashion (CS). That is, we fix
each variable as the response in turn and perform a variational spike-and-slab
regression. We obtain the final graphs using the same symmetrization and
thresholding scheme as used for W-PL. To incorporate $\mathbf{z}$, we apply
this estimation procedure for each of the covariate levels independently.
Because no information is shared between levels, this allows us to evaluate
the impact of the weighting scheme in W-PL. We use the implementation of the
variational spike-and-slab from Carbonetto et al. (2017), which employs a
hybrid hyperparameter specification scheme, wherein the $\pi$ candidates are
averaged and $\sigma^{2}$ and $\sigma^{2}_{\theta}$ are selected via Empirical
Bayes for each of the $\pi$ candidates.

In each of the experiments, we assign individuals $1,...,n_{1}$ to the first
level of a binary discrete covariate $\mathbf{z}_{i}=1$, and the remaining
individuals to the second level $\mathbf{z}_{i}=2$. We refer to the structure
of the sample as balanced when $n_{1}=n-n_{1}:=n_{2}$, and unbalanced when
$n_{1}\neq n_{2}$.

####  E.1 Covariate-Independent Setting

We first consider a covariate-independent setting where the ground truth
dependence structure is independent of $\mathbf{z}_{i}$ and set $n=100,p=10$.
To construct the precision matrices, we first define

$\displaystyle\lambda_{\mathbf{z}_{i}}$ $\displaystyle=$ $\displaystyle{[c{\bf
1}_{4}\leavevmode\nobreak\ \leavevmode\nobreak\ ,\leavevmode\nobreak\
\leavevmode\nobreak\ {\bf 0}_{p-3}]}^{\mathrm{T}},\text{ for both
}\mathbf{z}_{i}=1,2,$

where ${\bf 1}_{4}$ is a four-dimensional vector of ones, and ${\bf 0}_{p-3}$
is a $(p-3)$-dimensional vector of zeroes. We refer to $c=15$ as high signal,
and $c=3$ as reduced signal. We next define the precision matrix for the
$i$-th individual as
$\Omega_{i}=(\lambda_{\mathbf{z}_{i}}{\lambda_{\mathbf{z}_{i}}}^{\mathrm{T}}+10\mathbb{I}_{p+1})$.
The corresponding dependence structure we aim to estimate is

$\mathrm{G}^{*}=\begin{bmatrix}\mathbb{J}_{4}-\mathbb{I}_{4}&{\bf
0}_{4,p-3}\\\ {\bf 0}_{p-3,4}&{\bf 0}_{p-3,p-3}\end{bmatrix}$

where $\mathbb{J}_{k}$ is a $k\times k$ matrix where all entries are $1$.

We perform experiments in the covariate-independent setting on high signal
with balanced structure, reduced signal with balanced structure, and high
signal with unbalanced structure ($n_{1}=80,n_{2}=20$). We present results for
each of the experiments in Table 3. When the signal strength is high and the
sample is balanced, all three methods correctly detect all of the edges in the
ground truth structure. However, the performance of both competitors suffers
under the unbalanced sample structure, particularly for CS, while W-PL
correctly detects all edges. In the reduced signal setting, the differential
between W-PL and the competitors grows significantly.

$c$ | $n_{1}$ | $n_{2}$ | Method | Sensitivity$(\uparrow)$ | Specificity$(\uparrow)$
---|---|---|---|---|---
|  |  | W-PL | $\mathbf{0.8517}(0.1591)$ | $0.9843(0.0168)$
|  |  | mgm | $0.1417(0.1229)$ | $\mathbf{0.9988}(0.0033)$
3 | 50 | 50 | CS | $0.2950(0.1489)$ | $0.9929(0.0081)$
|  |  | W-PL | $\mathbf{1.0000}(0.0000)$ | $0.9941(0.0111)$
|  |  | mgm | $\mathbf{1.0000}(0.0000)$ | $\mathbf{0.9990}(0.0031)$
15 | 50 | 50 | CS | $\mathbf{1.0000}(0.0000)$ | $0.9955(0.0080)$
|  |  | W-PL | $\mathbf{1.0000}(0.0000)$ | $\mathbf{1.0000}(0.0000)$
|  |  | mgm | $0.9013(0.0602)$ | $0.9992(0.0033)$
15 | 80 | 20 | CS | $0.8147(0.0215)$ | $0.9980(0.0052)$
Table 3: Results for the covariate-independent setting, presented as
$\textit{mean}(\textit{standard deviation})$

####  E.2 Covariate-Free Setting

We next examine a setting identical to the covariate-independent one, again
with $n=100,p=10$, however, this time, assume that no information on the
covariates is available. In the absence of covariate information, W-PL selects
all weights to be equal to one. Thus, the graph estimates are identical for
all the individuals in this setting, akin to the usual graph selection
algorithms. Because mgm requires the timepoints to be specified, we omit it
from these experiments.

We present the results for experiments in the covariate-free setting with high
and low signal in Table 4. Unsurprisingly, results for both methods are
similar. The minor differences in performance may be attributed to the
differing hyperparameter specification schemes.

$c$ | Method | Sensitivity$(\uparrow)$ | Specificity$(\uparrow)$
---|---|---|---
| W-PL | $\mathbf{0.9533}(0.1168)$ | $\mathbf{0.9996}(0.0029)$
3 | CS | $0.9233(0.1313)$ | $0.9939(0.0103)$
| W-PL | $\mathbf{1.0000}(0.0000)$ | $\mathbf{1.0000}(0.0000)$
15 | CS | $\mathbf{1.0000}(0.0000)$ | $0.9963(0.0079)$
Table 4: Results for covariate-free setting

####  E.3 Covariate-Dependent Setting

We next consider the setting in which the precision matrix varies with the
covariate level. We define the relationship as

$\displaystyle\lambda_{\mathbf{z}_{i}}$ $\displaystyle=$ $\displaystyle{[c{\bf
1}_{4}\leavevmode\nobreak\ \leavevmode\nobreak\ ,\leavevmode\nobreak\
\leavevmode\nobreak\ {\bf 0}_{p-3}]}^{\mathrm{T}},\text{ if
}\mathbf{z}_{i}=1,\text{ and}$ $\displaystyle\lambda_{\mathbf{z}_{i}}$
$\displaystyle=$ $\displaystyle{[{\bf 0}_{p-3}\leavevmode\nobreak\
\leavevmode\nobreak\ ,\leavevmode\nobreak\ \leavevmode\nobreak\ c{\bf
1}_{4}]}^{\mathrm{T}},\text{ if }\mathbf{z}_{i}=2.$

As before, we define the precision matrices as
$\Omega_{i}=(\lambda_{\mathbf{z}_{i}}{\lambda_{\mathbf{z}_{i}}}^{\mathrm{T}}+10\mathbb{I}_{p+1})$,
and thus, the true graph structure $\mathrm{G}^{*}(\mathbf{z})$ for an
individual with covariate value $\mathbf{z}$ is

$\mathrm{G}^{*}(1)=\begin{bmatrix}\mathbb{J}_{4}-\mathbb{I}_{4}&{\bf
0}_{4,p-3}\\\ {\bf 0}_{p-3,4}&{\bf
0}_{p-3,p-3}\end{bmatrix},\quad\mathrm{G}^{*}(2)=\begin{bmatrix}{\bf
0}_{p-3,p-3}&{\bf 0}_{p-3,4}\\\ {\bf
0}_{4,p-3}&\mathbb{J}_{4}-\mathbb{I}_{4}\end{bmatrix}.$

We visualize these precision matrices and the corresponding dependence
structures for $p=10$ in Figure 4.

|  |
---|---|---
|  |
Figure 4: Left to right: True dependence structures, estimated inclusion
probabilities and estimated graphs for a sample simulation with two discrete
covariate levels with p=10.

In addition to varying signal strength and sample structure for $p=10$, we
additionally vary the dimension of the data to $p=30$ and $p=50$ with high
signal strength and balanced samples. In all experiments, we fix $n=100$. We
present results for these experiments in Table 5.

While the performance of W-PL and mgm are similar for $p=10$, as $p$
increases, the performance of mgm deteriorates. On the other hand, W-PL and CS
demonstrate robustness to the increased sample size. As in the covariate-
independent setting, the performance of both mgm and CS is significantly
harmed relative to W-PL when faced with reduced signal.

$p$ | $c$ | $n_{1}$ | $n_{2}$ | Method | Sensitivity$(\uparrow)$ | Specificity$(\uparrow)$
---|---|---|---|---|---|---
|  |  |  | W-PL | $\mathbf{0.5800}(0.1859)$ | $0.9900(0.0133)$
|  |  |  | mgm | $0.0867(0.0937)$ | $\mathbf{0.9980}(0.0046)$
10 | 3 | 50 | 50 | CS | $0.2950(0.1545)$ | $0.9927(0.0090)$
|  |  |  | W-PL | $\mathbf{1.0000}(0.0000)$ | $0.9849(0.0187)$
|  |  |  | mgm | $0.9950(0.0354)$ | $\mathbf{0.9982}(0.0053)$
10 | 15 | 50 | 50 | CS | $\mathbf{1.0000}(0.0000)$ | $0.9953(0.0072)$
|  |  |  | W-PL | $0.7973(0.0189)$ | $0.9867(0.0106)$
|  |  |  | mgm | $\mathbf{0.8127}(0.0212)$ | $\mathbf{0.9990}(0.0028)$
10 | 15 | 80 | 20 | CS | $0.8093(0.0166)$ | $0.9983(0.0046)$
|  |  |  | W-PL | $\mathbf{1.0000}(0.0000)$ | $0.9926(0.0036)$
|  |  |  | mgm | $0.7567(0.1812)$ | $\mathbf{0.9997}(0.0006)$
30 | 15 | 50 | 50 | CS | $\mathbf{1.0000}(0.0000)$ | $0.9976(0.0018)$
|  |  |  | W-PL | $0.9867(0.0425)$ | $0.9958(0.0016)$
|  |  |  | mgm | $0.4550(0.2022)$ | $\mathbf{0.9999}(0.0002)$
50 | 15 | 50 | 50 | CS | $\mathbf{0.9983}(0.0118)$ | $0.9982(0.0009)$
Table 5: Results for discrete covariate-dependent setting

####  E.4 High-Dimensional Setting

Our last series of experiments with the discrete covariate deals with a
challenging high-dimensional setting where $p\geq n$. To handle the increased
dimensionality, we found it necessary to modify our hyperparameter
specification scheme for this experiment. We use Carbonetto et al. (2012) to
obtain an Empirical Bayes estimate to the hyperparameter $\sigma^{2}$, and use
a grid search to optimize over the hyperparameters $\pi$ and
$\sigma_{\theta}^{2}$ using the ELBO as our objective function. For the
bandwidth hyperparameter, we consider an ad-hoc choice of $\tau=0.1$. We only
perform 20 trials per experiment in this setting.

We maintain the relationship between the covariates and the ground truth
structure as in Section E.3 and first consider an unbalanced setting with
$p=100,n_{1}=40,n_{2}=10$ and high signal. We present results from this
setting in Figure 5, and exclude mgm from this experiment due to its
deteriorating performance with large $p$ and high time-complexity. Note that
as only $10$ observations belong to level 2, separate estimation through CS
suffers significantly compared to W-PL when estimating the graph for level 2.

Figure 5: Results for covariate-dependent setting with $n=50,p=100$. Top row:
Specificity CS vs W-PL; Bottom row: sensitivity CS vs W-PL

Next, to demonstrate how the signal-to-noise ratio (SNR) influences the
performance of our approach, we study several further experiments in the high-
dimensional setting with $n=50,p=50$ and $n=50,p=100$, $n_{1}=20,n_{2}=30$,
$\lambda=(c{\bf 1}_{4},{\bf 0}_{p-3})^{\mathrm{T}}$ keeping other settings the
same. Note that the SNR is controlled by $c$.

To assess performance under sparsity and with weak signal strength, we analyze
the area under the receiver operating characteristic curve (AUC). By varying
the threshold for a posterior inclusion probability to indicate an edge, we
obtain a sequence of true and false positive ratios that we may use to
calculate the corresponding AUC. AUC can also be defined by the fraction of
pairs that the prediction ordered correctly: let $y_{1},...,y_{n}$ be the $0$
and $1$ responses and $p_{1},...,p_{n}$ be the corresponding predicted
probabilities. The AUC can then be calculated as

${\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbbm{1}{\\{y_{i}<y_{j}\\}}\mathbbm{1}{\\{p_{i}<p_{j}\\}}}/\sum_{i=1}^{n}\sum_{j=1}^{n}\mathbbm{1}_{\\{y_{i}<y_{j}\\}}$

We present results from these experiments in Figure 6. As the signal strength
(i.e., $c$) increases, the AUC also increases from around $0.5$ to $1$. When
the SNR is low, W-PL does not work well and produces an AUC close to $0.5$,
which essentially is a random guess as to the presence of an edge. However,
when there are sufficient observations and the SNR is high, the AUC for level
2 exceeds $0.9$.

Because of the low level of observations in the first level of the covariate
($n_{1}=20$), separate estimation with CS does not perform well. W-PL
consistently outperforms CS for both level 1 and level 2.

(a) $n=50$, $p=50$

(b) $n=50$, $p=100$

Figure 6: CS versus W-PL method, measured by the AUC scores between the true
graph structure and the calculated inclusion probabilities.

###  F Departure from Gaussian assumption

The method theoretically is built on the assumption that the true data
generation is Gaussian, while the pseudo-likelihood approach is used mostly as
a tool for estimation. To study the effects of departures from Gaussianity, we
have investigated two scenarios. Firstly, we consider the situation where the
data is contaminated, that is, the data comes from a Gaussian distribution
with an independent structure, $c\%$ of which is contaminated by data coming
from an unrelated independent Gaussian distribution. Figure 7 shows the
results with $5\%$ contamination. The results, however, get worse as the
amount of contamination increases. Secondly, we consider the $t$-distribution
with varying degrees of freedom. Figure 8 shows the sensitivity and
specificity for varying degrees of freedom. The results indicate that for
degrees of freedom greater than $6$, the results are stable, and naturally
shows improvement as the degrees of freedom increases. However, for degrees of
freedom less than 6, the performance suffers, as shown in the left panel.

|  |
---|---|---
|  |
Figure 7: Specificity and sensitivity comparisons between CS VS W-PL for data with $2\%$ (left), $5\%$ (middle) and $10\%$ (right) contamination.  |  |
---|---|---
|  |
Figure 8: Specificity and sensitivity comparisons between CS VS W-PL for t
distributed data with df 3 (left), 6(middle) and 12(right), with $n=50$ in
both groups and $p=11$.

###  G Comparison to Qiu et al. (2016)

Here, we provide a brief comparison of W-PL to the method of Qiu et al.
(2016). Although their method is applicable in the continuous covariate
setting, their work focuses on the case when there are time replicates per
subject. It is possible to extend their method to the case where there are not
replicates by estimating the subject-level covariance matrices using a kernel-
weighted average, however, the implementation provided by the authors does not
include this functionality. Thus, we mainly focused on loggle and mgm as
competitors for W-PL in our experiments, since the available implementations
of these models could directly handle data without replicates, and similar to
Qiu et al. (2016), both are kernel-based models that frame the precision
matrix as varying continuously with time.

We also note that although there are methods for modeling heterogeneous
dependence structures other than those that model the dependence structure as
varying with time, such as Ren et al. (2022), we are not aware of any that are
directly applicable to a continuous covariate. For example, Ren et al. (2022)
assumes that the data may be grouped into clusters such that the precision
matrix is homogeneous within each of the clusters, and that the means of the
clusters are sufficiently separable from one another. On the other hand, W-PL
can model the precision matrices as varying continuously and does not place
any restrictions on the mean structure of the data.

For comparison to Qiu et al. (2016), we used their simulation setting with 2
(the minimum allowed) time replicates for 100 individuals with 10 variables.
We consider estimating the graph for subject 1. For W-PL, we only used the
information at the first time point, whereas for Qiu et al. (2016), we used
the full information on both time points. In this experiment, we revert to the
hyperparameter specification scheme from Section E.4. We summarize results
from this experiment in Figure 9. Notably, W-PL obtains superior results even
when the method of Qiu et al. (2016) technically uses double the number of
observations.

|
---|---
|
Figure 9: Qiu’s HQ method versus W-PL method. Top left shows the average TPR
and FPR for HQ method (squares) versus our method (diamond) across 50
iterations. Top right shows the true dependence structure of the individual
for a particular iteration, and the bottom left shows our estimate. Bottom
right shows HQ estimate.

###  H Ground-Truth Dependence Structures

####  H.1 Unidimensional Covariate

In the unidimensional covariate setting, we can split the individuals in three
clusters based on their covariates:
$\mathcal{C}_{1}=\\{i:-3<\mathbf{z}_{i}<-1\\},\mathcal{C}_{2}=\\{i:-1<\mathbf{z}_{i}<1\\}$
and $\mathcal{C}_{3}=\\{i:1<\mathbf{z}<3\\}$. Then, the ground truth precision
structures for each of the clusters is given by:

$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{1})=\begin{bmatrix}0&1&0&\\\
1&0&1&{\bf 0}_{3,p-2}\\\ 0&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$
$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{2})=\begin{bmatrix}0&1&1&\\\
1&0&1&{\bf 0}_{3,p-2}\\\ 1&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$
$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{3})=\begin{bmatrix}0&0&1&\\\
0&0&1&{\bf 0}_{3,p-2}\\\ 1&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$

####  H.2 Multidimensional Covariate

In the multidimensional covariate setting, we can split the individuals in
nine clusters based on their covariates:

$\displaystyle\mathcal{C}_{1}=\\{i:\mathbf{z}_{i}\in(-3,-1)\times(-3,-1)\\}$
$\displaystyle\mathcal{C}_{2}=\\{i:\mathbf{z}_{i}\in(-3,-1)\times(-1,1)\\}$
$\displaystyle\mathcal{C}_{3}=\\{i:\mathbf{z}_{i}\in(-3,-1)\times(1,3)\\}$
$\displaystyle\mathcal{C}_{4}=\\{i:\mathbf{z}_{i}\in(-1,1)\times(-3,-1)\\}$
$\displaystyle\mathcal{C}_{5}=\\{i:\mathbf{z}_{i}\in(-1,1)\times(-1,1)\\}$
$\displaystyle\mathcal{C}_{6}=\\{i:\mathbf{z}_{i}\in(-1,1)\times(1,3)\\}$
$\displaystyle\mathcal{C}_{7}=\\{i:\mathbf{z}_{i}\in(1,3)\times(-3,-1)\\}$
$\displaystyle\mathcal{C}_{8}=\\{i:\mathbf{z}_{i}\in(1,3)\times(-1,1)\\}$
$\displaystyle\mathcal{C}_{9}=\\{i:\mathbf{z}_{i}\in(1,3)\times(1,3)\\}$

Then, the ground truth precision structures for each of the clusters is given
by:

$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{1}\cup\mathcal{C}_{4})=\begin{bmatrix}0&1&0&\\\
1&0&1&{\bf 0}_{3,p-2}\\\ 0&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$
$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{2}\cup\mathcal{C}_{3}\cup\mathcal{C}_{5}\cup\mathcal{C}_{6})=\begin{bmatrix}0&1&1&\\\
1&0&1&{\bf 0}_{3,p-2}\\\ 1&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$
$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{7})=\begin{bmatrix}0&0&0&\\\
0&0&1&{\bf 0}_{3,p-2}\\\ 0&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$
$\displaystyle\mathrm{G}^{*}(\mathcal{C}_{8}\cup\mathcal{C}_{9})=\begin{bmatrix}0&0&1&\\\
0&0&1&{\bf 0}_{3,p-2}\\\ 1&1&0&\\\ &{\bf 0}_{p-2,3}&&{\bf
0}_{p-2,p-2}\end{bmatrix}$

## References

  * Abramson (1982) Abramson, I. S. (1982).  On bandwidth variation in kernel estimates-a square root law.  The annals of Statistics , 1217–1223.
  * Atchadé (2019) Atchadé, Y. F. (2019).  Quasi-bayesian estimation of large gaussian graphical models.  Journal of Multivariate Analysis 173, 656–671.
  * Besag (1975) Besag, J. (1975).  Statistical analysis of non-lattice data.  Journal of the Royal Statistical Society: Series D (The Statistician) 24, 179–195.
  * Besag (1977) Besag, J. (1977).  Efficiency of pseudolikelihood estimation for simple gaussian fields.  Biometrika , 616–618.
  * Bhadra & Mallick (2013) Bhadra, A. & Mallick, B. K. (2013).  Joint high-dimensional bayesian variable and covariance selection with an application to eqtl analysis.  Biometrics 69, 447–457.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A. & McAuliffe, J. D. (2017).  Variational inference: A review for statisticians.  Journal of the American Statistical Association 112, 859–877.
  * Bolli et al. (2014) Bolli, N., Avet-Loiseau, H., Wedge, D. C., Van Loo, P., Alexandrov, L. B., Martincorena, I., Dawson, K. J., Iorio, F., Nik-Zainal, S., Bignell, G. R. et al. (2014).  Heterogeneity of genomic evolution and mutational profiles in multiple myeloma.  Nature communications 5, 1–13.
  * Cai et al. (2013) Cai, T. T., Li, H., Liu, W. & Xie, J. (2013).  Covariate-adjusted precision matrix estimation with an application in genetical genomics.  Biometrika 100, 139–156.
  * Carbonetto et al. (2012) Carbonetto, P., Stephens, M. et al. (2012).  Scalable variational inference for bayesian variable selection in regression, and its accuracy in genetic association studies.  Bayesian analysis 7, 73–108.
  * Carbonetto et al. (2017) Carbonetto, P., Zhou, X. & Stephens, M. (2017).  varbvs: Fast Variable Selection for Large-scale Regression.
  * Castillo et al. (2012) Castillo, I., van der Vaart, A. et al. (2012).  Needles and straw in a haystack: Posterior concentration for possibly sparse sequences.  The Annals of Statistics 40, 2069–2101.
  * ChuaHL et al. (2007) ChuaHL, B., NakshatriP, C. et al. (2007).  Nf-kappab represses e-cadherin expression and enhances epithelial to mesenchymal transition of mammary epithelial cells: potential involvement of zeb-1 and zeb-2.  Oncogene 26, 711.
  * Criswell & Arteaga (2007) Criswell, T. L. & Arteaga, C. L. (2007).  Modulation of nf$\kappa$b activity and e-cadherin by the type iii transforming growth factor $\beta$ receptor regulates cell growth and motility.  Journal of Biological Chemistry 282, 32491–32500.
  * Csiszár & Talata (2006) Csiszár, I. & Talata, Z. (2006).  Consistent estimation of the basic neighborhood of markov random fields.  The Annals of Statistics , 123–145.
  * Danaher et al. (2014) Danaher, P., Wang, P. & Witten, D. M. (2014).  The joint graphical lasso for inverse covariance estimation across multiple classes.  Journal of the Royal Statistical Society. Series B, Statistical methodology 76, 373.
  * Dasgupta et al. (2020) Dasgupta, S., Pati, D. & Srivastava, A. (2020).  A two-step geometric framework for density modeling.  Statistica Sinica 30, 2155–2177.
  * Fox & Dunson (2015) Fox, E. B. & Dunson, D. B. (2015).  Bayesian nonparametric covariance regression.  The Journal of Machine Learning Research 16, 2501–2542.
  * Freno et al. (2009) Freno, A., Trentin, E. & Gori, M. (2009).  Scalable pseudo-likelihood estimation in hybrid random fields.  In Proceedings of the 15th ACM SIGKDD international conference on Knowledge discovery and data mining.
  * Friedman et al. (2008) Friedman, J., Hastie, T. & Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9, 432–441.
  * Giudici & Green (1999) Giudici, P. & Green, P. (1999).  Decomposable graphical gaussian model determination.  Biometrika 86, 785–801.
  * Guha et al. (2020) Guha, N., Baladandayuthapani, V. & Mallick, B. K. (2020).  Quantile graphical models: a bayesian approach.  The Journal of Machine Learning Research 21, 3023–3069.
  * Guo et al. (2011) Guo, J., Levina, E., Michailidis, G. & Zhu, J. (2011).  Joint estimation of multiple graphical models.  Biometrika 98, 1–15.
  * Ha et al. (2015) Ha, M. J., Baladandayuthapani, V. & Do, K.-A. (2015).  Dingo: differential network analysis in genomics.  Bioinformatics 31, 3413–3420.
  * Han et al. (2016) Han, S. W., Chen, G., Cheon, M.-S. & Zhong, H. (2016).  Estimation of directed acyclic graphs through two-stage adaptive lasso for gene network inference.  Journal of the American Statistical Association 111, 1004–1019.
  * Haslbeck & Waldorp (2020) Haslbeck, J. M. B. & Waldorp, L. J. (2020).  mgm: Estimating Time-Varying Mixed Graphical Models in High-Dimensional Data.  Journal of Statistical Software 93, 1–46.
  * Heckerman et al. (1995) Heckerman, D., Geiger, D. & Chickering, D. M. (1995).  Learning bayesian networks: The combination of knowledge and statistical data.  Machine learning 20, 197–243.
  * Hoff & Niu (2012) Hoff, P. D. & Niu, X. (2012).  A covariance regression model.  Statistica Sinica , 729–753.
  * Huang et al. (2016) Huang, X., Wang, J. & Liang, F. (2016).  A variational algorithm for bayesian variable selection.  arXiv preprint arXiv:1602.07640 .
  * Ji et al. (1996) Ji, C., Seymour, L. et al. (1996).  A consistent model selection procedure for markov random fields based on penalized pseudolikelihood.  The annals of applied probability 6, 423–443.
  * Jordan et al. (1999) Jordan, M. I., Ghahramani, Z., Jaakkola, T. S. & Saul, L. K. (1999).  An introduction to variational methods for graphical models.  Machine learning 37, 183–233.
  * Kleijn et al. (2006) Kleijn, B. J., van der Vaart, A. W. et al. (2006).  Misspecification in infinite-dimensional bayesian statistics.  The Annals of Statistics 34, 837–877.
  * Kolar et al. (2010a) Kolar, M., P Parikh, A. & P Xing, E. (2010a).  On sparse nonparametric conditional covariance selection .
  * Kolar et al. (2010b) Kolar, M., Song, L., Ahmed, A., Xing, E. P. et al. (2010b).  Estimating time-varying networks.  The Annals of Applied Statistics 4, 94–123.
  * Lee & Liu (2012) Lee, W. & Liu, Y. (2012).  Simultaneous multiple response regression and inverse covariance matrix estimation via penalized gaussian maximum likelihood.  Journal of multivariate analysis 111, 241–255.
  * Liu et al. (2010) Liu, H., Chen, X., Wasserman, L. & Lafferty, J. D. (2010).  Graph-valued regression.  In Advances in Neural Information Processing Systems.
  * Lohr et al. (2014) Lohr, J. G., Stojanov, P., Carter, S. L., Cruz-Gordillo, P., Lawrence, M. S., Auclair, D., Sougnez, C., Knoechel, B., Gould, J., Saksena, G. et al. (2014).  Widespread genetic heterogeneity in multiple myeloma: implications for targeted therapy.  Cancer cell 25, 91–101.
  * Mani et al. (2007) Mani, S. A., Yang, J., Brooks, M., Schwaninger, G., Zhou, A., Miura, N., Kutok, J. L., Hartwell, K., Richardson, A. L. & Weinberg, R. A. (2007).  Mesenchyme forkhead 1 (foxc2) plays a key role in metastasis and is associated with aggressive basal-like breast cancers.  Proceedings of the National Academy of Sciences 104, 10069–10074.
  * Meinshausen et al. (2006) Meinshausen, N., Bühlmann, P. et al. (2006).  High-dimensional graphs and variable selection with the lasso.  Annals of statistics 34, 1436–1462.
  * Ni et al. (2019) Ni, Y., Stingo, F. C. & Baladandayuthapani, V. (2019).  Bayesian graphical regression.  Journal of the American Statistical Association 114, 184–197.
  * Nieman et al. (1999) Nieman, M. T., Prudoff, R. S., Johnson, K. R. & Wheelock, M. J. (1999).  N-cadherin promotes motility in human breast cancer cells regardless of their e-cadherin expression.  The Journal of cell biology 147, 631–644.
  * Ormerod & Wand (2010) Ormerod, J. T. & Wand, M. P. (2010).  Explaining variational approximations.  The American Statistician 64, 140–153.
  * Pearl et al. (2000) Pearl, J. et al. (2000).  Models, reasoning and inference.  Cambridge, UK: CambridgeUniversityPress .
  * Pensar et al. (2017) Pensar, J., Nyman, H., Niiranen, J., Corander, J. et al. (2017).  Marginal pseudo-likelihood learning of discrete markov network structures.  Bayesian analysis 12, 1195–1215.
  * Peterson et al. (2015) Peterson, C., Stingo, F. C. & Vannucci, M. (2015).  Bayesian inference of multiple gaussian graphical models.  Journal of the American Statistical Association 110, 159–174.
  * Pourahmadi (1999) Pourahmadi, M. (1999).  Joint mean-covariance models with applications to longitudinal data: Unconstrained parameterisation.  Biometrika 86, 677–690.
  * Pourahmadi (2000) Pourahmadi, M. (2000).  Maximum likelihood estimation of generalised linear models for multivariate normal covariance matrix.  Biometrika 87, 425–435.
  * Pourahmadi (2013) Pourahmadi, M. (2013).  High-dimensional covariance estimation: with high-dimensional data, vol. 882.  John Wiley & Sons.
  * Qiu et al. (2016) Qiu, H., Han, F., Liu, H. & Caffo, B. (2016).  Joint estimation of multiple graphical models from high dimensional time series.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 78, 487–504.
  * Raskutti et al. (2010) Raskutti, G., Wainwright, M. J. & Yu, B. (2010).  Restricted eigenvalue properties for correlated gaussian designs.  The Journal of Machine Learning Research 11, 2241–2259.
  * Ren et al. (2022) Ren, M., Zhang, S., Zhang, Q. & Ma, S. (2022).  Gaussian graphical model-based heterogeneity analysis via penalized fusion.  Biometrics 78, 524–535.
  * Rudelson & Vershynin (2013) Rudelson, M. & Vershynin, R. (2013).  Hanson-wright inequality and sub-gaussian concentration.  Electronic Communications in Probability 18, 1–9.
  * Van Kerm (2003) Van Kerm, P. (2003).  Adaptive kernel density estimation.  The Stata Journal 3, 148–156.
  * Wainwright & Jordan (2008) Wainwright, M. J. & Jordan, M. I. (2008).  Graphical models, exponential families, and variational inference.  Now Publishers Inc.
  * Wang & Kolar (2014) Wang, J. & Kolar, M. (2014).  Inference for sparse conditional precision matrices.  arXiv preprint arXiv:1412.7638 .
  * Yang & Peng (2020) Yang, J. & Peng, J. (2020).  Estimating Time-Varying Graphical Models.  Journal of Computational and Graphical Statistics 29, 191–202.
  * Yang et al. (2020) Yang, Y., Pati, D., Bhattacharya, A. et al. (2020).  $\alpha$-variational inference with statistical guarantees.  Annals of Statistics 48, 886–905.
  * Yin & Li (2011) Yin, J. & Li, H. (2011).  A sparse conditional gaussian graphical model for analysis of genetical genomics data.  The annals of applied statistics 5, 2630.
  * Yuan & Lin (2007) Yuan, M. & Lin, Y. (2007).  Model selection and estimation in the gaussian graphical model.  Biometrika 94, 19–35.
  * Zhang & Leng (2012) Zhang, W. & Leng, C. (2012).  A moving average cholesky factor model in covariance modelling for longitudinal data.  Biometrika 99, 141–150.
  * Zhou (2009) Zhou, S. (2009).  Restricted eigenvalue conditions on subgaussian random matrices.  arXiv preprint arXiv:0912.4045 .
  * Zhou et al. (2010) Zhou, S., Lafferty, J. & Wasserman, L. (2010).  Time varying undirected graphs.  Machine Learning 80, 295–319.


# An Empirical Determination of the Intergalactic Background Light from UV to
FIR Wavelengths Using FIR Deep Galaxy Surveys and the Gamma-ray Opacity of the
Universe

Floyd W. Stecker  Astrophysics Science Division, NASA/Goddard Space Flight
Center  Greenbelt, MD 20771<EMAIL_ADDRESS>Department of Physics
and Astronomy, University of California, Los Angeles  Sean T. Scully
Department of Physics, James Madison University  Harrisonburg, VA 22807;
<EMAIL_ADDRESS>Matthew A. Malkan  Department of Physics and Astronomy,
University of California, Los Angeles  Los Angeles, CA90095-1547;
<EMAIL_ADDRESS>
###### Abstract

We have previously calculated the intergalactic background light (IBL) as a
function of redshift from the Lyman limit in the far ultraviolet to a
wavelength of 5 $\mu$m near infrared range, based purely on data from deep
galaxy surveys. Here we utilize similar methods to determine the mid- and far-
infrared IBL from 5 $\mu$m to 850 $\mu$m. Our approach enables us to constrain
the range of photon densities, by determining the uncertainties in
observationally determined luminosity densities and spectral gradients. By
also including the effect of the 2.7 K cosmic background photons, we determine
upper and lower limits on the opacity of the universe to $\gamma$-rays up to
PeV energies within a 68% confidence band.

Our direct results on the IBL are consistent with those from complimentary
$\gamma$-ray analyses using observations from the Fermi $\gamma$-ray space
telescope and the H.E.S.S. air Čerenkov telescope. Thus, we find no evidence
of previously suggested processes for the modification of $\gamma$-ray spectra
other than that of absorption by pair production alone.

diffuse radiation – galaxies: observations – gamma-rays: theory

##  1 Introduction

We have previously employed a fully empirical approach to calculating the
intergalactic background light (IBL) to wavelengths up to 5 $\mu$m by using
observational data from deep galaxy surveys (Stecker, Malkan & Scully 2012
(hereafter SMS12); and Scully, Malkan & Stecker 2014 (hereafter SMS14)); for a
similar approach, see also Helgason & Kashlinsky (2012) and Khaire & Srianand
(2015). Here we extend our previous results from SMS12 and SMS14, determining
both the intergalactic background light (IBL) at longer wavelengths out to 850
$\mu$m and the subsequent $\gamma$-ray opacity of the Universe out to multi-
TeV energies. The spectra of astrophysical sources of such high energy
$\gamma$-rays are being studied by ground-based air Čerenkov telescopes and
linked air Čerenkov telescope arrays.

We accomplish the goals of this paper by using very recent deep galaxy survey
data at far infrared wavelengths where galaxy emission is produced by dust re-
radiation rather than starlight. We also include the $\gamma$-ray opacity from
photons of the 2.7 K cosmic background radiation (CBR). This enables us to
extend our calculations of the $\gamma$-ray opacity of the Universe to
energies much greater than a TeV.

Observations at wavelengths greater than 24 $\mu$m, have been covered by the
Multiband Imaging Photometer on the Spitzer space telescope (MIPS), now being
dramatically advanced by the availability of data from the Photoconductor
Array Camera and Spectrometer (PACS) and Spectral Photometric Imaging Receiver
(SPIRE) instruments on the Herschel space telescope, the Planck space
telescope, and by ground-based observations from the Atacama Large Millimeter
Array (ALMA) and observations from the Balloon-borne Large Aperture
Submillimeter Telescope (BLAST). By using only empirical data rather that
modeling, our approach is, by definition, model independent. We use published
luminosity functions and interpolations of luminosity densities between
observed wavebands, observations from high-redshift galaxy surveys being now
sufficient for this purpose.

Most past approaches to determining the IBL and its present-day spectral
energy distribution, referred to as the ”extralgalactic background light
(EBL)”, require assumptions about how the galaxy luminosity functions evolve
(e.g. Malkan & Stecker 1998, 2001; Kneiske, Mannheim & Hartmann, 2002:
Stecker, Malkan & Scully 2006 (SMS06); Franceschini et al. 2008; Finke,
Razzaque, & Dermer, 2010; Kneiske & Dole 2010). Other approaches make use of
semi-analytic models that require various assumptions regarding galaxy
evolution, stellar population synthesis modeling, or star formation rates the
properties of and dust attenuation, particularly for redshifts greater than 1
(e.g., Salamon & Stecker 1998; Gilmore et al. 2009; Somerville et al. 2012;
Dominguez et al. 2011; Inoue et al. 2013; Driver et al. 2016). In contrast,
our observationally-based approach is superior to model-based methods, since
it enables a determination of both the IBL and its observational
uncertainties. This is because we use observationally determined errors.
Approaches that rely on modeling cannot determine such uncertainties.

Observational studies of blazar $\gamma$-ray spectra have also been used to
probe the IBL (Ackermann et al. 2012; Abramowski et al. 2013; Biteau &
Williams 2015) through its opacity effect caused by electron-positron pair
production. This approach was originally suggested by Stecker, De Jager &
Salamon (1992) when the infrared EBL was unknown. Our present method of
determining the expected $\gamma$-ray opacity of the Universe by first using
observational data to determine the IBL and then calculating its effect on
$\gamma$-ray spectra of extragalactic sources complements the technique of
using direct $\gamma$-ray observations. This is because the later approach
requires a thorough knowledge of the intrinsic (unabsorbed) emission spectra
of the $\gamma$-ray sources. This requirement introduces unavoidable
uncertainties.

Our final results give the $\gamma$-ray opacity as a function of energy and
redshift to within a 68% confidence band that is based on observational data.
This opacity is solely caused by pair-production interactions of IBL photons
with extragalactic $\gamma$-rays . Thus, a direct comparison of this effect
with $\gamma$-ray spectra of extragalactic sources enables an assessment of
possible additional spectral modifications. One such possible spectral
modification has been suggested to be caused by secondary $\gamma$-ray
production from cosmic-ray interactions along the line of sight to the source
(Essey & Kusenko 2014). Another such modification might be caused by photon-
axion oscillations during propagation from the source to the Earth (e.g., De
Angelis et al. 2007; Mayer & Horns 2013).

In Section 2 we give the details of our calculations of upper and lower limits
on IR galaxy luminosity densities for wavelengths between 5 $\mu$m and 850
$\mu$m as a function of redshift, paying particular attention to the effects
of the emission spectra of polycyclic aromatic hydrocarbons (PAHs). In Section
3 we calculate the $z=0$ EBL based on the results from Section 2 and compare
it with both other work and observational limits. In Section 4 we compute the
resulting opacity of the universe to $\gamma$-rays out to a redshift of z = 5,
including the total opacity from $\gamma+\gamma\rightarrow e^{+}+e^{-}$
interactions with both IBL and thermal cosmic background radiation (CBR)
photons. In our conclusion section, Section 5, we compare our results with
those obtained by other methods in some of the papers mentioned above. We also
discuss the implications of our results.

##  2 Calculating the Infrared IBL

We have previously calculated the intergalactic background light (IBL) as a
function of redshift in the far ultraviolet to near infrared range, based
purely on data from deep galaxy surveys (SMS12; SMS14). Here we utilize
similar methods to extend the calculation into the mid- and far infrared IBL
out to a wavelength of 850 $\mu$m.

###  2.1 Determining the Luminosity Densities from Empirical Luminosity
Functions

In our previous work the observationally tolerated ranges of photon densities
were determined from the luminosity densities (LDs), $\rho_{L_{\nu}}$,
themselves, with errors provided by the various authors at wavelengths ranging
from the far UV to the 5 $\mu$m wavelength in the near IR. The LDs are
computed by integrating fits to the observationally determined luminosity
functions:

$\rho_{L_{\nu}}=\int_{L_{min}}^{L_{max}}dL_{\nu}\,L_{\nu}\Phi(L_{\nu};z)$ (1)

Cases where authors did not directly compute the LDs were excluded because
properly estimating their error requires knowledge of the covariance of the
errors in the fit parameters of the luminosity functions (LFs) and also
knowledge of any observational biases. To provide comprehensive redshift
coverage of the LDs, SMS12 and SMS14 made use of continuum colors between the
wavelength bands to fill in any gaps.

At wavelengths greater than 5 $\mu$m very few studies provide determinations
of luminosity densities (LDs) directly, as most authors are more concerned
with calculating the total IR density, integrated over all wavelengths. This
is because the total IR density is an observable that is correlated with the
star formation rate. Therefore, in this paper we used observer-given analytic
fits to the LFs at various wavelengths. When those fits were not provided, we
ourselves fit the observed infrared galaxy LFs at various wavelengths and
redshifts and then use equation (1) to obtain the LDs.

Galaxy LFs typically have characteristic shapes that are flatter at lower
luminosities but fall off more steeply at the highest luminosities. However,
at wavelengths greater than 5 $\mu$m the Schechter function, commonly used at
optical wavelengths, does not provide a good fit to the observed LFs, as its
exponential decrease falls off too quickly at the bright end compared with the
measured LFs. Most observers instead fit their LFs to a broken power-law
function that describes the data better (e.g., Marchetti et al. 2016). The
transition region between the power law that holds at low luminosities and
that holding at high luminosities is usually referred to as the knee of the
LF. The luminosity at the knee is usually designated as $L_{*}$. Most of the
luminosity density from an LF is produced by galaxies with luminosities in the
vicinity of the knee; much fainter galaxies do not contribute much to the LD
and the much brighter ones are quite rare. Thus, it is critical to determine
the location of the knee in order to estimate the LDs with any accuracy.

For the cases where we were required to determine our own fits, we chose a
double power-law fitting function of the form:

$\Phi(L)=\frac{c}{{L_{*}\left({{{\left({\frac{L}{L_{*}}}\right)}^{a}}+{{\left({\frac{L}{L_{*}}}\right)}^{b}}}\right)}}$
(2)

Here, the parameters $a$ and $b$ are the indices of the power-law fits in the
low luminosity and high luminosity ranges respectively. The overall
normalization of the LF is given by the parameter $c$, which has the dimension
of luminosity per Mpc3dex. To compute the LDs, we integrate equation (1) over
the galaxy luminosities between 4 $\times$ 107 L${}_{{}_{\mathord{\odot}}}$
and 1014 L${}_{{}_{\mathord{\odot}}}$.

In order to accurately compute the errors in the LDs as derived from the fit
parameters, which is essential for our calculation, we must do so in a way
that includes terms involving the off-diagonal elements of the error matrix of
the fits so as to account for covariance. In cases where the authors provided
there own fits, we have therefore re-derived the fits of the LFs to generate
the error matrix, retaining the same choice of fitting function and fixed
parameters, if given.

Our goal, as in SMS12 and SSM14, was to compute the observationally tolerated
ranges of luminosity densities. This requires that we represent the error on
these quantities as best we can. The statistical error is determined by
properly propagating the fit parameter errors accounting for covariance, lest
we overestimate the error. To compute the total error on each LD, we further
added in quadrature an additional systematic error that accounts for cosmic
variance. We compute cosmic variance based on the field sizes of the
individual studies. For the AKARI Wide Field IR Survey Explorer (WISE) and
GOODS fields this value is typically of order of 10%. In determining the LDs,
for consistency, all observational LFs are scaled to a Hubble parameter of
$h=0.7$

Our new additions to our observationally determined LDs extend our coverage of
rest frame galaxy photon production from the near IR to 850 $\mu$m in the far
IR, with enough determinations at each wavelength band to span the redshift
range $0\leq z\leq 2-3$. Using published results derived from observations by
the Spitzer and Herschel space telescopes, sufficient redshift coverage was
found for wavelength bands of $8,12,15,24,35,60,90,$ and $250$ $\mu$m .

There were two cases where we were required to combine LFs from different
observational studies in order to provide enough coverage to discriminate the
location of $L_{*}$ particularly for redshifts greater than $1$. At 12 $\mu$m
we combined the results of Perez-Gonzalez et al. (2005) and Rodighiero et al.
(2010) to compute LDs in redshift bins centered on redshifts of $1.2,1.6,$ and
$2.0$. In order to achieve sufficient redshift coverage at $90$ $\mu$m we
combined some of the higher redshift $100$ $\mu$m data from Lapi et al. (2011)
with that of the $90$ $\mu$m data from Gruppioni et al. (2013). This gains us
additional coverage at redshifts of $1.4,2.2,$ and $3.0$.

At 160, 350, 500, and 850 $\mu$m , LF data only exists for the very nearby
redshifts. We therefore assumed that their redshift evolution closely follows
that of the 250 $\mu$m band (Lapi et al. 2011; Marchetti et al. 2016). We use
local LDs calculated in those bands as a normalization to this evolution. This
assumption is justified because the emission in this wavelength region is
dominated by warm dust. At 160 $\mu$m we used the local LF given be Patel et
al. (2013). At 350 and 500 $\mu$m , we used the combined local LF data of
Marchetti et al. (2016) from Herschel/SPIRE, the Herschel/SPIRE estimate of
Vaccari et al. (2010), and the Planck satellite data from Negrello et al.
(2013). At 850 $\mu$m we computed the local LD from the LF provided by
Negrello et al. (2013). Figures 1 and 2 show the resulting derived values for
$\rho_{L_{\lambda}}(z)$ and their errors together with the observationally
determined $\pm$ 1 $\sigma$ confidence bands for the 8, 12, 15, 24 $\mu$m
bands and those of 35, 60, 90, and 250 $\mu$m respectively.

Our calculations of LDs at mid-IR to far-IR wavelengths presented here is an
extension of the work done in SMS12 and SSM14. Thus, the results given in
those papers for wavelengths less than 5 $\mu$m is almost unchanged from the
results presented here. However, we have updated the far UV calculations to
include the more recent work from Bouwens et al. (2015) for LDs in the
redshift range of $4$ to $7$. Even though the shape of the far UV band can
affect other bands when filling in for redshift gaps using colors, the overall
calculation yields results that are qualitatively the same as those presented
in SMS12 and SSM14, as the newer data do not significantly change the general
trend.

In order to place 68% upper and lower limits by using observational data on
$\rho_{L_{\nu}}$ we make as few assumptions about the luminosity density
evolution as possible. In SMS12 and SMS14 we utilized a robust rational
fitting function in the form of a broken power-law dependent on $(1+z)$ in
order to generate the confidence bands. We take the 68% confidence ranges of
these fits in each waveband as the $\pm$1 $\sigma$ confidence bands for the
LDs. At redshifts beyond the redshift at the peak of star formation, the LDs
decline with redshift. In accord with recent studies of the evolution of rest
frame LFs in the UV that trace the star formation rate (Finkelstein et al.
2015), we conservatively assume upper and lower limit power-law functions in
redshift to represent the rate of this decline as the highest redshifts. In
the upper limit case, we assume a decline proportional to $(1+z)^{-2}$. In the
lower limit case, we adopt a steeper decline proportional to $(1+z)^{-4}$.
These assumptions have almost no impact on the derivation of the opacity
confidence bands that we determined. Figures 1 and 2 show our results for LDs
as a function of redshift at various wavelengths.

###  2.2 Taking account of PAH emission

At wavelengths greater than $\sim$20 $\mu$m the LDs between the bands can be
determined by smoothly interpolating between our observationally based LDs at
specific wavelengths, since the spectral energy distributions (SEDs) from
galaxies in this wavelength range are smooth modified blackbody spectra
produced by dust re-radiation. However, in the 5 – 20 $\mu$m range the
situation is more complex.

In star-forming galaxies, the spectra between 7 and 13 $\mu$m are dominated
primarily by PAH emission. These PAH molecules are found in very small dust
grains in intergalactic media. They absorb the UV photons emitted by hot young
$O$ and $B$ stars and reemit them in molecular emission bands in the mid-IR.
They are thus a strong signature of active star formation in galaxies
(Peeters, Spoon & Tielens 2004). In this regard, we note that luminous star
forming galaxies at higher redshift have more prominent PAH emission features.
The importance of PAH features in the mid-IR at redshifts $0.5\leq z\leq 2.5$
has been shown by Lagache et al. (2005).

The average SED of nearby star-forming galaxies (from Spoon et al. 2007, see
also Smith et al. 2007) is shown in Figure 3 normalized to our best fit low
redshift LD confidence band. One can see that there is a relative ”valley”
between 9 and 11 $\mu$m. A simple direct interpolation between our 8 and 12
$\mu$m bands would therefore obtain an incorrectly high value for the LD in
this wavelength range. Since we do not have wavelength coverage in this
regime, we take this feature into account by lowering our interpolated LDs by
a factor of 3 for the upper limit and a factor of 5 for the lower limit at 10
$\mu$m . The factor of 5 is chosen as a lower limit based on the difference
between the peak at 8 $\mu$m and the depth of the valley near 10 $\mu$m from
the SED. For the upper limit, we relax the depth of this feature because,
while star forming galaxies make up the bulk of the IBL, contributions from
other galaxy types have a less pronounced PAH feature. Since, at high
redshifts, PAH emission correlates with the star formation rate (Shipley et
al. 2016), we assume that our relative PAH shape factors of 3 and 5 coevolve
with redshift.

##  3 The Extragalactic Background Light

At wavelengths above $\sim$ 100 $\mu$m the FIRAS and DIRBE instruments aboard
the Cosmic Background Explorer (COBE) have measured the EBL (Fixsen et al.
1998; Lagache et al. 1999). This diffuse background has now been in large part
resolved by recent galaxy count studies using the Herschel telescope (Berta et
al. 2010; Béthermin et al. 2012; Viero et al. 2015) along with ground based
studies using ALMA (Fujimoto et al. 2016) and BLAST (Marsden et al. 2011).
These more recent studies strongly support the COBE results.

However, there are no direct measurements of the EBL in the infrared range
below $\sim$ 100 $\mu$m owing to the predominance of the foreground radiation
of Zodiacal light from interplanetary dust re-radiation. The flux of Zodiacal
light is approximately two orders-of-magnitude larger than the EBL flux
(Spiesman et al. 1995). In this region only lower limits obtained from galaxy
counts exist.

Using our calculations of the LDs in the mid-IR and far-IR as a function of
wavelength, we have constructed a 68% confidence band for the spectral energy
distribution of the cosmic diffuse infrared background light (the IBL at
$z=0$). Figure 4 shows our new results, combined with our previous results at
shorter wavelengths, taken from SMS12 and SMS14. The light shaded band shows
the maximum effect of PAH emission. It can be seen that taking account of the
details of the PAH spectrum does not significantly affect our EBL results.
Figure 4 also shows the observational lower limits on the EBL obtained from
galaxy counts (in blue), extrapolations of mid-IR galaxy counts from Spitzer,
and direct measurements (in black).

##  4 The Optical Depth from $\gamma+\gamma\rightarrow e^{+}+e^{-}$
Interactions with IBL and 2.7 K CBR Photons

The co-moving radiation energy density for wavelength $\lambda$ at redshift
$z$, $u_{\nu}(z)$, where $\nu=c/\lambda$, is the time integral of the co-
moving luminosity density ${\rho}_{\nu}(z)$,

$u_{\nu}(z)=\int_{z}^{z_{\rm
max}}dz^{\prime}\,{\rho}_{\nu^{\prime}}(z^{\prime})\frac{dt}{dz}{(z^{\prime})},$
(3)

where $\nu^{\prime}=\nu(1+z^{\prime})/(1+z)$ and $z_{\rm max}$ is the redshift
corresponding to initial galaxy formation (Salamon & Stecker 1998), and

$\frac{dt}{dz}{(z)}={[H_{0}(1+z)\sqrt{\Omega_{\Lambda}+\Omega_{m}(1+z)^{3}}}]^{-1},$
(4)

with $\Omega_{\Lambda}=0.72$ and $\Omega_{m}=0.28$.

The upper and lower limits on our co-moving energy densities, derived using
equation (3), are shown in Figures 5 and 6.

In calculating the $\gamma$-ray opacities we use the relations for the photon
energy, $\epsilon_{\nu}=h\nu$, and the photon density,
$n_{\nu}=\rho_{\nu}/\epsilon$.

The cross section for photon-photon annihilation to electron-positron pairs
was first calculated by Breit & Wheeler (1934) as a solid result of quantum
electrodynamics. The threshold for this interaction follows from the Lorentz
invariance of the square of the four-momentum vector. This reduces to the
square of the threshold energy in the c.m.s., $s$, that is necessary to
produce twice the electron rest mass:

$s=2\epsilon E_{\gamma}(1-\cos\theta)=4m_{e}^{2}$ (5)

The Lorentz invariant quantity given by $s$ has been determined to hold to
better than one part in $10^{15}$ (Stecker & Glashow 2001; Jacobson, Liberati,
Mattingly & Stecker 2004).

The optical depth for $\gamma$-rays caused by electron-positron pair
production interactions with photons of the stellar radiation background is
given by

$\tau(E_{0},z_{e})=c\int_{0}^{z_{e}}dz\,\frac{dt}{dz}\int_{0}^{2}dx\,\frac{x}{2}\int_{0}^{\infty}d\nu\,(1+z)^{3}\left[\frac{u_{\nu}(z)}{h\nu}\right]\sigma_{\gamma\gamma}[s=2E_{0}hc/\lambda
x(1+z)],$ (6)

where $u_{\nu}(z)$ is the co-moving energy density of the photon field,
(Stecker, De Jager, & Salamon 1992).

In equations (5) and (6), $E_{0}$ is the observed $\gamma$-ray energy at
redshift zero, $\lambda$ is the wavelength at redshift $z$, $z_{e}$ is the
redshift of the $\gamma$-ray source at emission, $x=(1-\cos\theta)$,
$\theta$ being the angle between the $\gamma$-ray and the soft background
photon, and the pair production cross section $\sigma_{\gamma\gamma}$ is zero
for center-of-mass energy $\sqrt{s}<2m_{e}c^{2}$, $m_{e}$ being the electron
mass. Above this threshold, the pair production cross section is given by

$\sigma_{\gamma\gamma}(s)=\frac{3}{16}\sigma_{\rm
T}(1-\beta^{2})\left[2\beta(\beta^{2}-2)+(3-\beta^{4})\ln\left(\frac{1+\beta}{1-\beta}\right)\right],$
(7)

where $\sigma_{T}$ is the Thompson scattering cross section and
$\beta=(1-4m_{e}^{2}c^{4}/s)^{1/2}$ (Jauch & Rohrlich 1955).

It follows from equation (5) that the pair-production cross section has a
threshold at $\lambda=4.75\ \mu{\rm m}\cdot E_{\gamma}({\rm TeV})$.

The optical depth of the universe to the CBR is given by

$\tau_{CBR}=5.00\times
10^{5}\sqrt{{1.11PeV}\over{E_{\gamma}}}\int_{0}^{z}{dz^{\prime}~{}(1+z^{\prime})~{}{e^{-\left({1.11PeV}\over{E_{\gamma}(1+z^{\prime})^{2}}\right)}}\over{\sqrt{\Omega_{\Lambda}+\Omega_{m}(1+z^{\prime})^{3}}}}$
(8)

(Stecker 1969; SMS06).

Figure 7 shows the 68% confidence opacity bands for interactions with IBL
photons given for sources at $z=0.1,0.5,1,3~{}$and $5$, calculated using the
methods described above, along with the opacity produced by interactions of
$\gamma$-rays with photons of the 2.7 K cosmic background radiation. Note that
at the higher energies and redshifts where the opacity is dominated by
interactions with CRB photons, the uncertainty band becomes a very thin line,
since the CBR-dominated opacity is exactly determined by equation (8).

##  5 Discussion and Conclusions

###  5.1 Comparison with our previous backward evolution model

Ten years ago we made estimates the diffuse infrared background when hardly
any mid-IR and far-IR luminosity functions had been observationally determined
at high redshifts (SMS06). That work was therefore based on the assumptions of
a “backward evolution” model. Starting from the well-determined local (z=0) LF
at 60$\mu$m, we assumed that the average locally determined transformations
between different mid-IR and far-IR wavelengths applied, unchanged, at all
redshifts. The luminosity function used was a double power-law, similar to
that of equation (2) with the parameters $a=-1.35,b=-3.6$ and $L_{*}=8.5\times
10^{23}$ W/Hz at $z=0$ as determined at 60 $\mu$m, which was the wavelength
for which the most complete galaxy LF existed. We then assumed that the effect
of redshift evolution could be taken account purely by the evolution of
$L_{*}$, as described in SMS06.

We have compared the predictions of the SMS06 backwards evolution model with
the recently observed IR LFs as used in this paper. If we make relatively
small improvements in the assumed LF parameters, the backwards evolution model
agrees with the new observations surprisingly well. The data favor a slightly
flatter LF, with a low-luminosity slope of $a=1.25$ and a high-luminosity
slope of $b=3.25$, provided that we compensate by slightly decreasing the
normalization of the LF at $L_{*}$ and $z=0$ to $3.5\times 10^{-3}$ Mpc-3
dex-1. We then match the observed LFs at higher redshifts, by assuming that
$L_{*}$ increases as $(1+z)^{3.0}$ for $0<z<2.0$, with $L_{*}$ constant at
higher redshifts. In this way, we can obtain a good fit between the backward
evolution model and the observational data.

Although this slightly modified backwards evolution reproduces the observed
LFs at all redshifts well, there are a few discrepancies, in either direction.
The most substantial disagreement is with the 8$\mu$m observations of Huang et
al. (2007) at $z=0.15$. Below the knee of the LF the observed LF the data are
up to 0.3 dex higher than the LF of our backwards evolution model. On the
other hand, our best-fitting model overpredicts the 15 $\mu$m LFs by up to
0.25 dex at luminosities above the knee as compared with the data of Pozzi et
al. (2004), Le Floch et al (2005), and Rodighiero et al. (2010). The likely
explanation of both of these discrepancies is that our previously proposed
simple model SEDs are based on average 12$\mu$m luminosities as derived from
the 60 $\mu$m observations. Therefore, they do not include the strong
contributions from the PAH bands (See section 2.2). Thus our SEDs will
overpredict the dust continuum on either side of the PAH bump feature as shown
in Figure 3. Correspondingly, if we lower the normalization to better fit the
shorter and longer wavelengths, we under-estimate the broadband fluxes at
8$\mu$m rest wavelength.

Of course, our new observationally-based calculation presented in this paper
avoids the problems of backward evolution models, since here we use directly
observed LFs at 8, 15, and 24 $\mu$m. These data include the effect of the PAH
emission features and show their significance.

###  5.2 Our present results and comparison with other work

Figure 8 shows our 68% confidence band computed for the $z=0$ IBL (EBL) along
with EBL SEDs obtained from the models of Franceschini et al. (2008) and
Domínguez et al. (2011). Figure 9 shows our 68% opacity bands for
$z=0.1,0.5,1,3~{}$and $5$, in comparison with the opacity curves of
Franceschini et al. (2008) and Domínguez et al. (2011). We note that
Franceschini et al. (2008) used data only up to 8 $\mu$m that were available
at the time and did not include a PAH component in their model. The model of
Domínguez et al. (2011) assumes a redshift evolution at redshifts greater than
$\sim$ 1 that follows the evolution in the $K$-band given by Cirasuolo et al.
(2010).

###  5.3 Conclusion

In our previous papers (SMS12 and SMS14), we presented observationally based
results for the IBL as a function of wavelength and redshift for wavelengths
below of 5 $\mu$m . Based on those results, we computed the $\gamma$-ray
opacity of the Universe caused by electron-positron pair production up to a
$\gamma$-ray energy of $1.6/(1+z)$ TeV. In this paper we have extended our
determinations of the IBL within 68% confidence bands. This determination
defines upper and lower limits on the IBL to longer wavelengths that extend to
850 $\mu$m. Our model-independent results are based on observationally derived
luminosity functions from recent galaxy survey data from both local and high
redshift surveys. These data include results from Spitzer, Herschel and
Planck. We then use these results to calculate the opacity of the Universe to
$\gamma$-rays out to PeV energies. In doing so, we also take account of the
redshift dependence of interactions of $\gamma$-rays with photons of the 2.7 K
cosmic background radiation (CBR) (Stecker 1969), since the opacity from
interactions with CBR photons dominates over that from interactions with IBL
photons at the higher $\gamma$-ray energies and redshifts.

In figure 10 we give an energy-redshift plot showing the highest energy
photons from extragalactic sources as a function of redshifts as determined
from Fermi data (Abdo et al. 2010) plotted with our 68% confidence band for
$\tau=1$.

Our direct results on the IBL are consistent with those from complimentary
$\gamma$-ray analyses using observations from the Fermi-LAT $\gamma$-ray space
telescope and the H.E.S.S. air Čerenkov telescope. Figure 11 indicates how
well our opacity results for $z=1$ overlap with those obtained by the Fermi
collaboration (Ackermann et al. 2012). Our results are also compatible with
those obtained from higher energy $\gamma$-ray observations using H.E.S.S.
(Abramowski et al. 2013). This overlap of results from two completely
different methods strengthens confidence that both techniques are indeed
complimentary and supports the concept that the spectra of cosmic $\gamma$-ray
sources can be used to probe the IBL (Stecker et al. 1992).

Thus, we find no evidence for modifications of $\gamma$-ray spectra by
processes other than absorption by pair production, either by cosmic-ray
interactions along the line of sight to the source (Essey & Kusenko 2014) or
line-of-sight photon-axion oscillations during propagation (e.g., De Angelis
et al. 2007; Mayer & Horns 2013). In this regard, we note that the Fermi
Collaboration has very recently searched for irregularities in the
$\gamma$-ray spectrum of NGC 1275 that would be caused by photon-axion
oscillations and reported negative results (Ajello et al. 2016).

We conclude that modification of the high energy $\gamma$-ray spectra of
extragalactic sources occurs dominantly by pair production interactions of
these $\gamma$-rays with photons of the IBL. They therefore support the
concept of using the future Čerenkov Telescope Array instruments to probe the
cosmic background radiation fields at infrared wavelengths. This method can be
used in conjunction with future deep galaxy survey observations using the near
infrared and mid-infrared instruments aboard the James Webb Space Telescope.

## References

  * (1)
  * (2) Abdo, A., Ackermann, M., Ajello, M. et al. 2010, ApJ 723, 1082
  * (3)
  * (4) Abramowski, A., Acero, F., Aharonian, F. et al. 2013, A&A 550, A4
  * (5)
  * (6) Ackermann, M., Ajello, M., Allafort, A. et al. 2012, Science 338, 1190
  * (7)
  * (8) Ajello, M., Albert, A., Anderson, B. et al. 2016, Phys. Rev. Letters 116, 161101
  * (9)
  * (10) Babbedge, T.S.R., Rowan-Robinson, M., Vaccari, M. et al. 2006, MNRAS 370, 1159
  * (11)
  * (12) Béthermin, M., Le Floc’h, E., Ilbert, O. et al. 2012, A&A 542, A58
  * (13)
  * (14) Berta, S., Magnelli, B., Lutz, D. et al. 2010, A&A 518, L30
  * (15)
  * (16) Bouwens, R.J., Illingworth, G. D., Oesch, P. A. 2015, ApJ 803, 34
  * (17)
  * (18) Biteau, J. & Williams, D. A. 2015, ApJ 812:60
  * (19)
  * (20) Breit, G and Wheeler, J. A. 1934, Phys. Rev. 46, 1087
  * (21)
  * (22) Caputi, K.I., Lagache, G., Yan, L. et al. 2007, ApJ 660, 97
  * (23)
  * (24) Cirasuolo, M. McLure, R. J., Dunlop, J. S. et al. 2010, MNRAS401, 1166
  * (25)
  * (26) De Angelis, A., Roncadelli, M. & Mansutti, O. 2007, Phys. Rev. D 76, 121301
  * (27)
  * (28) Domínguez, A., Primack, J. R., Rosario, D. J. et al. 2011, MNRAS 410, 2556
  * (29)
  * (30) Driver, S.P., Andrews, S.K., Davies, L.J. et al. 2016, arXiv:1605.01523.
  * (31)
  * (32) Dye, S., Dunne, L., Eales, S. et al. 2010, A&A 518, L10
  * (33)
  * (34) Dwek, E. & Krennrich, F. 2013, Astropart. Phys. 43, 112
  * (35)
  * (36) Eales, S., Raymond, G.; Roseboom, I. G. et al. 2010, A&A 518, L23
  * (37)
  * (38) Essey, W. & Kusenko, A. 2014, Astropart. Phys. 57, 30.
  * (39)
  * (40) Fazio, G. G. & Stecker, F. W. 1970, Nature 226, 135
  * (41)
  * (42) Finke, J. D., Razzaque, S. & Dermer, C. D. 2010, ApJ 712, 238
  * (43)
  * (44) Finkelstein, S. L. Ryan, R. E., Jr., Papovich, C. et al. 2015, ApJ 810:71
  * (45)
  * (46) Fixsen, D.J. Dwek, E., Mather, J. C., Bennett, C. L. & Shafer, R. A., ApJ 508, 123 (1998)
  * (47)
  * (48) Franceschini, A., Rodighiero, G. & Vaccari, M. 2008, A&A, 487, 837
  * (49)
  * (50) Fujimoto, S., Ouchi, M., Ono, Y., et al. 2016, ApJ Supp. 222, 1
  * (51)
  * (52) Gilmore, R. C. and Madau, P., Primack, J. R. et al. 2009, MNRAS 399, 1694
  * (53)
  * (54) Goto, T., Oi, N., Ohyama, Y., et al. 2015, MNRAS 452, 1684
  * (55)
  * (56) Gruppioni, C., Pozzi, F., Rodighiero, G. et al. 2013, MNRAS 432, 23
  * (57)
  * (58) Huang, J.-S., Ashby, M. L. N., Barmby, P. et al. 2007, ApJ 664, 840
  * (59)
  * (60) Helgason, K. & Kashlinsky, A. 2012, ApJ 758, L13
  * (61)
  * (62) Inoue, Y., Inoue, S., Kobayashi, M. A. R. et al. 2013, ApJ 768:197
  * (63)
  * (64) Jacobson, T., Liberati, S. Mattingly, D. & Stecker, F.W. 2004, Phys. Rev. Letters 93, 021101
  * (65)
  * (66) Jauch, J. M. & Rohrlich, F. 1955, The Theory of Photons and Electrons (Cambridge, MA: Addison-Wesley)
  * (67)
  * (68) Khaire, V. & Srianand, R. 2015, ApJ 805:33
  * (69)
  * (70) Kneiske, T. M., Mannheim, K. & Hartmann, D. H. 2002, A&A 386, 1
  * (71)
  * (72) Kneiske, T. M. & Dole, H. 2010, A&A 515, A19
  * (73)
  * (74) Lagache, G., Abergel, A., Boulanger, F. et al. 1999, A&A 344, 322
  * (75)
  * (76) Lagache, G. Dole, H., Puget, J.-L. et al. 2004, ApJ Suppl. 154:112
  * (77)
  * (78) Lapi, A. et al. González-Nuevo, J., Fan, L. et al. 2011, ApJ 742:24
  * (79)
  * (80) Le Floc’h, E., Papovich, C., Dole, H. et al. 2005, ApJ 632, 169
  * (81)
  * (82) Malkan, M. A. & Stecker, F. W. 1998, ApJ 496, 13
  * (83)
  * (84) Malkan, M. A. & Stecker, F. W. 2001, ApJ 555, 641
  * (85)
  * (86) Marchetti, L., Vaccari, M., Franceschini, A.et al. 2016, MNRAS 456, 1999
  * (87)
  * (88) Marsden, G., Chapin, E. L., Halpern, M. et al. 2011, MNRAS 417, 1192
  * (89)
  * (90) Meyer, M. & Horns, D. 2013, Phys. Rev. D 87, 035027
  * (91)
  * (92) Negrello, M., Clemens, M., Gonzalez-Nuevo, J. et al. 2013, MNRAS 429, 1309
  * (93)
  * (94) Patel, H., Clemens, D.L., Vaccari, M., et al. 2013, MNRAS 428, 291
  * (95)
  * (96) Peeters, E., Spoon, H.W.W. & Tielens, A.G.G.M. 2004, ApJ 613, 986
  * (97)
  * (98) Pérez-González, P.G., Pablo G., Rieke, G. H. et al. 2005, ApJ 630, 82
  * (99)
  * (100) Pozzi, F., Gruppioni, C., Oliver, S. et al. 2004, ApJ 609, 122
  * (101)
  * (102) Rodighiero, G., Vaccari, M., Franceschini, A. et al. 2010, A&A 518, A8
  * (103)
  * (104) Salamon, M. H. & Stecker, F. W. 1998, ApJ 493, 547
  * (105)
  * (106) Sano, K., Kawara, K., Matsuura, S. et al. 2016, ApJ 818, 72
  * (107)
  * (108) Scully, S.T., Malkan, M.A. and Stecker, F.W. 2014 (SMS14), ApJ 784:138
  * (109)
  * (110) Shipley, H.V., Heath V.; Papovich, C. et al. 2016, ApJ 818:60
  * (111)
  * (112) Smith, J.D.T., Draine, B. T., Dale, D. A. et al. 2007, ApJ 656, 770
  * (113)
  * (114) Smith, D.J.B., Dunne, L., da Cunha, E. et al. 2012, MNRAS 427, 703
  * (115)
  * (116) Somerville, R. S., Gilmore, R. C., Primack, J. R. & Domínguez, A. 2012, MNRAS 423, 1992
  * (117)
  * (118) Spoon, H.W.W., Marshall, J. A., Houck, J. R. et al. 2007, ApJL 654, L49
  * (119)
  * (120) Spiesman, W.J., Hauser, M. G.; Kelsall, T. et al. 1995, ApJ 442, 662.
  * (121)
  * (122) Stecker, F.W. 1969, ApJ 157, 507
  * (123)
  * (124) Stecker, F. W., De Jager, O.C. & Salamon, M. H. 1992, ApJ 390, L49
  * (125)
  * (126) Stecker, F.W. & Glashow, S.L. 2001, Astropart. Phys. 16, 97
  * (127)
  * (128) Stecker, F.W., Malkan, M. A. & Scully, S. T. 2006 (SMS06), ApJ 648, 774
  * (129)
  * (130) Stecker, F.W., Malkan, M. A. & Scully, S. T. 2012 (SMS12), ApJ 761, 128
  * (131)
  * (132) Vaccari, M., Marchetti, L., Franceschini, A. et al. 2010, A&A 518, L20
  * (133)
  * (134) Viero, M.P., Moncelsi, L., Quadri, R. F. et al. 2015, ApJ 809, L22
  * (135)

Figure 1: The luminosity densities for 8, 12, 15, and 24 $\mu$m wavebands.
Data are from several sources. Some Spitzer data from all 4 wavebands are from
Rodighiero et al. (2010). AKARI 8 and 12 $\mu$m data are from Goto et al.
(2015), Spitzer data at 8 $\mu$m are from Huang et al. (2007) and Caputi et
al. (2007). Spitzer data at 8 and 24 $\mu$m are from Babbedge et al. (2006).
Spitzer data at 15 $\mu$m are from Le Floc’h et al. (2005). Spitzer data at 24
$\mu$m are from Pérez-González et al. (2005). The grey shading represents the
68% confidence bands (see text). Figure 2: The luminosity densities for the
35, 60, 90, and 250 $\mu$m wavebands.The grey shading represents the 68%
confidence bands. Herschel data at 35 60, and 90 $\mu$m are from Gruppioni et
al. (2013). Also 100 $\mu$m data from Lapi et al (2011) are plotted on the 90
$\mu$m graph. Herschel data at 250 $\mu$m are from Eales et al. (2010), Dye et
al. (2010) and Smith et al. (2012) (see text). Figure 3: Average low redshift
galaxy SED at MIR wavelengths based on the Class 1C SED of Spoon et al. (2007)
normalized to our low redshift LDs from Figures 1 and 2, as indicated by the
shaded region. Figure 4: Our spectral energy distribution of the EBL together
with empirical data based on our mid-IR LDs and far-IR LDs and the results of
SMS12 and SMS14. The light shaded area between $\sim$10 $\mu$m and $\sim$30
$\mu$m indicates the maximum effect of the PAH bands (see Sect. 2.2). The
lower limits from galaxy counts are shown in blue; direct measurements and
extrapolations from galaxy counts in the mid-IR are shown in black. References
for the empirical data before 2012 are given by Lagache et al. (2005) and Dwek
& Krennrich (2013). A 3.5 $\mu$m point is from Sano et al. (2016). The red
shaded area is based on the COBE-FIRAS results of Fixsen et al. (1998) with
limits described by modified black body spectra. Figure 5: Upper limit
envelope on the co-moving energy density as a function of energy and redshift.
Figure 6: Lower limit envelope on the co-moving energy density as a function
of energy and redshift. Figure 7: The optical depth of the universe from the
IBL and the CBR as well as the total optical depth as a function of energy,
given for redshifts of 0.1, 0.5, 1, 3, 5. It can be seen that the contribution
to the optical depth from the IBL dominates at lower $\gamma$-ray energies and
redshifts and that from the CBR photons dominates at the higher energies and
redshifts. The optical depth from CBR photons is an exact function of energy
as given by equation (8) and therefore the confidence band is becomes a thin
line. The dashed lines indicate the opacities $\tau=1$ and $\tau=3$.  Figure
8: A comparison of our confidence band with the models of Franceschini et al.
(2008) (solid black line) and Domínguez et al. (2011) (red dashed line).
Figure 9: Comparison of our opacity results with those obtained by the models
of Franceschini et al. (2008) (solid black line) and Domínguez et al. (2011)
(red dashed line). Figure 10: A $\tau=1$ energy-redshift plot (Fazio & Stecker
1970) showing our uncertainty band results compared with the Fermi plot of
their highest energy photons from FSRQs (red), BL Lacs (black) and and GRBs
(blue) vs. redshift (from Abdo et al. 2010). Figure 11: Comparison of our
results for $z=1$ with those obtained from an analysis of blazar $\gamma$-ray
spectra (Ackermann et al. 2012)