The probabilities of causation are commonly used to solve decision-making problems. Tian and Pearl derived sharp bounds for the probability of necessity and sufficiency (PNS), the probability of sufficiency (PS), and the probability of necessity (PN) using experimental and observational data. The assumption is that one is in possession of a large enough sample to permit an accurate estimation of the experimental and observational distributions. In this study, we present a method for determining the sample size needed for such estimation, when a given confidence interval (CI) is specified. We further show by simulation that the proposed sample size delivered stable estimations of the bounds of PNS.

§ INTRODUCTION

The probabilities of causation are widely used in many areas of industry, marketing, and health science, to solve decision-making problems. For example, Li and Pearl <cit.> proposed the “benefit function”, a linear combination of the probabilities of causation, which is the payoff/cost associated with selecting an individual with given characteristics to identify a set of individuals who are most likely to exhibit a desired mode of behavior. Mueller and Pearl <cit.> demonstrated, for example, that the probabilities of causation should be considered in personalized decision-making. Li et al. <cit.> showed that the probabilities of causation can improve the accuracy of machine learning algorithms.

Pearl <cit.> first used the structural causal model (SCM) to defined three binary probabilities of causation (i.e., PNS, PN, and PS) <cit.>. Tian and Pearl <cit.> then used experimental and observational data to bound those three probabilities of causation. Li and Pearl <cit.> established formal proof of those bounds. Mueller, Li, and Pearl <cit.> proposed narrowing the bounds of PNS using covariate information and the causal structure. Dawid et al. <cit.> also proposed using covariate information to narrow the bounds of PN. Li and Pearl <cit.> recently
established the theoretical bounds of nonbinary probabilities of causation.

All the abovementioned works are asymptotic (i.e., assuming a adequately large sample size to estimate the experimental and observational distributions). The proposed results in those works are relationships between the experimental and observational distributions and the probabilities of causation. However, the adequate sample size for obtaining those probabilities of causation remains unclear, thereby creating a barrier between the theoretical results and the real-world applications. Consider the following motivating example: a mobile carrier that wants to identify customers who are likely to discontinue their services within the next quarter based on customer characteristics (company management has access to user data, such as income, age, usage, and monthly payments). The carrier will then offer these customers a special renewal deal to dissuade them from discontinuing their services and to increase their service renewal rate. These offers provide considerable discounts to the customers, and the management prefers that these offers be made only to those customers who would continue to use the service if and only if they receive the offer. The manager decides to use Li and Pearl's unit selection model <cit.> but is unsure how many experimental and observational samples are required. Are $1000$ experimental and $1000$ observational samples adequate to bound the benefit function such that the error of the bounds are within $0.1$?

We present an assessment of the “adequate" of the sample size in the sense of CI in this study. We would then be able to answer the question, “How many samples are adequate to estimate the probability of causation?" as “This amount of samples is adequate to obtain the bounds of the probability of causation in $95\%$ CIs with margin of errors of $0.05$." The probabilities of causation in most cases are not identifiable; therefore, the CIs are for the bounds of the probabilities of causation in such cases.

§ PRELIMINARIES

We review the definitions for the three aspects of binary causation in this section, as defined in <cit.>. We use the language of counterfactuals in SCM, as defined in <cit.>. We use $Y_x=y$ to denote the counterfactual sentence “Variable $Y$ would have the value $y$, had $X$ been $x$". For the rest of the paper, we use $y_x$ to denote the event $Y_x=y$, $y_{x'}$ to denote the event $Y_{x'}=y$, $y'_x$ to denote the event $Y_x=y'$, and $y'_{x'}$ to denote the event $Y_{x'}=y'$. We assume that experimental distribution will be summarized in the form of the causal effects such as $P(y_x)$ and observational distribution will be summarized in the form of the joint probability function such as $P(x,y)$. If neither variable is specified, variable $X$ represents treatment and variable $Y$ represents effect.

The following are three prominent probabilities of causation:

Let $X$ and $Y$ be two binary variables in a causal model $M$, let $x$ and $y$ stand for the propositions $X=true$ and $Y=true$, respectively, and $x'$ and $y'$ for their complements. The probability of necessity is defined as the expression <cit.>

\begin{eqnarray}
\text{PN} \delequal P(Y_{x'}=false|X=true,Y=true)
\delequal  P(y'_{x'}|x,y) \nonumber
\label{pn}
\end{eqnarray}

\begin{eqnarray}
\text{PS} \delequal P(y_x|y',x') \nonumber
\label{ps}
\end{eqnarray}

\begin{eqnarray}
\text{PNS}\delequal P(y_x,y'_{x'}) \nonumber
\label{pns}
\end{eqnarray}

PNS denotes for the probability that $y$ would respond to $x$ both ways, and therefore measures both the sufficiency and necessity of $x$ to produce $y$.

Tian and Pearl <cit.> used Balke's program <cit.> to provide tight bounds for PNS, PN, and PS.

PNS has the following tight bounds:

\begin{eqnarray}
\max \left \{
\begin{array}{cc}
0, \\
P(y_x) - P(y_{x'}), \\
P(y) - P(y_{x'}), \\
P(y_x) - P(y)\\
\end{array}
\right \}
\le \text{PNS}
\le \min \left \{
\begin{array}{cc}
P(y_x), \\
P(y'_{x'}), \\
P(x,y) + P(x',y'), \\
P(y_x) - P(y_{x'}) + P(x, y') + P(x', y)
\end{array}
\right \}
\label{pnsub}
\end{eqnarray}

Note that we omitted the bounds of PN and PS because this study focuses primarily on the adequate sample size for estimating the bounds of PNS, it is simple to extend to other probabilities of causation.

§ MAIN RESULT

The bounds of PNS are the linear combination of the experimental distributions $P(y_x),P(y_{x'})$ and the observational distributions $P(x,y),P(x,y'),P(x',y),P(x',y')$ from Equation <ref>. Therefore, if we can obtain the CI of each of these distributions, then we can obtain the CIs of the bounds of the PNS. Let $R$ be a random variable such that $R=1$ if the event $y_x$ occurs and $R=0$ if the event ${y'}_x$ occurs, then it is clear that $R\sim Bernoulli(P(y_x))$. Therefore, if we use the frequentist to estimate the experimental and observational distributions, we have the following theorem and corollary (the detailed proof are in the appendix):

Given $m$ experimental samples and $n$ observational samples, if the frequentist is used to estimate the experimental and observational distributions, then the margin error of the bounds of PNS in $1-\alpha$ confidence interval is at most $z_{1-\alpha/2}(\sqrt{\frac{1}{m}} + \sqrt{\frac{1}{n}})$, where $z_{1-\alpha/2}$ can be found on z-table of standard normal distribution.

If the frequentist is used to estimate the experimental and observational distribution, to obtain the at most $0.05$ margin error of the bounds of PNS in $95\%$ confidence interval, we need $m$ experimental samples and $n$ observational samples, where $(\sqrt{\frac{1}{m}} + \sqrt{\frac{1}{n}}) \le 5/196.$ More specifically, if $m=n$, $6147$ experimental and $6147$ observational samples are adequate to obtain the at most $0.05$ margin error of the bounds of PNS in $95\%$ confidence interval.

This amount of samples ensures that the margin errors of the bounds of PNS in $95\%$ CI are no more than $0.05$. However, in practice, we usually do not need this amount of samples because there is only one term (i.e., $P(y_x)-P(y_{x'})+P(x,y')+P(x',y)$) in the PNS bounds of Equation <ref>, which consists of four distributions. The terms such as $P(y_x)$ only require $385$ experimental and $385$ observational samples to obtain the at most $0.05$ margin error of the $95\%$ CI, and the terms such as $P(y_x)-P(y_{x'})$ only require $1537$ experimental and $1537$ observational samples to obtain the at most $0.05$ margin error of the $95\%$ CI. We will illustrate the real errors of the estimations in simulated studies in the following section.

§ SIMULATION RESULTS

Here, we present simulated studies to show that the proposed number of experimental and observational samples are adequate to obtain the desired margin errors of the bounds of PNS using two SCMs.

§.§ Causal Model

To estimate the margin errors of the bounds, we must first understand the data generation process to have true experimental and observational distributions. The two models we are using are shown in Figure <ref> (two models have the same causal graph, but with different coefficients in SCMs; the generation method of the models is in the appendix), where $X$ is a binary treatment, $Y$ is a binary effect, and $Z$ is a set of $20$ independent confounders (say $Z_1,...,Z_{20}$). The structural equations are as follow (for simplicity reason, we let $x=1,x'=0$, and $y=1,y'=0$):

[->,>=stealth',node distance=2cm,
thick,main node/.style=circle,fill,inner sep=1.5pt]
[main node] (1) [label=above:$Z$];
[main node] (3) [below left =1cm of 1,label=left:$X$];
[main node] (4) [below right =1cm of 1,label=right:$Y$] ;
[every node/.style=font=]
(1) edge node  (3)
(1) edge node  (4)
(3) edge node  (4);

The Causal Model, where $X$ is a binary treatment, $Y$ is a binary effect, and $Z$ is a set of $20$ independent binary confounders.

\begin{eqnarray*}
Z_i &=& U_{Z_i} \text{ for } i \in \{1,...,20\},\\
%X &=& 1 \text{ if } M_X +U_X > 3,\\
X&=&f_X(M_X,U_X) =
\left \{
\begin{array}{cc}
1& \text{ if } M_X +U_X > 0.5, \\
0& \text{otherwize},
\end{array}
\right \}\\
%Y &=& 0.867991236792X+M_Y+U_Y,\\
Y&=&f_Y(X,M_Y,U_Y) =
\left \{
\begin{array}{cc}
1& \text{ if } 0< CX+M_Y+U_Y < 1 \text{ or } 1 < CX+M_Y+U_Y < 2,\\
% 1& \text{ if } 1 < CX+M_Y+U_Y < 2, \\
0& \text{otherwize},
\end{array}
\right \}\\
&&\text{where, } U_{Z_i}, U_X, U_Y \text{ are binary exogenous variables with Bernoulli distributions,}\\ &&\text{C is a constant, and }M_X,M_Y\text{ are linear combinatations of }Z_{1},...,Z_{20}.
% s.t., \\
% &&U_{Z_1} \sim \text{Bernoulli}(0.062053709067), U_{Z_2} \sim \text{Bernoulli}(0.060783320775),\\
% &&U_{Z_3} \sim \text{Bernoulli}(0.440283974413), U_{Z_4} \sim \text{Bernoulli}(0.062515746008),\\
% &&U_{Z_5} \sim \text{Bernoulli}(0.113914295859), U_{Z_6} \sim \text{Bernoulli}(0.153534224919),\\
% &&U_{Z_7} \sim \text{Bernoulli}(0.933699217108), U_{Z_8} \sim \text{Bernoulli}(0.602199323486),\\
% &&U_{Z_9} \sim \text{Bernoulli}(0.248045364151), U_{Z_{10}} \sim \text{Bernoulli}(0.577622106574),\\
% &&U_{Z_{11}} \sim \text{Bernoulli}(0.894331031858), U_{Z_{12}} \sim \text{Bernoulli}(0.612468910385),\\
% &&U_{Z_{13}} \sim \text{Bernoulli}(0.513732540844), U_{Z_{14}} \sim \text{Bernoulli}(0.206426380498),\\
% &&U_{Z_{15}} \sim \text{Bernoulli}(0.404311039089), U_{Z_{16}} \sim \text{Bernoulli}(0.195768108046),\\
% &&U_{Z_{17}} \sim \text{Bernoulli}(0.276887054651),\\
% &&U_{X} \sim \text{Bernoulli}(0.635925619967), U_{Y} \sim \text{Bernoulli}(0.112832229360),
% \end{eqnarray*}
% \begin{eqnarray*}
% &M_X& =
% \begin{bmatrix}
% Z_1~Z_2~...~Z_{17}
% \end{bmatrix}\times
% \begin{bmatrix}
% 0.135855821291\\0.748077631146\\0.566603758821\\0.829393512973\\	0.332734093899\\0.69349211629\\0.134668072184\\0.123198197011\\	0.38275550532\\0.54989366455\\0.829905428674\\0.917548456559\\	0.352152469922\\0.656879580936\\0.86086674095\\0.95100545228\\	0.352773184618
% %   x_{1}  \\
% %   x_{2}  \\
% %   \vdots \\
% %   x_{m}
% \end{bmatrix},
% M_Y =
% \begin{bmatrix}
% Z_1~Z_2~...~Z_{17}
% \end{bmatrix}\times
% \begin{bmatrix}
% 0.864339662391\\0.89355070787\\0.954829789303\\0.874456470217\\	0.936083805032\\0.562535633212\\0.421184877554\\0.129582468009\\	0.949521408339\\0.450139315368\\0.408096150159\\0.866196854316\\	0.656067001249\\0.770489996977\\0.42411338788\\0.697807293946\\	0.748871824623
% %   x_{1}  \\
% %   x_{2}  \\
% %   \vdots \\
% %   x_{m}
% \end{bmatrix}
\end{eqnarray*}

The value of $C,M_X,M_Y$ and the distributions of $U_X,U_Y,U_{Z_1},...,U_{Z_{20}}$ for the two models are provided in the appendix.

§.§ Informer Data

Based on the model from the previous section, individuals are determined by $22$ binary exogenous variables. From the viewpoint of the informer, we must know the actual experimental (i.e., $P(y_x),P(y_{x'})$) and observational distributions (i.e., $P(x,y),P(x,y'),P(x',y)$) to compute the true PNS bounds for the comparison purpose. Since the structural equations and the exogenous variable distributions are explicit, the experimental and observational distributions are as follows (we illustrated two terms, see the appendix for full details):

\begin{eqnarray*}
P(y_x) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})f_Y(1,M_Y,u_Y),\\
% P(y_{x'}) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})f_Y(0,M_Y,u_Y),\\
P(x,y) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times \\
&&\times f_X(M_X,u_X)f_Y(f_X(M_X,u_X),M_Y,u_Y),\\
% P(x,y') &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times\\
% && \times f_X(M_X,u_X)(1 - f_Y(f_X(M_X,u_X),M_Y,u_Y)),\\
% P(x',y) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times\\
% &&\times (1-f_X(M_X,u_X))f_Y(f_X(M_X,u_X),M_Y,u_Y).
% P(x',y') &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times \\
% &&\times (1-f_X(M_X,u_X))(1 - f_Y(f_X(M_X,u_X),M_Y,u_Y)).
\end{eqnarray*}

We can then obtain the informer view of the bounds of PNS using Equation <ref> and the above observational and experimental distributions.

§.§ Experimental Sample

Here is how we prepared the finite experimental samples. Again, an individual is determined by $U_X$, $U_Y$, and $U_{Z_1},...,U_{Z_{20}}$. Therefore, we first generated $(U_X,U_Y,U_{Z_1},...,U_{Z_{20}})$ at random using their distributions; then we generated $X$ at random using $Bernoulli(0.5)$; the value of $Y$ is then $f_Y(X,M_Y,U_Y)$. We then collect a experimental sample $(X,Y)$. $P(y_x)$ is then estimated as $\hat{P}(y_x)=\frac{\text{the number of experimental samples (1,1)}}{\text{the number of experimental samples (1,Y)}}$. $P(y_{x'})$ is then estimated as $\hat{P}(y_{x'})=\frac{\text{the number of experimental samples (0,1)}}{\text{the number of experimental samples (0,Y)}}$.

§.§ Observational Sample

Similarly to experimental sample, we first generated $(U_X,U_Y,U_{Z_1},...,U_{Z_{20}})$ at random using their distributions; the value of $X$ is then $f_X(M_X,U_X)$; the value of $Y$ is then $f_Y(X,M_Y,U_Y)$. Then, we collect a observational sample $(X,Y)$. $P(x,y),P(x,y'),P(x',y)$ estimates are then $\hat{P}(x,y)=\frac{\text{the number of observational samples (1,1)}}{\text{the number of total observational samples}}$, $\hat{P}(x,y')=\frac{\text{the number of observational samples (1,0)}}{\text{the number of total observational samples}}$, and $\hat{P}(x',y)=\frac{\text{the number of observational samples (0,1)}}{\text{the number of total observational samples}}$.

§.§ Simulation Results

For the two SCMs, we randomly generated $385, 1537,$ and $6147$ experimental and $385, 1537,$ and $6147$ observational samples respectively for $1000$ times as described in the previous sections. The true bounds of the PNS obtained from true distributions were then compared to  the estimated bounds obtained from finite experimental and observational samples. Figures <ref> and <ref> show the comparison results. We see that when the sample size is $385$ or even $1537$, the bounds are still in a relatively large unstable situation, with some of the estimated lower bounds being even larger than the true PNS value. When the sample size reached $6147$, the estimations appear very stable and always include the true PNS value. We also plot the sample size v.s. the average error of estimations (i.e., $|$estimated bounds $-$ true bounds$|/1000$) as shown in Figure <ref>. The average errors improved significantly before $1500$ sample sizes, therefore, $1537$ should be a good sample size to use in practice even though the adequate sample size from the theorem is $6147$. Besides, starting at roughly $300$ sample size, the average errors are less than the expected margin error of $0.05$. This is because the proposed theorem is of adequate size and considers the worst-case scenario.

$385$ experimental and $385$ observational samples.

$1537$ experimental and $1537$ observational samples.

$6147$ experimental and $6147$ observational samples.

Estimation of the bounds of PNS for the first model using different size of samples.

$385$ experimental and $385$ observational samples.

$1537$ experimental and $1537$ observational samples.

$6147$ experimental and $6147$ observational samples.

Estimation of the bounds of PNS for the second model using different size of samples.

§ DISCUSSION

We demonstrated how to evaluate the precision of the probability of causation estimation using PNS. It is simple to extend to nonbinary probabilities of causation as defined in <cit.>. In such cases, we can still define a random variable, $R$, that obeys the Bernoulli distribution and whose parameter is the experimental distribution $P(y_x)$. This work can be extended to unit selection problems <cit.>, because the unit selection problems are linear combination of the probabilities of causation. This work can be extended to estimate experimental distributions <cit.> using observational data.

§ CONCLUSION

We showed how to assess the accuracy of a PNS bounds estimation given by a finite number of experimental and observational samples, as well as how many experimental and observational samples are adequate to achieve a certain accuracy level of a PNS bounds estimation. Both theoretical and experimental results are provided. The proposed method can also be extended to general probabilities of causation.

First model.

Second model.

Average error of estimations using different size of data.

§ ACKNOWLEDGEMENTS

This research was supported in parts by grants from the National Science
Foundation [#IIS-2106908], Office of Naval Research [#N00014-17-S-12091
and #N00014-21-1-2351], and Toyota Research Institute of North America

§ APPENDIX

§.§ Proof of Theorem <ref>

Let $R$ be a random variable, such that $R=1$ if the event $y_x$ occurs and $R=0$ if the event $y'_x$ occrus, then $R\sim Bernoulli(P(y_x))$.

Let $R=(R_1,...,R_m)$ be $m$ random experimental samples of $R$, and $S_m=\sum_{i=1}^{m}R_i$.

Therefore, $S_m\sim binomial(m,P(y_x))$.

By Central limit theorem, we have,

$\frac{S_m-mp}{\sqrt{mP(y_x)(1-P(y_x))}}\rightarrow{}N(0,1)$ as $m \rightarrow{}\infty$.

We also have $\hat{P}(y_x)=S_m/m\xrightarrow{P} P(y_x)$,

Thus, we have,

$\sqrt{\frac{P(y_x)(1-P(y_x))}{\hat{P}(y_x)(1-\hat{P}(y_x))}}\xrightarrow{P}1$ and $\frac{\hat{P}(y_x)-P(y_x)}{\sqrt{P(y_x)(1-P(y_x))/m}}\xrightarrow{L}N(0,1)$.

Therefore, we have,

We know that $P(|\frac{\hat{P}(y_x)-P(y_x)}{\sqrt{\hat{P}(y_x)(1-\hat{P}(y_x))/m}}|\le z_{1-\alpha/2})=1-\alpha$,

we have the confidence interval of $\hat{P}(y_x)$ is

$[\hat{P}(y_x)-z_{1-\alpha/2}\cdot \sqrt{\hat{P}(y_x)(1-\hat{P}(y_x))/m},\hat{P}(y_x)+z_{1-\alpha/2}\cdot \sqrt{\hat{P}(y_x)(1-\hat{P}(y_x))/m}]$.

The margin of error of $\hat{P}(y_x)$ is $z_{1-\alpha/2}\cdot \sqrt{\hat{P}(y_x)(1-\hat{P}(y_x))/m}$, denoted as

$err(\hat{P}(y_x))=z_{1-\alpha/2}\cdot \sqrt{\hat{P}(y_x)(1-\hat{P}(y_x))/m}$.

We then have,

$err(\hat{P}(y_x))=z_{1-\alpha/2}\cdot \sqrt{\hat{P}(y_x)(1-\hat{P}(y_x))/m}\le z_{1-\alpha/2}\cdot \sqrt{\frac{1}{4m}}=\frac{z_{1-\alpha/2}}{2}\sqrt{\frac{1}{m}}$.

Similarly, we have,

We also have the fact that $err(A+B)\le err(A)+err(B)$ and $err(A-B)\le err(A)+err(B)$,

plug into Equation <ref>, we obtain that $err(\text{bounds of PNS})\le z_{1-\alpha/2}(\sqrt{\frac{1}{m}}+\sqrt{\frac{1}{n}})$.

§.§ Proof of Corollary <ref>

Let $\alpha = 0.05$.

We need $z_{0.975}(\sqrt{\frac{1}{m}}+\sqrt{\frac{1}{n}})\le 0.05$.

Simply plug $z_{0.975}=1.96$ in, we have,

$(\sqrt{\frac{1}{m}}+\sqrt{\frac{1}{n}}) \le 5/196$.

And if $m=n$, we have,

$\sqrt{\frac{1}{m}}\le 5/392$,

Therefore, $m\ge 6146.56$.

§.§ Two Causal Models

First, $M_X$ and $M_Y$ are both linear combinatation of $Z_{1},...,Z_{20}$. So for each model, we need to generate $20$ coefficients for $M_X$, $20$ coefficients for $M_Y$, a constant $C$, and $22$ Bernoulli distribution parameters for $U_X,U_Y,U_{Z_1},...,U_{Z_{20}}$. We generated the coeefficients for $M_X,M_Y$ and $C$ uniformly from $[-1,1]$, and generated the Bernoulli distribution parameters unoformly from $[0,1]$. The detailed two models are as follow:

§.§.§ Model 1

\begin{eqnarray*}
Z_i &=& U_{Z_i} \text{ for } i \in \{1,...,20\},\\
%X &=& 1 \text{ if } M_X +U_X > 3,\\
X&=&f_X(M_X,U_X) =
\left \{
\begin{array}{cc}
1& \text{ if } M_X +U_X > 0.5, \\
0& \text{otherwize},
\end{array}
\right \}\\
%Y &=& 0.867991236792X+M_Y+U_Y,\\
Y&=&f_Y(X,M_Y,U_Y) =
\left \{
\begin{array}{cc}
1& \text{ if } 0< CX+M_Y+U_Y < 1 \text{ or } 1 < CX+M_Y+U_Y < 2,\\
% 1& \text{ if } 1 < CX+M_Y+U_Y < 2, \\
0& \text{otherwize},
\end{array}
\right \}\\
&&\text{where, } U_{Z_i}, U_X, U_Y \text{ are binary exogenous variables with Bernoulli distributions.}
s.t., \\
&&U_{Z_1} \sim \text{Bernoulli}(0.352913861526), U_{Z_2} \sim \text{Bernoulli}(0.460995855543),\\
&&U_{Z_3} \sim \text{Bernoulli}(0.331702473392), U_{Z_4} \sim \text{Bernoulli}(0.885505026779),\\
&&U_{Z_5} \sim \text{Bernoulli}(0.017026872706), U_{Z_6} \sim \text{Bernoulli}(0.380772701708),\\
&&U_{Z_7} \sim \text{Bernoulli}(0.028092602705), U_{Z_8} \sim \text{Bernoulli}(0.220819399962),\\
&&U_{Z_9} \sim \text{Bernoulli}(0.617742227477), U_{Z_{10}} \sim \text{Bernoulli}(0.981975046713),\\
&&U_{Z_{11}} \sim \text{Bernoulli}(0.142042291381), U_{Z_{12}} \sim \text{Bernoulli}(0.833602592350),\\
&&U_{Z_{13}} \sim \text{Bernoulli}(0.882938907115), U_{Z_{14}} \sim \text{Bernoulli}(0.542143191999),\\
&&U_{Z_{15}} \sim \text{Bernoulli}(0.085023436884), U_{Z_{16}} \sim \text{Bernoulli}(0.645357252864),\\
&&U_{Z_{17}} \sim \text{Bernoulli}(0.863787135134), U_{Z_{18}} \sim \text{Bernoulli}(0.460539711624),\\
&&U_{Z_{19}} \sim \text{Bernoulli}(0.314014079207), U_{Z_{20}} \sim \text{Bernoulli}(0.685879388218),\\
&&U_{X} \sim \text{Bernoulli}(0.601680857267), U_{Y} \sim \text{Bernoulli}(0.497668975278),\\
\end{eqnarray*}

\begin{eqnarray*}
&M_X& =
\begin{bmatrix}
\end{bmatrix}\times
\begin{bmatrix}
0.259223510143\\ -0.658140989167\\ -0.75025831768\\ 0.162906462426\\ 0.652023463285\\ -0.0892939586541\\ 0.421469107769\\ -0.443129684766\\ 0.802624388789\\ -0.225740978499\\ 0.716621631717\\ 0.0650682260309\\ -0.220690334026\\ 0.156355773665\\ -0.50693672491\\ -0.707060278115\\ 0.418812816935\\ -0.0822118703986\\ 0.769299853833\\ -0.511585391002
\end{bmatrix},
M_Y =
\begin{bmatrix}
\end{bmatrix}\times
\begin{bmatrix}
-0.792867111918\\ 0.759967136147\\ 0.55437722369\\ 0.503970540409\\ -0.527187144651\\ 0.378619988091\\ 0.269255196301\\ 0.671597043594\\ 0.396010142274\\ 0.325228576643\\ 0.657808327574\\ 0.801655023993\\ 0.0907679484097\\ -0.0713852594543\\ -0.0691046005285\\ -0.222582013343\\ -0.848408031595\\ -0.584285069026\\ -0.324874831799\\ 0.625621583197
\end{bmatrix}
\end{eqnarray*}

§.§.§ Model 2

\begin{eqnarray*}
Z_i &=& U_{Z_i} \text{ for } i \in \{1,...,20\},\\
%X &=& 1 \text{ if } M_X +U_X > 3,\\
X&=&f_X(M_X,U_X) =
\left \{
\begin{array}{cc}
1& \text{ if } M_X +U_X > 0.5, \\
0& \text{otherwize},
\end{array}
\right \}\\
%Y &=& 0.867991236792X+M_Y+U_Y,\\
Y&=&f_Y(X,M_Y,U_Y) =
\left \{
\begin{array}{cc}
1& \text{ if } 0< CX+M_Y+U_Y < 1 \text{ or } 1 < CX+M_Y+U_Y < 2,\\
% 1& \text{ if } 1 < CX+M_Y+U_Y < 2, \\
0& \text{otherwize},
\end{array}
\right \}\\
&&\text{where, } U_{Z_i}, U_X, U_Y \text{ are binary exogenous variables with Bernoulli distributions.}
s.t., \\
&&U_{Z_1} \sim \text{Bernoulli}(0.524110233482), U_{Z_2} \sim \text{Bernoulli}(0.689566064108),\\
&&U_{Z_3} \sim \text{Bernoulli}(0.180145428970), U_{Z_4} \sim \text{Bernoulli}(0.317153536644),\\
&&U_{Z_5} \sim \text{Bernoulli}(0.046268153873), U_{Z_6} \sim \text{Bernoulli}(0.340145244411),\\
&&U_{Z_7} \sim \text{Bernoulli}(0.100912238566), U_{Z_8} \sim \text{Bernoulli}(0.772038172066),\\
&&U_{Z_9} \sim \text{Bernoulli}(0.913108434869), U_{Z_{10}} \sim \text{Bernoulli}(0.364272299067),\\
&&U_{Z_{11}} \sim \text{Bernoulli}(0.063667554704), U_{Z_{12}} \sim \text{Bernoulli}(0.454839320009),\\
&&U_{Z_{13}} \sim \text{Bernoulli}(0.586687215140), U_{Z_{14}} \sim \text{Bernoulli}(0.018824647595),\\
&&U_{Z_{15}} \sim \text{Bernoulli}(0.871017316787), U_{Z_{16}} \sim \text{Bernoulli}(0.164966968157),\\
&&U_{Z_{17}} \sim \text{Bernoulli}(0.578925020078), U_{Z_{18}} \sim \text{Bernoulli}(0.983082980658),\\
&&U_{Z_{19}} \sim \text{Bernoulli}(0.018033993991), U_{Z_{20}} \sim \text{Bernoulli}(0.074629121266),\\
&&U_{X} \sim \text{Bernoulli}(0.29908139311), U_{Y} \sim \text{Bernoulli}(0.9226108109253),\\
&&C= 0.975140894243,
\end{eqnarray*}

\begin{eqnarray*}
&M_X& =
\begin{bmatrix}
\end{bmatrix}\times
\begin{bmatrix}
0.843870221861\\ 0.178759296447\\ -0.372349746729\\ -0.950904544846\\ -0.439457721339\\ -0.725970103834\\ -0.791203963585\\ -0.843183562918\\ -0.68422616618\\ -0.782051030131\\ -0.434420454146\\ -0.445019418094\\ 0.751698021555\\ -0.185984172192\\ 0.191948271392\\ 0.401334543567\\ 0.331387702568\\ 0.522595634402\\ -0.928734581669\\ 0.203436441511
\end{bmatrix},
M_Y =
\begin{bmatrix}
\end{bmatrix}\times
\begin{bmatrix}
-0.453251661832\\ 0.424563325534\\ 0.0924810605305\\ 0.312680246141\\ 0.7676961338\\ 0.124337421843\\ -0.435341306455\\ 0.248957751703\\ -0.161303883519\\ -0.537653062121\\ -0.222087991408\\ 0.190167775134\\ -0.788147770713\\ -0.593030174012\\ -0.308066297974\\ 0.218776507777\\ -0.751253645088\\ -0.11151455376\\ 0.785227235182\\ -0.568046522383
\end{bmatrix}
\end{eqnarray*}

§.§ Informer Data

Detailed informer data can be obtained via the following equations:

\begin{eqnarray*}
P(y_x) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})f_Y(1,M_Y,u_Y),\\
P(y_{x'}) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})f_Y(0,M_Y,u_Y),\\
P(x,y) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times \\
&&\times f_X(M_X,u_X)f_Y(f_X(M_X,u_X),M_Y,u_Y),\\
P(x,y') &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times\\
&& \times f_X(M_X,u_X)(1 - f_Y(f_X(M_X,u_X),M_Y,u_Y)),\\
P(x',y) &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times\\
&&\times (1-f_X(M_X,u_X))f_Y(f_X(M_X,u_X),M_Y,u_Y).\\
P(x',y') &=& \sum_{U_X,U_Y,U_{Z_1},...,U_{Z_{20}}}P(u_X)P(u_Y)P(u_{Z_1})\times ...\times P(u_{Z_{20}})\times \\
&&\times (1-f_X(M_X,u_X))(1 - f_Y(f_X(M_X,u_X),M_Y,u_Y)).
\end{eqnarray*}


# Probing $f(R)$ gravity with PLANCK data on cluster pressure profiles

I. De Martino1,4 M. De Laurentis2,3,4 F. Atrio-Barandela1 S. Capozziello3,4 1
Física Teórica, Universidad de Salamanca, 37008 Salamanca, Spain;
2Department of Theoretical Physics, Tomsk State Pedagogical University (TSPU),
pr. Komsomolsky, 75, Tomsk, 634041, Russia;
3 Dipartimento di Fisica, Università di Napoli ”Federico II”;
4INFN sez. di Napoli Compl. Univ. di Monte S. Angelo, Edificio G, Via Cinthia,
I-80126 - Napoli, Italy<EMAIL_ADDRESS>
###### Abstract

Analytical $f(R)$-gravity models introduce Yukawa-like corrections to the
Newtonian potential in the weak field limit. These models can explain the
dynamics of galaxies and cluster of galaxies without requiring dark matter. To
test the model, we have computed the pressure profile of 579 X-ray galaxy
clusters assuming the gas is in hydrostatic equilibrium within the potential
well of the modified gravitational potential. We have compared those profiles
with the ones measured in the foreground cleaned SMICA released by the Planck
Collaboration. Our results show that Extended Theories of Gravity explain the
dynamics of self-gravitating systems at cluster scales and represent an
alternative to dark matter haloes.

##  1 Introduction

The Cosmic Microwave Background (CMB) temperature anisotropies measured by the
Wilkinson Microwave Anisotropy Probe (WMAP) [2] and the Planck satellites [3,
4, 5, 6] as well as luminosity distances measured of the Supernovae Type Ia
(SNeIa) [1], strongly favor the concordance $\Lambda$CDM model, where the
Universe is dominated by two unknown energy densities: Dark Energy (DE) and
Dark Matter (DM). As an alternative, analytical $f(R)$-models modify the
description of gravity instead of requiring extra matter components. In these
extensions of General Relativity (GR), the Hilbert-Einstein action is replaced
with a more general function of the Ricci scalar $R$ [7]. These models have
been found to introduce negligible corrections at the scale of the Solar
systems [8] so they verify the constraints imposed by the classical test of
GR, and the perturbations in self-gravitating systems behave similarly as in
GR [9]. Nevertheless, this modified gravity introduces a Yukawa-like
correction to the gravitational potential in the Post-Newtonian limit [7, 10].
This modified potential has been used to explain the dynamics of galaxies
without requiring DM [11, 12].

Since, the exact functional form of the Lagrangian is unknown, it is important
to test potential models using all the available data on self-gravitating
systems. Galaxy clusters are reservoirs of hot gas and when CMB photons cross
the cluster potential well, they are scattered off by the free electron of the
intracluster medium. The effect, known as Sunyaev-Zeldovich (SZ) [13, 14],
induces secondary anisotropies that have been extensively measured [15, 16, 5,
17]. Pressure profiles of cluster of galaxies based on the Navarro-Frenk-White
(NFW, [18]) DM density profile have been found to be in agreement with SZ [19]
and X-ray observations [20]. We have constructed the pressure profiles of
clusters of galaxies assuming the gas is in hydrostatic equilibrium within a
Newtonian potential with a Yukawa correction and compared them with Planck
data. This article is organized as follows: in Sec. 2, we describe the weak
field limit of $f(R)$ gravity; in Sec. 3, we summarize the pressure profiles
most commonly used in the literature and compare them with the pressure
profile for $f(R)$ models; in Sec. 4 we describe the data used in our
analysis; in Sec. 5 we present our results and in Sec. 6 we will summarize our
main conclusions.

##  2 Yukawa corrections to the Newtonian potential in $f(R)$-gravity

In Extended Theories of Gravity (ETGs), the field equations are obtained by
varying the action

${\cal A}=\frac{c^{4}}{16\pi G}\int d^{4}x\sqrt{-g}f(R)+{\cal L}_{m}\,,$ (1)

resulting in 111${\displaystyle f^{\prime}(R)={df(R)}/{dR}}$ is the first
derivative with respect to the Ricci scalar,
$\Box_{g}={{}_{;\sigma}}^{;\sigma}$ is the d’Alembertian with covariant
derivatives, $\displaystyle
T_{\mu\nu}=-2(-g)^{-1/2}{\delta(\sqrt{-g}\mathcal{L}_{m})/\delta g^{\mu\nu}}$
is the matter energy-momentum tensor, $T$ its trace, $g$ the determinant of
the metric tensor $g_{\mu\nu}$. Greek indices run from $0$ to $3$.

$f^{\prime}(R)R_{\mu\nu}-\frac{f(R)}{2}\,g_{\mu\nu}-f^{\prime}(R)_{;\mu\nu}+g_{\mu\nu}\Box_{g}f^{\prime}(R)\,=\,8\pi
GT_{\mu\nu}\,.$

To construct spherically symmetric solutions, we write the metric element as

$ds^{2}=g_{tt}c^{2}dt^{2}-g_{rr}dr^{2}-r^{2}d\Omega,$ (2)

where $d\Omega$ is the solid angle. We restrict our study to
$f(R)$-Lagrangians that are expandable in Taylor series around a fixed point
$R_{0}$

$f(R)=\sum_{n}\frac{f^{n}(R_{0})}{n!}(R-R_{0})^{n}\simeq
f_{0}+f^{\prime}_{0}R+\frac{f^{\prime\prime}_{0}}{2}R^{2}+...\,.$ (3)

In the Post-Newtonian limit of these models the gravitational potential can be
written as [7]

$\Phi_{grav}(r)=\frac{\Phi_{N}(r)}{(1+\delta)}\left(1+\delta
e^{-\frac{r}{L}}\right).$ (4)

In this expression, $\delta$ represents the deviation from GR at zero order
and the Newtonian limit of GR is recovered at $\delta\rightarrow 0$,
irrespective of the scale parameter $L$. The latter is the scale length of the
self-gravitating object [7, 10] and its effects are negligible at Solar System
scales where GR, is restored [10]. In terms of the coefficients in the Taylor
expansion they are given by $\delta=1-f^{\prime}_{0}=1$ and
$L=[-{f^{\prime}_{0}}/{(6f^{\prime\prime}_{0})}]^{1/2}$. In this context, the
dynamical effects of DM are now given by the modifications of the Newtonian
potential.

##  3 Cluster pressure profiles in $f(R)$ gravity

SZ effect induces two different type of secondary temperature anisotropies on
the CMB: a thermal contribution due to the motion of the electrons within the
cluster potential well and a kinematic one (KSZ) due to the peculiar motion of
the cluster. It is given by

$\frac{\Delta T}{T_{0}}=g(\nu)\frac{k_{B}\sigma_{T}}{m_{e}c^{2}}\int
n_{e}T_{e}dl,$ (5)

where $T_{e}$ is the electron temperature, $n_{e}$ the electron density and
the integration is carried out along the line of sight $l$. In Eq. (5) $k_{B}$
is the Boltzmann constant, $m_{e}c^{2}$ the electron annihilation temperature,
$c$ the speed of light, $\nu$ the frequency of observation, $\sigma_{T}$ the
Thomson cross section and $T_{0}$ the mean temperature of the CMB. Finally,
$g(\nu)=x\coth(x/2)-4$ is the frequency dependence of the TSZ effect, with the
reduced frequency $x=h\nu/KT_{0}$. For individual clusters, only the TSZ
anisotropy has been measured.

Eq. (5) shows that in order to compute the TSZ anisotropies we need to
estimate the cluster pressure profile, $n_{e}T_{e}$. Several cluster profiles,
based on X-ray data and numerical simulations, have appeared in the
literature:

$\bullet$ The isothermal $\beta$-model [21, 22] specifies the electron
density: $n_{e}(r)=n_{e,0}[1+(r/r_{c})^{2}]^{-3\beta/2}$. From X-ray surface
brightness data it has been found that $\beta=0.6-0.8$ [23].

$\bullet$ The universal profile [20] is a phenomenological model based on the
NFW DM-profile

$p(x)\equiv\frac{P_{0}}{(c_{500}x)^{\gamma_{a}}[1+(c_{500}x)^{\alpha_{a}}]^{(\beta_{a}-\gamma_{a})/\alpha_{a}}},$
(6)

where $x$ is the radial distance in units of $r_{500}$, the radius where the
average density is 500 times the critical density, and $c_{500}$ is the
concentration parameter at $r_{500}$. Different groups have fit the model
parameters $[c_{500},\alpha_{a},\beta_{a},\gamma_{a},P_{0}]$ to X-ray or CMB
data; their best fit values are given in Table 1.

Table 1: Parameters of the GNFW corresponding to [20, 15, 24], $\beta$ for Coma cluster and $f(R)$ is the best fit model to Planck data. _Model_ | $c_{500}$ | $\alpha_{a}$ | $\beta_{a}$ | $\gamma_{a}$ | $P_{0}$
---|---|---|---|---|---
_Arnaud_ | 1.177 | 1.051 | 5.4905 | 0.3081 | $8.403h_{70}^{3/2}$
_Planck_ | 1.81 | 1.33 | 4.13 | 0.31 | 6.41
_Sayers_ | 1.18 | 0.86 | 3.67 | 0.67 | 4.29
_$\beta$_ | $\beta$ | $n_{c,0}/m^{-3}$ |  $r_{c}$/Mpc |  $T_{e}$/keV |
| 2/3 | 3860. | 0.25 | 6.48 |
_$f(R)$_ | $\delta$ | L/Mpc | $\gamma$ |  |
| -0.98 | 0.1 | 1.2 |  |

To compute in $f(R)$ gravity the pressure profile $n_{e}T_{e}$ of Eq. (5) we
assume that the gas is in hydrostatic equilibrium within the modified
potential well of the cluster (without the DM contribution)

$\frac{dP(r)}{dr}=-\rho(r)\frac{d\Phi_{grav}(r)}{dr},$ (7)

and the physical state of the gas is described by a polytropic equation of
state

$P(r)\propto\rho^{\gamma}(r).$ (8)

Eqs. (7) and (8), together with the continuity equation

$\frac{dM(r)}{dr}=4\pi\rho(r),$ (9)

and the cluster gravitational potential given by eq. (4) form a close system
that can be solved numerically to obtain the pressure profiles of any given
cluster as a function of two gravitational parameters $(\delta,L)$ and the
polytropic index $\gamma$. Fig. 1 we represent the different profiles
integrated along the line of sight with the parameters given in Table 1. The
profiles are particularized for the parameters of the Coma cluster.

Figure 1:  Pressure profiles integrated along the line of sight for the Coma
cluster ($z=0.023$). We plot three GNFW profiles (dashed, solid and dash-
dotted lines), one $\beta=2/3$ model (long dashed line) and a $f(R)$ model
(red solid line). The parameters of each model are given in Table 1.

##  4 Data.

We use Planck data and a X-ray selected cluster catalog to constrain the
pressure profiles of clusters of galaxies.

$\bullet$ Our cluster catalog is fully described in [25] and lists position,
redshift, X-ray flux, X-ray luminosity. The X-ray electron temperature is
derived from the $L_{X}-T_{X}$ relation of [26]. A total of 579 clusters are
located outside the minimal Planck mask; the region where the CMB data is less
affected by galactic foreground residuals.

$\bullet$ The Planck Collaboration released in March 2013 nine maps of the CMB
temperature anisotropies spanning a frequency range from 30 to 845 GHz. Since
foreground cleaned maps at all frequencies were not made available, we
performed our analysis on a foreground cleaned map known as SMICA. This map
has been constructed using a component separation method by combining the data
at all frequencies [27]. The SMICA map has $5^{\prime}$ resolution.

For all except the brightest clusters, the underlying CMB anisotropies have
large amplitude than the TSZ contributions. Then, to compare the expected
profile with the data, we stack the signal of all the clusters in our sample.
The anisotropy is averaged in rings of width $\theta_{500}/2$, where
$\theta_{500}$ is the angular scale subtended by the $r_{500}$. For each data
point an error bar is obtained by evaluating 1,000 times the average profiles
at 579 random positions in the SMICA map. To avoid overlapping real and
simulated clusters we removed a disc of $80^{\prime}$ around each of the
clusters in our sample.

Using eqs. (4), (7), (8) and (9), we constructed the pressure profile of all
clusters in the data as a function of three parameters: $(\delta,L,\gamma)$.
We integrated the profiles along the line of sight and convolved them with a
Gaussian beam with the same resolution of the SMICA map. We considered two
parameterizations. In (A) $L=\zeta r_{500}$ differs from cluster to cluster
but scales linearly with $r_{500}$; in this parametrization we assume that
$\zeta=[0.1,4]$. In (B) $L$ has the same value for all the clusters; it varies
in the range $L=[0.1,20]$Mpc, from the scale of a typical core radius to the
mean cluster separation. Since if $\delta<-1$ the potential is repulsive and
if $\delta=-1$ the potential diverges, we took $\delta=[-0.99,1.0]$. Finally,
the polytropic index $\gamma$ was set to vary in the range $\gamma=[1.0,1.6]$,
that corresponds to an isothermal and adiabatic monoatomic gas, respectively.
We took 30 equally spaced steps in all intervals.

##  5 Results and discussion

We computed the likelihood function $\log{\cal L}=-\chi^{2}/2$ as

$\chi^{2}({p})=\Sigma_{i,j=0}^{N}(y({p},x_{i})-d(x_{i}))C_{ij}^{-1}(y({p},x_{j})-d(x_{j}))$
(10)

where $N=7$ is the number of data points, and ${p}=(\delta,L,\gamma)$. In eq.
(10), $d(x_{i})$ is the data and $C_{i,j}$ is the correlation function between
the average temperature anisotropy on discs and rings in bins of size
$\theta_{500}$. The process was computed on 579 random positions outside the
cluster locations and repeated 1,000 times.

In Figs. 2 and 3, we computed the 2D contours at the 68% and 95% confidence
levels of the marginalized likelihoods of pairs of parameters of Model A and
Model B, respectively. Since the contours are not closed, we can only quote
the following upper limits at the same level of confidence:
$\delta<-0.46,-0.10$, $\zeta<2.5,3.7$ and $\gamma>1.35,1.12$ for Model A and
$\delta<-0.43,-0.08$, $L<12,19$ Mpc and $\gamma>1.45,1.2$ for Model B. Notice
that $L$ is similar in both A and B and $\gamma$ dominated by the physical
boundary imposed.

Although model parameters are weakly constrained by the data, the value
$\delta=0$ is excluded at more than 95% confidence level. Since $\delta\simeq
0$ corresponds to the standard Newtonian potential without DM, our result
indicates that a self-gravitating gas, without DM or modified gravity, does
not fit the data.

Figure 2: Confidence contours for pairs of parameters of Model A. Contours are
at the 68% and 95% confidence level.
Figure 3: Same as in Fig. 2 for Model B.

##  6 Conclusions

We have constrained the analytical $f(R)$ model of gravity by comparing the
pressure profile of clusters of galaxies with CMB data. We have demonstrated
that these models can accurately fit cluster profiles without a potential well
dominated by DM. In this context, $f(R)$ gravity can be considered a viable
alternative to DM. Our results also showed that a profile made of gas, without
a modified potential or without significant DM contributions could not fit the
observed profiles.

Our constraints on model parameters are rather weak since we only used the
data on the cluster profile but we did not use the frequency dependence of the
TSZ effect, the reason being that the Planck Collaboration has not released
foreground cleaned maps at each frequency. Then, our results could be
substantially improved when including this frequency information, opening the
possibility of constraining - ruling out the modified Newtonian potential at
this scale.

## References

## References

  * [1] Suzuki N et al 2012 ApJ 746, 85
  * [2] Hinshaw G et al 2013 ApJS 208, 19
  * [3] Planck Collaboration. “Planck 2013 results. XV. CMB power spectra and likelihood”, 2014 A&A 571, 15
  * [4] Planck Collaboration. “Planck 2013 results. XVI. Cosmological parameters”, 2014 A&A 571, 16
  * [5] Planck Collaboration. “Planck 2013 results. XX. Cosmology from Sunyaev–Zeldovich cluster counts”, 2014 A&A 571, 20
  * [6] Planck Collaboration. “Planck 2013 results. XXI. Cosmology with the all-sky Planck Compton parameter $y$-map”, 2014 A&A 571, 21
  * [7] Capozziello S and De Laurentis M 2011 Physics Reports 509, 167
  * [8] Berry C P L and Gair J L 2011 Phys Rev D 83, 104022
  * [9] Capozziello S De Laurentis M de Martino I Formisano M and Odintsov S 2012 Phys Rev D 85, 044022
  * [10] Capozziello S and De Laurentis M 2012 Ann. Phys. 524, 1
  * [11] Cardone V F and Capozziello S 2011 MNRAS 414, 1301
  * [12] Napolitano N R Capozziello S Romanowsky A J Capaccioli M and Tortora C 2012 ApJ 748, 87
  * [13] Sunyaev R A and Zeldovich Y B 1972 Comments on Astrophys. Space Phys. 4 173
  * [14] Sunyaev R A and Zeldovich Y B 1980 MNRAS 190, 413
  * [15] Planck Collaboration. “Planck intermediate results V: Pressure profiles of galaxy clusters from the Planck Survey”, 2013 A&A 550, 131
  * [16] Planck Collaboration. “Planck intermediate results X: Physics of the hot gas in the Coma cluster”, 2013 A&A 554, 140
  * [17] Planck Collaboration. “Planck 2013 results. XXIX. Planck catalogue of Sunyaev-Zeldovich sources”, 2014 A&A 571, 29
  * [18] Navarro J F Frenk C S and White S D M 1997 ApJ 490, 493
  * [19] Atrio-Barandela F Kashlinsky A Kocevski D and Ebeling H 2008 ApJ 675, L57
  * [20] Arnaud M Pratt G W Piffaretti R Böhringer H Croston J H and Pointecouteau E 2010 A&A 517, 92
  * [21] Cavaliere A and Fusco-Femiano R 1976 A&A 49, 137
  * [22] Cavaliere A and Fusco-Femiano R 1978 A&A 70, 677
  * [23] Jones C and Forman W 1984 ApJ 276, 38
  * [24] Sayers J et al 2013 ApJ 768, 177
  * [25] Kocevski D and Ebeling H 2006 ApJ 645, 1043
  * [26] White D A Jones C and Forman W 1997 MNRAS 292, 419
  * [27] Planck Collaboration. “Planck 2013 results. XII. Component separation”, 2014 A&A 571, A12


####  C.2.2 Pure-dephasing processes

We consider pure-dephasing processes due to $1/f$ flux noise affecting the
coupler flux bias. However, modeling $1/f$ noise using a Lindblad master
equation is rather inconvenient because it involves time-dependent rates that
also depend on the details of the noise [59]. To study the gate fidelity using
a simpler metric, we derive a multi-level white-noise approximation to the
pure-dephasing rates $\gamma^{\varphi}_{\alpha}$, parametrized by the
coupler’s pure-dephasing time $T_{\varphi}$.

Multi-level white-noise approximation– We consider a split-transmon model for
the coupler and we account for flux noise via the effective Josephson energy

$E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})=E_{\mathrm{J}_{\Sigma}}\cos\left(\frac{\phi_{\mathrm{ext}}}{2}\right)\sqrt{1+d_{r}^{2}\tan^{2}\left(\frac{\phi_{\mathrm{ext}}}{2}\right)},$
(51)

where $E_{\mathrm{J}_{\Sigma}}=E_{\mathrm{J}_{1}}+E_{\mathrm{J}_{2}}$ adds the
Josephson energies of the junctions that form the SQUID loop,

$d_{r}=\frac{r-1}{r+1},$ (52)

and $r=E_{\mathrm{J}_{1}}/E_{\mathrm{J}_{2}}>1$ is the junction assymmetry
parameter [20].

Carrying out a derivation similar to the one outlined in Ref. [59], but
assuming a noise-spectral density $S[\omega]$ that is nonsingular at zero
frequency, we arrive at the multi-level dephasing rates

$\gamma^{\varphi}_{\alpha\beta}=\hbar^{-2}S[0]\langle\Phi_{\alpha}|\partial_{\phi_{\mathrm{ext}}}\hat{H}_{\mathrm{s}}|\Phi_{\alpha}\rangle\langle\Phi_{\beta}|\partial_{\phi_{\mathrm{ext}}}\hat{H}_{\mathrm{s}}|\Phi_{\beta}\rangle,$
(53)

where $\\{|\Phi_{\alpha}\rangle\\}$ are the eigenstates of the full device
Hamiltonian. In the two-state subspace spanned by
$\\{|\Phi_{\alpha}\rangle,|\Phi_{\beta}\rangle\\}$, the coherence
$\rho_{\alpha\beta}$ evolves in the interaction frame according to

$\dot{\rho}_{\alpha\beta}=-\frac{1}{2}\left(\gamma^{\varphi}_{\alpha\alpha}+\gamma^{\varphi}_{\beta\beta}-2\gamma^{\varphi}_{\alpha\beta}\right)\rho_{\alpha\beta}.$
(54)

Therefore, the coupler coherence time under pure-dephasing noise is given by

$T_{\varphi}=\frac{1}{\frac{\gamma^{\varphi}_{\alpha\alpha}+\gamma^{\varphi}_{\alpha\beta}}{2}-\gamma^{\varphi}_{\alpha\beta}},$
(55)

where $\alpha=000$ and $\beta=001$ correspond to the hybridized ground and
first-excited states of the coupler mode.

In our simulations, we consider $T_{\varphi}$ as a parameter and use Eq. 55 to
infer the multi-level pure-dephasing rates. We do so by first approximating
the matrix elements in Eq. 53 for the case of a transmon coupler, arriving at
the expression

$\frac{1}{T_{\varphi}}=S[0]\times\frac{E_{\mathrm{C}_{\mathrm{c}}}}{E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})}\times[\partial_{\phi_{\mathrm{ext}}}E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})/\hbar]^{2},$
(56)

which we rewrite as

$S[0]=\frac{1}{T_{\varphi}}\times\frac{E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})}{E_{\mathrm{C}_{\mathrm{c}}}}\times\frac{1}{[\partial_{\phi_{\mathrm{ext}}}E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})/\hbar]^{2}}.$
(57)

We use Eq. 57 to compute the rates in Eq. 53 for all pair of device
eigenstates in the model.

Comparison against 1/f noise estimations– We now consider a model of flux
noise where the spectral noise density takes the more realistic form

$S[f]=\frac{A_{\Phi}^{2}}{|f|}.$ (58)

Here, $A_{\Phi}$ quantifies the flux noise amplitude that is typically in the
range $10^{-6}-10^{-5}\,\mu\Phi_{0}$, with $\Phi_{0}=h/2e$ [60]. Next, we
follow the derivation in Ref. [59], arriving at an equation similar to Eq. 54.
However, because of the singular behavior of Eq. 58 at low frequency, the
pure-dephasing rates are now time-dependent. Taking this fact into
consideration, we define the coherence time $T_{\varphi}$ according to the
relation $\rho_{\alpha\beta}(T_{\varphi})=\rho_{\alpha\beta}(0)/e$, arriving
at the implicit equation

$2\times(2\pi)^{2}\times\left(\frac{A_{\Phi}}{\Phi_{0}}\right)^{2}\times
T_{\varphi}^{2}\left[\left(\frac{3}{2}-\gamma\right)-\log(\omega_{\mathrm{ir}}T_{\varphi})\right]\times\frac{E_{\mathrm{C}_{\mathrm{c}}}}{E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})}\times[\partial_{\phi_{\mathrm{ext}}}E_{\mathrm{J}_{\mathrm{c}}}(\phi_{\mathrm{ext}})/\hbar]^{2}=1,$
(59)

where we have introduced the infrared cutoff $\omega_{\mathrm{ir}}/2\pi\simeq
1$ Hz and the Euler constant $\gamma\approx 0.58$.

Figure 13b and c show the coherence times obtained by numerically solving Eq.
59 for a split-transmon qubit of frequency $5\,\mathrm{GHz}$ and
$-300\,\mathrm{MHz}$ anharmonicity. The coherence times are shown as a
function of the SQUID-junction assymmetry parameter $r$, for different values
of $A_{\Phi}$ that are typical in experiments. Panel b shows the result for
the external flux bias $\Phi_{\mathrm{ext}}=0.05\,\Phi_{0}$, while panel c
considers the case of $\Phi_{\mathrm{ext}}=0.25\,\Phi_{0}$. We observe that
coherence times as large as $100\,\mu$s are in principle possible for split-
transmon couplers with large junction assymmetry that opperate closer to their
upper sweet spot. We consider the SQUID-junction asymmetry ratio $r=7$ and
$\Phi_{\mathrm{ext}}=0.25$ for all simulations in the main text.

##  Appendix D Two-step frequency allocation

In this section, we describe additional details of the frequency allocation
method outlined in Sec. V.1.

###  D.1 Multi-qubit Hamiltonian

As in Sec. II.2, we model all circuit modes as Kerr nonlinear oscillators with
frequencies $\\{\omega_{{\mu}}\\}$ and anharmonicities $\\{\alpha_{{\mu}}\\}$.
For simplicity, we assume uniform first-neighbor (qubit-coupler) coupling
$J_{1}$ and next-neighbor (qubit-qubit and coupler-coupler) coupling $J_{2}$.
We also consider a single coupler drive at a time. The drive amplitude is
$\Omega$, and we work in a frame rotating at the drive frequency $\omega$. In
what follows, we assume that the couplings $J_{1}$ and $J_{2}$ are given, and
optimize over the mode parameters and the drive frequencies. To make sure that
the perturbative result is accurate, we work with moderate coupling parameters
$J_{1}/2\pi=25\,\mathrm{MHz}$, $J_{2}/2\pi=2\,\mathrm{MHz}$, and drive
strength $\Omega/2\pi=50\,\mathrm{MHz}$. We consider a total frequency
bandwidth for the qubit and coupler modes of $1.5\,\mathrm{GHz}$, and
anharmonicities in the range of $[-350,-250]\,\mathrm{MHz}$.

We write the unit-cell Hamiltonian $\hat{H}=\hat{H}^{0}+\eta\hat{V}$ as in
Sec. II.2.1, where $\hat{H}^{0}$ denotes the noninteracting part. Once the
interaction $\eta\hat{V}$ is turned on, the self-energies associated with the
computational eigenstates are approximated to second (leading) order in the
couplings (drive amplitude). In the limit where the self-energy
$\Sigma_{\alpha}$ of a computational state $|\Phi_{\alpha}\rangle$ is much
smaller than the detunings ${\epsilon}_{\alpha}^{0}-{\epsilon}_{\beta}^{0}$,
where $\beta$ labels a noncomputational state, the self-consistent nature of
Eq. 36 can be simplified by omitting the self-energy from all denominators in
that expression. In this approximation, a two-level truncation near the
resonance ${\epsilon}_{\alpha}^{0}-{\epsilon}_{\beta}^{0}\to 0$ leads to the
dispersive form of the self-energy
$\Sigma_{\alpha}\approx|G_{\alpha\beta}|^{2}/({\epsilon}_{\alpha}^{0}-{\epsilon}_{\beta}^{0})$,
where

$\begin{split}G_{\alpha\beta}=\sum_{k=1}^{\infty}\sum_{\alpha_{1},\dots,\alpha_{k}}&\langle\Phi_{\beta}^{0}|\eta\hat{V}|\Phi_{\alpha_{1}}^{0}\rangle\frac{\langle\Phi_{\alpha_{1}}^{0}|\eta\hat{V}|\Phi_{\alpha_{2}}^{0}\rangle}{{\epsilon}_{\alpha}^{0}-{\epsilon}_{\alpha_{1}}^{0}}\dots\\\
&\dots\frac{\langle\Phi_{\alpha_{k}}^{0}|\eta\hat{V}|\Phi_{\alpha}^{0}\rangle}{{\epsilon}_{\alpha}^{0}-{\epsilon}_{\alpha_{k}}^{0}},\end{split}$
(60)

is the virtual interaction rate between the eigenstates
$|\Phi_{\alpha}^{0}\rangle$ and $|\Phi_{\beta}^{0}\rangle$ of the bare
Hamiltonian $\hat{H}^{0}$. Here, we assume that $|G_{\alpha\beta}|$ is small
compared to $|{\epsilon}_{\alpha}^{0}-{\epsilon}_{\beta}^{0}|$. We use this
simplified expression to compute the self-energies of computational states in
the unit cell.

###  D.2 Cost function and optimization

Our two-step frequency optimization seeks to: i) minimize the static ZZ
interaction between all pairs of qubits, ii) minimize drive-activated ZZ
interactions on spectator qubits, iii) maximize desired two-qubit gate rates.
The cost function also incorporates the dispersive-coupling ratios
$\nu_{\alpha\beta}=|G_{\alpha\beta}/({\epsilon}_{\alpha}^{0}-{\epsilon}_{\beta}^{0})|$,
where $\alpha$ and $\beta$ are computational and noncomputational states,
respectively. Because the excitation number is constant in the absence of a
drive and under a rotating-wave approximation, we estimate the static ZZ
interaction taking into account states with up to two excitations. However, we
consider computational (noncomputational) states with up to two (three)
excitations in the presence of a drive.

Static ZZ interaction– We begin by generating several frequency layouts as
independent initial conditions for the optimizer. At this point, the qubits
and couplers anharmonicities are $-300\,\mathrm{MHz}$ and
$-350\,\mathrm{MHz}$, respectively. Next, we optimize these frequency layouts
independently using a least-squares algorithm. The optimization is done over
the mode frequencies and anharmonicities. The cost function incorporates
penalties to ensure that both the dispersive-coupling ratios and the the ZZ
interaction between any pair of qubits are below the chosen bounds of 0.05 and
$20\,\mathrm{kHz}$, respectively.

Two-qubit gate rate– Each pair of qubits offers many potential operating
points for the two-qubit gate, depending on the drive frequency. However, not
all resonances are equivalent when it comes to minimizing the impact of the
coupler drive on neighboring qubits. Thus, we take into account all possible
single- and two-photon resonances for a given pair of qubits, and computing
the driven ZZ interaction for any other qubit pair in the device.

More precisely, our cost function incorporates: i) the driven ZZ interaction
between the qubits that participate of the two-qubit gate, ii) the ratio
between the driven ZZ interaction for other pair of qubits and that calculated
in i), and iii) the ratios $\nu_{\alpha\beta}$ for qubits that do not
participate of the gate. The cost is evaluated for all coupler-drive
frequencies corresponding to a transition between a computational and a
noncomputational state of a specific two-qubit pair, and then added up for all
pair of neighboring qubits in the device. While information in i-iii) is
somewhat redundant, we find that partitioning the cost function in this way
gives us enough flexibility for the optimization. Only the coupler-drive
frequencies are optimized during this second step.

##  Appendix E Expanded-space Hamiltonian for multiple drives

With the definitions introduced in Sec. V.2, the system Hamiltonian [see Eq.
7] takes the form

$\hat{H}(\hat{\bm{\phi}},\hat{\bm{n}},t)=\hat{H}_{\mathrm{s}}(\hat{\bm{\phi}},\hat{\bm{n}})+\hat{H}_{\mathrm{drive}}[\hat{\bm{\phi}},\hat{\bm{n}},\bm{\Omega}(t),\bm{\theta}(t)].$
(61)

Next, we promote the drive-phase variables to a vector operator
$\bm{\theta}(t)\to\hat{\bm{\vartheta}}$, with conjugate momenta
$\hat{\bm{m}}\to-i{\partial}_{\bm{\vartheta}}$ in the phase representation.
Under this transformation, the effective Hamiltonian that extends Eq. 9 to $d$
microwave drives takes the form

$\hat{H}_{\mathrm{eff}}(\hat{\bm{\vartheta}},\hat{\bm{m}},\hat{\bm{\phi}},\hat{\bm{n}},t)=\hat{H}(\hat{\bm{\phi}},\hat{\bm{n}},t)+\hbar\bm{\omega}_{\mathrm{eff}}(t)\cdot\hat{\bm{m}},$
(62)

where $\bm{\omega}_{\mathrm{eff}}(t)=\dot{\bm{\theta}}(t)$ groups the drive
frequencies.

## References

  * Takita  _et al._ [2016] M. Takita, A. D. Córcoles, E. Magesan, B. Abdo, M. Brink, A. Cross, J. M. Chow, and J. M. Gambetta, Demonstration of weight-four parity measurements in the surface code architecture, Physical Review Letters 117, 210505 (2016).
  * Brink  _et al._ [2018] M. Brink, J. M. Chow, J. Hertzberg, E. Magesan, and S. Rosenblatt, Device challenges for near term superconducting quantum processors: frequency collisions, in  _2018 IEEE International Electron Devices Meeting (IEDM)_ (IEEE, 2018) pp. 6–1.
  * Ku  _et al._ [2020] J. Ku, X. Xu, M. Brink, D. C. McKay, J. B. Hertzberg, M. H. Ansari, and B. Plourde, Suppression of unwanted z z interactions in a hybrid two-qubit system, Physical Review Letters 125, 200504 (2020).
  * Krinner  _et al._ [2020a] S. Krinner, S. Lazar, A. Remm, C. K. Andersen, N. Lacroix, G. J. Norris, C. Hellings, M. Gabureac, C. Eichler, and A. Wallraff, Benchmarking coherent errors in controlled-phase gates due to spectator qubits, Physical Review Applied 14, 024042 (2020a).
  * Chow  _et al._ [2012] J. M. Chow, J. M. Gambetta, A. D. Corcoles, S. T. Merkel, J. A. Smolin, C. Rigetti, S. Poletto, G. A. Keefe, M. B. Rothwell, J. R. Rozen, _et al._ , Universal quantum gate set approaching fault-tolerant thresholds with superconducting qubits, Physical Review Letters 109, 060501 (2012).
  * Kandala  _et al._ [2021] A. Kandala, K. Wei, S. Srinivasan, E. Magesan, S. Carnevale, G. Keefe, D. Klaus, O. Dial, and D. McKay, Demonstration of a high-fidelity cnot gate for fixed-frequency transmons with engineered z z suppression, Physical Review Letters 127, 130501 (2021).
  * Chow  _et al._ [2013] J. M. Chow, J. M. Gambetta, A. W. Cross, S. T. Merkel, C. Rigetti, and M. Steffen, Microwave-activated conditional-phase gate for superconducting qubits, New Journal of Physics 15, 115012 (2013).
  * Krinner  _et al._ [2020b] S. Krinner, P. Kurpiers, B. Royer, P. Magnard, I. Tsitsilin, J.-C. Besse, A. Remm, A. Blais, and A. Wallraff, Demonstration of an all-microwave controlled-phase gate between far-detuned qubits, Physical Review Applied 14, 044039 (2020b).
  * Mitchell  _et al._ [2021] B. K. Mitchell, R. K. Naik, A. Morvan, A. Hashim, J. M. Kreikebaum, B. Marinelli, W. Lavrijsen, K. Nowrouzi, D. I. Santiago, and I. Siddiqi, Hardware-efficient microwave-activated tunable coupling between superconducting qubits, Physical Review Letters 127, 200502 (2021).
  * Berke  _et al._ [2020] C. Berke, E. Varvelis, S. Trebst, A. Altland, and D. P. DiVincenzo, Transmon platform for quantum computing challenged by chaotic fluctuations, arXiv preprint arXiv:2012.05923 (2020).
  * Hertzberg  _et al._ [2021] J. B. Hertzberg, E. J. Zhang, S. Rosenblatt, E. Magesan, J. A. Smolin, J.-B. Yau, V. P. Adiga, M. Sandberg, M. Brink, J. M. Chow, _et al._ , Laser-annealing josephson junctions for yielding scaled-up superconducting quantum processors, npj Quantum Information 7, 1 (2021).
  * DiCarlo  _et al._ [2009] L. DiCarlo, J. M. Chow, J. M. Gambetta, L. S. Bishop, B. R. Johnson, D. Schuster, J. Majer, A. Blais, L. Frunzio, S. Girvin, _et al._ , Demonstration of two-qubit algorithms with a superconducting quantum processor, Nature 460, 240 (2009).
  * Neeley  _et al._ [2010] M. Neeley, R. C. Bialczak, M. Lenander, E. Lucero, M. Mariantoni, A. O’connell, D. Sank, H. Wang, M. Weides, J. Wenner, _et al._ , Generation of three-qubit entangled states using superconducting phase qubits, Nature 467, 570 (2010).
  * Chen  _et al._ [2014] Y. Chen, C. Neill, P. Roushan, N. Leung, M. Fang, R. Barends, J. Kelly, B. Campbell, Z. Chen, B. Chiaro, _et al._ , Qubit architecture with high coherence and fast tunable coupling, Physical Review Letters 113, 220502 (2014).
  * Barends  _et al._ [2019] R. Barends, C. Quintana, A. Petukhov, Y. Chen, D. Kafri, K. Kechedzhi, R. Collins, O. Naaman, S. Boixo, F. Arute, _et al._ , Diabatic gates for frequency-tunable superconducting qubits, Physical Review Letters 123, 210501 (2019).
  * McKay  _et al._ [2016] D. C. McKay, S. Filipp, A. Mezzacapo, E. Magesan, J. M. Chow, and J. M. Gambetta, Universal gate for fixed-frequency qubits via a tunable bus, Physical Review Applied 6, 064007 (2016).
  * Caldwell  _et al._ [2018] S. Caldwell, N. Didier, C. Ryan, E. Sete, A. Hudson, P. Karalekas, R. Manenti, M. da Silva, R. Sinclair, E. Acala, _et al._ , Parametrically activated entangling gates using transmon qubits, Physical Review Applied 10, 034050 (2018).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. Brandao, D. A. Buell, _et al._ , Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019).
  * Rol  _et al._ [2020] M. A. Rol, L. Ciorciaro, F. K. Malinowski, B. M. Tarasinski, R. E. Sagastizabal, C. C. Bultink, Y. Salathe, N. Haandbæk, J. Sedivy, and L. DiCarlo, Time-domain characterization and correction of on-chip distortion of control pulses in a quantum processor, Applied Physics Letters 116, 054001 (2020).
  * Hutchings  _et al._ [2017] M. Hutchings, J. B. Hertzberg, Y. Liu, N. T. Bronn, G. A. Keefe, M. Brink, J. M. Chow, and B. Plourde, Tunable superconducting qubits with flux-independent coherence, Physical Review Applied 8, 044003 (2017).
  * Grifoni and Hänggi [1998] M. Grifoni and P. Hänggi, Driven quantum tunneling, Physics Reports 304, 229 (1998).
  * Mundada  _et al._ [2020] P. S. Mundada, A. Gyenis, Z. Huang, J. Koch, and A. A. Houck, Floquet-engineered enhancement of coherence times in a driven fluxonium qubit, Physical Review Applied 14, 054033 (2020).
  * Huang  _et al._ [2021] Z. Huang, P. S. Mundada, A. Gyenis, D. I. Schuster, A. A. Houck, and J. Koch, Engineering dynamical sweet spots to protect qubits from 1/f noise, Physical Review Applied 15, 034065 (2021).
  * Petrescu  _et al._ [2021] A. Petrescu, C. L. Calonnec, C. Leroux, A. Di Paolo, P. Mundada, S. Sussman, A. Vrajitoarea, A. A. Houck, and A. Blais, Accurate methods for the analysis of strong-drive effects in parametric gates, arXiv preprint arXiv:2107.02343 (2021).
  * Peskin and Moiseyev [1993] U. Peskin and N. Moiseyev, The solution of the time-dependent schrödinger equation by the (t, t’) method: Theory, computational algorithm and applications, The Journal of Chemical Physics 99, 4590 (1993).
  * Breuer  _et al._ [1988a] H. P. Breuer, K. Dietz, and M. Holthaus, The role of avoided crossings in the dynamics of strong laser field—matter interactions, Zeitschrift für Physik D Atoms, Molecules and Clusters 8, 349 (1988a).
  * Breuer  _et al._ [1988b] H. P. Breuer, K. Dietz, and M. Holthaus, Strong laser fields interacting with matter i, Zeitschrift für Physik D Atoms, Molecules and Clusters 10, 13 (1988b).
  * Breuer and Holthaus [1989] H. Breuer and M. Holthaus, Adiabatic processes in the ionization of highly excited hydrogen atoms, Zeitschrift für Physik D Atoms, Molecules and Clusters 11, 1 (1989).
  * Breuer  _et al._ [1990] H. P. Breuer, K. Dietz, and M. Holthaus, Transport of quantum states of periodically driven systems, Journal de Physique 51, 709 (1990).
  * Drese and Holthaus [1999] K. Drese and M. Holthaus, Floquet theory for short laser pulses, The European Physical Journal D-Atomic, Molecular, Optical and Plasma Physics 5, 119 (1999).
  * Guérin [1997] S. Guérin, Complete dissociation by chirped laser pulses designed by adiabatic floquet analysis, Physical Review A 56, 1458 (1997).
  * Weinberg  _et al._ [2017] P. Weinberg, M. Bukov, L. D’Alessio, A. Polkovnikov, S. Vajna, and M. Kolodrubetz, Adiabatic perturbation theory and geometry of periodically-driven systems, Physics Reports 688, 1 (2017).
  * Hone  _et al._ [1997] D. W. Hone, R. Ketzmerick, and W. Kohn, Time-dependent floquet theory and absence of an adiabatic limit, Physical Review A 56, 4045 (1997).
  * Young and Deal Jr [1970] R. H. Young and W. J. Deal Jr, Adiabatic response to an oscillatory field, Journal of Mathematical Physics 11, 3298 (1970).
  * Martínez-Garaot  _et al._ [2015] S. Martínez-Garaot, A. Ruschhaupt, J. Gillet, T. Busch, and J. G. Muga, Fast quasiadiabatic dynamics, Physical Review A 92, 043406 (2015).
  * García-Ripoll  _et al._ [2020] J. J. García-Ripoll, A. Ruiz-Chamorro, and E. Torrontegui, Quantum control of frequency-tunable transmon superconducting qubits, Physical Review Applied 14, 044035 (2020).
  * Yatsenko  _et al._ [2004] L. Yatsenko, S. Guérin, and H.-R. Jauslin, Pulse-driven near-resonant quantum adiabatic dynamics: Lifting of quasidegeneracy, Physical Review A 70, 043402 (2004).
  * McKay  _et al._ [2017] D. C. McKay, C. J. Wood, S. Sheldon, J. M. Chow, and J. M. Gambetta, Efficient z gates for quantum computing, Physical Review A 96, 022330 (2017).
  * Holthaus and Just [1994] M. Holthaus and B. Just, Generalized $\pi$ pulses, Physical Review A 49, 1950 (1994).
  * Zhang and Dykman [2017] Y. Zhang and M. Dykman, Preparing quasienergy states on demand: A parametric oscillator, Physical Review A 95, 053841 (2017).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Wood and Gambetta [2018] C. J. Wood and J. M. Gambetta, Quantification and characterization of leakage errors, Physical Review A 97, 032306 (2018).
  * Motzoi  _et al._ [2009] F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm, Simple pulses for elimination of leakage in weakly nonlinear qubits, Physical Review Letters 103, 110501 (2009).
  * Ficheux  _et al._ [2021] Q. Ficheux, L. B. Nguyen, A. Somoroff, H. Xiong, K. N. Nesterov, M. G. Vavilov, and V. E. Manucharyan, Fast logic with slow qubits: microwave-activated controlled-z gate on low-frequency fluxoniums, Physical Review X 11, 021026 (2021).
  * Paik  _et al._ [2016] H. Paik, A. Mezzacapo, M. Sandberg, D. McClure, B. Abdo, A. Córcoles, O. Dial, D. Bogorin, B. Plourde, M. Steffen, _et al._ , Experimental demonstration of a resonator-induced phase gate in a multiqubit circuit-qed system, Physical Review Letters 117, 250502 (2016).
  * Ni  _et al._ [2021] Z. Ni, S. Li, L. Zhang, J. Chu, J. Niu, T. Yan, X. Deng, L. Hu, J. Li, Y. Zhong, _et al._ , Scalable method for eliminating residual $zz$ interaction between superconducting qubits, arXiv preprint arXiv:2111.13292 (2021).
  * Didier  _et al._ [2019] N. Didier, E. A. Sete, J. Combes, and M. P. da Silva, Ac flux sweet spots in parametrically modulated superconducting qubits, Physical Review Applied 12, 054015 (2019).
  * Naaman  _et al._ [2008] O. Naaman, J. Aumentado, L. Friedland, J. Wurtele, and I. Siddiqi, Phase-locking transition in a chirped superconducting josephson resonator, Physical Review Letters 101, 117005 (2008).
  * Murch  _et al._ [2011] K. Murch, R. Vijay, I. Barth, O. Naaman, J. Aumentado, L. Friedland, and I. Siddiqi, Quantum fluctuations in the chirped pendulum, Nature Physics 7, 105 (2011).
  * Murch  _et al._ [2012] K. Murch, E. Ginossar, S. Weber, R. Vijay, S. Girvin, and I. Siddiqi, Quantum state sensitivity of an autoresonant superconducting circuit, Physical Review B 86, 220503 (2012).
  * Shalibo  _et al._ [2012] Y. Shalibo, Y. Rofe, I. Barth, L. Friedland, R. Bialczack, J. M. Martinis, and N. Katz, Quantum and classical chirps in an anharmonic oscillator, Physical Review Letters 108, 037701 (2012).
  * Salis  _et al._ [2020] G. Salis, N. Moll, M. Roth, M. Ganzhorn, and S. Filipp, Time-resolved tomography of a driven adiabatic quantum simulation, Physical Review A 102, 062611 (2020).
  * Pechal  _et al._ [2021] M. Pechal, F. Roy, S. A. Wilkinson, G. Salis, M. Werninghaus, M. J. Hartmann, and S. Filipp, Direct implementation of a perceptron in superconducting circuit quantum hardware, arXiv preprint arXiv:2111.12669 (2021).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annual Review of Condensed Matter Physics 11, 369 (2020).
  * Gyenis  _et al._ [2021] A. Gyenis, A. Di Paolo, J. Koch, A. Blais, A. A. Houck, and D. I. Schuster, Moving beyond the transmon: Noise-protected superconducting quantum circuits, PRX Quantum 2, 030101 (2021).
  * Malekakhlagh  _et al._ [2020] M. Malekakhlagh, E. Magesan, and D. C. McKay, First-principles analysis of cross-resonance gate operation, Physical Review A 102, 042605 (2020).
  * Morvan  _et al._ [2021] A. Morvan, L. Chen, J. M. Larson, D. I. Santiago, and I. Siddiqi, Optimizing frequency allocation for fixed-frequency superconducting quantum processors, arXiv preprint arXiv:2112.01634 (2021).
  * Pop  _et al._ [2014] I. M. Pop, K. Geerlings, G. Catelani, R. J. Schoelkopf, L. I. Glazman, and M. H. Devoret, Coherent suppression of electromagnetic dissipation due to superconducting quasiparticles, Nature 508, 369 (2014).
  * Di Paolo  _et al._ [2021] A. Di Paolo, T. E. Baker, A. Foley, D. Sénéchal, and A. Blais, Efficient modeling of superconducting quantum circuits with tensor networks, npj Quantum Information 7, 1 (2021).
  * Koch  _et al._ [2007] J. Koch, M. Y. Terri, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the cooper pair box, Physical Review A 76, 042319 (2007).