Raghavendran Balu
email: rbalu@inria.fr

Teddy Furon
email: teddy.furon@inria.fr

Sébastien Gambs
email: sgambs@irisa.fr

Challenging differential privacy: the case of non-interactive mechanisms

Keywords:

In this paper, we consider personalized recommendation systems in which before publication, the profile of a user is sanitized by a non-interactive mechanism compliant with the concept of differential privacy. We consider two existing schemes offering a differentially private representation of profiles: BLIP (BLoom-and-flIP) and JLT (Johnson-Lindenstrauss Transform). For assessing their security levels, we play the role of an adversary aiming at reconstructing a user profile. We compare two inference attacks, namely single and joint decoding. The first one decides of the presence of a single item in the profile, and sequentially explores all the item set. The latter strategy decides whether a subset of items is likely to be the user profile, and considers all the possible subsets. Our contributions are a theoretical analysis as well as a practical implementation of both attacks, which were evaluated on datasets of real user profiles. The results obtained clearly demonstrates that joint decoding is the most powerful attack, while also giving useful insights on how to set the differential privacy parameter ǫ.

Introduction

Most of the social applications, like recommender systems or private matching, require computing some kind of pairwise similarity between the profiles of different users. Some of the challenges that such systems face include privacy and scalability issues. For instance, privacy concerns arise naturally due to the potentially sensitive nature of profiles, and some users may even refuse to participate if they have no guarantees on the privacy of their profiles.

To address these concerns, the concept of differential privacy [START_REF] Dwork | Differential privacy[END_REF] has been introduced by Dwork in the context of private analysis on statistical databases and has known a widespread adoption in the privacy community. In a nutshell, the main privacy guarantee provided by differential privacy is that for any computation that will be performed on the database, adding or removing a single row from the database will not significantly change the probability of a particular output. Usually, the database is composed of the collection of the individuals' data, and differential privacy protects the privacy of a particular individual, which corresponds to a row of the database. In contrast in our context, the "database" is actually the profile of an individual (e.g., composed of the items he has liked) and therefore the guarantees provided by differential privacy applies to the protection of the items contained in the profile.

One of the usual limits of differential privacy is that each time a differentially private computation takes place, the user loses a little bit of privacy (as measured by the value of the privacy parameter ǫ). Therefore, if this computation takes place too many times, the user may spend all his privacy budget and remains with no privacy left. The adversary is then able to reconstruct almost entirely the user's profile. One possible approach to solve this problem is to sanitize the profile of a user with a non-interactive mechanism compliant with the concept of differential privacy before his publication. In particular, this paper investigates the privacy guarantees offered by two non-interactive mechanisms offering a differentially private representation of profiles: BLIP (BLoom-and-flIP) [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] and JLT (Johnson-Lindenstrauss Transform) [START_REF] Kenthapadi | Privacy via the johnsonlindenstrauss transform[END_REF].

In this paper, we propose two inference attacks that help to assess the privacy guarantee provided by the BLIP and JLT mechanisms. We provide an analysis of the utility and the protection offered by BLIP and JLT against these attacks, by deriving theoretical bounds on the resulting approximation error generated by a specific value of the privacy parameter. Furthermore, we evaluate experimentally the trade-off between privacy and utility achieved by these mechanisms. These attacks helps to better understand the privacy guarantees offered by a differentially-private mechanism, while also enabling the privacy practitioner to tune ǫ experimentally.

A detailed survey on inference attacks on sanitized data can be found in [START_REF] Liu | A survey of attack techniques on privacypreserving data perturbation methods[END_REF] and [START_REF] Chen | A survey of multiplicative perturbation for privacy-preserving data mining[END_REF]. Common attacks include eigen-analysis [START_REF] Guo | On the use of spectral filtering for privacy preserving data mining[END_REF][START_REF] Huang | Deriving private information from randomized data[END_REF], MAP estimation [START_REF] Huang | Deriving private information from randomized data[END_REF], Independent Component Analysis (ICA) [START_REF] Guo | Deriving private information from arbitrarily projected data[END_REF] and distribution analysis [START_REF] Agrawal | On the design and quantification of privacy preserving data mining algorithms[END_REF]. MAP estimation and ICA make direct assumptions on the distribution of the original data, whereas distribution analysis and our approach estimate it from publicly available information. In addition, eigen-analysis makes even stronger assumptions on the representation of data and thus is not generic enough to apply to representations studied in this paper. Furthermore, the possibility of using probabilistic inference techniques to attack sanitized histogram data has been illustrated in [START_REF] Diaconis | Algebraic algorithms for sampling from conditional distributions[END_REF] and [START_REF] Dobra | Measuring the disclosure risk for multi-way tables with fixed marginals corresponding to decomposable log-linear models[END_REF]. In these works, bounds of records count are estimated from histogram of attributes coming from a Markov Chain Monte Carlo (MCMC) simulation. This line of work is different from our approach aiming at reconstructing a user profile from perturbed data. Application of probabilistic inference techniques for parameter estimation on differentially private data is illustrated in [START_REF] Williams | Probabilistic inference and differential privacy[END_REF]. In this work, the authors have also experimentally validated their approach using MCMC on parameter estimation of logistic regression and probabilistic inference of principal components. Although their objective was not directly the reconstruction of data, their approach demonstrates that probabilistic inference is possible on differentially private data.

The outline of the paper is the following. First in Section 2, we give an overview of the concept of differential privacy. Then, we describe two non-interactive differentially private mechanisms in Sections 3 and 4 that have been recently proposed: BLIP (BLoom-and-flIP ) [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] and one based on the Johnson-Lindenstrauss Transform (JLT in short) [START_REF] Kenthapadi | Privacy via the johnsonlindenstrauss transform[END_REF]. These mechanisms transform the profile of a user into a compact representation that estimate the similarity between profiles while hiding the presence or absence of a particular item in the profile (in the sense of differential privacy). Afterwards in Section 5, we provide a theoretical analysis showing that the joint decoding strategy is more powerful than the single decoding strategy in reconstructing the profile of a user. Finally in Section 6, we propose a tractable implementation of this strategy based on the MCMC algorithm before reporting in Section 7 the results on two real datasets.

Differential privacy

In this paper, we are interested in a strong privacy notion called differential privacy [START_REF] Dwork | Differential privacy[END_REF]. Differential privacy aims at providing strong privacy guarantees with respect to the input of some computation by randomizing the output of this computation, and this independently of the auxiliary information that the adversary might have gathered. In our setting, the input of the computation is the profile of a user and the randomized output will be a perturbed version of a compact representation of this profile (e.g., a Bloom filter or a random projection).

Two profiles x and x ′ are said to differ in at most one element or said to be neighbors if they are equal except for possibly one entry.

Definition 1 (Differential privacy [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF]). A randomized function F : D n → D n is ǫ-differentially private, if for all neighboring profiles x, x ′ ∈ D n and for all t ∈ D n :

P[F(x) = t] e ǫ • P[F(x ′ ) = t] .
This probability is taken over all the coin tosses of F and e is the base of the natural logarithm.

The parameter ǫ is public and may take different values depending on the application (for instance it could be 0.1, 0.25, 1.5 or even 10). The smaller the value of ǫ, the higher the privacy but also as a consequence the higher the impact might be on the utility of the resulting output. A relaxed notion differential privacy called (ǫ,δ)-differential privacy [START_REF] Dwork | Our data, ourselves: Privacy via distributed noise generation[END_REF], can be seen as a probabilistic variant in which the guarantees of differential privacy hold with probability of 1δ.

Originally, differential privacy was developed within the context of private data analysis and the main guarantee is that if a differentially private mechanism is applied on a dataset composed of the personal data of individuals, no output would become significantly more (or less) probable whether or not a single participant contributes to the dataset. This means that observing the output of the mechanism only gains negligible information about the presence (or absence) of a particular individual in the database. This statement is a statistical property about the behavior of the mechanism (i.e., function) and holds independently of the auxiliary knowledge that the adversary might have gathered. More specifically, even if the adversary knows the whole database but one individual row, a mechanism satisfying differential privacy still protects the privacy of this individual row. In our setting, the database that we want to protect is the profile of a user and the objective of a differentially private mechanism is to hide the presence or absence of a particular item in the profile.

Dwork, McSherry, Nissim and Smith have designed a generic technique, called the Laplacian mechanism [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF], that achieves ǫ-differential privacy for a function f by adding random noise to the true answer of f before releasing it. Subsequently, McSherry and Talwar have proposed the exponential mechanism [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF] which unlike the Laplacian mechanism that works only for functions with numerical output, provides differential privacy for functions whose output is more structured (e.g., graphs or trees). Both previous mechanisms (i.e., Laplacian and Exponential mechanisms) are interactive as they require a two-way communication protocol between the curator (the entity in charge of the database) and the client performing the query. Therefore, the curator has to be online in order to receive the query and prepare the associate response to this query.

On the other hand, a non-interactive mechanism computes some function from the original database and releases it once and for all, which corresponds to a one-way communication protocol. The output released by the non-interactive mechanism can later be used by anyone to compute the answer to a particular class of queries (usually not just a single specific query), without requiring any further interactions with the curator. It is important to understand that the answer is computed from the output released by the non-interactive mechanism, thus after publishing this output the curator can go offline. One particular type of non-interactive mechanism is the generation of a synthetic dataset that allows the answer to certain class of queries (but not necessarily all) to be approximated. Examples of non-interactive mechanisms for differential privacy include [START_REF] Beimel | Distributed private data analysis: on simultaneously solving how and what[END_REF][START_REF] Li | Compressive mechanism: utilizing sparse representation in differential privacy[END_REF].

In the next sections, we describe two non-interactive mechanisms that have recently been proposed. The first mechanism is based on randomizing a Bloom filter representation of the profile [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] while the second relies on the application of the Johnson-Lindenstrauss transform and the addition of noise [START_REF] Kenthapadi | Privacy via the johnsonlindenstrauss transform[END_REF]. Both mechanisms preserve some global properties such as the ability to compute a distance between two profiles while hiding the details of the profiles themselves.

BLIP

The main objective of BLIP [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] is to prevent the adversary from learning the presence (or absence) of an item in the profile of a user by observing the Bloom filter representation of this profile. Our theoretical analysis provided in Section 5 is based on the model of profiles and the BLIP sanitization described thereafter.

Setup of BLIP

The setup that we consider for the theoretical analysis is the following. We assume that a profile P is a list of c items randomly picked from a set of N ∈ N ⋆ possible items: P = {j 1 , . . . , j c }. We denote the set of items by [N ], with

[N ] {1, . . . , N } and the set of all possible profiles by P. This set is a subset of the power set of [N ] and we have |P| = N c . For the moment, we make the assumption that c is publicly known, but this hypothesis will be lifted later by inferring this value directly from the Bloom filter.

The profile is first encoded in the form of a Bloom filter, which is a binary string of L bits. Each item j ∈ P is hashed through K different hash functions (h 1 , . . . , h K ). Each hash function yields a position h k (j) in the Bloom filter, pseudo-randomly selected based on the identifier of the item j. One simple technique to implement this is to rely on K cryptographic hash functions modulo L. We call the codeword X j associated to item j the following string of L bits:

X j (ℓ) = 1 if ∃k ∈ [K] such that h k (j) = ℓ, 0 otherwise. ( 1
)
The Bloom filter associated to the profile P = {j 1 , . . . , j c } is denoted by B P and computed as the aggregation of the codewords:

B P = X j1 ∨ . . . ∨ X jc , (2)
in which ∨ denotes the bit-wise (inclusive) OR operator. Our presentation of Bloom filters is different than usual to stress the link with our general problem.

The BLIP mechanism adds noise to the Bloom filter representation of a profile before publishing it. We denote the output of BLIP by BP :

BP = B P ⊕ N, (3)
in which ⊕ corresponds to the bit-wise logical (exclusive) XOR operator and N ∈ {0, 1} L is a random binary string of size L, whose symbols are i.i.d. (independent and identically distributed) as a Bernoulli distribution B(p ǫ ) (i.e. , N (ℓ) ∈ {0, 1} and

P[N (ℓ) = 1] = p ǫ , ∀ℓ ∈ [L]
). Alaggan, Gambs and Kermarrec [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] proved that the BLIP mechanism ensures ǫ-differential privacy for the items of the profile if

p ǫ = 1/(1 + e ǫ/K ). ( 4
)

The simple model

We assume that the hash functions produce independently random outputs, which means that the probability that h k (j) "points" to a given index is 1/L. This assumption implies that the bits of the codewords can be modeled as independent Bernoulli random variables:

X j (ℓ) ∼ B(p), ∀(j, ℓ) ∈ [N ] × [L] with p P[X j (ℓ) = 1] = 1 -1 - 1 L K . ( 5
)
For a random P composed of c items, we have

B P (ℓ) ∼ B(π c ), ∀ℓ ∈ [L],
with

π c P[B P (ℓ) = 1] = 1 -(1 -p) c = 1 -1 - 1 L cK . (6)
As for the BLIP, BP contains i.i.d. random symbols BP (ℓ) ∼ B(π c ) with

πc P[ BP (ℓ) = 1] = (1 -p ǫ )π c + p ǫ (1 -π c ). ( 7
)

More complex models

This subsection presents two possible extensions of the simple model, in which we no longer assume that c is fixed in advance and publicly known.

To account for this, we introduce the probability P[|P | = c], in which |P | denotes the number of items in P . Then, we have to replace π c by:

π c → π = c>0 π c P[|P | = c]. (8)
This new expression leads to π = (1p ǫ )π + p ǫ (1π). Not knowing c may not be a big challenge for the adversary because he can easily infer the number of items in a profile. The quantity ω( BP )/L, in which ω(.) is the Hamming weight of a binary string (the number of bits set to one), is an unbiased estimator of πc . Inverting ( 7) is possible when p ǫ = 1/2 (i.e. , ǫ > 0) since p ǫ is public:

πc = ω( BP )/L -p ǫ 1 -2p ǫ , (9)
which in turn gives an estimator ĉ by inverting [START_REF] Guo | On the use of spectral filtering for privacy preserving data mining[END_REF]. In the same way, a confidence interval for πc based on ω( BP )/L yields a confidence interval [c min , c max ] on c.

An even more refined model consists in taking into account the popularity of the items. Indeed, popular items impact the Bloom filter by ensuring that some of its bits are more likely to be set to one. To tackle this issue, we still pretend that the bits are independent but distributed according their own Bernoulli law:

B P (ℓ) ∼ B(π(ℓ)), ∀ℓ ∈ [L]. The same model holds for the BLIP symbols: BP (ℓ) ∼ B(π(ℓ)), with π(ℓ) = (1 -p ǫ )π(ℓ) + p ǫ (1 -π(ℓ)).

JLT

Kenthapadi and co-authors [START_REF] Kenthapadi | Privacy via the johnsonlindenstrauss transform[END_REF] proposed another mechanism to prevent the adversary from learning the presence (or absence) of an item in the profile, although their scheme tackles a different data type (i.e., real vector). In the sequel, we denote this proposal by JLT because it is based on the Johnson-Lindenstrauss Transform.

Description

The profile is encoded in the form of a real vector of length L as follows. A codeword X j associated to item j is a real vector. Its L components have been independently and identically drawn such that X j (i)

i.i.d. ∼ N (0, 1/L), ∀(i, j) ∈ [L] × N . The codebook (X 1 , X 2 , • • • , X N ) is generated
once for all and is public. Profile P is encoded into vector Y P = j∈P X j , then the user adds a noise N (private data) before publishing ỸP = Y P + N. The authors of [START_REF] Kenthapadi | Privacy via the johnsonlindenstrauss transform[END_REF] recommend a white Gaussian noise:

N (i) i.i.d. ∼ N (0, σ 2 ). According to [2, Lemma 2], if L ≥ 2(log(N ) + log(2/δ)), σ ≥ 4 ǫ log(1/δ) and ǫ < log(1/δ) (10)
then this mechanism complies with (ǫ, δ)-differential privacy (for 0 < δ < 1).

A simple probabilistic model

The adversary does not know the profile P and therefore he models the observation ỸP as a white Gaussian noise since ỸP is the sum of c + 1 white Gaussian noises. As these patterns are statistically independent, their powers sum up so that ỸP (i)

i.i.d.
∼ N (0, σ 2 + c/L). We assume now that σ 2 is a recommended noise power, and thus that it is a public parameter. This allows the adversary to estimate the number of items in profile P in the following manner:

ĉ = L L -1 L i=1 ỸP (i) 2 -Lσ 2 . ( 11
)
Consider now the case in which the adversary knows that the item j is in the profile. This knowledge stems into a refined statistical model of the observation: ỸP (i)

i.i.d.
∼ N (X j (i), σ 2 + (c -1)/L). In the same way, knowing the profile P ends up with ỸP (i)

i.i.d. ∼ N j∈P X j (i), σ 2 .

Theoretical analysis

In this section, we propose two decoders that can be used by an adversary to reconstruct the profile of a given user out of his public representation. This analysis is detailed for the BLIP mechanism, but similar concepts hold for the JLT scheme. The expressions of the information theoretical quantities are given in Appendix 9 for BLIP and Appendix 10 for JLT.

Single decoder

From the observation of one BLIPed representation b, the adversary would like to infer which item belongs to the original profile. The adversary can conduct this inference by analyzing the L symbols of b and making an hypothesis test about the presence of item j in the underlying profile.

-H 0 : Item j is not in the profile, which means that the observed BLIP symbols are statistically independent from the symbols of codeword X j :

P[ BP (ℓ), X j (ℓ)] = P[ BP (ℓ)]P[X j (ℓ)], ∀ℓ ∈ [L].
-H 1 : Item j belongs to P , and thus there is a slight dependency between the symbols of the observed BLIP and that of codeword

X j : P[ BP (ℓ), X i (ℓ)] = P[ BP (ℓ)|X i (ℓ)]P[X i (ℓ)], ∀ℓ ∈ [L].
For a given item, this test may make two types of error: 1) False positive rate α 1 : The probability of detecting the presence of an item that does not belong to the profile; 2) False negative rate α 2 : The probability of missing the presence of an item that belongs to the profile. Information theory gives an upper bound on the performance of the test thanks to the Stein's lemma. More precisely, for a given α 2 , the probability of false positive cannot be lower than

α 1 ≥ e -(I( BP ;X)+1)/(1-α2) , (12)
in which I( BP ; X) is the mutual information between a BLIPed filter and the codeword of an item of the profile. This test concerns a particular item, but an adversary that wants to reconstruct the whole profile needs to repeat it for the whole ensemble of size N . This repetition increases the global probability of false positive η 1 :

η 1 = 1 -(1 -α 1 ) N -c N α 1 , (13)
in which we assume that N α 1 ≪ 1 and c ≪ N . η 1 is the probability that at least one item not in the profile is detected as belonging to the profile. At the end, for targeted error probabilities (α 2 , η 1 ), inequality [START_REF] Williams | Probabilistic inference and differential privacy[END_REF] constraints the size of the item ensemble the adversary can deal with:

log(N ) ≤ I( BP ; X) 1 -α 2 + log η 1 . (14)
The last inequality stresses the important role of I( BP ; X). Appendices 9 and 10 provide expressions of this quantity for the BLIP and JLT mechanisms.

Joint decoder

Let us consider another strategy. From the observation b, the adversary would like to test whether P was the original profile that gave birth to this BLIPed representation. The difference with the previous approach is that the presence of items are not tested independently but jointly, hence the name "joint decoder".

Basically, the analysis is the same as previously except that the information theoretic quantity is now I( BP ; P ) = I( BP ; (X j1 , . . . , X jc )) and that the ensemble of profiles is much bigger. Roughly, log(|P|) ≈ c log N , thus we have:

log(N ) ≤ I( BP ; P ) c(1 -α 2 ) + log η 1 . (15)
Stated differently, the performance of this approach is driven by the quantity I( BP ; P )/c. Theorem [18, Eq. (3.4)] states that I( BP ; (X j1 , . . . , X jc ))/c ≥ I( BP ; X j ), which means that considering the items jointly yields better performances. Appendices 9 and 10 provide expressions of this quantity for respectively the BLIP and JLT mechanisms. For this first scheme, subsection 9.2 shows that the difference I( BP ; (X j1 , . . . , X jc ))/c -I( BP ; X j ) can be be substantial for practical setups. We also provide upper bounds simply depending on ǫ.

Practical decoders

The previous section can be summarized as follows: joint decoding is theoretically more powerful than single decoding. However, no complexity argument has been so far taken into account. This section deals with this issue by proposing practical implementations of a single and a joint decoder. Again, we take the example of BLIP but our approach is more generic as it works also with JLT.

Single decoders

In practice, a single decoder computes from the observed BLIPed profile a score s j for any item j ∈ [N ], which reflects the likelihood of belonging to the profile (i.e., the most likely item has the highest score). The score is compared to a threshold to decide whether or not the item should be included in the reconstructed profile. The complexity of this single decoder is O(N ) since it is exhaustive and goes through all the possible items.

As a practical implementation, we propose the Maximum Likelihood decoder in which the score s j = log P[ BP = b|j∈P ]

P[ BP = b]
equals, by independence of the symbols:

s j = n 11 log 1 -p ǫ π + n 01 log p ǫ 1 - π , with: (16
)

n 11 = |{ℓ ∈ [L]| b(ℓ) = 1 AND X j (ℓ) = 1}|, (17)
n 01 = |{ℓ ∈ [L]| b(ℓ) = 0 AND X j (ℓ) = 1}|. ( 18
)
This decoder is derived from models that are more realistic in which π c ≈ π c-1 ≈ π, so that the score of item j only takes into account the (n 11 + n 01 ) symbols in which X j (ℓ) = 1 (i.e., at most K symbols over L).

Joint decoder

In practice, a joint decoder computes from the observed BLIPed filter a score for any profile P ′ ∈ P, which reflects the likelihood that P ′ is the true profile. This score is computed by taking into account L symbols but the complexity of a joint decoder is proportional to |P| (i.e., O(N c )), which is computationally expensive. Yet, there exists at least three possible approaches that approximate joint decoding with a reasonable complexity: 1) Monte Carlo Markov Chain (MCMC) [START_REF] Knill | Interpretation of pooling experiments using the Markov chain Monte Carlo method[END_REF][START_REF] Furon | Decoding fingerprints using the Markov Chain Monte Carlo method[END_REF], 2) Belief Propagation Decoder [START_REF] Sejdinovic | Note on noisy group testing: asymptotic bounds and belief propagation reconstruction[END_REF] and 3) Joint Iterative Decoder [START_REF] Meerwald | Toward practical joint decoding of binary Tardos fingerprinting codes[END_REF].

In this paper, we investigate the first approach. The MCMC decoder is based on two key ideas. First, it receives as input an observed BLIPed filter b and then creates a Markov Chain that will be used to sample profiles according to the posterior distribution P[P | b]. This sampling requires a burn-in period after which the Markov Chain has converged. Once this convergence has occurred, it samples profiles with the targeted posterior distribution. During a second phase, some profiles are sampled and statistics are computed such as the marginal a posteriori distribution P[j ∈ P | b] that item j belongs to the true profile.

Posterior distribution. The objective is to sample profiles according to the posterior distribution P[P | b], which can be written as:

P[P | b] = P[ BP = b|P ]P[P ] P[ BP = b] . (19)
In this equation, P[P ] is the a priori probability of P . To simplify our presentation, we consider only the simple model exposed in Section 3.2. We denote by |P | the size of profile P (i.e., the number of his items), and we set by

P[P ] = 0 if |P | = c
, and 1/|P| otherwise. Any profile is equally likely provided it has exactly c items. When we use more realistic models in our experimental work, the prior will be substantially different. We denote by ω(B) the Hamming weight of a binary vector B (i.e., the number of bits set to 1). The probability P[ BP = b|P ] = P[N = B P ⊕ b] has the following expression

P[ BP = b|P ] = p ω(B P ⊕ b) ǫ (1 -p ǫ ) L-ω(B P ⊕ b) . ( 20
)
The evaluation of the last quantity P[ BP = b] in ( 19) is more involved:

P[ BP = b] = P ∈P P[ BP = b|P ]P[P ]. (21)
It requires a screening of P, which is intractable for large c and N , which is why we will rely on the Markov chain.

Markov Chain. A Markov Chain is an iterative process with an internal state (i.e., a profile in our case) taking value P (t) at iteration t. The next iteration draws a new state P (t+1) according to a transition probability distribution P[P (t+1) |P (t) ]. The Markov Chain is initialized randomly at state P (0) . The probability distribution of transitions is crafted with care to enforce a convergence of the distribution of sampled profile P (t) to the posterior P[P | b] of [START_REF] Knill | Interpretation of pooling experiments using the Markov chain Monte Carlo method[END_REF] as t → ∞ (see Section 6.3). In practice, the convergence occurs after the first T iterations, the so-called burn-in period. Once this period has passed, it means that the Markov Chain has forgotten its starting point (i.e., the samples are now independent of P (0) ) and that the distribution of the sample profiles has converged.

Monte Carlo method. After the burn-in period, the Markov Chain keeps on sampling for M more iterations. The marginal a posteriori probabilities are then estimated by a Monte Carlo method, which computes the empirical frequency that item j is present in sample P (t) :

P[j ∈ P | b] = |{t ∈ [T + 1, T + M ]|j ∈ P (t) }|/M. ( 22
)
From these estimations, several post-processing are possible such as:

-inferring the most likely items of the true profile by ranking them in decreasing marginal probabilities, -reconstructing the profile as the set of items whose marginal probability is above a given threshold, -reconstructing the profile as the set of items with highest marginal.

Transition probabilities

Algorithmic coding of a profile. Section 3.3 describes how to infer from the observed BLIP a maximum number c max of items of the corresponding profile. In this algorithm, we code a profile as a vector of c max components taking values in [N ] ∪ {0}. Some of these components may take the value "0" meaning an "empty item", while the others have different values (i.e., there is no pair of non-zero components with the same value). For instance, for c max = 5, P = (0, 3, 2, 0, 4) represents the profile of 3 items: #2, #3 and #4.

We define V(P 0 , i) as the neighborhood of profile P 0 in the following manner:

V(P 0 , i) = {P ∈ P|P (k) = P 0 (k) ∀k = i}. ( 23
)
This neighborhood profile is the set of all profiles whose coding differs at most from the i-th component. Note that P 0 ∈ V(P 0 , i). If P 0 (i) = 0, this neighborhood comprises profiles having at most one more item. Otherwise if P 0 (i) > 0, this neighborhood contains profiles having at most one different item (i.e., P 0 (i) is substituted by another item) and one profile having one less item (i.e., item P 0 (i) is substituted by 0, the "empty item").

Multi-stage Gibbs sampling. Instead of computing the transition probabilities for all the possible profiles, we restrict the transitions to the neighborhood of the actual state. At the iteration t + 1, an integer i is first uniformly drawn in [c max ] that indicates the subset V(P (t) , i). Then, the following transition probability distribution is computed: ∀P ∈ V(P (t) , i)

P[P (t+1) = P |P (t) ] = P[ BP = b|P ]P[P ] P ′ ∈V(P (t) ,i) P[ BP ′ = b|P ′ ]P[P ′ ] (24)
Iteration t + 1 ends by randomly drawing state P (t+1) from this distribution. This choice of probabilistic transitions is called a multi-stage Gibbs sampler with random scan [START_REF] Robert | Monte Carlo statistical methods[END_REF]Alg. A.42]. It guarantees that the law of sampled profiles converges to the stationary distribution P[P | b], which legitimates our approach [START_REF] Robert | Monte Carlo statistical methods[END_REF]Sect. 10.2.1]. The unknown multiplicative constant P[ BP = b] in [START_REF] Knill | Interpretation of pooling experiments using the Markov chain Monte Carlo method[END_REF] has disappeared in the ratio. This transition probability distribution only depends on the priors P[P ] (which depends on the mathematical model of a profile), and the conditional probabilities P[ BP = b|P ] (which depends on the privacy-preserving mechanism). For instance, for the JLT mechanism, P[ ỸP = ỹ|P ] ∝ exp(-ỹ -j∈P X j 2 /2σ 2 ).

Experiments

Setup

In this section, we test the inference attacks designed on two real datasets: Digg and MovieLens. The Digg dataset has been collected on a social news aggregator and the profile of a user is composed of the news he has read. The MovieLens dataset is a snapshot from a movie recommendation site and in this dataset the profile of a user is composed of the movies he likes. For the experiments, we split both datasets into two parts : the training set and the testing set. The characteristics of these datasets are summarized in Table 1, in which c avg is the average number of items per profile and sparsity is the average occupancy of items among the user profiles.

During the experiments, we assume that the adversary has access to some raw profiles of users to estimate the item priors (i.e., popularities of items). This is similar to assuming that the adversary has access to some global information about the general distribution of items in the population. We rely on the training dataset for computing the frequencies of items while the testing dataset is used solely for evaluating the performance of the attacks. In terms of parameters, for BLIP we set the number of hash functions K = 20 and the number of bits of the representation to L = 5, 000. The values of ǫ are from the set {59, 28, 17, 8, 6, 5, 3, 2, 0}, which equivalently translate to the corresponding flipping p ǫ from the range {0.05, 0.2, 0.3, 0.4, 0.42, 0.44, 0.46, 0.48, 0.5}. For the JLT scheme, we set the size of the representation L to 1, 000. L is set to a lower value as the representation, a dense real vector, is richer than the binary version of BLIP. The privacy parameter ǫ takes value in the set {600, 6, 3, 2, 1, 0.75, 0.5, 0.25, 0.1}, which translates into a noise level σ in {0, 1, 2, 3, 6, 8, 12, 24, 61}.

For MCMC, we used a burn-in period of T = 1, 000 samples and estimation sample size of M = 19, 000 for all the experiments. In practice, we observed that the performance is not very sensitive to the burn-in period length. As with other MCMC based approaches proper initialization for sampling is highly desirable for a faster convergence to the stationary distribution. We used the input public representation of the profile to estimate ĉ and started with ĉ random items. A poor estimation of ĉ has to be traded-off with a longer burn-in period. We also prefilter items that are to be tested against the public profile for joint decoder, to reduce the search space. To realize this, we first predict the f × ĉ most probable items for a given profile (f ∈ [2, 6]) using single decoder and then run the joint decoder on the filtered items to return ĉ items. This prefiltering decreases significantly the running time of the algorithm without impacting the prediction as only unlikely items will not be considered by the joint decoder.

Reconstruction attacks

We benchmark four attacks that produce a score per item:

-The single decoder described in [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF].

-The popularity-based attack in which the score of an item is its prior estimated from the training data, independent of the given public representation.

-Our MCMC joint decoder with and without priors (i.e., with flat priors) in which the scores are the estimated marginal a posteriori probabilities.

Reconstruction P is then the list of the top ĉ items ranked based on their scores. We measure the performance of a reconstruction attack by computing the cosine similarity between the reconstruction P and the true profile P as expressed in [START_REF] Alvim | On the relation between differential privacy and quantitative information flow[END_REF] for all the profiles of the testing set.

Afterwards, we compute the following statistics: average, the 10% and the 90% quantiles of the cosine similarities. The plots in Figure 1 show that the performance of the reconstruction attack is better for high values of ǫ while it degrades as ǫ → 0. In this case, p ǫ → 0.5 and every profile becomes equiprobable so that inferring the original profile becomes impossible. In addition, ĉ depends on ǫ and low value results in a poor estimation of ĉ, which impacts the similarity measure as only top ĉ items of the prediction is considered in the reconstructed profile. As the estimation of ĉ is performed similarly for all the four attacks, the performance drop is common to all of them. Overall the performance of our MCMC attack is better than the single decoder of [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] for almost all ǫ values over the two datasets. Another way to see this is to find the range of ǫ in which a given attack performs worse than the baseline (i.e., the popularity-based attack). For instance, by setting ǫ = 8, the designer is sure that the single attack is no longer a threat. However, a skilled adversary can reconstruct almost 50% of the profile thanks to our MCMC attack.

Taking into account the prior of items improves the efficiency in the reconstruction significantly, provided that the estimation is reliable. This improvement is clearly observed on the MovieLens dataset. As for the Digg setup, priors of the training set do not generalized to the test set, hence they do not help much. We conducted the same experiment with the JLT scheme. The figure is not included in the paper due to a lack of space, but the results that we obtained are very close from the one of BLIP and thus we can draw the same conclusions.

Identifying the presence of an item

When ǫ is very small, Figure 1 clearly shows that the adversary cannot hope to reconstruct the full profile. In this situation, we evaluate the prediction of top R items, with R ≪ c, as another assessment of the privacy guarantees. The success is measured in terms of the mean Average Precision at R (mAP@R) given in (26), which is the mean over the Q profiles in the test dataset of the average of the precisions at rank 1 ≤ r ≤ R. The precision(r) refers to the fraction of correct items out of the top r predicted items. The mAP is sensitive to the order of the correct results and is a better gauge of the quality of a ranking.

mAP@K = 1 Q Q q=1 1 R R r=1
precision q (r) . The characteristics of mAP@R depicted in Figure 2 are almost similar to the exact reconstruction measurement. Even if the exact reconstruction of profile is hardly possible for a given ǫ, predicting the top R items work. For instance, the maximum reconstruction for ǫ = 0 for Movielens is 0.23 whereas the mean average precision is close to 0.5. The same conclusion holds for the Digg dataset.

Utility-privacy trade-off

Finally, we also studied the achievable trade-off between privacy and utility. Since BLIP and JLT are used for similarity estimation, we quantify the utility in terms of the recall, which is defined as the probability of identifying the knearest neighbors (we set k = 10 in our experiments). In this experiment, we measure privacy as 1cos(P, P ) (see [START_REF] Alvim | On the relation between differential privacy and quantitative information flow[END_REF]) based on the joint decoder. Figure 3 illustrates the utility-privacy trade-off obtained for various ǫ. The trade-off is almost similar on the two datasets. The privacy preserving properties of JLT transform is slightly better than BLIP, at least for the parameters we used in our simulation. This difference in performance is due partially to the representation superiority of dense real vector over binary vector. However, BLIP offers a more compact representation of the profile (5, 000 bits versus 1, 000 scalars). The plot is helpful in fixing ǫ giving good utility without compromising much on privacy.

Conclusion

In differential privacy, the trade-off between utility and privacy is set by the parameter ǫ. However, being able to choose an appropriate value for this parameter is still an open research question, which has not been deeply investigated, with a few exceptions [START_REF] Lee | How much is enough? Choosing ǫ for differential privacy[END_REF][START_REF] Alvim | On the relation between differential privacy and quantitative information flow[END_REF]. In this paper, we have made a step forward to answer this question by proposing two generic inference attacks, namely single and joint decoding, whose objective is to reconstruct the profile of a user out of a differentially-private representation produced through a non-interactive mechanism. The first inference attack decides of the presence of a single item and sequentially explores all the item set, while the latter strategy decides whether a subset of items is likely to be the user profile and considers all possible subsets.

We have evaluated the effectiveness of the attack on two schemes producing differentially private representations: BLIP (BLoom-and-flIP) [START_REF] Alaggan | BLIP: Non-interactive Differentially-Private Similarity Computation on Bloom Filters[END_REF] and JLT (Johnson-Lindenstrauss Transform) [START_REF] Kenthapadi | Privacy via the johnsonlindenstrauss transform[END_REF]. Our theoretical analysis as well as the experimental results clearly shows that joint decoding is more powerful than single decoding. Overall, we believe that this attack helps better understanding the privacy guarantees offered by a wide class of differentially-private mechanisms (interactive or not) as well as for the privacy practitioner to tune experimentally ǫ to ensure the maximum utility without compromising much on privacy. 9 Appendix A: BLIP mechanism

The probabilities πc and πc-1 appear in (27) because we assume that the profiles are of identical size c. When considering more complex but also more practical models, this difference vanishes as πc and πc-1 are replaced by π:

I( BP ; X) ≈ Lp(h b (π) -h b (p ǫ )). (29)
As for the joint decoding, Bloom filter being a deterministic process, we write:

I
K = log(2) L c , (31)
which is the recommended number of hash functions in Bloom filter design, and we apply the model of Section 3.2 with ǫ = 20. For a given c, too small L means too few observed symbols for reliably estimating the profile. Too large L implies a big K and therefore, p ǫ tends to 1/2 according to (4). Figure 4 (right) shows that I( BP ; P )/c -I( BP ; X) can be substantial: a joint decoding allows the adversary to tackle up to 3.5 (i.e. e 1.25 ) times more items.

Upper bounds. As ǫ → 0, p ǫ → 1/2 as well as π, so that I( BP ; X) → 0 and also I( BP ; P ) → 0. When ǫ = 0, observing the BLIP is useless since it brings no information. In this situation, neither the single nor the joint decoding can do anything. We can bound the quantity in common in both expressions as follows: Typical Bloom filter setup. Figure 4 shows that estimating an important number of items is possible provided that L grows linearly with c. Indeed, it is also common practice in the design of Bloom filter to set:

L = -c log(P f p ) (log 2) 2 , (33)
in which P f p is the probability of false positive of the Bloom filter (i.e., to detect the presence of an item not belonging to P ). Inserting (31) and (33) in the expression of the mutual informations, we get quantities independent of c: This shows that if the Bloom filter is properly designed, the power of the attack does not depend on c but solely on the values oflog(P f p ) and ǫ. Moreover, the joint decoder is 1/ log(2) ∼ 1.44 more "powerful" than the single decoder.

10 Appendix B: JLT mechanism

The same analysis holds for the JLT representation described in Section 4. The main difference lies in the fact that we manipulate differential entropies because the JLT representation is a real vector. The quantities at stake respectively for the single and joint decoders are upper bounded, thanks to conditions (10)

I( ỸP ; X) = L 2 log 1 + 1 (c -1) + Lσ 2 ≤ ǫ 32 + 2ǫ(c -1)L , (36)
I( ỸP ; P ) c = L 2c log 1 + c Lσ 2 ≤ ǫ 32 , (37)

  cos(P, P ) = |P. P | |P || P |

Fig. 1 .

 1 Fig.1. Values of the cosine similarity (average, 10% quantile and 90% quantile) of BLIP for MCMC with prior, with no prior and single decoding for various ǫ on Movielens (left) and Digg (right) dataset.

Fig. 2 .

 2 Fig. 2. Mean Average Precision for R = 10 for BLIP for MCMC with prior, with no prior and single decoding for various ǫ on Movielens (left) and Digg (right) dataset.

Fig. 3 .

 3 Fig. 3. Utility against privacy for BLIP and JLT for various ǫ on Movielens (left) and Digg (right) datasets.

9. 1

 1 Single and joint decodingWe have I( BP ; X) = H( BP ) -H( BP |X), in which H is the (Shannon) entropy of a random variable. With the simple model detailed in Section 3.2, we get thatI( BP ; X) = L(h b (π c ) -(1p)h b (π c-1 )ph b (p ǫ )),(27)with h b (p) the entropy of a Bernoulli distribution B(p) (in hats): h b (p) -p log(p) -(1p) log(1p) = h b (1p).

  ( BP ; P ) = I( BP ; B P ) = H( BP ) -H( BP |B P ) = H( BP ) -H(N) = L(h b (π c )h b (p ǫ )).Figure 4 (left) shows I( BP ; P )/c as a function of c and L. From a particular (c, L), we set

hFig. 4 .

 4 Fig. 4. (Left): Mutual information of the joint decoder I( BP ; P )/c in nats as a function of (c, L). (Right): Difference I( BP ; P )/c -I( BP ; X) in nats as a function of (c, L).

Table 1 .

 1 Datasets characteristics Nb of users Training set size Testing set size N cavg Sparsity %

	Digg	531	331	200	1237 317 25.63%
	MovieLens	943	600	343	1682 106	6.30%


A Grankin

E Brion

E Bimbard

R Boddeda

I Usmani

A Ourjoumtsev

P Grangier

Quantum statistics of light transmitted through an intracavity Rydberg medium

We theoretically investigate the quantum statistical properties of light transmitted through an atomic medium with strong optical non-linearity induced by Rydberg-Rydberg van der Waals interactions. In our setup, atoms are located in a cavity and non-resonantly driven on a two-photon transition from their ground state to a Rydberg level via an intermediate state by the combination of the weak signal field and a strong control beam. To characterize the transmitted light we compute the second-order correlation function g (2) (τ ). The simulations we obtained on the specific case of rubidium atoms suggest that the bunched or antibunched nature of the outgoing beam can be chosen at will by appropriately tuning the physical parameters.

I. INTRODUCTION

In an optically non-linear atomic medium, dispersion and absorption of a classical light beam depend on powers of its amplitude [START_REF] Boyd | Nonlinear Optics[END_REF]. At the quantum level, dispersive optical nonlinearities translate into effective interactions between photons. The ability to achieve such strong quantum optical non-linearities is of prominent importance in quantum communication and computation for it would allow to implement photonic conditional two-qubit gates.

The standard Kerr dispersive non linearities obtained in non-interacting atomic ensembles, either in off-resonant two-level or resonant three-level configurations involving Electromagnetically Induced Transparency (EIT), are too small to allow for quantum non-linear optical manipulations. To further enhance such non-linearities, EIT protocols were put forward in which the upper level of the ladder is a Rydberg level. In such schemes, the strong van der Waals interactions between Rydberg atoms result in a cooperative Rybderg blockade phenomenon [START_REF] Lukin | Dipole Blockade and Quantum Information Processing in Mesoscopic Atomic Ensembles[END_REF][START_REF] Saffman | Quantum information with Rydberg atoms[END_REF][START_REF] Comparat | Dipole blockade in a cold Rydberg atomic sample[END_REF], where each Rydberg atom prevents the excitation of its neighbors inside a "blockade sphere". This Rydberg blockade deeply changes the EIT profile and leads to magnified non-linear susceptibilities [START_REF] Pritchard | Cooperative Atom-Light Interaction in a Blockaded Rydberg Ensemble[END_REF][START_REF] Dudin | Strongly Interacting Rydberg Excitations of a Cold[END_REF][START_REF] Peyronel | Quantum nonlinear optics with single photons enabled by strongly interacting atoms[END_REF][START_REF] Maxwell | Storage and Control of Optical Photons Using Rydberg Polaritons[END_REF]. In particular, giant dispersive non-linear effects were experimentally obtained in an off-resonant Rydberg-EIT scheme using cold rubidium atoms placed in an optical cavity [START_REF] Parigi | Observation and Measurement of Interaction-Induced Dispersive Optical Nonlinearities in an Ensemble of Cold Rydberg Atoms[END_REF][START_REF] Stanojevic | Dispersive optical nonlinearities in an EIT-Rydberg medium[END_REF]. In this paper, we theoretically investigate the quantum statistical properties of the light generated in the latter protocol. Note that, contrary to other theoretical works, e.g. [START_REF] Gorshkov | Photon-Photon Interactions via Rydberg Blockade[END_REF][START_REF] Gorshkov | Dissipative Many-Body Quantum Optics in Rydberg Media[END_REF], here, we are interested in the dispersive regime.

Moreover, since we place the atoms in a cavity rather than in free-space, the theoretical framework and calculations we perform also differ from [START_REF] Gorshkov | Photon-Photon Interactions via Rydberg Blockade[END_REF][START_REF] Gorshkov | Dissipative Many-Body Quantum Optics in Rydberg Media[END_REF]. In particular, a technical benefit of our approach is that we are not restricted to considering only photon pairs but could, in principle, investigate higher-order correlations.

We first write the dynamical equations for the system of interacting three-level atoms coupled to the strong control field and the non-resonant cavity mode, fed by the probe beam. We show that, under some assumptions, the system effectively behaves as a large spin coupled to the cavity mode [START_REF] Guerlin | Cavity quantum electrodynamics with a Rydberg-blocked atomic ensemble[END_REF]. We then compute the steady-state second-order correlation function to characterize the emission of photons out of the cavity. Our numerical simulations suggest that the bunched or antibunched nature of the outgoing light as well as its coherence time may be controlled through adjusting the detuning between the cavity mode and probe field frequencies.

The paper is structured as follows. In Sec. II, we present our setup and the assumptions we make to compute its dynamics. We also explain the analytical and numerical methods we employ to calculate the second-order g (2) correlation function of the outgoing light beam.

In Sec. III, we present and interpret the results of the simulations we obtained for g (2) (0) and g (2) (τ > 0) on the specific experimental case considered in [START_REF] Parigi | Observation and Measurement of Interaction-Induced Dispersive Optical Nonlinearities in an Ensemble of Cold Rydberg Atoms[END_REF]. Finally, we conclude in Sec. IV by evoking open questions and perspectives of our work. Appendices address supplementary technical details which are omitted in the text for readability.

II. MODEL AND METHODS

The system we consider comprises N atoms which present a three-level ladder structure The transitions |g ↔ |e and |e ↔ |r are respectively driven by a weak probe field of frequency ω p and a strong control field of frequency ω cf . To limit absorption, both fields are off-resonant, the respective detunings are given by ∆ e ≡ (ω pω e ) and ∆ r ≡ (ω p + ω cfω r ).

Moreover, to enhance dispersive effects while keeping a high input-output coupling efficiency, the atoms are placed in an optical low-finesse cavity. The transition |g ↔ |e is supposed in the neighbourhood of a cavity resonance. The frequency and annihilation operator of the corresponding mode are denoted by ω c and a, respectively ; the detuning of this mode with the probe laser is defined by ∆ c ≡ (ω pω c ) and α denotes the feeding rate of the cavity mode with the probe field, which is supposed real for simplicity. Finally, we introduce g and Ω cf which are the single-atom coupling constant of the transition |g ↔ |e with the cavity mode and the Rabi frequency of the control field on the transition |e ↔ |r , respectively. In the following paragraphs, we study the dynamics of the system which, under some assumptions, is equivalent to a damped harmonic oscillator, i.e. the cavity mode, coupled to an assembly of spins three-level atoms is neglected, and the leading non-linear effect comes from the Rydberg-Rydberg collisional effects. The system therefore consists of N effective two-level atoms {|g , |r }, with an effective power-broadened dipole decay rate from the Rydberg level

γ r = γ r + Ω 2 cf γ e 4 (∆ 2 e + γ 2 e )
, coupled to the cavity mode of effective decay rate

γ c = γ c + g 2 Nγ e ∆ 2 e + γ 2
e increased by the coupling to the atomic ensemble. The Hamiltonian reads

H = -∆r N n=1 σ (n) rr + N m<n=1 κ mn σ (m) rr σ (n) rr -∆c a † a + α a + a † + g eff a N n=1 σ (n) rg + h.c.
In this expression, we introduced the atomic operators σ At this point, following [START_REF] Guerlin | Cavity quantum electrodynamics with a Rydberg-blocked atomic ensemble[END_REF], we introduce the Rydberg bubble approximation. In this approach, the strong Rydberg interactions are assumed to effectively split the sample into

N b bubbles {B α=1,...,N b } each of which contains n b = N N b
atoms but can only accomodate a single Rydberg excitation, delocalized over the bubble. Note that the number of atoms per bubble n b is approximately given by [START_REF] Parigi | Observation and Measurement of Interaction-Induced Dispersive Optical Nonlinearities in an Ensemble of Cold Rydberg Atoms[END_REF]

n b = 2π 2 ρ at 3 |C 6 | ∆ r -Ω 2 cf /(4∆ e )
where ρ at is the atomic density. Each bubble can therefore be viewed as an effective spin 1 2 whose Hilbert space is spanned by

|-α = |G α ≡ iα∈Bα |g iα |+ α = |R α ≡ 1 √ n b {|rg . . . g + . . . + |g . . . gr }
the ground state of the bubble B α and its symmetric singly Rydberg excited state, respectively. Introducing the bubble spin-1 2 operators s (α)

-= |-α + α | -the operator s (α)
-corresponds to the lowering operator of the spin and the annihilation of a Rydberg excitation, one can write the Hamiltonian under the approximate form (see Appendix A)

H ≈ -∆c a † a + α a + a † -∆r N b 2 + J z +g eff √ n b aJ + + a † J -
where we introduced the collective angular momentum

J -≡ N b α=1 s (α)
-. The system is therefore equivalent to a large spin, i.e. the assembly of spin- 1 2 Rydberg bubbles, coupled to a harmonic oscillator. Its density matrix satisfies the master equation

∂ t ρ = Lρ (II.1) = 1 i H, ρ + γ c 2aρa † -a † aρ -ρa † a + γ r N b α=1 2s (α) -ρs (α) + -s (α) + s (α)
-ρρs

(α) + s (α) -
One can also write the Heisenberg-Langevin equations for the time-dependent operators

a (t) , J -(t) ∂ t a = i ∆c -γc a -iα + ig eff √ n b J -+ ãin (II.2) ∂ t J -= i ∆r -γr J -+ i g eff NN b a + Jin (II.3)
where ãin , Jin ≡ N n=1 F (n) gr are the Langevin forces associated to a and J -, respectively. Note that we neglected the effect of extra dephasing due to, e.g., collisions or laser fluctuations.

To study the quantum properties of the light transmitted through the cavity, we shall compute the function g

(2)
out , which characterizes the two-photon correlations. In the inputoutput formalism [START_REF] Walls | Quantum Optics[END_REF], one shows that this function simply equals the function g (2) for the intra-cavity field (see Appendix B for details) given by

g (2) (τ ) = Tr a † ae Lτ aρ ss a † Tr [a † aρ ss ] 2 (II.4)
where ρ ss denotes the steady state of the system defined by Lρ ss = 0, see Eq. (II.1).

In the regime of small feeding parameter α, one can compute ρ ss numerically by propagating in time the initial state ρ 0 ≡ |N r = 0 N r = 0| ⊗ |n c = 0 n c = 0| (here |N r = 0, 1, . . . , N b represents the symmetric state in which N r ≡ N b 2 + Jz bubbles are excited, and |n c = 0, 1, . . . are the Fock states of the cavity mode). To this end, one applies the Liouvillian evolution operator e Lt in a truncated basis, restricted to states of low numbers of excitations (typically with n c + N r ≤ 6). The steady state is reached in the limit of large times -ideally when t → ∞. The denominator of the ratio Eq.(II.4) is directly obtained from ρ ss . To compute its numerator, one first propagates in time aρ ss a † from t = 0 to τ , using the same procedure as above, then applies the operator a † a and takes the trace. In the regime of weak feeding, it is also possible to get a perturbative expression for g (2) (0) by computing the expansion of a † a † aa ss and a † a ss in powers of α. To this end, one uses the Heisenberg equations of the system Eqs.(II.2,II.3) to derive the hierarchy of equations relating the different mean values and correlations . . . ss up to the fourth order in α. After straightforward though lengthy algebra, one gets an expression for g (2) (0) which is too cumbersome to be reproduced here but allows for faster calculations than the numerical approach. Such a fully analytical treatment, however, cannot, to our knowledge, be extended to g (2) (τ > 0); for τ > 0 we therefore entirely rely on numerical simulations.

To conclude this section, we consider the regime of large number of bubbles and low number of excitations, i.e. N b ≫ 1 and Jz ≪ N b . As shown in Appendix A, the operator b ≡ J - √ N b is then approximately bosonic, and the term N b 2 + Jz can be put under the form

N b 2 + J z ≈ J + J - 2 (N b + 1) + (J + J -) 2 4 (N b + 1) 3 ≈ N b (N b + 1) b † b + N 2 b (N b + 1) 3 b † bb † b ≈ b † b + 1 N b b † b † bb
Finally, we get the following approximate expression for the effective Hamiltonian

H ≈ -∆c a † a + α a + a † -∆r b † b - κ 2 b † b † bb + g eff √ N ab † + a † b
where κ ≡ 2 ∆r /N b . In this regime, the system therefore behaves as two coupled oscillators:

one is harmonic, the cavity field, the other is anharmonic, the Rydberg bubble field.

In the following section, we present and discuss the results we obtained with the specific system used in [START_REF] Parigi | Observation and Measurement of Interaction-Induced Dispersive Optical Nonlinearities in an Ensemble of Cold Rydberg Atoms[END_REF]. It appears that one can choose the bunched or antibunched behaviour of the light transmitted through the cavity by adjusting the detuning ∆ c . We also show that the time behaviour of the function g (2) (τ ) depends on the regime considered, and can be roughly understood as resulting from the damped exchange of a single excitation between atoms and field.

III. NUMERICAL RESULTS AND DISCUSSION

We consider the physical setup presented in [START_REF] Parigi | Observation and Measurement of Interaction-Induced Dispersive Optical Nonlinearities in an Ensemble of Cold Rydberg Atoms[END_REF], i.e. an ensemble of 87 Rb atoms, whose state space is restricted to the levels |g = 5s 1 /γ e . The numerical and theoretical curves are in such a good agreement for the regime considered that the corresponding curves cannot be distinguished. One notes a strong bunching peak (B) θ B = -4.9 and a deep antibunching area centered on (A) θ A = 0. This suggests that around (A), photons are preferably emitted one by one, while around (B) they are preferably emitted by pairs. Note, however, that, as a ratio, g (2) (0) gives only information on the relative importance of pair and single-photon emissions. Its peaks therefore do not correspond to maxima of photon pair emission, but to the best compromises between a † a † aa ss and a † a We now investigate the behaviour of g (2) (τ > 0) for two different values of the detuning, i.e. θ B = -4.9 and θ A = 0 which respectively correspond to the peak (B) and minimum (A) of g (2) (0). The numerical simulations we obtained are given in Fig. III.2. The plot relative to (B) exhibits damped oscillations, alternatively showing a bunched g (2) (τ ) > 1 or antibunched g (2) (τ ) < 1 behaviour. The plot corresponding to (A) always remains on the antibunched side, though asymptotically tending to 1.

The features observed can be understood and satisfactorily accounted for by a simple three-level model. Indeed, due to the weakness of α, the system, in its steady state, is expected to contain at most two excitations (either photonic or atomic). After a photon detection at t = 0, it contains at most one excitation which can be exchanged between the cavity field and atoms, as it has been known for long [START_REF] Brecha | Observation of oscillatory energy exchange in a coupled atom cavity system[END_REF][START_REF] Brune | Quantum Rabi Oscillation: A Direct Test of Field Quantization in a Cavity[END_REF]. In other words, the operator aρ ss a † can be expanded in the space restricted to the three states

{|00 ≡ |N r = 0, n c = 0 , |01 ≡ |N r = 0, n c = 1 , |10 ≡ |N r = 1
, n c = 0 } and the effective non-Hermitian Hamiltonian for the system, in this subspace, takes the following form:

H 3 =      0 α 0 α -∆c -iγ c g eff √ N 0 g eff √ N -∆r -iγ r     
The oscillatory dynamical behaviour observed for g (2) (t) in the specific cases (A,B) is correctly recovered by this Hamiltonian, which validates the schematic model we used and suggests it comprises the main physical processes at work.

To conclude this section, it is worth mentioning that the two-boson approximation, though strictly speaking not applicable here -the parameters considered in this section indeed correspond to a number of bubbles N b ≃ 2, yields, however, the qualitative behaviour for g (2) (0). The minimum is correctly located, though slightly higher than in the spin model; the antibunching peak is slightly shifted towards positive detunings and is weaker than in the previous treatment. These discrepancies result from too low a value of the non-linearity parameter κ ; they can be corrected through replacing κ = 2 ∆/N b by κ′ = 2 ∆/ (N b -1)

in the two-boson Hamiltonian. We first note that κ and κ′ coincide in the regime of large number of bubbles. Moreover, κ′ makes sense in the regime of low number of bubbles: in particular, when N b → 1, i.e. when only one bubble is available, the non-linearity, proportional to κ′ , diverges accordingly, therefore forbidding the boson field to contain more than one excitation. Finally, let us mention that κ′ can also be recovered via a perturbative treatment of the full model which will be presented in a future paper.

IV. CONCLUSION

In this work, we studied how the strong Rydberg-Rydberg van der Waals interactions in an atomic medium may affect the quantum statistical properties of an incoming light beam.

In our model, atoms are located in a low finesse cavity and subject to a weak signal beam and a strong control field. These two fields non-resonantly drive the transition from the ground to a Rydberg level. The system was shown to effectively behave as a large spin coupled to a damped harmonic oscillator, i.e. the assembly of Rydberg bubbles and the cavity mode, respectively. The strong anharmonicity of the atomic spin affects the quantum statistics of the outgoing light beam. To demonstrate this effect, we performed analytical and numerical calculations of the second-order correlation function g (2) (τ ≥ 0). The results we obtained on a specific physical example with rubidium atoms show indeed that the transmitted light presents either bunched or antibunched characters, depending on the detuning between the cavity mode and the probe field. This suggests that, in such a setup, one could design light of arbitrary quantum statistics through appropriately adjusting the physical parameters.

In this work, we performed the Rydberg bubble approximation, which allowed us to derive a tractable effective Hamiltonian. This scheme is, however, questionable: interactions between bubbles are indeed neglected, and the different spatial arrangements of the bubbles in the sample are not considered. Though challenging, it would be interesting to run full simulations of the system, rejecting those states which are too far off-resonant due to Rydberg-Rydberg interactions. Besides validating the assumption of the present work, this would indeed enable us to consider other regimes, such as, for instance, the case of resonant transition towards the Rydberg level. We also implicitly made the assumption that the cavity mode and control beam were homogeneous. Spatial variations should be included in the model and their potential influence studied in a future work. Finally, due to the very weak probe field regime considered in this paper, we only presented results on the function g (2) (τ ): the production of n = 3, 4, . . . correlated photons is indeed very unlikely. In principle, we can, however, numerically compute g (n) (τ ) for any n > 2, which might be relevant in a future work, if addressing stronger probe fields.

Ṽa-c , where

Ha = -∆ e N n=1 σ (n) ee -∆ r N n=1 σ (n) rr + Ω cf 2 N n=1 σ (n) re + σ (n) er + N m<n=1 κ mn σ (m) rr σ (n) rr Hc = -∆ c a † a + α a + a † Ṽa-c = N n=1 g aσ (n) eg + a † σ (n)
ge with the detunings ∆ c ≡ (ω pω c ), ∆ e ≡ (ω pω e ), and ∆ r ≡ (ω p + ω cfω r ).

The corresponding Heisenberg-Langevin equations are:

d dt a = (i∆ c -γ c ) a -iα -ig N i σ (i) ge + a in (A.1) d dt σ (i) ge = (i∆ e -γ e ) σ (i) ge -i Ω cf 2 σ (i) gr + iga σ (i) ee -σ (i) gg + F (i) ge (A.2) d dt σ (i) gr = (i∆ r -γ r ) σ (i) gr -i Ω cf 2 σ (i) ge + igaσ (i) er (A.3) -iσ (i) gr N j =i κ ij σ (j) rr + F (i) gr d dt σ (i) er = {i (∆ r -∆ e ) -γ er } σ (i) er + i Ω cf 2 σ (i) rr -σ (i) ee (A.4) +iga † σ (i) gr -iσ (i) er N j =i κ ij σ (j) rr + F (i)
er where a in and

F (i)
αβ denote Langevin forces.

Elimination of the intermediate state

Let us now simplify the system. First, one deduces from Eq.(A.4) that σ er is of second order in the small feeding constant α. The term aσ

(i)
er can therefore be neglected in Eq.(A.3). Moreover, since the ground state population remains dominant during the evolution of the system we can write σ (i) eeσ (i) gg ≃ -I ; from Eq.(A.2), the steady-state solution for σ (i) ge in the far detuned regime is therefore

σ (i) ge ≃ Ω cf 2 (∆ e +
F (i) gr = F (i) gr + Ω cf 2 (∆ e + iγ e ) F (i) ge ≈ F (i) gr + Ω cf 2∆ e F (i)
ge Note that, in the absence of collisional terms, one simply recovers the standard three-level EIT susceptibility in the far-detuned regime

da dt =   i∆ c -γ c - g 2 N Ω 2 cf 4(γr -i∆r) -i∆ e   a -iα + ãin
Finally, we get the effective Hamiltonian

H = -∆r N n=1 σ (n) rr + N m<n=1 κ mn σ (m) rr σ (n) rr -∆c a † a + α a + a † + g eff a N n=1 σ (n) rg + h.c.

Rybderg bubble approximation

As described in the main text, we introduce the Rydberg bubble approximation. In this approach, the strong Rydberg interactions are assumed to effectively split the sample into

n b = 2π 2 ρ at 3 |C 6 | ∆ r -Ω 2 cf /(4∆ e )
where ρ at is the atomic density. Each bubble can therefore be viewed as an effective spin 1 2 whose Hilbert space is spanned by

|-α = |G α ≡ iα∈Bα |g iα |+ α = |R α ≡ 1 √
-= |-α + α | -the operator s (α)
-corresponds to the lowering operator of the spin and the annihilation of a Rydberg excitation, one can write

N n=1 σ (n) gr = N b α=1 iα∈Bα σ (iα) gr ≈ N b α=1 s (α) - -α iα∈Bα σ (iα) gr + α ≈ √ n b N b α=1 s (α) - = √ n b J -
where we introduced the collective angular momentum J -≡ N b α=1 s (α)

-. In the same way,

N n=1 σ (n) rr = N b α=1 iα∈Bα σ (iα) rr ≈ N b α=1 |+ α + α | + α iα∈Bα σ (iα) rr + α ≈ N b α=1 1 2 + s (α) z ≈ N b 2 + J z
where we used

|+ α + α | ≡ 1 2 + s (α) z .
Finally, the Hamiltonian of the system takes the approximate form

H ≈ -∆c a † a + α a + a † -∆r N b 2 + J z +g eff N b aJ + + a † J -
which represents the interaction of the large spin J -with the cavity mode a.

Regime of large number of bubbles and low number of excitations

From the well-known relation J + J -= J 2 -J 2 z + J z we deduce the second-order operator equation

J 2 z -J z -2 N b 2 N b + 2 2 + J + J -= 0
In the regime of large number of bubbles N b ≫ 1 and for low excitation numbers, i.e. eigenstates of the total angular momentum j = N b 2 ; m = -N b 2 + k with k ≪ N b , the solution of this equation is approximately given by

J z ≈ - N b 2 + J + J - 2 (N b + 1) + (J + J -) 2 4 (N b + 1) 3
whence, at the lowest order in the excitation number,

N b 2 + J z ≈ J + J - 2 (N b + 1) + (J + J -) 2 4 (N b + 1) 3 (A.7) [J + , J -] ≈ -2 N b (A.8)
Injecting Eq.(A.7) into the previous form of the Hamiltonian we get

H ≈ -∆c a † a + α a + a † -∆r J + J - 2 (N b + 1) + (J + J -) 2 4 (N b + 1) 3 + g eff √ N a J + √ N b + a † J - √ N b
Moreover, from Eq.(A.8) we deduce that the operator b ≡ J - √ N b is approximately bosonic and therefore, the Hamiltonian can finally be put under the form

H ≈ -∆c a † a + α a + a † -∆r b † b - κ 2 b † b † bb + g eff √ N ab † + a † b
where κ ≡ 2 ∆r /N b .

Appendix B: Calculation of g

(2) out

By definition, the second-order correlation function for the outgoing field is g

out (t Here we used the relation

[X (t 1 ) , a in (t 2 )] = 2γ c θ (t 1 -t 2 ) [X, a]
where X is any system operator [START_REF] Walls | Quantum Optics[END_REF] and where θ (τ ) is the Heaviside step-function (with θ (0) =

  with a ground |g , intermediate |e and Rydberg states |r (see Fig. II.1). The energy of the atomic level |k = g, e, r is denoted by ω k (by convention ω g = 0) and the dipole decay rates from the intermediate and Rydberg states are denoted by γ e and γ r , respectively.

Figure II. 1 :

 1 Figure II.1: a) The setup consists of N cold atoms placed in an optical cavity which is fed by a weak (classical) laser beam of frequency ω p and a strong control laser field of frequency ω cf . b) The atoms present a three-level ladder structure {|g , |e , |r }. The transitions |g ↔ |e and |e ↔ |r are non-resonantly driven by the injected probe and control laser fields, respectively, with the respective coupling strength and Rabi frequency g and Ω cf (see the text for the definitions of the different detunings represented here).

2 e + γ 2 e

 22 kl ≡ I(1) ⊗ . . . ⊗ I (n-1) ⊗ |k l| ⊗ I (n+1) ⊗ . . . ⊗ I (N ) for (k, l) = g, e, r as well as the effective detunings∆r ≡ ∆ r -Ω 2 cf ∆ e 4 (∆ 2 e + γ 2 e ) and ∆c ≡ ∆ c -g 2 N∆ e ∆respectively shifted from ∆ r and ∆ c by the AC Stark shift of the control beam and by the linear atomic susceptibility. The quantity κ mn ≡ C 6 / r mr n 6 is the van der Waals interaction between atoms (m, n) in their Rydberg level -when atoms are in the ground or intermediate states, their interactions are neglected, whileg eff = gΩ cf 2∆ eis the effective coupling strength of the two-photon transition |g → |r driven by the cavity mode and the control laser.

2 ; F = 2 , |e = 5p 3 2 ; F = 3 and |r = 95d 5 2 ;

 2222 F = 4 with the decay rates γ e = 2π × 3 MHz, and γ r = 2π × 0.03 MHz. The other physical parameters must be designed so that strong non-linearities may be observed at the single-photon level. In the specific system considered here, we find this is achieved for a cavity decay rate γ c = 2π × 1 MHz, a volume of the sample V = 40π × 15 × 15µm 3 , a sample density n at = 0.4µm -3 , a control laser Rabi frequency Ω cf = 10γ e , a cooperativity C = 1000, a detuning of the intermediate level ∆ e = -35γ e , a detuning of the Rydberg level ∆ r = 0.4γ e , a cavity feeding rate α = 0.01γ e . For these parameters, the cavity detuning ∆ (0) c = -6.1γ e corresponds to the maximal average number of photons in the cavity. Note that these physical parameters are experimentally realistic and feasible.Let us first focus on the second-order correlation function at zero time g(2) (0), represented on Fig. III.1 a) as a function of the reduced detuning θ ≡ ∆ c -∆ (0) c

2 ss

 2 , as can be checked by comparison of Fig. III.1 a) and b). Hence, pair emission might dominate in a regime where the number of photons coming out from the cavity is actually very small.

N

  b bubbles {B α=1,...,N b } each of which contains n b = N N batoms but can only accomodate a single Rydberg excitation, delocalized over the bubble. Note that the number of atoms per bubble n b is approximately given by[START_REF] Parigi | Observation and Measurement of Interaction-Induced Dispersive Optical Nonlinearities in an Ensemble of Cold Rydberg Atoms[END_REF]

Figure 2 ss = n 2 Figure III. 2 :

 222 Figure III.1: a) Second-order correlation function at zero time g (2) (0) (numerical and analytical plots coincide), as a function of the reduced detuning θ ≡ ∆ c -∆ (0) c/γ e . In the neighbourhood of the minimum (A) θ A = 0, a strong antibunching region is observed (see inset); a strong bunching area is obtained around the peak (B) θ B = -4.9. b) Average number of pairs a † a † aa ss = n (n -1) ss (thin line) and square of the average number of photons a † a 2 ss = n 2 ss in the steady state (thick line). The position of the peak of the correlation function g(2) (0) is signaled by a vertical line.

  n b {|rg . . . g + . . . + |g . . . gr } the ground state of the bubble B α and its symmetric singly Rydberg excited state, respectively. Introducing the bubble Pauli operators s

	(α)

  and keeping only non-zero terms (all terms like a † in .... and ....a in equal zero), one obtains in the numerator four non-zero termsa † (t 1 )a † (t 2 )a in (t 2 )a(t 1 ) a † (t 1 )a † (t 2 )a(t 2 )a(t 1 ) a † (t 1 )a † in (t 2 )a(t 2 )a(t 1 ) a † (t 1 )a † in (t 2 )a in (t 2 )a(t 1 )Let us consider the first term. Using the standard commutation relations between a and a in operators we have:a † (t 1 )a † (t 2 )a in (t 2 )a(t 1 ) = a † (t 1 )a † (t 2 )a(t 1 )a in (t 2 ) + a † (t 1 )a † (t 2 ) [a in (t 2 ), a(t 1 )] = 2γ c θ(t 1t 2 ) a † (t 1 )a † (t 2 ) [a(t 2 ), a(t 1 )]

1 , t 2 ) = a † out (t 1 ) a † out (t 2 ) a out (t 2 ) a out (t 1 ) a † out (t 2 ) a out (t 2 ) a † out (t 1 ) a out (t 1 ) Using the relations [14] a † out (t) a out (t) = 2γ c a † (t) a (t) a out (t) = 2γ c a (t)a in (t)

  1 2 ). Evaluating the other terms in the same way one finally obtainsg (2) out (t 1 , t 2 ) = a † (t m ) a † (t M ) a (t M ) a (t m ) a † (t 1 ) a (t 1 ) a † (t 2 ) a (t 2 )where t m ≡ min (t 1 , t 2 ) and t M ≡ max (t 1 , t 2 ).

Acknowledgments

This work was supported by the EU through the ERC Advanced Grant 246669 DELPHI and the Collaborative Project 600645 SIQS.

Appendix A: Derivation of the effective Hamiltonian

Rotating Wave Approximation

The full Hamiltonian of the system can be written under the form

ge where σ

(n) ) , ω α is the energy of the atomic level |α for α = e, r (with the convention ω g = 0), and κ mn ≡ C 6 rm-rn 6 denotes the van der Waals interaction between atoms in the Rydberg level -when atoms are in the ground or intermediate states, their interactions are neglected.

We switch to the rotating frame defined by |ψ → ψ = exp -it H 0 where

rr and perform the Rotating Wave Approximation to get the new Hamiltonian H = Ha + Hc +