footmiscCommand

11institutetext: School of Information, University of Texas at Austin, TX, USA
11email<EMAIL_ADDRESS>
22institutetext: Dept. of Computer Eng., TOBB Economy and Tech. University,
Ankara, Turkey
22email<EMAIL_ADDRESS>
# Understanding and Predicting the Characteristics of Test Collections

Md Mustafizur Rahman 11 Mucahid Kutlu 22 Matthew Lease 11

###### Abstract

Shared-task campaigns such as NIST TREC select documents to judge by pooling
rankings from many participant systems. Therefore, the quality of the test
collection greatly depends on the number of participants and the quality of
submitted runs. In this work, we investigate i) how the number of
participants, coupled with other factors, affects the quality of a test
collection; and ii) whether the quality of a test collection can be inferred
prior to collecting relevance judgments. Experiments on six TREC collections
demonstrate that the required number of participants to construct a high-
quality test collection varies significantly across different test collections
due to a variety of factors. Furthermore, results suggest that the quality of
test collections can be predicted.

###### Keywords:

Evaluation; Test Collections; Pooling; Reusability

##  1 Introduction

The system-based IR evaluation relies on having a test collection, which
consists of a document collection, a set of search topics, and relevance
judgments for selected document-topic pairs. Typically a test collection is
constructed by organizing a shared-task where the organizers (e.g., NIST TREC)
usually provide a collection of documents and a set of topics developed by
experts, and ask the participating groups to submit their ranked list of
documents (i.e., runs) for each topic. Subsequently, the documents to be
judged are selected from the ranked lists of participants using various
methods such as pooling (Sparck Jones and Van Rijsbergen, 1975) and bandit
techniques (Voorhees, 2018).

While shared-tasks are extremely useful to build test collections, the cost of
test collections is still very high due to the huge amount of required human
judgments. Besides, organizing the shared task itself is also difficult, slow,
and expensive. Furthermore, documents selected for human annotation completely
depend on the participating groups and their submitted runs. Therefore, if
very few groups participate in a shared-task, the resultant test collection
might not be of high quality as desired. Then, what should task organizers do
when there are few participants? They might a) cancel the shared-task, or b)
accept the risk of building a low-quality test collection and proceed as it is
with all the costs of organizing the shared-task and collecting human
annotation, or c) increase the publicity of the shared-task to attract more
participants, or d) ask the existing participants to submit manual runs to
increase the diversity of the results.

In this work, we conduct experiments to shed light on the impact of the number
of participants on the quality of pooling-based test collections, particularly
the reusability of a test collection. We investigate the following two key
research questions. RQ-1) i) how does the number of participants coupled with
other factors affect the quality of a test collection? RQ-2) Can we predict
the quality of a test collection prior to collecting relevance judgments?

For RQ-1, our experiments vary the number of shared-task participants by down-
sampling participants from past TREC tasks in order to construct simulated
test collections. Then we analyze the impact of the interaction between the
number of participants and other factors (e.g., the number of topics,
collection size, etc.) on resultant test collection quality. For RQ-2, we
develop a model to predict test collection quality using only the number of
participants, the number of topics, collection size, and pool depth. We
analyze the generalization performance of the prediction model using our
designed “Leave-one-test-collection-out” setup111In order to ensure the
reproducibility of our findings, we will share our source code in the final
version of the manuscript..

Our analysis leads to the following recommendations for shared-tasks:

  * •

With few participants, increasing the pool depth increases test collection
quality more than increasing the number of topics.

  * •

With few participants, organizers should maximally encourage manual runs.

  * •

With larger document collections, it is particularly important that organizers
attract strong participation in order to achieve test collection quality.

  * •

To help guide track planning, organizers should consider using a predictive
model such as ours to infer the expected quality of the test collection based
on the planned pool depth and expected number of participants, so that they
can adjust course as needed if the predicted quality is low.

##  2 Factors Impacting the Qualities of Test Collections

The quality of a test collection depends on a variety of factors, such as: i)
the number of topics; ii) pool depth (assuming pool-based judging); iii) the
number of participants; iv) the collection size; and v) the types of runs
(e.g., manual runs and automatic runs) and their quality. One might also
consider vi) the target evaluation metric (e.g., MAP@1000, NDCG@10) in
assessing how well a test collection supports reliable evaluation for a given
retrieval task, as measured by a particular metric. Since constructing a test
collection is expensive, rusability is desirable. Reusability is often
measured by how a given run contributing to the pool would have been assessed
if excluded from the pool. In this study, we focus on reusability as a key
measure of the test collection quality.

While considerable work has investigated how the quality of a test collection
is impacted by the above-mentioned factors, prior studies have not explored
how the number of participants interacts with other factors.

Number of topics. Sparck Jones and Van Rijsbergen (1976) suggest that 250
topics can be acceptable, but 1000 topics are needed for reliable evaluation.
Voorhees (2000) performs an empirical analysis on the TREC6 test collection
and shows that system rankings computed based on 5 or 10 topics are relatively
unstable, whereas a set of 25 or more topics produces a stable ranking of IR
systems. Buckley and Voorhees (2000) calculate the error rate of various
evaluation measures and find that for reliable evaluation the required number
of topics should be at least 25. Webber et al. (2008) recommend that a set of
150 topics is required to statistically distinguish the performance of one IR
system from other IR systems. Zobel (1998) finds that a set of 25 topics can
reasonably predict the performance of IR systems on a separate set of 25
topics.

Pool depth. Prior work has also studied the trade-off between collecting fewer
judgments with more topics (i.e., Wide and Shallow (WaS) judging) vs. more
judgments with fewer topics (i.e., Narrow and Deep (NaD) judging). Carterette
et al. (2008) report on TREC Million Query track and conclude that WaS judging
produces a more reliable evaluation of IR systems than NaD judging. Kutlu et
al. (2018) find that NaD judging is preferred to WaS judging if we consider
intelligent topic selection or other hidden costs of shallow judging, such as
topic creation time and noisier judgments. Voorhees (2018) investigates the
impact of varying pool depth on the reusability of a test collection.

Run types. To see how different types of runs (e.g., manual and automatic)
impact collection quality, Büttcher et al. (2007) adapt the “leave-one-group-
out”(Zobel, 1998) experiment by removing all unique documents contributed by
the manual runs from the pool. The authors find that their setup ranks the
runs differently than found in the original TREC 2006 Terabyte task.

Collection size. Hawking and Robertson (2003) study the effect of the document
collection size on Very Large Collection track (Hawking et al., 1998) and
observe a higher value of P@20 for runs in the larger collection.
Interestingly, they find no change in value between the large and small
document collection when the runs are ranked based on MAP.

##  3 The Impact of Varying the Number of Participants

Our experimental design for analyzing the impact of the number of groups is
shown in Algorithm 1. First, we construct the original qrels ($Q_{o}$) using
all participants from set $G$ [Line 2]. Then we evaluate all runs using this
original qrels ($Q_{o}$) in terms of a ranking metric (e.g., MAP@1000,
NDCG@10, etc.) and store the ranking of runs in $E$ [Line 3]. Next, we change
the group number, $g$, from $1$ to $|G|$ with an increment of $1$ at each
iteration to create test collection with varying number of participants [Line
5]. At each iteration, we randomly sample $g$ number of groups ($\hat{G_{i}}$)
from the set of groups $G$ [Line 7] and construct the simulated test
collection (i.e., qrels) $Q_{g}$ using only the participants in $\hat{G_{i}}$
[Line 8]. Then, we evaluate all participating runs in set $G$ by using
simulated test collection $Q_{g}$ in terms of a ranking metric (e.g.,
MAP@1000, NDCG@10, etc.) and store these new ranking of runs in $E_{g}$ [Line
9]. We calculate the performance difference in terms of $\tau_{ap}$ (Yilmaz et
al., 2008) and Max Drop (Voorhees, 2018) (i.e., the maximum drop in a run’s
rank, between the original ranking of runs $E$ and the ranking of runs
obtained via the respective simulated test collection $E_{g}$) [Line 11]. Note
that we calculate average scores across different group samples for a
particular parameter setup. In addition, we can also utilize Algorithm 1 to
experiments with a varying number of topics and varying pool depth because it
takes the set of topics $T$ and pool depth $P$ as inputs.

1

Input : Set of groups $G$ $\bullet$ Number of samples for groups $N$ $\bullet$
Set of topics $T$ $\bullet$ Pool depth $P$

Output : E, A set of performance score indexed by group number

2

3$R\leftarrow\text{Total number of runs from all groups in G}$

4 $Q_{o}\leftarrow Construct\\_Qrels(G,T,P)$ $\triangleright$ official qrels

5 $E\leftarrow Evaluate\\_runs(R,Q_{o})$ $\triangleright$ Evaluate all runs
with $Q_{o}$

6 $\hat{E}\leftarrow\emptyset$ $\triangleright$ keeps scores of systems with
reduced qrels

7for  _group no $g\leftarrow 1$ to $|G|$_ do

8 for  _sample number $i\leftarrow 1$ to $N$_ do

9 $\hat{G_{i}}\leftarrow\text{randomly sample g groups from }G$

10 $Q_{g}\leftarrow\text{Construct\\_Qrels(}\hat{G_{i}},T,P)$

11 $E_{g}\leftarrow Evaluate\\_runs(R,Q_{g})$ $\triangleright$ Evaluate all
runs using qrels $Q_{g}$

12 $\hat{E}\leftarrow\hat{E}\cup E_{g}$

return _Evaluate_Performance_Difference(E, $\hat{E}$)_

Algorithm 1 Experimental Design

###  3.1 Datasets

We conduct our experiments on six TREC tracks and datasets: the 2013-2014 Web
Tracks on ClueWeb12222lemurproject.org/clueweb12, the 2006 Terabyte track on
Gov2333ir.dcs.gla.ac.uk/test_collections/gov2-summary.htm, and the 2004 Robust
Retrieval Task (Robust’04), the 1999 TREC-8 ad hoc track, the 1998 TREC-7 ad
hoc track on TIPSTER disks 4-5444trec.nist.gov/data/docs_eng.html (excluding
the congressional record). Table 1 provides statistics about test collections
we use. Later tracks have fewer participants than earlier tracks both in terms
of the number of groups and the submitted runs. Later tracks also tend to use
larger document collections, without commensurate increase in pool depth,
leading to an increasing prevalence of relevant documents in judged pools,
from $\sim$5% to $\sim$40%.

Table 1: Statistics about the test collections used in this study. Track |  #Groups | #Manual | #Auto |  Pool | Collection | #Topics | #Docs |  #Judged | %Rel
---|---|---|---|---|---|---|---|---|---
|  | Runs | Runs |  Depth |  |  |  |  |
WT’14  | 9  | 4  | 26  | 25 | ClueWeb12 | 50 | 52,343,021 | 14,432 | 39.2%
WT’13  | 13  | 3  | 31  | 10 and 20 | ClueWeb12 | 50 | 52,343,021 | 14,474 | 28.7%
TB’06  | 20  | 19  | 61  | 50 | Gov2 | 50 | 25,205,179 | 31,984 | 18.4%
Adhoc’99  | 40  | 9  | 62  | 100 | Disks45-CR | 50 | 528,155 | 86,830 | 5.4%
Adhoc’98  | 41  | 16  | 68  | 100 | Disks45-CR | 50 | 528,155 | 80,345 | 5.8%
Robust’04  | 14  | 0  | 110  | 100 | Disks45-CR | 249 | 528,155 | 311,410 | 5.6%

Figure 1: $\tau_{ap}$ (first column), and Max Drop (second column) obtained by
simulated test collections with a varying number of topics on the Robust’04
dataset. The x-axis represents the number of groups and the y-axis shows
results when runs are ranked using MAP (top row) and NDCG@10 (bottom row).

###  3.2 Results and Discussion

####  3.2.1 Impact of Number of Topics.

We first consider how the number of topics interacts with the number of groups
in relation to test collection quality. To explore this, we study Robust’04
using various subsets of its 249 topics, randomly sampling $m$ topics, with
$m\in\\{50,100,150,200,249\\}$. Robust’04 topics can be categorized into 5
sets: 301-350 (Adhoc’97), 350-400 (Adhoc’98), 401-450 (Adhoc’99), 601-650
(Robust’03-Hard), and 651-700 (Robust’04). To ensure coverage over the
different sets in our sampling, we apply stratified sampling. Algorithm 1
implements our experimental setup for a given topic subset (we vary the topic
subset outside of the algorithm). The pool depth is set to $100$, and we
evaluate MAP@1000 and NDCG@10.

Figure 2: $\tau_{ap}$ (first row), and Max Drop (second row) obtained by
simulated test collections with a varying number of topics along with a
varying pool depth on the Robust04 dataset. The x-axis represents the number
of groups and the y-axis shows results when runs are ranked using MAP.

Figure 1 shows $\tau_{ap}$ (larger is better) and Max Drop (smaller is better)
on the Robust’04 test collection using all runs. Each line shows an average
value of the computed performance metric (e.g., $\tau_{ap}$) across $4$
different random samples (i.e., N is set to 4 in Algorithm 1) of groups when
$m$ topics are randomly sampled from 249 topics. The first row presents
$\tau_{ap}$ and Max Drop when the runs are ranked using MAP (i.e., MAP@1000)
while the second row reports the same metric when runs are ranked using
NDCG@10. We report the Area under the Curve (AUC) for each of the line plots
and the Pearson correlation ($\rho$) between the number of groups and the
corresponding performance metrics.

Let us first consider when runs are ranked by MAP (Figure 1, first row). For
$\tau_{ap}$ and any number of groups $g$, we do not see any significant
difference in AUC when we down-sample to 50, 100, 150, or 200 topics vs. all
249 topics. Furthermore, we achieve a $\tau_{ap}$ correlation of $0.9$ using
only $3$ participating groups, irrespective of the number of topics. We
observe the same outcome when using NDCG@10 (Figure 1, bottom row): there is
no significant difference in AUC, and we only need $2$ participating groups to
achieve a $\tau_{ap}$ correlation of 0.9 or above. Since NDCG@10 is far
shallower than MAP@1000, it is reasonable to observe that we might need fewer
groups to achieve $\tau_{ap}>=0.9$ using NDCG@10. On the other hand, Max Drop
for MAP and NDCG shows a significant difference in AUC, when we down-sample
topics. We do not see any decreasing pattern in AUC for MAP or NDCG as we
increase the number of topics.

In these experiments, we see that for any given number of participating
groups, increasing the number of topics does not improve test collection
reusability. For both $\tau_{ap}$ and Max Drop metrics, this observation
holds. However, we must note the extent of our experiments. First, since we
experiment using only Robust’04 (given its 249 topics), further analysis on
other test collections is needed. Second, our experiments only vary the number
of topics above 50, so results with more spartan topic sizes might vary.
Third, we assume a fixed pool depth of $100$. Fourth, Robust’04 contains only
automatic runs; manual runs often find other unique relevant documents. Fifth,
Robust’04 dataset is relatively small compared to a modern collection such as
ClueWeb’12 (Table 1). Since larger collections tend to contain more relevant
documents, increasing the number of topics for these larger document
collections may be more valuable.

####  3.2.2 Impact of Pool Depth.

In the previous experiment, we observe how varying the number of topics
interacts with the varying number of participating groups to build a reusable
test collection while keeping the pool depth fixed. In this experiment, we
also change the pool depth along with the number of topics and the number of
groups. The experimental setup for this experiment is the same as discussed in
the previous experiment except we vary the pool depth $p$ where $p$ takes
values from the set $\\{20,40,60,80,100\\}$.

The results for this experiment are reported in Figure 2. The top row
represents results for $\tau_{ap}$ whereas the bottom rows present results for
Max Drop. The columns of Figure 2 indicate the number of topics sampled for
each of those pool depth variations. Due to space constraints, we report the
results where the runs are ranked using MAP. Please note that from the
previous experiments, considering it’s limitation discussed above, we already
know that increasing the number of topics does not improve the reusability of
a test collection when the pool depth is 100. The conclusion still holds at
each varying pool depth reported in Figure 2.

By observing Figure 2, we find that for a fixed number of topics, increasing
the pool depth improves the AUC in terms of $\tau_{ap}$ and lowers the value
of AUC, which indicates a better quality test collection. The greatest
improvement in AUC happens when we increase the pool depth from 20 to 40 in
all topic sets we investigate, suggesting that it has a high return on
investment. Based on our results, if we have enough evaluation budget, using a
pool depth of at least 40 seems a reasonable choice.

Another interesting observation from Figure 2 is that how the number of groups
interacts with the pool depth. For example, in the plot for 50 topics, when we
have a pool depth of 20 only, we need at least 7 groups (half of the total
number of groups in Robust’04 test collection), to achieve a $\tau_{ap}$
correlation of 0.9 or above. However, if the number of participants goes down
to 3 groups (one-fifth of the total number of groups in Robust’04 test
collection), we need a pool depth of $80$ to achieve the same $\tau_{ap}$
correlation. The observation holds for all other varying numbers of topics
sampled in Figure 2.

Based on the above discussion, we can conclude that the number of participants
is the most important factor for the quality of test collections. Therefore,
shared-task organizers should pay attention to the publicity of shared tasks.
However, if we have a few participating groups, for a given judging budget,
rather than increasing the number of topics, we should increase pool depth in
order to have a reusable test collection. It should be noted this conclusion
is subjected to all the limitations discussed in the previous experiment
except for the pool depth.

Figure 3:  $\tau_{ap}$, and Max Drop obtained by simulated test collections
with a varying number of groups with manual runs (top row) and without manual
runs (bottom row) on the five TREC datasets. The x-axis represents the number
of groups. Runs are ranked using MAP.

Figure 4:  The number of unique relevant documents obtained by simulated test
collections on six different test collections. Left plot does not include
Robust’04 (RB’04) dataset as it contains only automatic runs. In the right
plot, Robust’04 test collection has a higher number of unique relevant
documents because it has 249 topics whereas the other datasets have 50 topics
only.

####  3.2.3 Impact of the Document Collection Size.

We conduct experiments on five different test collections, namely Adhoc’99,
Adhoc’98, TB’06, WT’13, and WT’14, which contain both types of runs (i.e.,
manual and automatic runs), as seen in Table 1. We vary the number of groups
and report respective $\tau_{ap}$, Max Drop scores (See Figure 3), and the
number of relevant documents (See Figure 4), following Algorithm 1. We do not
report AUC since the AUC is not comparable across different collections. The
top row of Figure 3 considers both manual and automatic runs in the
construction of a test collection and the bottom row only considers the
automatic runs. We present results here using only MAP due to space
constraints.

From the Table 1, we can see that the size of document collection for Adhoc’99
and Adhoc’98 test collection is very small ($\approx 0.5$ Million) compared to
the size of document collection ($\approx 52$ Million) for WT’13 and WT’14
test collections. Although the number of unique relevant documents found
depends on the number of topics, pool depth, and the number of participating
groups, it is intuitive that for a fixed number of topics and pool depth, a
larger document collection usually has a higher number of unique relevant
documents. Therefore, the size of the document collection coupled with a
varying number of participating groups will affect the reusability of a test
collection. By observing Figure 3 (top row), we can see that for Adhoc’98 and
Adhoc’99 collections with both manual and automatic runs, we can achieve a
$\tau_{ap}$ of 0.9 or above when there are 8 and 9 participating groups,
respectively (i.e., approximately $20\%$ of the original number of
participating groups for the respective test collections). However, for WT’13
and WT’14 datasets, we need 12 and 8 groups (Figure 3, top row), respectively,
which is around on average $90\%$ of the original number of participating
groups for those respective test collections.

Although there is no acceptable range of values for Max Drop, if we consider
that the $90\%$ of the original number of groups participate for WT’13 and
WT’14 test collections, we find the Max Drop is in between 6 and 10. Regarding
Pearson correlation ($\rho$) scores computed between the number of groups and
the corresponding $\tau_{ap}$ and Max Drop, we observe higher Pearson
correlation values for WT’13 and WT’14 collections than other collections.
This is because $\tau_{ap}$ keeps increasing and Max Drop keeps decreasing
with an increasing number of groups in WT’13 and WT’14 test collections. On
the other hand, both Adhoc’99 and Adhoc’98 collections have low Pearson
correlation scores because $\tau_{ap}$ increases and Max Drop decreases when
the number of groups is increased from 2 to 10. However, after having 10
groups, $\tau_{ap}$ and Max Drop scores become almost stable for these two
collections.

In summary, we can conclude that we would need a higher number of
participating groups to have a reusable collection if the underlying document
collection is very large. This is probably due to the fact that runs in these
larger document collections would return more unique-relevant documents than
runs in other test collections (Figure 3). Therefore, a run might be highly
affected if it does not contribute to the pool. As another hypothesis, the
groups participating in the recent shared-tasks are able to develop a more
diverse set of IR runs than the groups participated in earlier shared-tasks
due to the progress in the field of IR.

However, it should be noted that the finding reported in this experiment has
certain limitations. Firstly, all five test collections used in this analysis
have only 50 topics. Secondly, for WT’13 and WT’14 collections, the employed
pool depth is very shallow (Table 1). A deeper pool depth and a higher number
of topics might provide us a different conclusion than the one stated here.

####  3.2.4 Impact of Manual Runs.

In order to see the impact of manual runs on test collection reusability, we
remove all manual runs from the simulated test collections and conduct the
same experiments as described above. Figure 3, bottom row presents the results
for this particular setup. Before going into further details, we should notice
from Figure 4 that the number of unique-relevant documents noticeably reduces
when we do not use any manual run in AdHoc’98, AdHoc’99, and TB’06 test
collections, confirming that manual runs usually provide more unique-relevant
documents than automatic runs. We also observe that WT’13 and WT’14 test
collections are less affected in terms of the number of unique relevant
documents by not having as many manual runs as other collections (See Table
1).

Comparing results between the top and bottom rows of Figure 3, we find that
not any using manual run increases the required number of participating groups
to achieve a $\tau_{ap}$ correlation of 0.9 or above. For example, for
Adhoc’98, and Adhoc’99 collections, we need at least $50\%$ of the original
number of participants (Table 3, bottom row) to achieve a $\tau_{ap}$
correlation of 0.9 or above which is $20\%$ (Table 3, top row) when we include
manual runs in the test collections. For TB’06 and WT’13 collections, we
actually need $100\%$ of the original number of participants. Our observations
are also similar for Max Drop scores. This suggests that the unique-relevant
documents detected by manual runs could not have been detected by any of the
automatic runs, thereby affecting the ranking of manual runs during the
evaluation.

The experimental evidence informs us that if we have a very few participating
groups, it is always better to have as much as possible manual runs to develop
a reusable test collection. Again, we should note that this finding is limited
by the facts that we have only 50 topics and the pool depth is very shallow
for the WT’13 and WT’14 test collections. Further experiments on other test
collections might suggest differently.

Table 2: Performance of the multiple linear regression model using Leave-one-test-collection-out setup on five different test collections. Training | Testing | Intersection of | MSE
---|---|---|---
Set | Set | Document Collection | $\tau_{ap}$
|  | between Train & Test Set |
{Adhoc’99 , TB’06 , WT’13 , WT’14} | Adhoc’98 | Yes | 0.005
{Adhoc’98 , TB’06 , WT’13 , WT’14} | Adhoc’99 | Yes | 0.003
{Adhoc’98 , Adhoc’99 , WT’13 , WT’14} | TB’06 | No | 0.171
{Adhoc’98 , Adhoc’99 , TB’06 , WT’14} | WT’13 | Yes | 0.060
{Adhoc’98 , Adhoc’99 , TB’06 , WT’13} | WT’14 | Yes | 0.055

##  4 Predicting the Qualities of Test Collections

We investigate whether it is possible to forecast the quality of a test
collection even before gathering the ranked lists of participants. Our
rationale is that the shared-task organizers can act accordingly based on the
predicted quality of a test collection even before spending budget on
collecting relevance judgments. In this study, we focus on predicting
$\tau_{AP}$ as a measure of reusability.

To generate data for our model, we use the same simulated test collections
constructed from Algorithm 1 and employ MAP to compute $\tau_{ap}$. We utilize
the following features for the prediction model: i) the number of
participating groups ($G$), ii) the number of topics ($T$), iii) the pool
depth ($P$) and iv) the size of the document collection ($C$). Then we fit a
Multiple Linear Regression model:
$\hat{y}=W_{0}+W_{1}*G+W_{2}*T+W_{3}*P+W_{3}*C$ on the training data to
predict $\tau_{ap}$. Here W’s are the learned weights for the features of the
model. As a performance measure of our prediction model, we report Mean
Squared Error (MSE)
$=\sum_{i=1}^{n}\frac{(\hat{y_{i}}-y_{i})*(\hat{y_{i}}-y_{i})}{n}$, where
$y_{i}$, and $\hat{y_{i}}$ are the predicted and true target value of
$\tau_{ap}$ respectively and $n$ is the total number of data points. A lower
value of MSE indicates a better model.

To understand the generalization performance of the prediction model, we
employ “Leave-one-test-collection-out” (LOTO) strategy. In this LOTO setup, in
turn, we hold out one test collection and utilize the remaining test
collections from the set of test collections to train the prediction model,
and then we test the predictive performance of the trained model on the held-
out test collection. For example, in the first row of Table 2, we utilize
{Adhoc’99, TB’06, WT’13 , WT’14} test collections as training set and Adhoc’98
test collection as testing set for our model. Since Adhoc’98 in the test set
shares the same document collection (Table 1) with Adhoc’99 from the train
set, the 3rd column of 1st row indicates “Yes”. In contrast, TB’06 in the
testing set (3rd row) does not share the document collection with any of the
test collections utilized in the training set and thus the 3rd column of the
3rd row indicates “No”. MSE for predicting $\tau_{ap}$ is reported in 4th
column of Table 2.

From Table 2, we observe that the model can predict $\tau_{ap}$ with a very
high accuracy (MSE $\leq 0.06$) for all five test collections except for
TB’06. This can be because a) runs in Adhoc’98 and Adhoc’99 might be similar
due to closeness of two shared-tasks in terms of time, and b) both have the
same document collection (Column 3 of Table 2). We can apply the same
reasoning for WT’13 and WT’14 test collections.

In summary, our prediction performance improves when we use test collections
from the same document collection for both training and testing sets of the
model. In practice, our prediction model can be especially useful after its
first year as TREC usually continues a track for more than one year. Thus, the
results in the first year can be used to forecast the quality of the test
collections in the following years for the same track. Test collection
construction parameters (i.e., topic set size and pool depth) can be set based
on the predictions.

##  5 Conclusion and Future Work

In this work, we investigate how varying the number of participating groups
coupled with other factors in a shared-task affects the reusability of a test
collection. Our main findings based on our extensive experiments conducted on
six TREC test collections are as follows. Firstly, when we have a very few
participating groups in a shared-task, increasing the pool depth provides a
more reusable test collection than increasing the number of topics. Secondly,
the size of the document collection and the types of runs play an important
role when the number of participants is very few. Thirdly, the reusability of
a test collection can be predicted with high accuracy when the test
collections used for training and testing the model have the same underlying
document collection. We believe that our findings will be useful for shared-
task organizers to take necessary actions while organizing the event.

There are many possible future directions for this work. For example, our
experimental analysis is conducted on pooling-based test collections. However,
we plan to extend our work on analyzing test collections constructed using
other techniques. For example, the TREC 2017 Common Core track (Allan et al.,
2017) test collection is constructed using a bandit-based technique (Losada et
al., 2016). In addition, we also plan to address the mentioned limitations of
our analysis by using more test collections, preferable recent ones in which
deep learning based IR models are used. Future work might also explore the
prediction of test collection quality when the ranked lists of documents are
submitted but relevance judgments are not collected yet. This scenario will
enable using different features such as the number of unique documents
retrieved by groups, rank correlation of ranked lists, and others with more
sophisticated machine learning models.

## References

  * (1)
  * Allan et al. (2017) James Allan, Donna Harman, Evangelos Kanoulas, Dan Li, Christophe Van Gysel, and Ellen M Voorhees. 2017.  TREC 2017 Common Core Track Overview.. In _TREC_.
  * Buckley and Voorhees (2000) Chris Buckley and Voorhees. 2000.  Evaluating Evaluation Measure Stability.  (2000), 33–40.  https://doi.org/10.1145/345508.345543
  * Büttcher et al. (2007) Stefan Büttcher, Charles LA Clarke, Peter CK Yeung, and Ian Soboroff. 2007.  Reliable information retrieval evaluation with incomplete and biased judgements. In _Proceedings of the 30th annual international ACM SIGIR conference on Research and development in information retrieval_. ACM, 63–70.
  * Carterette et al. (2008) Ben Carterette, Virgil Pavlu, Evangelos Kanoulas, Javed A. Aslam, and James Allan. 2008\.  Evaluation over thousands of queries. In _Proceedings of the 31st annual international ACM SIGIR conference on Research and development in information retrieval - SIGIR ’08_. ACM Press, Singapore, Singapore, 651.  https://doi.org/10.1145/1390334.1390445
  * Hawking et al. (1998) David Hawking, Nick Craswell, and Paul Thistlewaite. 1998\.  Overview of TREC-7 very large collection track.  _NIST SPECIAL PUBLICATION SP_ (1998), 93–106.
  * Hawking and Robertson (2003) David Hawking and Stephen Robertson. 2003.  On collection size and retrieval effectiveness.  _Information retrieval_ 6, 1 (2003), 99–105.
  * Kutlu et al. (2018) Mucahid Kutlu, Tamer Elsayed, and Matthew Lease. 2018\.  Intelligent topic selection for low-cost information retrieval evaluation: A New perspective on deep vs. shallow judging.  _Information Processing & Management_ 54, 1 (Jan. 2018), 37–59.  https://doi.org/10.1016/j.ipm.2017.09.002
  * Losada et al. (2016) David E. Losada, Javier Parapar, and Álvaro Barreiro. 2016\.  Feeling lucky?: multi-armed bandits for ordering judgements in pooling-based evaluation. In _Proceedings of the 31st Annual ACM Symposium on Applied Computing - SAC ’16_. ACM Press, Pisa, Italy, 1027–1034.  https://doi.org/10.1145/2851613.2851692
  * Sparck Jones and Van Rijsbergen (1975) K Sparck Jones and C Van Rijsbergen. 1975.  Report on the Need for and Provision of an” Ideal.  _Information Retrieval Test Collection_ (1975).
  * Sparck Jones and Van Rijsbergen (1976) Karen Sparck Jones and Cornelis Joost Van Rijsbergen. 1976\.  Information retrieval test collections.  _Journal of documentation_ 32, 1 (1976), 59–75.
  * Voorhees (2000) Voorhees. 2000.  Variations in relevance judgments and the measurement of retrieval effectiveness.  _Information processing & management_ 36, 5 (2000), 697–716.
  * Voorhees (2018) Voorhees. 2018.  On Building Fair and Reusable Test Collections Using Bandit Techniques. In _Proceedings of the 27th ACM International Conference on Information and Knowledge Management_ _(CIKM ’18)_. ACM, New York, NY, USA, 407–416.  https://doi.org/10.1145/3269206.3271766
  * Webber et al. (2008) William Webber, Alistair Moffat, and Justin Zobel. 2008\.  Statistical power in retrieval experimentation. In _Proceeding of the 17th ACM conference on Information and knowledge mining - CIKM ’08_. ACM Press, Napa Valley, California, USA, 571\.  https://doi.org/10.1145/1458082.1458158
  * Yilmaz et al. (2008) Emine Yilmaz, Javed A Aslam, and Stephen Robertson. 2008\.  A new rank correlation coefficient for information retrieval. In _Proceedings of the 31st annual international ACM SIGIR conference on Research and development in information retrieval_. ACM, 587–594.
  * Zobel (1998) Justin Zobel. 1998\.  How reliable are the results of large-scale information retrieval experiments?. In _Proceedings of the 21st annual international ACM SIGIR conference on Research and development in information retrieval_. ACM, 307–314.


[b]Carsten Schneider

#  DESY 22–117, DO-TH 22/18, RISC Report Series 22-08, SAGEX-22-28-E
Computer Algebra and Hypergeometric Structures for Feynman Integrals

Johannes Blümlein  Marco Saragnese

###### Abstract

We present recent computer algebra methods that support the calculations of
(multivariate) series solutions for (certain coupled systems of partial)
linear differential equations. The summand of the series solutions may be
built by hypergeometric products and more generally by indefinite nested sums
defined over such products. Special cases are hypergeometric structures such
as Appell-functions or generalizations of them that arise frequently when
dealing with parameter Feynman integrals.

##  1 Introduction

For certain classes of topologies Feynman parameter integrals in terms of the
dimensional parameter $\varepsilon=D-4$, with $D$ the dimension of space–time
can be represented in terms of higher transcendental functions of the
hypergeometric type, cf. e.g. [1, 2, 3, 4]. This includes, e.g., general
off–shell representations at the one–loop level for multi–leg diagrams, cf.
e.g. [5, 6, 7, 8, 9] and at the two- and three–loop level for various
scattering processes, cf. e.g. [12, 13, 10, 14]. Such representations of
higher transcendental function representations, covering as special cases
hypergeometric functions or Appell functions [11, 4, 16] and their variations,
are often preferential: they provide further inside in the classification of
the arising functions. In particular, one can expand the summands in
$\varepsilon$ and can apply symbolic methods [15] and in particular symbolic
summation tools, see [17, 16], to obtain the coefficients of special functions
in terms of iterative sums [18, 19, 20, 21, 22, 23, 15, 24, 25] or integrals
[26, 20, 21, 22, 23, 24, 25]

In most cases such representations cannot be extracted directly from the given
Feynman parameter integrals due to their complicated mathematical structure.
In such instances one often succeeds in computing linear difference equations
in a free discrete parameter (often the Mellin variable $n$) or linear
differential equation in a continuous variable $z$ (often in its $z$-space
representation). One successful tool is the Almkvist-Zeilberger algorithm [27,
28] to compute from a given integral (with an hyperexponential integrand) a
linear difference or differential equation. Another tactic is to apply
successively Newton’s binomial theorem and Mellin–Barnes decompositions on the
integrand, implemented in different packages [29]. Applying afterwards the
residue theorem produces definite multiple sums and one can utilize symbolic
summation tools [30, 17] in order to compute linear differential or difference
equations in the free parameters. Finally, a highly successful approach are
IBP methods [31] that enable one to represent such integrals in terms of
master integrals which themselves can be represented as solutions of coupled
systems of linear differential equations. Utilizing then uncoupling methods
[32] lead again to linear differential equations that contain the desired
parameter Feynman integrals as solutions. For further details on these and
related technologies we refer, e.g., to [33] and references therein.

In Section 2 we summarize the available tools to find solutions in terms of
iterative sums or integrals for a given linear difference or differential
equation. Utilizing these tools we will illustrate useful methods to find
hypergeometric series solutions and generalizations of them for ordinary
linear differential equations (Section 3.1), partial linear differential
equations (Section 3.2) and special coupled systems of partial linear
differential equations (Section 3.3). We conclude the article in Section 4.

##  2 Solving difference or differential equations in terms of iterative sums
or integrals

As described in the introduction one can use, e.g., symbolic integration and
summation tools or IBP methods to compute linear difference or differential
equations that contain the parameter Feynman integrals (or parts of them) as
solutions.

Given such a linear difference equation, one may use the summation package
Sigma [35, 17, 34] that enables one to solve the following problem.

GIVEN a recurrence

$a_{0}(n)F(n)+\dots+a_{\delta}(n)F(n+\delta)=h(n)$

where the coefficients $a_{0}(n),\dots,a_{\delta}(n)$ and the inhomogeneous
part $h(n)$ are expressions in terms of indefinite nested sums defined over
hypergeometric products, together with initial values, say111In the following
we assume that $\mathbb{K}$ is an appropriate computable field of
characteristic $0$ that contains all the relevant constants for our concrete
calculations. $F(0),\dots,F(\delta-1)\in\mathbb{K}$.
DECIDE constructively if $F(n)$ can be expressed in terms of indefinite nested
sums defined over hypergeometric products.

We note that this class of indefinite nested sums defined over hypergeometric
products contains as special cases all iterative sums that arose so far in QCD
calculations. In particular, it contains harmonic sums [18, 19] such as

$S_{2,1}(n)=\sum_{i=1}^{n}\frac{1}{i^{2}}\sum_{j=1}^{i}\frac{1}{j},$

generalized harmonic sums [20, 21] such as

$S_{1,1}(2,\tfrac{1}{2},n)=\sum_{k=1}^{n}\frac{2^{k}}{k}\sum_{i=1}^{k}\frac{\displaystyle(\tfrac{1}{2})^{i}}{i}\sum_{j=1}^{i}\frac{\displaystyle
S_{1}(j)}{j},$

cyclotomic harmonic sums [22] such as

$CS_{(2,1,2),(1,0,2),(2,1,1)}(n)=\sum_{k=1}^{n}{\frac{1}{(2k+1)^{2}}}\sum_{j=1}^{k}\frac{1}{j^{2}}\sum_{i=1}^{j}\frac{1}{2i+1},$

binomial sums [15, 23] such as

$\sum_{j=1}^{n}\frac{4^{j}S_{1}({j-1})}{\binom{2j}{j}j^{2}}\hskip 28.45274pt$

or generalized binomial sums [24, 25] such as

$\sum_{k=1}^{n}(\tfrac{1}{4})^{k}\big{(}1-{\eta}\big{)}^{k}\binom{2k}{k}\sum_{j=1}^{k}\frac{4^{j}}{j^{2}\binom{2j}{j}}$

where the extra parameter $\eta$ encodes the ratio of two masses. More
generally, sums of the form

$\sum_{k=1}^{n}\left(\prod_{i=1}^{k}\frac{1+i+i^{2}}{i+1}\right)\sum_{j=1}^{k}\frac{1}{j\binom{4j}{3j}^{2}}$

are covered where $f(k)=\prod_{i=1}^{k}\frac{1+i+i^{2}}{i+1}$ is
hypergeometric in $k$, i.e., its shift-quotient is a rational function in $k$:

$\frac{f(k+1)}{f(k)}=\frac{k^{2}+3k+3}{k+2}\in\mathbb{Q}(k).$

Similarly, one may activate the package HarmonicSums based on [37, 35, 38, 36]
to tackle the following problem.

GIVEN a linear differential equation

$b_{0}(x)f(x)+\dots+b_{\lambda}(x)D_{x}^{\lambda}f(x)=0$ (1)

with polynomials $b_{0}(x),\dots,b_{\lambda}(x)\in\mathbb{K}[x]$ and initial
values, say $f(0),\dots,D_{x}^{\lambda-1}f(x)|_{x=0}\in\mathbb{K}$.
DECIDE constructively if $f(x)$ can be expressed in terms of iterated
integrals defined over hyperexponential functions.

As in the difference equation case, this solver contains as special cases all
solutions in terms of iterated integrals that arose so far in QCQ-
calculations. Namely, it covers harmonic polylogarithms [26] such as

$H_{1,-1}(x)=\int_{0}^{x}\frac{1}{1-\tau_{1}}\int_{0}^{\tau_{1}}\frac{1}{1+\tau_{2}}d\tau_{2}d\tau_{1},$

generalize harmonic polylogarithms [20, 21] such as

$H_{2,-2}(x)=\int_{0}^{x}\frac{1}{2-\tau_{1}}\int_{0}^{\tau_{1}}\frac{1}{2+\tau_{2}}d\tau_{2}d\tau_{1},$

cyclotomic harmonic polylogarithms [22] such as

$\int_{0}^{x}\frac{1}{1+\tau_{1}+\tau_{1}^{2}}\int_{0}^{\tau_{1}}\frac{1}{1+\tau_{2}^{2}}d\tau_{2}d\tau_{1},$

radical integrals [23] such as

$\int_{0}^{x}\frac{1}{\sqrt{1+\tau_{1}}}\int_{0}^{\tau_{1}}\frac{1}{1+\tau_{2}}d\tau_{2}d\tau_{1},$

or generalized radical integrals [24, 25] such as

$\int_{0}^{x}\frac{1}{1-\tau_{1}+\eta\tau_{1}}\int_{0}^{\tau_{1}}\sqrt{1-\tau_{2}}\sqrt{1-\tau_{2}+\eta\tau_{2}}d\tau_{2}d\tau_{1}$

where $\eta$ is the ration of two masses. A more general example is

$\int_{0}^{x}e^{\int_{1}^{\tau_{1}}\frac{1}{1+y+y^{2}}dy}\int_{0}^{\tau_{1}}\frac{1}{1+\tau_{2}}d\tau_{2}d\tau_{1}$

where $f(\tau_{1})=e^{\int_{1}^{\tau_{1}}\frac{1}{1+y+y^{2}}dy}$ is
hyperexponential, i.e. its logarithmic derivative is a rational function in
$\tau_{1}$:

$\frac{D_{\tau_{1}}f(\tau_{1})}{f(\tau_{1})}=\frac{1}{1+\tau_{1}+\tau_{1}^{2}}\in\mathbb{Q}(\tau_{1}).$

We remark that HarmonicSums contains also Kovacic’s algorithm [39, 36] that
enables one to find Liouvillian solutions of second order linear differential
equations. As a consequence one can find also expressions in terms of
iterative integrals over double routed radicals such as
$\sqrt[3]{1-\sqrt{1+\tau_{1}}}$.

##  3 Extracting hypergeometric structures from differential equations

Given a linear differential equation, one may succeed in finding sufficiently
many linearly independent solutions in terms of iterative integrals over
hyperexponential functions by using the algorithms given in [37, 39, 35, 36].
But in most cases such a representation is not possible. If the underlying
parameter integral, say $f(x,\varepsilon)$, depends on the dimensional
parameter $\varepsilon$, one may use refined methods described in [28] (and
inspired by [40]) to find closed form representations of the
$\varepsilon$-expansion

$f(x,\varepsilon)=\phi_{l}(x)\varepsilon^{l}+\phi_{l+1}(x)\varepsilon^{l+1}+\dots+\phi_{r}(x)\varepsilon^{r}+O(\varepsilon^{r+1})$

for some $l,r\in\mathbb{Z}$ with $l\leq r$. More precisely, given the
differential equation and appropriate initial values, one can decide
algorithmically if the coefficients
$\phi_{l}(x),\phi_{l+1}(x),\dots,\phi_{r}(x)$ can be represented in terms of
iterative integrals over hyperexponential functions. But in the most general
case (in particular, for more complicated Feynman integrals), also here the
available methods will fail.

However, if one searches for representations in terms of hypergeometric
structures, the holonomic machinery [41] in combination with the available
recurrence solver of Sigma might lead to the desired result.

###  3.1 Ordinary linear differential equations

In the ordinary case, the holonomic approach can be stated as follows. Let

$f(x)=\sum_{n=0}^{\infty}F(n)x^{n}$ (2)

be a (formal) power series. Then there exist
$b_{0}(x),\dots,b_{\lambda}(x)\in\mathbb{K}[x]$ (not all zero) with

$b_{0}(x)f(x)+\dots+b_{\lambda}(x)D_{x}^{\lambda}f(x)=0$ (3)

if and only if there exist $a_{0}(x),\dots,a_{\delta}(x)\in\mathbb{K}[x]$ (not
all zero) with

$a_{0}(n)F(n)+\dots+a_{\delta}(n)F(n+\delta)=0.$ (4)

We emphasize that this correspondence is fully algorithmic. If the linear
difference equation or the linear differential equation is given, one can
compute the other equation explicitly. Here will exploit only one direction:
we are given a linear differential equation of an unknown power series (2) and
compute the corresponding linear difference equation for the coefficients
$F(n)$; this can be accomplished by plugging the ansatz (2) into the equation
and performing coefficient comparison w.r.t. $x^{n}$. If one succeeds in
finding a solution of the recurrence (3) with Sigma (see Section 2), one
obtains automatically a (formal) solution of the linear differential equation
(4).

Example 1: In order to find a power series solution

$f(x)=\sum_{n=0}^{\infty}F(n)x^{n}$

for

$-\left(x^{4}-64x^{3}\right)f^{(4)}(x)-2\left(5x^{3}-144x^{2}\right)f^{(3)}(x)\\\
-\left(25x^{2}-208x\right)f^{\prime\prime}(x)-(15x-8)f^{\prime}(x)-f(x)=0,$

we compute the recurrence

$8(n+1)(2n+1)^{3}F(n+1)-(n+1)^{4}F(n)=0$

for the coefficients $F(n)$ of the series $f(x)$. Using Sigma we obtain the
solution

$F(n)=\frac{1}{\binom{2n}{n}^{3}}=\frac{(1)_{n}(1)_{n}(1)_{n}(1)_{n}}{\big{(}\tfrac{1}{2}\big{)}_{n}\big{(}\tfrac{1}{2}\big{)}_{n}\big{(}\tfrac{1}{2}\big{)}_{n}n!}\frac{1}{64^{n}}.$

Thus the full set of power series solutions for the given linear differential
equation is given by

$f(x)=c\cdot\sum_{n=0}^{\infty}\frac{x^{n}}{\binom{2n}{n}^{3}}=c\cdot{}_{4}F_{3}\biggl{[}\genfrac{.}{.}{0.0pt}{}{1\mskip
8.0mu1\mskip 8.0mu1\mskip 8.0mu1}{\tfrac{1}{2}\mskip 8.0mu\tfrac{1}{2}\mskip
8.0mu\tfrac{1}{2}};\frac{x}{64}\biggr{]},\quad c\in\mathbb{K}.$

Example 2: To find power series solutions

$f(x)=\sum_{n=0}^{\infty}F(n)x^{n}$

for

$\left(x^{6}-32x^{5}+256x^{4}\right)f^{(6)}(x)+\left(23x^{5}-528x^{4}+2560x^{3}\right)f^{(5)}(x)\\\
+\left(171x^{4}-2552x^{3}+6272x^{2}\right)f^{(4)}(x)+2\left(245x^{3}-2002x^{2}+1728x\right)f^{(3)}(x)\\\
+2\left(253x^{2}-786x+72\right)f^{\prime\prime}(x)+4(35x-12)f^{\prime}(x)+4f(x)=0$

we compute the underlying recurrence

$(n+2)(n+1)^{3}F(n)-4(n+2)(2n+1)^{2}(2n+3)F(n+1)+16(2n+1)^{2}(2n+3)^{2}F(n+2)=0.$

Using Sigma we obtain the solutions

$F(n)=\frac{1}{\binom{2n}{n}^{2}}\Big{(}c_{1}+c_{2}S_{1}({n})\Big{)}=\frac{(1)_{n}(1)_{n}(1)_{n}}{\big{(}\tfrac{1}{2}\big{)}_{n}\big{(}\tfrac{1}{2}\big{)}_{n}n!}\frac{1}{16^{n}}\Big{(}c_{1}+c_{2}S_{1}({n})\Big{)},\quad
c_{1},c_{2}\in\mathbb{K};$

note that we obtained two linearly independent solutions and thus the solution
space of the recurrence is completely determined. Thus the full set of power
series solutions is given by

$f(x)=c_{1}\cdot{}_{3}F_{2}\biggl{[}\genfrac{.}{.}{0.0pt}{}{1\mskip
8.0mu1\mskip 8.0mu1}{\tfrac{1}{2}\mskip
8.0mu\tfrac{1}{2}};\frac{x}{16}\biggr{]}+c_{2}\sum_{n=0}^{\infty}\frac{S_{1}(n)}{\binom{2n}{n}^{2}}x^{n},\quad
c_{1},c_{2}\in\mathbb{K}.$

In the two examples above we obtained the full solution set of the underlying
linear recurrence and thus could provide all power series solutions of the
given linear differential equation. In general, utilizing the algorithms from
[35, 17, 34] we can find all power series solutions whose coefficients can be
given in terms of indefinite nested sums over hypergeometric products. This
toolbox has been utilized in various concrete QCD calculations such as [44].
In particular, we can find all ${}_{p}F_{q}$-solutions of the form
${}_{p}F_{q}[\dots,h(x)]$ with $h(x)=c\cdot x$ where $c\in\mathbb{K}^{*}$. We
note further that one can search also for solutions for more complicated
arguments, like $h(x)=\frac{x^{2}(x^{2}-9)^{2}}{(x^{2}+3)^{2}}$ as elaborated
in [43]. Furthermore, one can hunt also for Puiseux series solutions by
variants of this method.

###  3.2 Partial linear differential equations

For the multivariate case the problem can be stated as follows: Find a power
series solution

$f(x_{1},\dots,x_{r})=\sum_{n_{1}=0}^{\infty}\dots\sum_{n_{r}=0}^{\infty}F(n_{1},\dots,n_{r})x_{1}^{n_{1}}\dots
x_{r}^{n_{r}}$ (5)

for a partial linear differential equation

$\sum_{(s_{1},\dots,s_{r})\in
T}b_{(s_{1},\dots,s_{r})}(x_{1},\dots,x_{r})D^{s_{1}}_{x_{1}}\dots
D^{s_{r}}_{x_{r}}f(x_{1},\dots,x_{r})=0$

where the finite structure set $T\subset\mathbb{N}^{r}$ and the polynomial
coefficients
$b_{(s_{1},\dots,s_{r})}(x_{1},\dots,x_{r})\in\mathbb{K}[x_{1},\dots,x_{r}]$
are given.

Already the problem to decide if there is a polynomial solution
$f(x_{1},\dots,x_{r})\in\mathbb{K}[x_{1},\dots,x_{r}]$ is algorithmically
unsolvable. More precisely, it has been shown in [42] that a solution to this
problem is equivalent to solve Hilbert’s 10th problem (which is unsolvable).
Applying the holonomic translation mechanism [41], i.e., plugging in the
ansatz (5) into the equation and performing coefficient comparison w.r.t.
$x_{1}^{n_{1}}\dots x_{r}^{n_{r}}$, one gets a partial linear difference
equation of the form

$\sum_{(s_{1},\dots,s_{r})\in
S}a_{(s_{1},\dots,s_{r})}(n_{1},\dots,n_{r})F(n_{1}+s_{1},\dots,n_{r}+s_{r})=0$

with a finite structure set $S\subset\mathbb{Z}^{r}$ and polynomial
coefficients
$a_{(s_{1},\dots,s_{r})}(n_{1},\dots,n_{r})\in\mathbb{K}[n_{1},\dots,n_{r}]$.
However, as for the differential case, it has been shown in [42] that in
general one cannot decide algorithmically if there exists a polynomial
solution $F(n_{1},\dots,n_{r})\in\mathbb{K}[n_{1},\dots,n_{r}]$. But in the
discrete case there are at least methods available [45] that enable one to
hunt for rational solutions. Recently, these techniques have been extended in
[16] to search for solutions in terms of iterative sums defined over
hypergeometric products provided that one can predict the set of sums that
arises.

For instance, given the partial linear difference equation

$(n+1)^{2}\left(k+n^{2}+2\right)\left(3kn^{2}-4k^{2}-5kn-12k+2n^{3}+2n^{2}-8n-8\right){\color[rgb]{0,0,1}F(n,k+1)}\\\
+(n+1)^{2}\left(k+n^{2}+3\right)\left(2k^{2}-2kn^{2}+2kn+6k-n^{3}-n^{2}+4n+4\right){\color[rgb]{0,0,1}F(n,k+2)}\\\
+(n+1)^{2}(k+n+1)\left(2k-n^{2}+n+4\right)\left(k+n^{2}+1\right){\color[rgb]{0,0,1}F(n,k)}\\\
-(k+1)n^{2}(n+2)^{2}\left(k+n^{2}+2n+2\right){\color[rgb]{0,0,1}F(n+1,k)}\\\
+kn^{2}(n+2)^{2}\left(k+n^{2}+2n+3\right){\color[rgb]{0,0,1}F(n+1,k+1)}=0$

and given the predicted set of sums $W=\\{S_{1}(k),S_{1}(n+k),S_{2,1}(n+k)\\}$
together with a degree bound $b=5$ for the maximal total degree of the arising
objects in the numerator, one can compute $37$ solutions

$\frac{p}{(1+n)^{2}(1+k+n^{2})}$

with

$\displaystyle p\in\Big{\\{}$ $\displaystyle
1+\frac{1}{2}nS_{1}({k+n}),k,n,kn,kn^{2},kn^{3},kn^{4},kS_{1}({n}),knS_{1}({n}),kn^{2}S_{1}({n}),kn^{3}S_{1}({n}),kS_{1}({n})^{2},$
$\displaystyle
knS_{1}({n})^{2},kn^{2}S_{1}({n})^{2},knS_{1}({n})^{3},kS_{1}({n})^{4},kS_{2,1}({n}),knS_{2,1}({n}),kn^{2}S_{2,1}({n}),kn^{3}S_{2,1}({n}),$
$\displaystyle
kS_{1}({n})S_{2,1}({n}),knS_{1}({n})S_{2,1}({n}),kn^{2}S_{1}({n})S_{2,1}({n}),kS_{1}({n})^{2}S_{2,1}({n}),knS_{1}({n})^{2}S_{2,1}({n}),$
$\displaystyle
kS_{1}({n})^{3}S_{2,1}({n}),kS_{2,1}({n})^{2},knS_{2,1}({n})^{2},kn^{2}S_{2,1}({n})^{2},kS_{1}({n})S_{2,1}({n})^{2},knS_{1}({n})S_{2,1}({n})^{2},$
$\displaystyle
kS_{1}({n})^{3},kS_{1}({n})^{2}S_{2,1}({n})^{2},kS_{2,1}({n})^{3},knS_{2,1}({n})^{3},kS_{1}({n})S_{2,1}({n})^{3},kS_{2,1}({n})^{4}\Big{\\}}$

by using the package SolvePLDE [16]. In particular, one can search for
solutions with hypergeometric contributions. In particular, the techniques
from [40] have been incorporated in this package to search for the
coefficients of $\varepsilon$-expansions.

###  3.3 Coupled systems of partial linear differential equations

Solving coupled system of partial linear differential (and difference)
equations is a widely open research topic. However, as elaborated in [16] one
succeeds in finding hypergeometric structures for at least some interesting
special cases. For instance, given the system

$\displaystyle(x-1)y{\color[rgb]{0,0,1}D_{xy}f(x,y)}+(x(2\varepsilon+\tfrac{7}{2})-\varepsilon+1){\color[rgb]{0,0,1}D_{x}f(x,y)}$
$\displaystyle+(x-1)x{\color[rgb]{0,0,1}D_{x}^{2}f(x,y)}+y(2\varepsilon+1){\color[rgb]{0,0,1}D_{y}f(x,y)}+\tfrac{3}{2}(2\varepsilon+1){\color[rgb]{0,0,1}f(x,y)}$
$\displaystyle=0,$ $\displaystyle
x(y-1){\color[rgb]{0,0,1}D_{xy}f(x,y)}+x(4-\varepsilon){\color[rgb]{0,0,1}D_{x}f(x,y)}+(y-1)y{\color[rgb]{0,0,1}D_{y}^{2}f(x,y)}$
$\displaystyle+(y(\tfrac{13}{2}-\varepsilon)-\varepsilon+1){\color[rgb]{0,0,1}D_{y}f(x,y)}+\tfrac{3(4-\varepsilon)}{2}{\color[rgb]{0,0,1}f(x,y)}$
$\displaystyle=0,$

we can find a power series solution

$f(x,y)=\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}F(n,m)x^{n}y^{m}$

as follows. Plugging in this ansatz and comparing coefficients w.r.t.
$x^{n}y^{n}$ yields the first-order coupled system

$\displaystyle\tfrac{3}{2}(2\varepsilon+1)F(n,m)-n(\varepsilon-1)F(n+1,m)$
$\displaystyle=0,$
$\displaystyle-\tfrac{3}{2}(\varepsilon-4)F(n,m)-m(\varepsilon-1)F(n,m+1)$
$\displaystyle=0.$

In general, the Ore–Sato theorem [46] provides a criterion for such first-
order systems whether a solution can be given in terms of Gamma-functions and
Pochhammer symbols. In [16] a surprisingly simple method has been elaborated
and implemented in the package HypSeries that enables one to produce such
solutions in terms of hypergeometric products. In our concrete example we
compute the solution

$\displaystyle F(n,m)=$
$\displaystyle\Big{(}\prod_{i=1}^{n}\frac{(1+2i)(3+i-\varepsilon)}{2i(-2+i+\varepsilon)}\Big{)}\prod_{i=1}^{m}\frac{(1+2i+2n)(i+2\varepsilon)}{2i(-2+i+n+\varepsilon)}=\frac{\big{(}\frac{3}{2}\big{)}_{m+n}(4-\varepsilon)_{n}(1+2\varepsilon)_{m}}{m!n!(-1+\varepsilon)_{m+n}}$

and thus the hypergeometric series solution

$\displaystyle
f(x,y)=\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\frac{\big{(}\frac{3}{2}\big{)}_{m+n}(4-\varepsilon)_{n}(1+2\varepsilon)_{m}}{m!n!(-1+\varepsilon)_{m+n}}$

which is closely related to Appell-like structures. Expanding the summand in
$\varepsilon$ gives

$f(x,y)=\varepsilon^{-1}\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}F_{-1}(n,m)+\varepsilon^{0}\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}F_{0}(n,m)+\dots$

with

$\displaystyle F_{-1}(n,m)$
$\displaystyle=-\frac{1}{6}\frac{x^{m}y^{n}(3+n)!\big{(}\frac{3}{2}\big{)}_{m+n}}{n!(-2+m+n)!}$
$\displaystyle F_{0}(n,m)$ $\displaystyle=\bigg{[}\dots
6S_{1}({n})+6S_{1}({m+n})-12S_{1}({m})\bigg{]}\frac{x^{m}y^{n}(3+n)!\big{(}\frac{3}{2}\big{)}_{m+n}}{n!(-2+m+n)!}.$

Finally, using the summation package Sigma with its summation algorithms [17]
one obtains

$\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}F_{-1}(n,m)=-\frac{15x^{6}}{4(x-y)^{4}(1-x)^{7/2}}-\frac{15y^{3}Q(x,y)}{64(x-y)^{4}(1-y)^{13/2}}$

for some polynomial $Q(x,y)$. For further details, in particular for the
result of the constant coefficient
$\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}F_{0}(n,m)$ we refer to [16].

##  4 Conclusion

We illustrated up-to-date computer algebra methods and implementations that
assist the user to find hypergeometric structures for various classes of
linear differential equations. While the ordinary case (see Sections 2 and
3.1) has been pushed forward non-trivially in various directions, the partial
case (see Sections 3.2 and 3.3) has been neglected so far. We hope that the
recent results from [16] and the supplementary material of this article will
open up new developments for improved differential equation solvers that are
instrumental for future explorations and calculuations of parameter Feynman
integrals.

Acknowledgment. This work was supported by the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska–Curie grant
agreement No. 764850, SAGEX and the Austrian Science Fund (FWF) grants SFB F50
(F5009-N15) and P33530.

## References

  * [1] R. Hamberg, Second order gluonic contributions to physical quantities, Ph.D. Thesis, Leiden University, 1991.
  * [2] A.I. Davydychev and M.Y. Kalmykov, Nucl. Phys. B 699 (2004) 3–64 [arXiv:hep-th/0303162 [hep-th]].
  * [3] I. Bierenbaum, J. Blümlein and S. Klein, Nucl. Phys. B 780 (2007) 40–75 [arXiv:hep-ph/0703285 [hep-ph]].
  * [4] M. Kalmykov, V. Bytev, B.A. Kniehl, S.O. Moch, B.F.L. Ward and S.A. Yost, Hypergeometric Functions and Feynman Diagrams, [arXiv:2012.14492 [hep-th]], in: Antidifferentiation and the Calculation of Feynman Amplitudes, (Springer, Heidelberg, 2021), J. Blümlein and C. Schneider, eds.
  * [5] E.E. Boos and A.I. Davydychev, Theor. Math. Phys. 89 (1991) 1052–1063.
  * [6] J. Fleischer, F. Jegerlehner and O.V. Tarasov, Nucl. Phys. B 672 (2003) 303–328 [arXiv:hep-ph/0307113 [hep-ph]].
  * [7] N. Watanabe and T. Kaneko, J. Phys. Conf. Ser. 523 (2014) 012063 [arXiv:1309.3118 [hep-ph]].
  * [8] J. Blümlein, K.H. Phan and T. Riemann, Acta Phys. Polon. B 48 (2017) 2313 [arXiv:1711.05510 [hep-ph]].
  * [9] K.H. Phan and T. Riemann, Phys. Lett. B 791 (2019) 257–264 [arXiv:1812.10975 [hep-ph]].
  * [10] S. Bauberger, M. Böhm, G. Weiglein, F.A. Berends and M. Buza, Nucl. Phys. B Proc. Suppl. 37 (1994) no.2, 95–114 [arXiv:hep-ph/9406404 [hep-ph]].
  * [11] P. Appell and J. Kampé de Fériet, Fonctions Hypergéométriques et Hypersphériques, Polynomes D’ Hermite, (Gauthier-Villars, Paris, 1926);
P. Appell, Les Fonctions Hypergëométriques de Plusieur Variables, (Gauthier-
Villars, Paris, 1925);
M.J. Schlosser, in: Computer Algebra in Quantum Field Theory: Integration,
Summation and Special Functions, C. Schneider, J. Blümlein, Eds., p. 305–324,
(Springer, Wien, 2013) [arXiv:1305.1966]

  * [12] C. Anastasiou, E.W.N. Glover, and C. Oleari, Nucl. Phys. B572 (2000) 307–360 [hep-ph/9907494].
  * [13] C. Anastasiou, E.W.N. Glover, and C. Oleari, Nucl. Phys. B565 (2000) 445–467, [hep-ph/9907523].
  * [14] J. Ablinger, J. Blümlein, A. Hasselhuhn, S. Klein, C. Schneider and F. Wißbrock, Nucl. Phys. B 864 (2012) 52–84 [arXiv:1206.2252 [hep-ph]].
  * [15] A.I. Davydychev and M. Yu Kalmykov, Nucl. Phys. B 699 (2004) 3–64;
S. Weinzierl, J. Math. Phys. 45 (2004) 2656–2673

  * [16] J. Blümlein, M. Saragnese, and C. Schneider, Hypergeometric Structures in Feynman Integrals [arXiv:2111.15501]
  * [17] C. Schneider, Sém. Lothar. Combin. 56 (2007) 1–36 article B56b;
C. Schneider, Simplifying Multiple Sums in Difference Fields, in: Computer
Algebra in Quantum Field Theory: Integration, Summation and Special Functions
Texts and Monographs in Symbolic Computation eds. C. Schneider and J. Blümlein
(Springer, Wien, 2013) 325–360 [arXiv:1304.4134]

  * [18] J.A.M. Vermaseren, Int. J. Mod. Phys. A 14 (1999) 2037–2076 [hep-ph/9806280].
  * [19] J. Blümlein and S. Kurth, Phys. Rev. D 60 (1999) 014018 [hep-ph/9810241].
  * [20] S. Moch, P. Uwer and S. Weinzierl, J. Math. Phys. 6 (2002) 3363–3386
  * [21] J. Ablinger, J. Blümlein and C. Schneider, J. Math. Phys. 54 (2013) 082301 [arXiv:1302.0378]
  * [22] J. Ablinger, J. Blümlein and C. Schneider, J. Math. Phys. 52 (2011) 102301 [arXiv:1105.6063]
  * [23] J. Ablinger, J. Blümlein, C.G.Raab and C. Schneider, J. Math. Phys. 55 (2014) 112301 [arXiv:1407.1822]
  * [24] J. Ablinger, J. Blümlein, A. De Freitas, A. Goedicke, C. Schneider and K. Schönwald, Nucl. Phys. B 932 (2018) 129–240 [arXiv:1804.02226]
  * [25] J. Ablinger, J. Blümlein, A. De Freitas, A. Goedicke, M. Saragnese, C. Schneider and K Schönwald, Nucl. Phys. B 955 (2020) 115059 [arXiv:2004.08916]
  * [26] E. Remiddi and J.A.M. Vermaseren, Int. J. Mod. Phys. A 15 (2000) 725–754 [hep-ph/9905237].
  * [27] M. Apagodu and D. Zeilberger, 2006 Adv. Appl. Math. 37 139–152;
J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C.
Schneider, Comput. Phys. Commun. 202 (2016) 33–112 [arXiv:1509.08324]

  * [28] J. Ablinger, 2021 Extensions of the AZ-algorithm and the Package MultiIntegrate Anti-Differentiation and the Calculation of Feynman Amplitudes ed Blümlein J and Schneider C (Cham: Springer International Publishing) pp 35–61
  * [29] M. Czakon, Comput. Phys. Commun. 175 (2006) 559–571;
A. Smirnov and V. Smirnov, Eur. Phys. J. C62 (2009) 445–449;
J. Gluza, K. Kajda and T. Riemann Comput. Phys. Commun. 177 (2007) 879–893;
B. Ananthanarayan, S. Banik, S. Friot and S. Ghosh, Phys. Rev. Lett. 127
(2021) 151601 [arXiv:hep-th/2012.15108]

  * [30] H. Wilf and D. Zeilberger, Invent. Math. 108 (1992) 575–633;
K. Wegschaider, Computer generated proofs of binomial multi-sum identities
(1997) Master’s thesis RISC, J. Kepler University;
R. Lyons, P. Paule and A. Riese, Appl. Algebra Engrg. Comm. Comput. 13 (2002)
327–333

  * [31] K.G. Chetyrkin and F.V. Tkachov, Nucl. Phys. B 192 (1981) 159–204;
S. Laporta, Int. J. Mod. Phys. A 15 (2000) 5087–5159 [hep-ph/0102033];
P. Marquard and D. Seidel, The Crusher algorithm, unpublished;
C. Studerus, Comput. Phys. Commun. 181 (2010) 1293–1300 [arXiv:0912.2546
[physics.comp-ph]];

  * [32] B. Zürcher, Rationale Normalformen von pseudo-linearen Abbildungen (1994) Ph.D. thesis Mathematik, ETH Zürich;
A. Bostan, F. Chyzak and É. de Panafieu, Complexity estimates for two
uncoupling algorithms, in: Proc. ISSAC’13 (2013, Boston);
S. Gerhold, Uncoupling Systems of Linear Ore Operator Equations (2002)
Master’s thesis RISC, J. Kepler University Linz.

  * [33] J. Blümlein and C. Schneider, The SAGEX Review on Scattering Amplitudes, Chapter 4: Multi-loop Feynman Integrals, J. Phys. A: Math. Theor., in press, 2022 [arXiv:2203.13015 [hep-th]
  * [34] C. Schneider, Symbolic Summation in Difference Fields, Ph.D. Thesis, RISC, Johannes Kepler University Linz, 2001;
C. Schneider, J. Differ. Equations Appl. 11 (2005) 799–821;
C. Schneider, J. Symb. Comput. 72 (2016) 82–127; S. Abramov, M. Bronstein, M.
Petkovšek and C. Schneider, J. Symb. Comput. 107 (2021) 23–66;
J. Ablinger and C. Schneider, Solving linear difference equations with
coefficients in rings with idempotent representations, in: Proc. ISSAC 21, ed
Mezzarobba M. (2021) 27–34.

  * [35] S. Abramov and M. Petkovšek, 1994 D’Alembertian solutions of linear differential and difference equations, in: Proc. ISSAC’94 ed von zur Gathen J (ACM Press) pp 169–174
  * [36] J. Ablinger, Computing the inverse Mellin transform of holonomic sequences using Kovacic’s algorithm PoS (RADCOR2017) 069 ed Hoang A and Schneider C (2017) pp 1–8
  * [37] M. Bronstein, Linear ordinary differential equations: Breaking through the order 2 barrier, in: Proc. ISSAC ’92 (1992) 42–48
  * [38] J. Blümlein, Comput. Phys. Commun. 159 (2004) 19–54 [arXiv:hep-ph/0311046].
J. Ablinger and C. Schneider, Annals of Combinatorics 22 (2018) 213–244
[arXiv:1510.03692].

  * [39] J.J. Kovacic, J. Symb. Computation 2 (1986) 3–43
  * [40] J. Blümlein, S. Klein, C. Schneider and F. Stan J. Symb. Comput. 47 2012 1267–1289 [arXiv:1011.2656]
  * [41] B. Salvy and P. Zimmermann, ACM Trans. Math. Software 20 (1994) 163–177;
C. Mallinger, Algorithmic Manipulations and Transformations of Univariate
Holonomic Functions and Sequences Master’s thesis RISC, J. Kepler University
Linz, 1996;
M. Kauers and P. Paule, The Concrete Tetrahedron, Text and Monographs in
Symbolic Computation (Springer Wien), 2011.

  * [42] S.A. Abramov and M. Petkovšek, On Polynomial Solutions of Linear Partial Differential and (q-)Difference Equations, in: Proc. CASC (2012) 1–11.
  * [43] E. Imamoglu and M. van Hoeij, J. Symbolic Comput. 83 (2017) 245–271;
J. Ablinger, J. Blümlein, A. De Freitas, M. van Hoeij, E. Imamoglu, C.G. Raab,
C.S. Radu, and C. Schneider J. Math. Phys. 59 (2018) 062305 [arXiv:1706.01299]

  * [44] J. Blümlein,, A. Hasselhuhn, S. Klein, and C. Schneider, Nucl. Phys. B 866 (2013) 196–211 [arXiv:1205.4184];
J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C.
Schneider, Nucl. Phys. B 890 (2014) 48–151 [arXiv:409.1135];
J. Ablinger, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel, M.
Round, C. Schneider and F. Wißbrock, Nucl. Phys. B 882 (2014) 263–288
[arXiv:1402.0359];
A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C. Schneider,
Nucl. Phys. B 897 (2015) 612–644 [arXiv:1504.08217];
A. Behring, J. Blümlein, A. De Freitas, A. Hasselhuhn, A. von Manteuffel and
C. Schneider, Phys. Rev. D 92 (2015) 114005 [arXiv:1508.01449];
J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C.
Schneider, Comput. Phys. Commun. 202 (2016) 33–112 [arXiv:1509.08324];
A. Behring, J. Blümlein, G. Falcioni, A. De Freitas, A. von Manteuffe and C.
Schneider, Phys. Rev. D 94 (2016) 114006 [arXiv:1609.06255];
J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel and C.
Schneider, Nucl. Phys. B 922 (2017) 1–40 [arXiv:1705.01508]

  * [45] M. Kauers and C. Schneider, Partial denominator bounds for partial linear difference equations, in: Proc. ISSAC’10 (2010) 211–218;
M. Kauers and C. Schneider, A refined denominator bounding algorithm for
multivariate linear difference equations, in: Proc. ISSAC’11 (2011) 201–208

  * [46] S.A. Abramov and M. Petkovšek, Adv. in Appl. Math. 29 (2002) 386–411.


# Ultrasonic wave transport in concentrated disordered resonant emulsions

Benoit Tallon  Department of Physics and Astronomy, University of Manitoba,
Winnipeg, Manitoba R3T 2N2, Canada  Thomas Brunet  I2M, Université de Bordeaux
- CNRS - Bordeaux INP, Talence, France F-33405  John H. Page
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Manitoba, Winnipeg, Manitoba Canada R3T 2N2, Canada

###### Abstract

We show how resonant (near-field) coupling affects wave transport in
disordered media through ultrasonic experiments in concentrated suspensions.
The samples consist of resonant emulsions in which oil droplets are suspended
in a liquid gel. By varying the droplet concentration, the limits of the
Independent Scattering Approximation are experimentally demonstrated. For the
most concentrated samples, the proximity of resonant scatterers induces a
renormalization of the surrounding medium, leading to a reducing of scattering
strength. We point out an optimal volume fraction $\phi$ of oil droplets for
which non-diffusive wave transport is experimentally demonstrated. Our
demonstration of maximum scattering at an intermediate droplet concentration
is very relevant for designing materials for the study of wave transport
phenomena such as Anderson Localization.

Anderson localization of classical waves in three dimensions has been
unambiguously demonstrated more than ten years ago Hu  _et al._ (2008). Such
an unusual effect (the halt of wave propagation due to wave interferences)
only occurs in strongly scattering heterogeneous media. When the scattering
mean free path $\ell_{s}$ (the wave attenuation characteristic length due to
scattering events) is larger than the incident wavelength $\lambda$,
interferences are negligible and descriptions based on the independent
scattering and ladder approximations are valid to describe the wave intensity
transport Sheng (2006). Under these approximations, the average intensity
transport is well described by the diffusion equation Page  _et al._ (1995).
However, when $\ell_{s}\sim\lambda$, mesoscopic wave interferences are
significant leading to “sub-diffusive” transport or localized signatures such
as non-Rayleigh speckle statistics Chabanov  _et al._ (2000); Hu  _et al._
(2008), wavefunction multifractality Faez  _et al._ (2009), infinite range
correlations Hildebrand  _et al._ (2014) and spatial confinement of intensity
Hu  _et al._ (2008); Cobus  _et al._ (2018), which is also revealed in the
dynamics of the coherent backscattering cone Cobus  _et al._ (2016).
Observations of such mesoscopic phenomena have led to much activity designing
new strongly-scattering samples for both optics and acoustics experiments
Lagendijk  _et al._ (2009); Skipetrov, S. E. and Page, J. H. (2016); Tallon
_et al._ (2017); Haberko  _et al._ (2020); Goïcoechea  _et al._ (2020); Vynck
_et al._ (2021); Monsarrat  _et al._ (2021); Matis  _et al._ (2022). In this
context, one way to decrease the scattering mean free path is to take
advantage of scattering resonances. In acoustics, the mechanical contrast
between heterogeneities and the surrounding medium leads to shape resonances
Brunet  _et al._ (2012) (analogous to optical Mie resonances Mie (1908)),
characterized by large deformations of the scattering inclusions (particles,
or droplets in the case of two-fluid media). The wave intensity generated by
these deformations is characterized by the cross section $\sigma_{s}$, which
is maximized around resonant frequencies. To link the scattering properties of
an isolated particle to the global scattering strength of the medium, one can
invoke the Independent Scattering Approximation (ISA) Lax (1951). The basic
idea of the ISA is to assume that in a diluted medium, a particle will scatter
the wave at most once in a scattering sequence (recurrent scattering events
such as loops are neglected). Under this assumption, the scattering can be
expressed easily in terms of the cross section
$\sigma_{s}\propto\textrm{Im}\\{f(0)\\}$, with $f(0)$ being the forward
scattering amplitude for a single scattering event. Thus, the scattering mean
free path can be written as $\ell_{s}=1/\eta\sigma_{s}$, with $\eta$ being the
number of scatterers per unit volume.

In this paper, we show that an increase of the scatterers’ concentration
$\eta$ is not necessarily a condition that will lead to a decrease in
$\ell_{s}$. By probing both coherent and diffusing ultrasonic wave transport
in resonant emulsions, we observe that coupling between near-by resonant
scatterers may result in a weakening of the scattering strength of the medium.
Our all-fluid suspension of fluorinated (FC40) oil droplets randomly dispersed
in a water-based gel matrix constitutes an excellent model system for this
study. In diluted emulsions (volume fraction $\phi=$ 5%), we have already
observed Tallon  _et al._ (2017) a strong resonant scattering regime due to
high sound speed contrast between oil droplets ($v_{1}\equiv
v_{\textrm{oil}}=0.64$ mm/$\mu$s) and the surrounding gel ($v_{0}\equiv
v_{\textrm{gel}}=1.48$ mm/$\mu$s). By increasing the concentration, we observe
initially that the scattering strength increases but at a slower rate than
predicted by the ISA. As the concentration is further increased, we
demonstrate here that an optimal concentration exists for maximizing the
scattering strength. At this point, the scattering is sufficiently strong that
the transport becomes subdiffusive and is well described by the self-
consistent theory of localization. Beyond this concentration, however, the
scattering becomes weaker, and simple diffusive transport is again found, with
the diffusion coefficient at $\phi=$ 40% becoming comparable to its value at
$\phi=$ 5%. Since the velocity contrast in this emulsion ($v_{1}<<v_{0}$) is
similar to that found in some strongly scattering optical systems, our
observation of maximum scattering at a concentration that is less than half
that of random close packing may be relevant to one of the challenges that has
been encountered in experimental observations of the localization of light
Skipetrov, S. E. and Page, J. H. (2016).

The samples studied here are relatively concentrated all-fluid suspensions (5%
$<\phi\leq$ 40 %) with a droplet mean radius $\bar{a}=$ 0.17 mm. This droplet
size allows the first set of shape resonances to be observed in the MHz range,
as these resonances occur when the wavelength is comparable to the droplet
size. Very low droplet polydispersity ($\sim 3$ %), crucial for studying the
influence of resonances on wave propagation Mascaro  _et al._ (2013), was
ensured thanks to micro-fluidic techniques controlled by robotics Tallon  _et
al._ (2020). Another advantage of this emulsion system is its almost
negligible acoustic dissipation, due to the low viscosity of both FC40 oil and
water-based gel. Thus the effects of resonances on wave transport were not
obscured by dissipation effects, allowing diffusive wave transport in pulsed
experiments to be observed over a long range of times and consequences of
mesoscopic interferences to be probed.

Ultrasonic characterization of the samples’ basic acoustic properties was
performed via ballistic measurements of ensemble-averaged transmitted acoustic
wave pulses, as described in detail in Refs. Tallon  _et al._ (2017, 2020).
Hence, accurate measurements of both the amplitude attenuation coefficient
$\alpha$ and the phase velocity $v_{\textrm{ph}}$ were performed, quantifying
the impact of droplet shape resonances on wave propagation as illustrated in
Fig. 1 for $\phi=25\%$. Note that because dissipation is negligibly small,
$\alpha$ provides a measure of the scattering mean free path
($\ell_{s}=\frac{1}{2}\alpha^{-1}$ to an excellent approximation Tallon  _et
al._ (2017, 2020)).

Comparison of our experimental data with theory in Fig. 1 reveals that
predictions based on both the independent scattering approximation (ISA) Lax
(1951) and the Lloyd and Berry model (L&B) Lloyd and Berry (1967) do not
describe well the average wave field in these concentrated resonant emulsions.
In these theories, the wave number $k$ of the ballistic wave can be written
as:

$k^{2}=(k^{\prime}+\textrm{j}k^{\prime\prime})^{2}=k_{0}^{2}+\eta\delta_{k_{0}},$
(1)

where $k^{\prime}=\omega/v_{\textrm{ph}}$, $k^{\prime\prime}=\alpha$ and
$\omega$ is the angular frequency of the incident plane wave. Here,
$\delta_{k_{0}}$ is the modification of the wave number in the pure matrix
$k_{0}$ due to scattering and can be expressed as a function of the amplitude
of a wave scattered by an isolated droplet $f(\theta)$:

$\begin{split}\delta_{k_{0}}=4\pi f(0)&\qquad\textrm{(ISA)}\\\
\delta_{k_{0}}=4\pi
f(0)&-\eta\frac{4\pi^{2}}{k_{0}^{2}}\bigg{[}f^{2}(0)-f^{2}(\pi)\\\
+&\int_{0}^{\pi}\frac{\textrm{d}f^{2}(\theta)}{\textrm{d}\theta}\frac{\textrm{d}\theta}{\textrm{sin}(\theta/2)}\bigg{]}\qquad\textrm{(L\&B)},\end{split}$
(2)

where $\theta$ is the angle between the directions of the incident and
scattered waves. Thus, $f(0)$ and $f(\pi)$ represent the amplitudes of forward
and backward scattered waves respectively. Uncorrelated point scatterers are
assumed in the ISA model, making this theory relevant only for dilute media
($\eta\ll 1$). When $\eta$ increases, the finite size, $a$, of scatterers
induces spatial correlations (the “hole correction” Fikioris and Waterman
(1964)) included in L&B theory.

For a concentrated heterogeneous medium in the intermediate frequency regime
($\lambda\sim a$), waves incident on a given scatterer may include substantial
contributions from the waves scattered by near-by scatterers. In this case,
the surrounding matrix may be viewed as an effective medium that depends on
scatterers’ properties Sheng (2006) and the correction $\delta$ depends on the
wave number $k$. Hence, the implicit equation $k^{2}=k_{0}^{2}+\eta\delta_{k}$
needs to be solved. A theory using these concepts is the spectral function
approach (SpFA) Jing  _et al._ (1992); Sheng (2006). As shown in Fig. 1,
better agreement is found between our experiments and SpFA predictions than
between these experiments and the ISA or L&B theories.

Figure 1:  Attenuation (a) and phase velocity (b) measurements (open circles),
for $\phi=25\%$, compared to predictions from Lloyd and Berry theory (L&B,
black dash-dotted line), independent scattering approximation (ISA, red dashed
line) and spectral function approach (SpFA, black solid line).

The scattering concept that is incorporated in the SpFA is illustrated in Fig.
2(a). We calculate the scattering amplitude $f_{\kappa_{e}}(\theta)$ of the
wave scattered by an oil droplet that is coated with the water-based gel and
immersed in an effective medium, with wave number $\kappa_{e}$, whose value is
to be determined at each frequency $\omega$. The radius of the coating $b$ is
expressed as a function of the volume fraction of scatterers $b=a/\phi^{1/3}$
while the presence of the coating enables a multiple scattering contribution
from structural correlations between droplets and the matrix to be included.
Physically, the properties of the effective medium take into account the
effect of this scattering from neighboring particles, resulting in a weakening
of the effective contrast between a scatterer and the surrounding matrix. The
spectral function may be approximated as

$S(\omega,\kappa_{e})=-\textrm{Im}\left\langle G\right\rangle\hskip
5.69046pt\approx-\textrm{Im}\frac{1}{4\pi\eta f_{\kappa e}(0)},$ (3)

where $\left\langle G\right\rangle$ is the average Green function of the
effective medium. For each frequency $\omega$, the algorithm scans trial
values of $\kappa_{e}$ to find the optimum value $\kappa_{e}=k^{\prime}$ that
gives a maximum of the spectral function $S(\omega,\kappa_{e})$ (i.e., the
$\kappa_{e}$ value corresponding to the least scattering) Sheng (2006)111Note
that $S(\omega,\kappa_{e})$ has to be scanned at constant $\omega$ to find the
peaks that correspond to an ultrasonic propagation experiment. The maxima of
$S$ for all frequencies $\omega$ identify the propagating “quasimodes”, which
have a finite width due to radiative damping. These solutions of Eq. (3) allow
the calculation of the “renormalized” forward scattering amplitude
$f_{\kappa_{e}=k^{\prime}}(0)$ from which we can deduce $\ell_{s}(\omega)$. As
shown in Fig. 2, the locations of the peaks of $S$ determine the dispersion
curve $\omega$ vs. $k^{\prime}$, from which $v_{\text{ph}}$ is obtained. For
dilute media, $b\rightarrow\infty$, and the solutions obtained from the peaks
of $S$ approach the ISA predictions.

Figure 2:  (a) Illustration of the scattering problem that has to be solved to
find the dispersion relation of the quasimodes. A droplet (pale orange) is
coated with a layer of gel (blue) and embedded in a uniform medium (textured
green) with wavenumber $\kappa_{e}$. (b) Map of the calculated spectral
function $S(\omega,\kappa_{e})$ at $\phi=25\%$. The color scale [from smallest
to largest $S(\omega,\kappa_{e})$] ranges from purple through dark and light
blue, green and yellow to light and dark red. The peak values of
$S(\omega,\kappa_{e})$ yield the dispersion curve (solid white line, with
black edge), which is compared with the experimental data (solid grey circles,
with black edge). Predictions of the ISA (dot-dashed pink and black curve), as
well as the dispersion line in pure gel (dashed yellow and black line), are
also shown.

Very good agreement is found between the calculated dispersion relation of the
quasimodes and the experimental data (Fig. 2(b)). This good agreement persists
across the entire frequency range, even though it is somewhat challenging to
resolve the quasimode peaks between about 2.5 and 3 MHz where the scattering
is especially strong ($\ell_{s}\approx\lambda/2$).

Figure 3:  (a) Temporal evolution of the average transmitted intensity
$I(\rho,t)$ for different values of the transverse distance $\rho$ and for a
frequency $\omega/2\pi=$ 2.43 MHz. The transmitted intensities were measured
at the surface of the sample, and are normalized so that
$I(0,t)_{\text{max}}=1$. (b) Corresponding $w_{\rho}^{2}(t)/L^{2}$, where $L=$
12 mm is the thickness of the slab. Experimental data are represented by
symbols and SCT predictions by solid lines. The black dashed line represents
the linear behavior of $w^{2}(t)/L^{2}$ expected for classical diffusion.

We characterize next the transport behavior of the multiply scattered waves by
using the transverse confinement method Hu  _et al._ (2008); Cobus  _et al._
(2018). Using a point-like source (here a focusing ultrasonic transducer), the
acoustic transmission across a slab of emulsion was measured with a needle
hydrophone in order to avoid spatial averaging in the detector, and resolve
the position dependence of the multiply scattered wave field along the output
surface of the sample Tallon  _et al._ (2020). When $\ell_{s}\gg\lambda$, wave
interferences can be neglected, the scattered wave transport is analogous to a
random walk process, and the diffusion approximation Page  _et al._ (1995) can
be reliably used to describe the transport of energy by the multiply scattered
waves. In this case, the normalized transmitted average wave intensity
$I(\rho,t)$ has a Gaussian spatial shape (the diffusive halo)

$I(\rho,t)/I(0,t)=\textrm{exp}(-\rho^{2}/w^{2}(t)),$ (4)

where $\rho$ is the transverse distance with respect to the source (for
$\rho=$ 0, the source and receiver are aligned on opposite sides of the slab).
The squared width of the halo, $w^{2}(t)$, increases linearly with time, with
$w^{2}(t)=4D_{B}t$ depending on the (Boltzmann) diffusion coefficient $D_{B}$.
By contrast, when $\ell_{s}\approx\lambda$ or smaller, destructive
interferences slow down the diffusive wave transport, leading to a non-linear
evolution of $w^{2}(t)$ Hu  _et al._ (2008). The ultimate limit is Anderson
localisation, for which $w^{2}(t)$ tends to a constant value Hu  _et al._
(2008); Cherroret  _et al._ (2010) at long times, indicating a trapping of
wave energy in the vicinity of the source.

Our experimental results for the average transmitted intensity at $\phi=25\%$
are presented for several $\rho$ values in Fig. 3(a), from which the temporal
evolution of $w_{\rho}^{2}(t)$ in Fig. 3(b) is determined. Note that the width
squared $w_{\rho}^{2}(t)$, as defined by Eq. (4), now varies slightly with
$\rho$. The clear difference between our experimental data and a linear
diffusive evolution of $w^{2}(t)$ (black dashed line) emphasizes the sub-
diffusive wave transport in this resonant emulsion. In addition to its non-
linear temporal evolution, which was observed over the concentration range
$20\%\leq\phi\leq 30\%$, the dependence of $w_{\rho}^{2}(t)$ on the transverse
distance from the source $\rho$ is another characteristic of sub-diffusive
wave transport, for which the halo deviates from a Gaussian spatial shape Hu
_et al._ (2008); Cobus  _et al._ (2018). These experimental observations are
well described by the “Self-Consistent Theory” (SCT) of localization Skipetrov
and van Tiggelen (2006), which accounts for the renormalization of the
diffusion coefficient by interferences and for the resulting spatial
dependence of $D(\mathbf{r})$ in finite samples. From the best fit of SCT to
the measured $w^{2}(t)$, we find that the scattering strength $k^{\prime}\ell$
is only $4\%$ above its critical value $k^{\prime}\ell_{c}$, below which
Anderson localization is reached. Thus even though the data are quite close to
the mobility edge, the localization regime has not yet been reached.

An advantage of resonant emulsions is their all-fluid nature. This can
simplify both measurements and their interpretations by (i) probing the
wavefield directly inside the medium and (ii) considering only scalar waves in
the models. This in situ measurement could allow the spatial dependence of the
diffusivity $D$ to be checked experimentally in the sub-diffusive regime
Skipetrov and van Tiggelen (2006). These advantages motivate an examination of
the optimal concentration of droplets for maximizing the scattering and
possible localization effects. This effect is assessed in Fig. 4 for the
multiply scattered waves by estimating an effective diffusion coefficient $D$
from the initial slope of the width squared as the concentration is increased
222At $\phi=5\%$ and $40\%$, $D=D_{B}$, whereas at intermediate
concentrations, $D<D_{B}$ due to the interference effects associated with
stronger scattering that lead to the bending over of $w_{\rho}^{2}(t)$ at long
times [3(b)].. Experimentally, we find a minimum in $D$ at $\phi=25\%$, where
the greatest deviation from linearity in $w_{\rho}^{2}(t)$ is also observed,
providing convincing evidence that the scattering strength is maximized at
this concentration. This behavior contrasts with that observed for suspensions
of hard particles (e.g., glass spheres), where the resonances are weaker and
the strongest scattering was found for random close packing Page  _et al._
(1995); Cowan  _et al._ (1998). In Fig. 4(a), we also compare our experimental
emulsion data for $D$ with the predictions of the ISA and the spectral
function approach, using calculations similar to those described in Ref.
Tallon  _et al._ (2020). The relatively good agreement between experiment and
the SpFA model for the diffusion coefficient points to the mechanism
underlying the observation of maximum scattering around $\phi=25\%$. When
there are very strong scattering resonances, the scattering of waves from
nearby scatterers effectively weakens the overall scattering strength relative
to what would be predicted in the ISA, an effect that is encapsulated in the
model _via_ an approximate effective medium whose properties become closer to
those of individual scatterers. As the concentration of scatterers increases,
the scattering strength initially increases simply because there are more
scatterers, but it then passes through a maximum as the renormalization of the
effective medium around each scatterer becomes more pronounced, before
decreasing at higher concentrations where $D$ becomes larger again. While the
SpFA model captures the overall behavior well, it underpredicts the
experimentally observed increase in $D$ at $\phi=40\%$, likely reflecting a
limitation due to approximations in the model.

Figure 4: (a) Effective diffusion coefficient $D$ versus volume fraction
$\phi$ for the frequency 2.5 MHz. Predictions based on ISA and spectral
function are represented in dashed blue and solid red lines respectively. The
experimental values of $D$ are extracted from the linear fits to $w^{2}(t)$ at
short times shown in (b) for $\rho=25$ mm.

Having found the optimal droplet concentration for maximizing the effective
scattering in resonant emulsions, and showing that Anderson localization
cannot be realized in this system, it is of interest to consider alternatives
for reaching this goal in fluid-like suspensions. Soft matter techniques have
great potential for creating suspensions with higher scattering proprieties,
as they allow the design of monodisperse inclusions with controlled sound
speeds and even stronger scattering resonances Brunet  _et al._ (2013); Raffy
_et al._ (2016). For example, soft porous silicone rubber has a very low sound
speed ($v_{1}\approx$ 100 m/s) Brunet  _et al._ (2015), resulting in
suspensions with scattering strength as high as
$k^{\prime}\ell_{\textrm{ext}}\approx$ 0.05 Ba  _et al._ (2017). Such
suspensions are therefore promising candidates once the challenge of
fabricating a large number of these spherical particles has been overcome.

In conclusion, the acoustic scattering strength of fluid suspensions increases
with droplet concentration up to a certain threshold. Beyond this threshold,
due to coupling between near-by resonant particles, the effective scattering
becomes so weakened that the diffusivity actually increases with
concentration. This coupling effect manifests itself as a deviation from the
ISA or L&B theories. However, by taking this effect into account _via_ a
renormalized effective medium surrounding the scatterers, the spectral
function approach provides a good overall description of this reduction of
scattering strength as concentration increases. For our resonant emulsions, we
find an optimal volume fraction of oil droplets of $25\%$ that maximizes the
scattering strength, leading to multiply scattered wave transport with a
significant deviation from diffusive behavior that is the most subdiffusive
possible in this system. Despite the renormalization effects causing the
reduction in effective wave velocity contrast between the surrounding medium
and the scatterers, it should be possible to employ soft matter techniques to
design new resonant suspensions for which this limitation does not prevent
observations of extreme behavior such as Anderson localization. Since the
difference in velocities between scattering inclusions and their surroundings
are similar for soft matter suspensions in acoustics and typical strongly
scattering media in optics, our findings should be relevant to wave transport
in a quite wide range of heterogeneous materials.

This work was supported by the LabEx AMADEus (ANR10-LABX-42) within IdEx
Bordeaux (ANR-10-IDEX-03-02), i.e., the Investissements d’Avenir programme of
the French government managed by the Agence Nationale de la Recherche, and the
Natural Sciences and Engineering Research Council of Canada’s Discovery Grant
Program (RGPIN-2016-06042 and RGPIN-2022-04130).

## References

  * Hu  _et al._ (2008) H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Nat. Phys. 4, 945 (2008).
  * Sheng (2006) P. Sheng,  _Introduction to Wave Scattering, Localization, and Mesoscopic Phenomena_ (Springer, Heidelberg, 2006).
  * Page  _et al._ (1995) J. Page, H. Schriemer, A. Bailey, and D. Weitz, Phys. Rev. E 52, 3106 (1995).
  * Chabanov  _et al._ (2000) A. A. Chabanov, M. Stoytchev, and A. Z. Genack, Nature 404, 850 (2000).
  * Faez  _et al._ (2009) S. Faez, A. Strybulevych, J. H. Page, A. Lagendijk, and B. A. Van Tiggelen, Phys. Rev. Lett. 103, 155703 (2009).
  * Hildebrand  _et al._ (2014) W. K. Hildebrand, A. Strybulevych, S. E. Skipetrov, B. A. Van Tiggelen, and J. H. Page, Phys. Rev. Lett. 112, 073902 (2014).
  * Cobus  _et al._ (2018) L. A. Cobus, W. K. Hildebrand, S. E. Skipetrov, B. A. van Tiggelen, and J. H. Page, Phys. Rev. B. 98, 193901 (2018).
  * Cobus  _et al._ (2016) L. A. Cobus, S. E. Skipetrov, A. Aubry, B. A. van Tiggelen, A. Derode, and J. H. Page, Phys. Rev. Lett. 116, 214201 (2016).
  * Lagendijk  _et al._ (2009) A. Lagendijk, B. A. Van Tiggelen, and D. S. Wiersma, Phys. Today 62, 24 (2009).
  * Skipetrov, S. E. and Page, J. H. (2016) Skipetrov, S. E. and Page, J. H., New J. Phys. 18, 021001 (2016).
  * Tallon  _et al._ (2017) B. Tallon, T. Brunet, and J. H. Page, Phys. Rev. Lett. 119, 164301 (2017).
  * Haberko  _et al._ (2020) J. Haberko, L. S. Froufe-Pérez, and F. Scheffold, Nat. Commun. 11, 4867 (2020).
  * Goïcoechea  _et al._ (2020) A. Goïcoechea, S. E. Skipetrov, and J. H. Page, Phys. Rev. B 102, 220201 (2020).
  * Vynck  _et al._ (2021) K. Vynck, R. Pierrat, R. Carminati, L. S. Froufe-Pérez, F. Scheffold, R. Sapienza, S. Vignolini, and J. J. Sáenz, arXiv:2106.13892v1 (2021), 10.48550/ARXIV.2106.13892.
  * Monsarrat  _et al._ (2021) R. Monsarrat, R. Pierrat, A. Tourin, and A. Goetschy, arXiv:2110.12034 (2021), 10.48550/arXiv.2110.12034.
  * Matis  _et al._ (2022) B. R. Matis, S. W. Liskey, N. T. Gangemi, A. D. Edmunds, W. B. Wilson, V. D. Wheeler, B. H. Houston, J. W. Baldwin, and D. M. Photiadis, Commun. Phys. 5, 21 (2022).
  * Brunet  _et al._ (2012) T. Brunet, S. Raffy, B. Mascaro, J. Leng, and R. Wunenburger, Appl. Phys. Lett 101, 011913 (2012).
  * Mie (1908) G. Mie, Ann. Phys.-Berlin 25, 377 (1908).
  * Lax (1951) M. Lax, Rev. Mod. Phys. 23, 287 (1951).
  * Mascaro  _et al._ (2013) B. Mascaro, T. Brunet, O. Poncelet, C. Aristégui, S. Raffy, O. Mondain-Monval, and J. Leng, J. Acoust. Soc. Am. 133, 1996 (2013).
  * Tallon  _et al._ (2020) B. Tallon, T. Brunet, J. Leng, and J. H. Page, Phys. Rev. B 101, 054202 (2020).
  * Lloyd and Berry (1967) P. Lloyd and M. V. Berry, Proc. Phys. Soc. 91, 678 (1967).
  * Fikioris and Waterman (1964) J. G. Fikioris and P. C. Waterman, J. Math. Phys. 5, 1413 (1964).
  * Jing  _et al._ (1992) X. Jing, P. Sheng, and M. Zhou, Phys. Rev. A 46, 6513 (1992).
  * Note (1) Note that $S(\omega,\kappa_{e})$ has to be scanned at constant $\omega$ to find the peaks that correspond to an ultrasonic propagation experiment.
  * Cherroret  _et al._ (2010) N. Cherroret, S. E. Skipetrov, and B. A. van Tiggelen, Phys. Rev. E 82, 056603 (2010).
  * Skipetrov and van Tiggelen (2006) S. E. Skipetrov and B. A. van Tiggelen, Phys. Rev. Lett. 96, 043902 (2006).
  * Note (2) At $\phi=5\%$ and $40\%$, $D=D_{B}$, whereas at intermediate concentrations, $D<D_{B}$ due to the interference effects associated with stronger scattering that lead to the bending over of $w_{\rho}^{2}(t)$ at long times [3(b)].
  * Cowan  _et al._ (1998) M. L. Cowan, K. Beaty, J. Page, Z. Liu, and P. Sheng, Phys. Rev. E 58, 6626 (1998).
  * Brunet  _et al._ (2013) T. Brunet, J. Leng, and O. Mondain-Monval, Science 342, 323 (2013).
  * Raffy  _et al._ (2016) S. Raffy, B. Mascaro, T. Brunet, O. Mondain-Monval, and J. Leng, Adv. Mater. 28, 1760 (2016).
  * Brunet  _et al._ (2015) T. Brunet, A. Merlin, B. Mascaro, K. Zimny, J. Leng, O. Poncelet, C. Aristégui, and O. Mondain-Monval, Nat. Mater. 14, 384 (2015).
  * Ba  _et al._ (2017) A. Ba, A. Kovalenko, C. Aristégui, O. Mondain-Monval, and T. Brunet, Sci. Rep. 7, 40106 (2017).