Gaël Poux-Médard
email: gael.poux-medard@univ-lyon2.fr

Julien Velcin
email: julien.velcin@univ-lyon2.fr

Sabine Loudcher
email: sabine.loudcher@univ-lyon2.fr

Dynamic Mixed Membership Stochastic Block Model for Weighted Labeled Networks

Keywords: Information systems → Clustering, Social networks, Social recommendation, Data mining, Collaborative filtering, Network data models, • Mathematics of computing → Probabilistic inference problems, Bayesian networks Block models, Dynamic networks, Clustering, Recommender systems

Most real-world networks evolve over time. Existing literature proposes models for dynamic networks that are either unlabeled or assumed to have a single membership structure. On the other hand, a new family of Mixed Membership Stochastic Block Models (MMSBM) allows to model static labeled networks under the assumption of mixed-membership clustering. In this work, we propose to extend this later class of models to infer dynamic labeled networks under a mixed membership assumption. Our approach takes the form of a temporal prior on the model's parameters. It relies on the single assumption that dynamics are not abrupt. We show that our method significantly differs from existing approaches, and allows to model more complex systems -dynamic labeled networks. We demonstrate the robustness of our method with several experiments on both synthetic and real-world datasets. A key interest of our approach is that it needs very few training data to yield good results. The performance gain under challenging conditions broadens the variety of possible applications of automated learning tools -as in social sciences, which comprise many fields where small datasets are a major obstacle to the introduction of machine learning methods.

 (top), or listen to either of these without any specific pattern (bottom). For 200 epochs containing only 5 observations each, our approach (in red) infers any smooth dynamic membership pattern, and does it more accurately than static models (in green [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF]) and models that consider each time slice independently (in blue [START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF]).

INTRODUCTION

Dynamic networks are powerful tools to visualize and model interactions between different entities that can evolve over time. The network's nodes represent the interacting entities, and ties between these nodes represent an interaction. In many real-world situations, ties strength can vary over time -on music streaming websites for instance, users' affinity with various musical genres vary over time [START_REF] Kumar | Predicting Dynamic Embedding Trajectory in Temporal Interaction Networks[END_REF][START_REF] Villermet | Follow the Guides: Disentangling Human and Algorithmic Curation in Online Music Consumption[END_REF] -see Fig. 1 as an illustration. The network is said dynamic. Now, every interaction does not have the same significance. A music listener might like both Rock and Jazz, but favor one over the other. This user's ties to musical genres do not have the same intensity; each tie is associated to a number, representing the strength of the interaction. The network is said to be dynamic and weighted.

Finally, simple dynamic weighted networks might not be fit to grasp the complexity of a given situation. A music listener may have different opinions on musical genres; they can like it, dislike it, be bored of it, prefer to listen to some only in the morning, at night, etc. Each of these relations can be represented by their own tie in the network, each associated to their own strength. The network is said to be dynamic, weighted and labeled. Such networks are high-dimensional objects, whose direct inference is a difficult problem. Several ways to achieve this task in a static setting have been proposed; the Stochastic Block Models (SBM) family is one of the most popular approaches [START_REF] Cobo-López | Optimal prediction of decisions and model selection in social dilemmas using block models[END_REF][START_REF] Guimerà | A Network Inference Method for Large-Scale Unsupervised Identification of Novel Drug-Drug Interactions[END_REF][START_REF] Holland | Stochastic blockmodels: First steps[END_REF]. The underlying assumption is that some sets of nodes behave in a similar way regarding other sets of nodes. They can be grouped in a cluster; instead of modeling every edge for every node in a network, only the edges between clusters are modeled, making the inference task much more tractable. Each cluster is associated to labeled edges, and each node is associated to a cluster. A variant of SBM that allows more expressive power is called the Mixed-Membership SBM (MMSBM), where each node can belong to several clusters with different proportions [START_REF] Edoardo | Mixed Membership Stochastic Blockmodels[END_REF][START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF][START_REF] Wang | Stochastic Blockmodels for Directed Graphs[END_REF]. The inferred network is then weighted and labeled, but not dynamic. A major advantage of the SBM family is that they yield readily interpretable results (interpretable clusters), unlike most existing neural-network approaches for labeled networks [START_REF] Fan | On Interpretability of Artificial Neural Networks: A Survey[END_REF].

The goal of this paper is to infer networks that are dynamic, weighted and labeled by using a mixed-membership SBM. Our contribution consists in the extension of a broad class of MMSBM by using a temporal Dirichlet prior. We first conduct a careful review of dynamic network inference literature, and we show that no prior work addresses the task at hand. Although some previous works tackle similar problems, they are not fit for inferring dynamic, weighted and labeled networks. Besides, our approach is conceptually much simpler than those proposed in the literature. Last but not least, our dynamic extension is readily plugable into most existing MMSBM for labeled and weighted networks.

We develop an EM optimization algorithm that scales linearly with the size of the dataset, demonstrate the effectiveness of our approach on both synthetic and real-world datasets, and further detail a possible application scenario.

BACKGROUND 2.1 Notations

We consider a network of 𝐼 nodes and 𝑂 labels. All the clustering models discussed in this section can be represented by a restricted set of parameters: a matrix 𝜃 (𝑡) ∈ R 𝐼 ×𝐾 accounts for the membership of each of 𝐼 nodes to each of 𝐾 possible clusters, and a block-interaction matrix 𝑝 (𝑡) ∈ R 𝐾 ×𝐾 ×𝑂 represents the probability for each of 𝐾 clusters to yield each of 𝑂 labels. Both 𝜃 (𝑡) and 𝑝 (𝑡) can vary over time; in the following, this temporal dependence is implicit (𝑥 (𝑡) := 𝑥) unless specified otherwise. The network is said to be unlabeled when 𝑂 = 1, and binary (as opposed to weighted) if an edge can only exist or not exist.

Dynamic unlabeled networks -Single-membership

Single-membership SBMs consider a membership matrix such as 𝜃 ∈ {0; 1} 𝐼 ×𝐾 : each membership vector equals 1 for one cluster, and 0 everywhere else -"hard" clustering in the literature. In [START_REF] Xu | Dynamic Stochastic Blockmodels: Statistical Models for Time-Evolving Networks[END_REF][START_REF] Xu | Dynamic Stochastic Blockmodels for Time-Evolving Social Networks[END_REF], the authors proposed to model a binary unlabeled dynamic network using a label-switching inference method, optimized by a Sequential Monte Carlo algorithm [START_REF] Jin | A Survey of Community Detection Approaches: From Statistical Modeling to Deep Learning[END_REF]. Both the membership and the interaction matrices can vary over time, thus supposing two independent underlying Markov processes [START_REF] Jin | A Survey of Community Detection Approaches: From Statistical Modeling to Deep Learning[END_REF].

In [START_REF] Tang | Detecting Social Media Hidden Communities Using Dynamic Stochastic Blockmodel with Temporal Dirichlet Process[END_REF][START_REF] Yang | Detecting communities and their evolutions in dynamic social networks -a Bayesian approach[END_REF], the authors propose to model a binary dynamic and unlabeled network. The cluster interaction matrix 𝑝 can vary over time while keeping the memberships 𝜃 static. The entries of 𝑝 are drawn from a Dirichlet distribution and expressed as a Chinese Restaurant Process. This process converges to a Dirichlet distribution, and allows to infer a potentially infinite number of clusters. This model is therefore non-parametric and inferred using an MCMC algorithm.

A conceptually novel way of modeling dynamic unlabeled networks under the single membership assumption has been proposed in [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF][START_REF] Matias | A semiparametric extension of the stochastic block model for longitudinal networks[END_REF]. The authors propose to model the cluster interaction and membership matrices as Poisson processes, that explicitly represent the continuous temporal dependency without slicing the dataset into episodes. The method allows to infer varying membership and interaction matrices for dynamic binary or weighted networks, but their results have shown that allowing both to vary simultaneously leads to identifiability and label switching issues [START_REF] Funke | Stochastic block models: A comparison of variants and inference methods[END_REF]. This conclusion seems reasonable, given none of these SBM algorithms can reach a global optimum of the likelihood function. A model where both 𝜃 and 𝑝 can vary over time is more likely to get stuck into a local optimum.

Finally, we mention the existence of SBM variants that account for dynamic degree-correction [START_REF] Wilson | Modeling and detecting change in temporal networks via the degree corrected stochastic block model[END_REF] or that enforce a scale-free characteristic [START_REF] Wu | Dynamic Stochastic Block Model with Scale-Free Characteristic for Temporal Complex Networks[END_REF].

All these methods consider unlabeled networks, and consider a hard clustering which does not allow for as much expressive power as the Mixed-Membership approaches.

Dynamic unlabeled networks -Mixed-membership

Mixed-membership SBMs consider a membership matrix such as 𝜃 ∈ R 𝐼 ×𝐾 , where each membership vector 𝜃 𝑖 has a 𝐿 1 norm of 1.

Literature also refer to it as "soft" clustering. Similar to [START_REF] Tang | Detecting Social Media Hidden Communities Using Dynamic Stochastic Blockmodel with Temporal Dirichlet Process[END_REF][START_REF] Yang | Detecting communities and their evolutions in dynamic social networks -a Bayesian approach[END_REF], a method for inferring dynamical binary unlabeled networks has been proposed in [START_REF] Fan | Dynamic Infinite Mixed-Membership Stochastic Blockmodel[END_REF]. The membership vector of each piece of information is drawn from a Chinese Restaurant Process (CRP) according to the number of times a node has already been associated to each cluster before. The process yields a distribution over an infinity of available clusters. The formulation as a CRP arises naturally, because the prior on membership vectors is typically a Dirichlet distribution -a CRP naturally converges to a draw from a Dirichlet distribution [START_REF] Richard Arratia | Poisson Process Approximations for the Ewens Sampling Formula[END_REF]. The block-interaction matrix 𝑝 does not vary over time. The article shows a complexity analysis that suggest the methods runs with a complexity of O (𝑁 2 ) which makes it unfit for large-scale real-world applications.

The work the most closely related to ours is [START_REF] Xing | A state-space mixed membership blockmodel for dynamic network tomography[END_REF]. This seminal work proposed the dMMSB as a way to model dynamic binary unlabeled networks using a variational algorithm [START_REF] Lee | A review of stochastic block models and extensions for graph clustering[END_REF]. To do so, the authors modify the original MMSBM [START_REF] Edoardo | Mixed Membership Stochastic Blockmodels[END_REF] to consider a logistic normal distribution as a prior on the membership vectors 𝜃 𝑖 . This choice allows to model correlations between membership vectors' evolution [START_REF] Ahmed | Seeking The Truly Correlated Topic Posterior -on tight approximate inference of logistic-normal admixture model[END_REF]. The membership vectors are then embedded in a state space model, that is a space where one can define a linear transition between two time points for a given variable. The authors define such trajectory for the membership vectors as a linear function of the previous time point. The trajectory is estimated and smoothed using a Kalman Filter. This approach is the most closely related to ours, as it models the temporal dependency using a prior distribution over memberships, noted 𝑃 (𝜃 ). [START_REF] Xing | A state-space mixed membership blockmodel for dynamic network tomography[END_REF] has been extended to consider 𝑃 (𝜃 ) as a logistic normal mixture prior [START_REF] Ho | Evolving Cluster Mixed-Membership Blockmodel for Time-Evolving Networks[END_REF], which improves its expressive power.

However, this model and related extensions are not fit for the task at hand. First, it considers unlabeled and binary networks [START_REF] Lee | A review of stochastic block models and extensions for graph clustering[END_REF], and extension to labeled and weighted networks, if possible, is not trivial. The proposed optimization algorithm requires a loop until convergence at each EM iteration, making it unable to handle large datasets. Besides, the clusters interaction matrix 𝑝 must remain static for the approach to work, which we alleviate here. And most importantly, is assumes a linear transition between time steps in the state space, while we do not assume any kernel function in our proposed approach.

Static labeled networks -Mixed-membership

Recent years saw a rise of Bayesian methods for inferring static weighted labeled networks using MMSBM variants [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Complex decision-making strategies in a stock market experiment explained as the combination of few simple strategies[END_REF][START_REF] Poux-Médard | Information Interactions in Outcome Prediction: Quantification and Interpretation Using Stochastic Block Models[END_REF][START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF].

Here, nodes are associated a type 𝑎; each type of node has its own layer in a multipartite network, and its associated set of available clusters through a membership matrix 𝜃 𝑎 . Nodes from one layer belong to every cluster corresponding to this layer with mixed proportions. The interaction between clusters of all layers then yields a distribution over 𝑂 possible output labels. This results in a static MMSBM for labeled networks. Our contribution extends all of [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Complex decision-making strategies in a stock market experiment explained as the combination of few simple strategies[END_REF][START_REF] Poux-Médard | Information Interactions in Outcome Prediction: Quantification and Interpretation Using Stochastic Block Models[END_REF][START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF] to the dynamic case. These works will be further detailed in Section 3. Note that in [START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF], the authors consider a temporal slicing of the data and consider each slice as independent from the others; a time slice is considered as a node in a tripartite network. We will compare our approach to this modeling later.

We use these works as a base model, that we couple to a dynamic prior distribution on parameters. Our work focuses on making the prior probability of both 𝜃 and 𝑝 time-dependent, in order to model these parameters dynamics. We provide a ready-to-use temporal plug-in for each of the works presented in this section. It applies to dynamical, weighted and labeled networks in a mixed-membership context, and inference is conducted with a scalable variational EM algorithm.

DYNAMIC LABELED MMSBM 3.1 Base model

For clarity, we choose to consider the simplest form of a labeled MMSBM: each of 𝐼 nodes is associated to each of 𝐾 clusters in mixed proportions, and each of 𝐾 clusters is in turn associated to 𝑂 labels (SIMSBM(1) in [START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF]). As shown in [START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF], our demonstration on SIMSBM(1) trivially extends to [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF] (SIMSBM(1,1)) [START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF] (SIMSBM(2,1)) and [START_REF] Poux-Médard | Information Interactions in Outcome Prediction: Quantification and Interpretation Using Stochastic Block Models[END_REF] (SIMSBM(2)) -and in general any model formulated as SIMSBM(x,y,...).

We consider a set of 𝐼 nodes that can be associated to 𝑂 possible labels on a discrete time interval, or epoch, written 𝑡. We assume that, at each time step, each of 𝐼 nodes belongs to a mixture of 𝐾 available clusters, each of which are in turn yield a probability distribution over 𝑂 labels. The membership of each of 𝐼 nodes to each of the 𝐾 possible clusters at time 𝑡 is encoded in the membership matrix 𝜃 (𝑡 ) ∈ R 𝐼 ×𝐾 . One vector 𝜃 (𝑡 ) 𝑖 represents the probability that 𝑖 belongs to any of the 𝐾 clusters at time 𝑡, and is normalized as:

∑ 𝑘 ∈𝐾 𝜃 (𝑡 ) 𝑖,𝑘 = 1 ∀𝑖, 𝑡 (1)
The probability for each of 𝐾 clusters to be associated to each of 𝑂 labels at time 𝑡 is encoded in the matrix 𝑝 (𝑡 ) ∈ R 𝐾 ×𝑂 . An entry 𝑝 (𝑡 ) 𝑘 (𝑜) represents the probability that cluster 𝑘 is associated to label 𝑜 at time 𝑡, and thus is normalized as:

∑ 𝑜 ∈𝑂 𝑝 (𝑡 ) 𝑘 (𝑜) = 1 ∀𝑘, 𝑡 (2)
Finally, the probability that a node 𝑖 is associated to label 𝑜 at time 𝑡 (i.e. the probability of an edge between 𝑖 and 𝑜 at time 𝑡) is written:

𝑃 (𝑖 → 𝑜 |𝑡) = ∑ 𝑘 ∈𝐾 𝜃 (𝑡 ) 𝑖,𝑘 𝑝 (𝑡 ) 𝑘 (𝑜) (3)
Given a set 𝑅 • of observed triplets (𝑖, 𝑜, 𝑡), the model's posterior distribution can be written [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF]:

𝑃 (𝜃, 𝑝 |𝑅 • ) ∝ 𝑃 (𝑅 • |𝜃, 𝑝) 𝑡 𝑃 (𝜃 (𝑡 ) )𝑃 (𝑝 (𝑡 ) ) (4) = (𝑖,𝑜,𝑡 ) ∈𝑅 • ∑ 𝑘 ∈𝐾 𝜃 (𝑡 ) 𝑖,𝑘 𝑝 (𝑡 ) 𝑘 (𝑜) 𝑡 𝑖 𝑃 (𝜃 (𝑡 ) 𝑖 ) 𝑘 𝑃 (𝑝 (𝑡 ) 𝑘 )
Now, before we describe the optimization procedure, we must choose the priors 𝑃 (𝜃 (𝑡 ) ) and 𝑃 (𝑝 (𝑡 ) ).

Simple Dynamic SBM (SDSBM) prior

We formulate the prior distribution over 𝜃 (𝑡 ) and 𝑝 (𝑡 ) under a single assumption: the parameters at a given time do not vary abruptly at small time scales. It means an entry 𝜃

(𝑡 1 )
𝑖,𝑘 should not differ significantly from 𝜃

(𝑡 2 )
𝑖,𝑘 for every 𝑡 2 close enough to 𝑡 1 . The entries close to a reference time are called temporal neighbours.

Our a priori knowledge on each entry 𝜃 (𝑡 ) 𝑖 and 𝑝

(𝑡 )

𝑘 is that they should not differ significantly from their temporal neighbours. This is a fundamental difference with [START_REF] Xing | A state-space mixed membership blockmodel for dynamic network tomography[END_REF], where the next parameters values are estimated using from a Kalman Filter that only considers the previous time step. Moreover, the authors assume a linear transition function, while we do not make such hypothesis. An illustration of the proposed approach is given Fig. 2, where the prior probability of a membership over time depends on its temporal neighbours.

Dirichlet distribution. Since each entry 𝜃 (𝑡 ) 𝑖 and 𝑝 (𝑡 )

𝑘 is normalized to 1, we consider a Dirichlet distribution as a prior, which naturally yields normalized vectors such that 𝑛 𝑥 𝑛 = 1. It reads:

𝐷𝑖𝑟 (𝑥 |𝛼) = 1 𝐵(𝛼) 𝑛 𝑥 𝛼 𝑛 -1 𝑛 (5)
where 𝐵(•) is the multivariate beta function. In Eq. 5, the vector 𝛼 is called the concentration parameter and must be provided to the model. This parameter defines the mode and the variance of the Dirichlet distribution.

We consider a concentration parameter as 𝛼 = 1 + 𝛽𝛼 0 , so that when 𝛽 = 0 we recover a uniform prior over the simplex, and 𝛽 ≥ 0 |𝑡 -𝑡 ′ | . When inferring a parameter 𝑥 (𝑡 ) at a time 𝑡, the variance of its prior probability 𝑃 (𝑥 (𝑡 ) ) depends on 𝑡 relative to the temporal neighbours. Here for instance, the variance is null at 𝑡 = 2 because 𝜅 (𝑡, 𝑡 ′ ) diverges, and so does 𝛼 (𝑡 ) , hence the variance collapsing to 0. so that the prior has a unique mode. The most frequent value drawn from Eq.5 (or mode) is

M(𝑥 𝑛 ) = 𝛼 𝑛 -1 ′ 𝑛 (𝛼 ′ 𝑛 -1) = 𝛼 0,𝑛 ′ 𝑛 𝛼 0,𝑛
. We recover a uniform prior for 𝛽 = 0; the variance vanishes with 𝛽 ≫ 1 as 1 𝛽 . The effect of various values of 𝛽 on the prior distribution is illustrated in Fig. 3. 𝑘 do not vary abruptly over small time scales. To do so, we define their prior probability's mode with respect to close temporal neighbours. The hyper-parameter 𝛽 controls the variance of the prior -that is, how much our constraint should impact the inference procedure. We express the Simple Dynamic prior parameters for 𝜃 (𝑡 ) 𝑖 as:

Simple

𝛼 (𝑡,𝜃 ) 𝑖,𝑘 = 1 + 𝛽 𝑡 ′ ≠𝑡 𝜅 (𝑡, 𝑡 ′ )𝜃 (𝑡 ′ ) 𝑖,𝑘 𝑡 ′ ≠𝑡 𝜅 (𝑡, 𝑡 ′ ) ⟨𝜃 (𝑡 ) 𝑖,𝑘 ⟩ (6)
where 𝜅 (𝑡, 𝑡 ′ ) is a weight function, and 𝛼 (𝑡,𝜃 ) corresponds to the concentration parameter for 𝜃 at time 𝑡. In following experiments, we define the weight function as 𝜅 (𝑡, 𝑡 ′ ) =

𝑁 𝑡 ′ |𝑡 -𝑡 ′ |
, where 𝑁 𝑡 ′ is the number of observations made at time 𝑡 ′ . This way, temporal neighbours' influence decrease as the inverse of temporal distance. We illustrate the influence of this particular kernel function on the prior probability on membership at all times in Fig. 2. In particular we see that with this expression, the prior probability variance goes to 0 when the considered time it very close to a temporal neighbour.

𝑃 (𝜃 (𝑡 ) 𝑖 |{𝜃 (𝑡 ′ ) 𝑖,𝑘 } 𝑡 ′ ≠𝑡 ) ∝ 𝑘 𝜃 (𝑡 ) 𝑖,𝑘 𝛽 ⟨𝜃 (𝑡 ) 𝑖,𝑘 ⟩ (7) 𝑃 (𝑝 (𝑡 ) 𝑘 (𝑜)|{𝑝 (𝑡 ′ ) 𝑘 (𝑜)} 𝑡 ′ ≠𝑡 ) ∝ 𝑜 𝑝 (𝑡 ) 𝑘 (𝑜) 𝛽 ⟨𝑝 (𝑡 ) 𝑘 (𝑜 ) ⟩
We omitted the normalisation factor for clarity -it does not influence the inference procedure.

Inference

3.3.1 E step. We develop an EM inference procedure for maximizing the log-posterior distribution defined Eq.8. The expectation step computes the expected probability of a latent variable (here a cluster 𝑘) being chosen given each entry of 𝑅 • . Since such latent variables do not appear in the priors expressions, the expectation step remains unchanged by the introduction of the Simple Dynamic Priors; in general, prior distributions do not intervene in the computation of the expectation step [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. The E step for such labeled networks has already been discussed on several occasions [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Information Interactions in Outcome Prediction: Quantification and Interpretation Using Stochastic Block Models[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF] In Eq.8, we introduced a proposal distribution 𝜔

(𝑡 )
𝑖,𝑜 (𝑘), that represents the probability of one cluster allocation 𝑘 given the observation (𝑖, 𝑜, 𝑡). The last line followed from Jensen's inequality assuming that 𝑘 𝜔 (𝑡 ) 𝑖,𝑜 (𝑘) = 1. We notice that Jensen's inequality holds as an equality when:

𝜔 (𝑡 ) 𝑖,𝑜 (𝑘) = 𝜃 (𝑡 ) 𝑖,𝑘 𝑝 (𝑡 ) 𝑘 (𝑜) 𝑘 ′ 𝜃 (𝑡 ) 𝑖,𝑘 ′ 𝑝 (𝑡 ) 𝑘 ′ (𝑜) (9)
which provides us with the expectation of the latent variable 𝑘 given an observation (𝑖, 𝑜, 𝑡) ∈ 𝑅 [START_REF] Ho | Evolving Cluster Mixed-Membership Blockmodel for Time-Evolving Networks[END_REF] where 𝜕(𝑖, 𝑡) = {𝑜 |(𝑖, •, 𝑡) ∈ 𝑅 • } is the subset of labels associated to both 𝑖 and 𝑡, and 𝑁 𝑖,𝑡 = |𝜕(𝑖, 𝑡)| is the size of this set. Note that for 𝛽 = 0 we recover the M-step of standard static MMSBM models [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Information Interactions in Outcome Prediction: Quantification and Interpretation Using Stochastic Block Models[END_REF][START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF].

The derivation of the M-step for the entries 𝑝

Discussion

3.4.1 Plug-in for existing models. Using SDSBM (for Simple Dynamic SBM) as a temporal extension in existing MMSBM models requires very few changes. In [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Information Interactions in Outcome Prediction: Quantification and Interpretation Using Stochastic Block Models[END_REF][START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF], its introduction boils down to adding a term 𝛽 ⟨𝑥⟩ to the numerator of maximization equations, and the corresponding normalizing term 𝛽 to the denominator. This way, our approach is ready-to-use to make these models for for modelling dynamic networks.

Flexible dynamic modeling.

The prior allows to consider that some parameters are dynamic and that others are not. For instance, when several membership matrices are involved, as in [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF][START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF][START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF]), setting 𝛽 = 0 for some makes them time-invariant (or universal). The Simple Dynamic prior also allows to choose whether the blockinteraction tensor 𝑝 is dynamic. Moreover, 𝛽 does not have to be the same for every membership matrix, or even every entry 𝑖 of each of them. Finally, 𝛽 itself can vary over time. To summarize, 𝛽 allows to control the temporal scale over which parameters may vary. This allows to jointly model universal parameters (𝛽 = 0) and dynamical ones (𝛽 ≠ 0).

Tunable temporal dependence.

Finally, the choice of the averaging kernel function 𝜅 (𝑡, 𝑡 ′ ) is important. It allows to choose the range over which the inference of a variable should rely on its temporal neighbours. A formulation as the inverse of time difference seems relevant: the weight of a neighbour appearing at a time 𝛿𝑡 later should diverge as 𝛿𝑡 → 0, so that continuity is ensured. Besides, one could control the smoothness of the curve with respect to time by tuning the weight function as 𝜅 (𝑡, 𝑡 ′ ) =

EXPERIMENTS 4.1 Synthetic data

In this section, we develop several situations in which our method (abbreviated SDSBM for Simple Dynamic MMSBM) could be useful 1 . Experiments are run for 𝐼 = 100, 𝐾 = 3 and 𝑂 = 3, which are standard testing parameters in the literature of dynamic networks 1 Code and data available at https://anonymous.4open.science/r/SDSBM-26DC inference [START_REF] Fan | Dynamic Infinite Mixed-Membership Stochastic Blockmodel[END_REF][START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF]. We choose to infer a dynamic membership matrix 𝜃 (𝑡 ) and to provide a universal block-interaction matrix 𝑝 ∀𝑡. Note that the model yields good performances when 𝑝 also has to be inferred, but due to identificability and label-switching issues raised in [START_REF] Matias | Statistical clustering of temporal networks through a dynamic stochastic block model[END_REF], there is not unbiased way to assess the correctness of the inferred memberships values. Additional experiments showed that inferring jointly 𝑝 and 𝜃 yields the same prediction accuracy.The expression we use for this matrix 𝑝 is given Eq.13. We systematically test two variation patterns for 𝜃 : a sinusoidal pattern (Fig. 1-top) and a broken-line pattern (Fig. 1-bottom). Each pattern is generated with different coefficient for each item; the memberships still sum to 1 at all times.

𝑝 =       1 -𝑠 𝑠 0 0 1 -𝑠 𝑠 𝑠 0 1 -𝑠       (13)
To the best of our knowledge, the only attempt to model dynamic parameters in labeled weighted and dynamic networks using a MMSBM is [START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF]. In this work, each epoch is modeled independently from the others. We refer to this baseline as the "No coupling" or "NC" baseline. For reference, we also compare to a baseline that does not consider the temporal dimension and infers a single universal value for each variable (SIMSBM(1)) [START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF]. It has already been assessed that MMSBMs typically outperform state of the art models (NMF, TF, KNN, etc.) in [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF] and [START_REF] Poux-Médard | Serialized Interacting Mixed Membership Stochastic Block Model[END_REF]. Such study is not repeated in the present article, for being out of the scope of our demonstration. We recall that we purposely chose the simplest form of MMSBM for clarity; the reported results may therefore be poor in absolute, but underline an improvement over equally simple models (NC and SIMSBM(1)) brought by our Simple Dynamic prior approach.

We systematically perform a 5-folds cross validation. The model is trained on 80% of the data, 𝛽 is tuned using 10% as a validation set, and the model is evaluated on the 10% left. We choose as metrics the AUC-ROC and the RMSE on the real values of 𝜃 (black line in Fig. 1). The procedure is repeated 5 times; the error bars reported in the experimental results represent the standard error over these folds.

4.1.1 SDSBM unveils complex temporal patterns. In Fig. 4a., we consider 1 000 observations for each item 𝑖 ∈ 𝐼 and vary the number of epochs from 10 to 1 000. In the expression of 𝑝, 𝑠 is set to 0.05. For both the sinusoidal and line-broken memberships, the model shows better predictive performances (in terms of AUC ROC) than the proposed baselines. Interestingly, the SDSBM performances remain stable as the number of epochs increases unlike the NC baseline, which means it alleviates a bias of the temporal modeling proposed in [START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF]. The RMSE with respect to the true parameters remains low over the whole range of tested number of epochs. The RMSE increases as the number of epochs grows because the number of parameters to estimate increases with it; this makes the inference more subject to local variations, which in turns mechanically increases the RMSE. Overall, SDSBM recovers dynamic variations of the membership vectors with a good performance; a sample of the inferred dynamic memberships is shown Fig. 1.

SDSBM works with little data.

A major problem that arises when considering temporal data is the scarcity of observations, because slicing implies reducing the number of observations in each slice. This concern largely arises in social sciences, where data retrieval cannot be automated and requires tedious human labor. Here, we demonstrate that our method works in challenging conditions, when data is scarse. In Fig. 4b., we vary the number of observations available for each item from 100 to 10.000, distributed over 100 epochs. Thus, in the most challenging situation, there is only one observation per epoch used to determine 𝐼 dynamic memberships over 3 clusters. In the expression of 𝑝, 𝑠 is set to 0.05. We see Fig. 4b. that for both patterns, the predictive power of SDSBM remains high in such conditions. Moreover, the RMSE on the true dynamic memberships in this case is fairly low, and decreases rapidly as the number of observations increases. This is due to SDSBM linking the time slices together during the training phase: the smoothing constraint makes it so that every time slices indirectly benefits from the training data of its neighbours.

When the number of observations is high, the "no coupling" baseline [START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF] reaches the performances of SDSBM. This is because as the number of observations in each slice goes to infinity, the models need less to rely on temporal neighbours. However, even for 10.000 observations per item (100 observations per epoch), SDSBM yields better results. As an illustration, the results in Fig. 1 have been obtained using only 5 observations per epoch.

4.1.

3 SDSBM handles highly stochastic interaction patterns. Finally, we control the deterministic character of the block-interaction matrix 𝑝 by varying 𝑠. We express such character as the mean entropy of 𝑝 ⟨𝑆 (𝑝)⟩ with respect to its possible outputs: ⟨𝑆 (𝑝)⟩ = 1 𝐾 𝑘 ∈𝐾 𝑜 𝑝 𝑘 (𝑜) log 𝑝 𝑘 (𝑜). The maximum entropy for the proposed expression of 𝑝 is reached for 𝑠 = 0.5. We consider 1 000 observations spread over 100 epochs. We show in Fig. 4c. that the predictive performance of all three methods drops as the entropy increases. This is expected, as observations are generated from the true model with a higher variance; each observation becomes less informative about the generative process as 𝑠 grows. However, the RMSE on the real parameters inferred using SDSBM remains low at the maximum entropy, meaning the model recovers the correct membership parameters.

Real-world data

4.2.1 Experimental setup. Finally, we demonstrate the validity of our approach on real-world data to argue for its usefulness and scalability. SDSBM builds on previous works on labeled MMSBM and shares the same linear complexity O (|𝑅 • |) with |𝑅 • | the size of the dataset [START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF]. For our experiments, we consider the recent and documented datasets from [START_REF] Kumar | Predicting Dynamic Embedding Trajectory in Temporal Interaction Networks[END_REF]. The Reddit dataset (10.000 users, 984 subreddits they contribute to, ∼670k observations), the LastFm dataset (980 users, 1000 songs users listened to, ∼1.3M observations) and the Wikipedia (Wiki) dataset (8227 users, 1000 pages users edited, ∼157k observations). The goal is to predict over time what subreddit a user will contribute to, what songs they will listen to, and what Wikipedia pages they will edit. The Reddit and Wikipedia datasets contain 1 month of data; we slice them in 1 day long temporal intervals. The LastFm dataset spans over approximately 5 years; we slice it into periods of 3 days each. In addition, we build an additional dataset (Epi) about historical Table 1: Numerical results on real-world datasets -Metrics abbreviations stand for the area under the ROC curve (ROC), the Average Precision (AP), the Normalized Coverage Error (NCE). Metrics for models stand for Simple Dynamic SDM (SDSBM), No Coupling baseline (NC) and the classical static mixed membership SBM (SIMSBM(1)). Overall, our approach allows for a higher predictive power.

ROC AP NCE

Epi SDSBM 0.9025 [START_REF] Ho | Evolving Cluster Mixed-Membership Blockmodel for Time-Evolving Networks[END_REF] 2 . The dataset is made of 117.000 latin inscriptions comprising one or several of 18 social status (slave, soldier, senator, etc.) and its location as one of 62 possible regions, along with an estimated datation spanning from 100BC to 400AD. The goal is to guess the region where a status has been found, with respect to time. The goal is to recover statuses diffusion in roman territories. We slice this dataset in epochs of one year each.

Evaluation is again conducted using a 5-folds cross validation with 80% of training data, 10% of validation data and 10% of testing data for each fold. For each pair (𝑖, 𝑜 𝑡𝑟𝑢𝑒 ) in the test set, we query the probability for every output 𝑜 given 𝑖 and build the confusion matrix by comparing them to 𝑜 𝑡𝑟𝑢𝑒 . In Table 1, we present the results of our method compared to the proposed baselines for various metrics: AUC-ROC (ROC), Average Precision (AP) and Normalized Coverage Error (NCE). The first two metrics evaluate how well models assign probabilities to observations, and the latter evaluates the order in which possible outputs are ranked.

Results

. Overall, we see that our method exhibits a greater predictive power, except for the Reddit dataset where the static SBM performs as well as SDSBM. We explain this by the lack of significant temporal variation over the considered interval. This could be expected, since the dataset comprises roughly 80% of repeated actions [START_REF] Kumar | Predicting Dynamic Embedding Trajectory in Temporal Interaction Networks[END_REF], meaning that users do not significantly explore new communities over a month. This result shows that SDSBM also works well in the static case. On the other datasets, SDSBM performs better often by a large margin, especially for the ROC-AUC, meaning that SDSBM is efficient at distinguishing classes from each other. We also checked that our initial assumption on the smooth temporal variation of the parameters holds. For every real-world dataset, the absolute average change in membership between two consecutive time slices is systematically less than 0.005±0.01. The Figure 5: Geographic evolution of status distribution from latin graves (100BC -500AD) -We applied the SDSBM to the Epigraphy dataset. We recall that our goal is to predict a roman region (e.g. Illyria, Hispania, etc.) given a status (e.g. Slave, Senator, etc.) and a year. We plot the temporal evolution of statuses membership to the five manually labeled clusters (in gray).

For clarity, we removed small membership transfers from the data, which explains why the total cluster's population may vary from one time to another. This plot allows to visualize some global historical trends about the evolution of the Roman Empire (e.g. 3rd century crisis, spread of military presence in Europe, Italy demilitarization, etc.).

absolute variation after 10 time slices is less than 0.03±0.04. Therefore the inferred parameters do not vary abruptly over time. We recall that the model used here is deliberately simplistic for demonstration; low metrics do not means the Simple Dynamic prior does not work, but instead that it should be coupled to a more complex model (any of SIMSBM(x,y,...)).

Case study.

As an illustration of what SDSBM has to offer, we plot in Fig. 5 a possible visualization of the memberships evolution in time for the epigraphy dataset. On the left and on the right, we show the items that are considered in the visualization. The time goes from left (100BC) to right (500AD), and the flows represent the membership transfers between epochs. The grey bars represent the clusters. We manually annotated them by looking at their composition. From this figure, we can recover several historical facts: military presence in Rome was scarce for most of the times considered; Italy concentrates less military presence as time goes (due to its spread over the now extended empire), until the 3rd century crisis that led to its re-militarization; most of the slaves that have been accorded an inscription are located in Italia throughout time; the religious functions (Augusti) are evenly spread on the territory at all times; the libertii (freed slaves) inscriptions are essentially present in Rome and in Italy, etc. Obviously, dedicated works are needed in order to support these illustrative claims, and we believe SDSBM can provide such extended comprehension of these processes.

CONCLUSION

We introduced a simple way to model time in dynamic weighted and labeled networks by assuming a dynamic Mixed-Membership SBM. Our method consists in defining the Simple Temporal prior, ready to plug into any of a whole class of existing static MMSBMs. Time is considered under the single assumption that network's ties do not vary abruptly.

We assessed the performance of the proposed method by defining the SDSBM and testing it in several controlled situations on synthetic datasets. In particular, we show that our prior allows stable performances with respect to the dataset slicing, and that it works well under challenging conditions (small amounts of data or high entropy blocks interaction matrix). We tested SDSBM on large scale real-world datasets and showed it yields better results than two proposed baselines. Finally, we illustrate an application interest on a dataset of Latin inscriptions that indirectly narrates the evolution of the Roman Empire.

We discussed the advantages of using our approach: uneven slicing of observations in time, heterogeneous dynamic time-scales, time-dependent blocks-interaction matrix, and informativeness of the prior; exploring these directions on real-world data may help retrieving meaningful clusters on useful applications. On a further note, we believe that a key interest of our approach is the small amount of data needed to get satisfactory performances. This point is fundamental to a number of social sciences, and we believe our approach could ease the incorporation of automated learning methods in these fields.

Figure 1 :

 1 Figure1: Users' attachment to groups can vary over time -A music listener could cyclically prefer Rock, Jazz or Pop music (top), or listen to either of these without any specific pattern (bottom). For 200 epochs containing only 5 observations each, our approach (in red) infers any smooth dynamic membership pattern, and does it more accurately than static models (in green[START_REF] Godoy-Lorite | Accurate and scalable social recommendation using mixed-membership stochastic block models[END_REF]) and models that consider each time slice independently (in blue[START_REF] Tarrés-Deulofeu | Tensorial and bipartite block models for link prediction in layered networks and temporal networks[END_REF]).

  times Intermediate memberships' prior variance

Figure 2 :

 2 Figure2: Prior probability's variance on memberships at all times according to the temporal neighbourhood -Variance of the prior over a membership entry (filled curves, we represented 3 such entries in different colors as illustration) as a function of time, given some temporal neighbours (black dots). This illustration considers an averaging kernel as 𝜅 (𝑡, 𝑡 ′ ) = 1 |𝑡 -𝑡 ′ | . When inferring a parameter 𝑥 (𝑡 ) at a time 𝑡, the variance of its prior probability 𝑃 (𝑥(𝑡 ) ) depends on 𝑡 relative to the temporal neighbours. Here for instance, the variance is null at 𝑡 = 2 because 𝜅 (𝑡, 𝑡 ′ ) diverges, and so does 𝛼(𝑡 ) , hence the variance collapsing to 0.

  Dynamic prior's mode. Our single constraint is that 𝜃 (𝑡 ) 𝑖 and 𝑝 (𝑡 )

Figure 3 :

 3 Figure3: Prior probability on a membership vector for various values of 𝛽 according to temporal neighbourhood -Darker means higher probability. Projected on a simplex tri-space (each of 3 axes ranges from 0 to 1). The white dots represent the temporal neighbours of the considered 3D vector. Their average is given as ⟨𝜃 ⟩ using a uniform weight function 𝜅 (𝑡, 𝑡 ′ ) for illustration purpose. 𝛽 controls the variable's prior variance around its neighbours.

  . It can be derived as: log 𝑃 (𝜃, 𝑝 |𝑅 • ) ∝ log 𝑃 (𝑅 • |𝜃, 𝑝)

									𝑃 (𝜃	(𝑡 ) 𝑖 )	𝑃 (𝑝	(𝑡 ) 𝑘 )
								𝑡	𝑖	𝑘
	=		∑	log	∑	𝜃	(𝑡 ) 𝑖,𝑘 𝑝	(𝑡 ) 𝑘 (𝑜)	(8)
		(𝑖,𝑜,𝑡 ) ∈𝑅 •	𝑘 ∈𝐾
		+	∑	∑	log 𝑃 (𝜃	(𝑡 ) 𝑖 )	∑	log 𝑃 (𝑝	(𝑡 ) 𝑘 )
			𝑡	𝑖					𝑘
	≥	∑ (𝑖,𝑜,𝑡 ) ∈𝑅 •	∑ 𝑘 ∈𝐾	𝜔	(𝑡 ) 𝑖,𝑜 (𝑘) log	𝜃	(𝑡 ) 𝑖,𝑘 𝑝 𝜔 (𝑡 ) (𝑡 ) 𝑘 (𝑜) 𝑖,𝑜 (𝑘)
		+	∑	∑	log 𝑃 (𝜃	(𝑡 ) 𝑖 )	∑	log 𝑃 (𝑝	(𝑡 ) 𝑘 )
			𝑡	𝑖					𝑘

  𝜕 log 𝑃 (𝜃, 𝑝 |𝑅 • ) -𝑖 ′ ,𝑡 ′ 𝜙

								parameters. Starting with the membership matrix entries 𝜃	(𝑡 ) 𝑖,𝑘 :
												𝜕𝜃	(𝑡 )	(𝑡 ′ ) 𝑖 ′ ( 𝑘 ′ 𝜃	(𝑡 ′ ) 𝑖 ′ ,𝑘 ′ -1)	= 0
												𝑖,𝑘
								↔	∑ 𝑜 ∈𝜕 (𝑖,𝑡 )	𝜔	𝑖,𝑘 (𝑡 ) 𝑖,𝑜 (𝑘) 𝜃 (𝑡 )	+	𝑖,𝑘 𝛽 ⟨𝜃 (𝑡 ) 𝑖,𝑘 ⟩ 𝜃 (𝑡 )	-𝜙	(𝑡 ) 𝑖	= 0
								↔	∑	𝜔	(𝑡 ) 𝑖,𝑜 (𝑘) + 𝛽 ⟨𝜃	(𝑡 ) 𝑖,𝑘 ⟩ = 𝜙	(𝑡 ) 𝑖 𝜃	(𝑡 ) 𝑖,𝑘
									𝑜 ∈𝜕 (𝑖,𝑡 )
								↔	∑	∑	𝜔	(𝑡 ) 𝑖,𝑜 (𝑘)	+𝛽	∑	⟨𝜃	(𝑡 ) 𝑖,𝑘 ⟩	= 𝜙	(𝑡 ) 𝑖
									𝑜 ∈𝜕 (𝑖,𝑡 )	𝑘	𝑘
												=1 (Eq.9)	=1 (Eq.1)
								↔	𝑜 ∈𝜕 (𝑖,𝑡 ) 𝜔 𝑁 𝑖,𝑡 + 𝛽 (𝑡 ) 𝑖,𝑜 (𝑘) + 𝛽 ⟨𝜃	(𝑡 ) 𝑖,𝑘 ⟩	= 𝜃	(𝑡 ) 𝑖,𝑘
					(𝑡 ) 𝑖,𝑜 (𝑘).
	Using this expression, we can rewrite the log-likelihood log 𝑃 (𝑅 • |𝜃, 𝑝)
	as [4, 9, 20]:
	log 𝑃 (𝑅 • |𝜃, 𝑝) =	∑ (𝑖,𝑜,𝑡 ) ∈𝑅 •	∑ 𝑘 ∈𝐾	𝜔	(𝑡 ) 𝑖,𝑜 (𝑘) log	𝜃	(𝑡 ) 𝑖,𝑘 𝑝 𝜔 𝑖,𝑜 (𝑘) (𝑡 ) 𝑘 (𝑜) (𝑡 )	(10)
	3.3.2 M step. Taking back the first line of Eq.8 and substituting
							(𝑡 ) 𝑖 ), and finally
	compute the maximization equations with respect to the model's

• , written 𝜔 with Eq.7 and Eq.10, we get an unconstrained expression of the posterior distribution. We introduce Lagrange multipliers to account for the constraints of Eq.1 (𝜙 (𝑡 ) 𝑖 ) and Eq.2 (𝜓

  𝑎 = 1, 2, ... for instance, where 𝑁 𝑡 ′ is the number of observations in the time slice 𝑡 ′ . Overall, the Simple Dynamic prior works by inferring the variables using both microscopic and mesoscopic temporal scales. If a time slice 𝑡 has few observations but some of its neighbours have a greater number of them, learning the parameters at 𝑡 is helped mostly by the population of its closest (1|Δ𝑡 | ) and most populated (𝑁 𝑡 ′ ) neighbours, and less influenced by further and less populated epochs. This is what is illustrated in Fig.2.

	a.	0.80 0.85			Sinusoids				No coupling SDSBM Broken lines
		0.75											No temporal dependence
		0.60 0.65 0.70 AUC ROC				SDSBM No coupling No temporal dependence
		0.55
		0.50
			0.4		SDSBM No coupling							SDSBM No coupling
					No temporal dependence					No temporal dependence
			0.3
		RMSE	0.2
			0.1
			0.0	0	200	400 Number of epochs 600	800		1000	0	200	400 Number of epochs 600	800	1000
	b.	0.80 0.85											SDSBM No coupling
		0.75											No temporal dependence
		0.60 0.65 0.70 AUC ROC				SDSBM No coupling No temporal dependence
		0.55
		0.50
			0.4				SDSBM No coupling						SDSBM No coupling
							No temporal						No temporal dependence
			0.3
		RMSE	0.2
			0.1
			0.0	0	2000 Number of observations per item 4000 6000 8000	10000	0	2000 Number of observations per item 4000 6000 8000	10000
	c.	0.80 0.85				No coupling SDSBM						SDSBM No coupling
		0.75				No temporal dependence				No temporal dependence
		0.60 0.65 0.70 AUC ROC
		0.55
		0.50
			0.4				SDSBM No coupling						SDSBM No coupling
			0.3				No temporal dependence				No temporal dependence
		RMSE	0.2
			0.1
			0.0	0.0	0.1	0.2	0.3 Entropy of p 0.4	0.5	0.6	0.7	0.0	0.1	0.2	0.3 Entropy of p 0.4	0.5	0.6	0.7
	Figure 4: Results on synthetic data -(a.) SDSBM retrieves
	the correct dynamic memberships and is little influenced
	by the data slicing. (b.) SDSBM works well on tiny datasets.
	(c.) SDSBM retrieves dynamic memberships in challenging
	situations.
																𝑁 𝑡 ′ |𝑡 -𝑡 ′ | 𝑎 where

Clauss-Slaby repository, http://www.manfredclauss.de/fr/index.html


Asemare Mengistie Taddese
email: asemare-mengistie.taddese@uha.fr

Mohamed Lo

Nicolas Verrier
email: nicolas.verrier@uha.fr

Matthieu Debailleul

Olivier Haeberlé

Jones Tomographic Diffractive Microscopy with a Polarized Array Sensor

Tomographic Diffractive Microscopy (TDM) based on scalar light-field approximation is widely implemented. Samples exhibiting anisotropic structures however necessitate to account for the vectorial nature of light, leading to 3-D quantitative polarimetric imaging. In this work, we have developed a high-numerical aperture (at both illumination and detection) Jones TDM system, with detection multiplexing via a Polarized Array Sensor (PAS), for imaging optically birefringent samples at high resolution. The method is first studied through image simulations. To validate our setup, experiment using a sample containing both birefringent and non-birefringent objects is performed. Araneus diadematus spider silk fiber and Pinna Nobilis oyster shell crystals are finally studied, allowing to assess both birefringence and fast-axis orientation maps.

Introduction

Polarization imaging is a promising technique in several fields, and hence gaining much attention, especially due to its capability of identifying structures that are unseen by conventional imaging.

Applications of polarization imaging range from macroscopic to microscopic scales, and from the simplest and robust to unique and comprehensive techniques [START_REF] Tuchin | Tissue optics and photonics: light-tissue interaction[END_REF]. Single cell study and tissue polarization imaging [START_REF] Stender | Single cell optical imaging and spectroscopy[END_REF], as well as clinical applications of such imaging technologies, mainly for qualitative analyzing of biological samples in large Field of View (FOV) [START_REF] Nan | Linear polarization difference imaging and its potential applications[END_REF], are also of great interest.

In recent years, several research groups have shown interest towards quantitative 2-D Jones matrix polarization sensitive holography [START_REF] Shin | Reference-free polarization-sensitive quantitative phase imaging using single-point optical phase conjugation[END_REF][START_REF] Wang | Jones phase microscopy of transparent and anisotropic samples[END_REF][START_REF] Colomb | Polarization microscopy by use of digital holography: application to optical-fiber birefringence measurements[END_REF][START_REF] Aknoun | Quantitative retardance imaging by means of quadri-wave lateral shearing interferometry for label-free fiber imaging in tissues[END_REF][START_REF] Wang | Optical polarization in biomedical applications[END_REF][START_REF] Yang | Single-shot digital holographic microscopy for quantifying a spatially-resolved jones matrix of biological specimens[END_REF][START_REF] Park | Jones matrix microscopy for living eukaryotic cells[END_REF]. For instance, such an approach has been used to study anisotropic samples in 2-D holographic microscopy for live cell imaging [START_REF] Wang | Jones phase microscopy of transparent and anisotropic samples[END_REF], and more recently, synthetic aperture holographic tomography has been developed for noninvasive imaging of cancer cells, with better sensitivity compared to standard microholography [START_REF] Park | Jones matrix microscopy for living eukaryotic cells[END_REF].

For 3-D imaging, polarization sensitive optical coherence tomography [START_REF] De Boer | Polarization sensitive optical coherence tomography-a review[END_REF], 3-D polarized light imaging [START_REF] Menzel | A jones matrix formalism for simulating three-dimensional polarized light imaging of brain tissue[END_REF] and fluorescence confocal microscopy [START_REF] Smalyukh | Three-dimensional imaging of orientational order by fluorescence confocal polarizing microscopy[END_REF] are the well established techniques.

While these approaches have shown their invaluable importance for clinical application along

with their sub micrometer resolutions, quantification of birefringence in 3-D is difficult, and only qualitative anisotropic contrast is revealed. Recently, large scale polarization contrast TDM in transmission configuration has been implemented for improving image contrast of anisotropic specimens such as zebrafish embryos [START_REF] Van Rooij | Polarization contrast optical diffraction tomography[END_REF].

Three-dimensional and truly quantitative polarimetric imaging modalities are less common.

Proposed by E. Wolf in 1969 [START_REF] Wolf | Three-dimensional structure determination of semi-transparent objects from holographic data[END_REF], TDM under scalar approximation has been used in several research areas [START_REF] Devaney | A filtered backpropagation algorithm for diffraction tomography[END_REF][START_REF] Haeberlé | Tomographic diffractive microscopy: basics, techniques and perspectives[END_REF][START_REF] Park | Quantitative phase imaging in biomedicine[END_REF][START_REF] Jin | Tomographic phase microscopy: principles and applications in bioimaging[END_REF][START_REF] Debailleul | Holographic microscopy and diffractive microtomography of transparent samples[END_REF][START_REF] Ossowski | Near-infrared, wavelength, and illumination scanning holographic tomography[END_REF][START_REF] Denneulin | Gsure criterion for unsupervised regularized reconstruction in tomographic diffractive microscopy[END_REF]. However, such configurations fail at quantifying samples with anisotropic nature. Although V. Lauer reformulated Wolf's work in 2002, providing full vectorial formulation by accounting for polarization of light [START_REF] Lauer | New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope[END_REF], it is only very recently that polarimetric TDM regained attention.

Advances towards polarization sensitive TDM have confirmed the importance of quantitative 3-D polarimetric approaches. Polarization was first taken into account in reflection TDM by Zhang in 2013, [START_REF] Zhang | Full-polarized tomographic diffraction microscopy achieves a resolution about one-fourth of the wavelength[END_REF], and recently polarization sensitive transmission TDM configurations have also been implemented [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF][START_REF] Shin | Tomographic measurement of dielectric tensors at optical frequency[END_REF].

In general, experimental implementations and image reconstruction processes of polarization sensitive TDM to image weakly birefringent specimens, such as biological sample (cells, cell structure, and tissue), with high-resolution and good sensitivity will need to be well explored and understood before becoming a routine technique. In this work, we have developed a highresolution polarization sensitive TDM (PS-TDM) based on a Polarized Array Sensor (PAS) to image weakly birefringent samples. The use of a PAS paves the way for implementing a simplified Jones and/or Mueller matrix PS-TDM by multiplexing data acquisition. Image simulations of a synthetic sample as well as experiments using a sample containing both birefringent and non-birefringent objects is performed to validate the system.

Theory

Vectorial Born Approximation

The vectorial Helmholtz equation expresses a monochromatic light-wave propagating through a non-magnetic and weakly scattering anisotropic object embedded in a homogeneous and isotropic background medium, and writes as:

[∇ 2 + 𝑘 2 0 ]E 𝑡 (r) = 𝑘 2 0 [𝜖 𝑟 (r) -I]E 𝑡 (r). (1)
E 𝑡 (r) is the total propagating vectorial electric field detected at position r = (𝑥, 𝑦, 𝑧). 𝜖 𝑟 (r)

and I are the relative electric permittivity and identity matrix, respectively. The modulus of the wave vector is 𝑘 0 = 2𝜋𝑛 0 /𝜆, with 𝜆 and 𝑛 0 being the wavelength in vacuum and background medium refractive index (RI), respectively. Considering 𝑛 as the RI tensor of the birefringent sample and a weakly birefringent sample, |𝛿𝑛| I, 𝜖 𝑟 (r) is given as:

𝜖 𝑟 (r) = (𝑛 0 I + 𝛿𝑛(r)) 2 ≈ 𝑛 2 0 I + 2𝑛 0 𝛿𝑛(r), (2)
where 𝑛(r) = 𝑛 0 I + 𝛿𝑛, 𝛿𝑛 being the RI tensor contrast of the anisotropic object.

The total field vector is given as the sum of scattered E 𝑠 (r) and incident field E 𝑖 (r). The simplified solution for Eq. ( 1) writes as [START_REF] Lauer | New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope[END_REF][START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF]:

E 𝑡 (r) = E 𝑖 (r) + E 𝑠 (r) = E 𝑖 (r) + ∫ G(r, r ) × V(r ) × E 𝑡 (r, r )𝑑r , (3)
where V(r) = 𝑘 2 0 (𝑛(r) 2 /𝑛 2 0 -I) is the object's scattering potential tensor (Born's scattering field vector). For a general birefringent sample, V(r) is a third order tensor. The Green's tensor G(r, r ) in the homogeneous and isotropic background medium is diagonal :

G(r, r ) = I𝑔(r, r ),
with 𝑔(𝑟, 𝑟 ) being the scalar solution, given as:

𝑔(𝑟, 𝑟 ) = exp(𝑖𝑘 0 𝑛 0 |r -r |) 4𝜋|r -r | . (4)
Eq. ( 3) is further simplified by considering a weakly scattering birefringent sample, the total field vector inside this sample can be regarded as the incident field vector (which is the consequence of first order Born approximation [START_REF] Born | Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF]), i.e. E 𝑠 (r)

E 𝑖 (r) => E 𝑡 (r) ≈ E 𝑖 (r). E 𝐵 (r) = ∫ G(r, r ) × V(r ) × E 𝑖 (r, r )𝑑r , (5)
where E 𝐵 (r) is Born's scattered field vector. Considering a plane wave illumination, E 𝑖 (r) = E 𝑖 (r) exp(𝑖k 𝑖 .r), we have:

E 𝐵 (r) = ∫ G(r, r ) × V(r ) × E 𝑖 (r) exp(𝑖k 𝑖 .r )𝑑r , (6)
where E 𝑖 (r) is the illumination beam's vectorial complex amplitude.

It should be noted that, under scalar approximation, it is possible to easily adapt these developments to Rytov or Born approximation, as studied in literature. The Rytov complex field is given by E 𝑟 (r) = ln(E 𝐵 (r)). Imaginary and real parts are computed as

{E 𝑟 (r)} = arg [E 𝑡 /E 𝑖 ],
and

{E 𝑟 (r)} = ln(|E 𝑡 |/|E 𝑖 |)
, respectively. In that case, the phase difference

𝜙 𝑡 -𝜙 𝑖 = arg [E 𝑡 /E 𝑖 ]
has to be unwrapped.

Jones Formalism

Quantitative polarimetric imaging techniques relying on a coherent light commonly utilize the Jones formalism to relate the input and output fields [START_REF] Wang | Jones phase microscopy of transparent and anisotropic samples[END_REF][START_REF] Park | Jones matrix microscopy for living eukaryotic cells[END_REF][START_REF] Ferrand | Ptychography in anisotropic media[END_REF][START_REF] Jones | A new calculus for the treatment of optical systemsv. a more general formulation, and description of another calculus[END_REF]. Aiming at retrieving the scattered field amplitude, we start by developing approximated mathematical expression for PS-TDM imaging using Jones calculus. The Jones tensor J obj of the birefringent sample under examination relates the scattered field vector E 𝑠 (r) and the incident field vector E 𝑖 (r) as:

E 𝑠 (r) = J obj (r)E 𝑖 (r). (7)
Recalling Born approximation given in Eq. ( 6), one gets:

E 𝐵 (r) = J obj (r)E 𝑖 (r) exp(𝑖k 𝑖 .r) = ∫ G(r, r ) [3×3] × V(r ) [3×3] × E [3×3] 𝑖 exp(𝑖k 𝑖 .r )𝑑𝑟 , (8)
where:

J obj =          𝐽 𝑥 𝑥 𝐽 𝑥 𝑦 𝐽 𝑥𝑧 𝐽 𝑦 𝑥 𝐽 𝑦𝑦 𝐽 𝑦𝑧 𝐽 𝑧 𝑥 𝐽 𝑧 𝑦 𝐽 𝑧𝑧          . ( 9
)
Following vectorial first order Born approximation, in the same way as for scalar approximation [START_REF] Wolf | Three-dimensional structure determination of semi-transparent objects from holographic data[END_REF], one gets a linear relationship between the Jones tensor components and the sample's permittivity tensor of the anisotropic sample:

J[3×3] obj (k -k 𝑖 ) = 2𝜋𝑖 𝑘 𝑧 Ṽ(k -k 𝑖 ) [3×3] , (10)
where 𝑘 𝑧 is the wave-vector's 𝑧 component magnitude. Ṽ(kk 𝑖 ) [3×3] and

J[3×3] obj (k -k 𝑖 ) are
the Fourier transforms of V(r) and J obj (r).

From now on, we approximate the RI tensor by neglecting its third component, as in [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF].

Such approximation is valid for anisotropic objects exhibiting weak polarization response along the optical axis. Doing so makes it possible to compute all the 9 components of Ṽ(kk 𝑖 ) [3×3] with only two independent polarization components, as one is able to measure only the x and y components of the scattered field using a PS-TDM. As a consequence, Eq. ( 10) now reduces to a second order tensor:

Japprox obj (k -k 𝑖 ) = 2𝜋𝑖 𝑘 𝑧 Ṽ(k -k 𝑖 ) [2×2] , (11)
where the birefringent sample's Jones tensor 𝐽 𝑠 (r) := J approx obj (r) is given as:

𝐽 𝑠 (r) =       𝐽 𝑥 𝑥 𝐽 𝑥 𝑦 𝐽 𝑦 𝑥 𝐽 𝑦𝑦       . ( 12
)
Jones matrix analysis of birefringent samples necessitates illumination along two independent incident polarization states (𝑥 and 𝑦), as well as analysis of the scattered light along these two axes. Hence, for full Jones matrix retrieval, four measurement are usually required [START_REF] Van Rooij | Polarization contrast optical diffraction tomography[END_REF][START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF][START_REF] Shin | Tomographic measurement of dielectric tensors at optical frequency[END_REF],

which reduces speed, and can also decrease sensitivity due to measurement artifacts.

A first advantage of using a PAS as detector, providing mosaicked images in a single acquisition for various polarization states, is that data acquisition can be simplified, by multiplexing detection.

Also, for interferometric systems, linearly cross-polarized configurations lead to zero background intensity, making phase estimation of the reference beam ambiguous. This results in a random background phase, which makes phase unwrapping challenging [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF][START_REF] Zhang | Far-field diffraction microscopy at 𝜆/10 resolution[END_REF]. Conversely, the use of circularly polarized beam greatly alleviates this problem, since it results in a non-zero background intensity, allowing for easier phase unwrapping. In such configuration, the sample is sequentially illuminated with right-, then left-circular polarization. We therefore now extend the formalism proposed in [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF] to the case of circular polarization illumination.

Let us symbolize the output field at each micro-analyzer's orientation:

• 𝐸 𝑜 1𝑥 : right-circular illumination and detection at 0 The polarization sensitive hologram to be recorded using the polarization array sensor results from interference between the output field 𝑂 𝑖 𝑗 and the reference 𝑅 𝑖 𝑗 beam,

𝐼 𝑖 𝑗 = |𝑂 𝑖 𝑗 + 𝑅 𝑖 𝑗 | 2 = |𝑂 𝑖 𝑗 | 2 + |𝑅 𝑖 𝑗 | 2 + 𝑂 𝑖 𝑗 𝑅 * 𝑖 𝑗 + 𝑂 * 𝑖 𝑗 𝑅 𝑖 𝑗 , (13)
where 𝑖 and 𝑗 represents the illumination and analysis axes, respectively. Extracting the first diffraction order leads to:

𝑂 𝑖 𝑗 𝑅 * 𝑖 𝑗 =       𝐸 𝑜 1𝑥 𝐸 𝑜 2𝑥 𝐸 𝑜 1𝑦 𝐸 𝑜 2𝑦       , (14)
which permits to write relationship between the measured complex scattered field along the specific analyzer axes and the birefringent sample's Jones tensor components (see Appendix):

      𝐽 𝑥 𝑥 𝐽 𝑥 𝑦 𝐽 𝑦 𝑥 𝐽 𝑦 𝑦       =       𝐸 𝑜 1𝑥 + 𝐸 𝑜 2𝑥 𝑖(𝐸 𝑜 1𝑥 -𝐸 𝑜 2𝑥 ) 𝑖(𝐸 𝑜 2𝑦 -𝐸 𝑜 1𝑦 ) 𝐸 𝑜 1𝑦 + 𝐸 𝑜 2𝑦       . ( 15
)
The above expression provides a linear relationship between the scattered field and the birefringent sample's Jones tensor components. For light illuminating the specimen at normal incidence, and propagating in 𝑧 direction, the illumination beam has either 𝑥 and/or 𝑦 polarization components, and can be written as:

𝐸 𝑖 = 𝐸 𝑖 𝑥 ê𝑥 + 𝐸 𝑖 𝑦 ê𝑦 , (16)
ê𝑥 and ê𝑦 being unit vectors along 𝑥 and 𝑦 axis. 𝐸 𝑖 𝑥 and 𝐸 𝑖 𝑦 are complex input fields, with amplitude and phase. But in illumination rotation TDM configuration, the object is illuminated with various angular directions. Hence, to correctly extract the polarization information for each illumination beam, one has to decompose it into its respective polarization contributions. To do so, we use the approach proposed in [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF], introducing the rotation matrix R 𝑚 as an orthogonal 3

× 3 matrix satisfying R -1 𝑚 = R 𝑇 𝑚 and R 𝑚 R 𝑇 𝑚 = R 𝑇 𝑚 R 𝑚 = I: R 𝑚 =          cos 2 (𝜙) cos(𝜃) + sin 2 (𝜙) sin(𝜙) cos(𝜙) [cos(𝜃) -1] sin(𝜃) cos(𝜙) sin(𝜙) cos(𝜙) [cos(𝜃) -1] sin 2 (𝜙) cos(𝜃) + cos 2 (𝜙) sin(𝜃) sin(𝜙) -sin(𝜃) cos(𝜙) -sin(𝜃) sin(𝜙) cos(𝜃)          , (17)
with 𝜙 and 𝜃 being the azimuthal and polar angles, respectively, as depicted in Fig. 1.

So, the illumination beam, of amplitude 1, can now be specified by its wave vector 𝑘 𝑖 =

[cos(𝜙) sin(𝜃), sin(𝜙) sin(𝜃), cos(𝜃)]. Hence, for a general polarized beam entering the condenser, the illumination beam can be computed from the input beam 𝐸 𝑖𝑛 𝑝 at any polar angle by using R 𝑚 , as given in [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF]:

E 𝑖 (𝑟) = R 𝑚 [𝐸 𝑖𝑛 𝑝 𝑥 𝐸 𝑖𝑛 𝑝 𝑦 𝐸 𝑖𝑛 𝑝 𝑧 ] 𝑇 , (18)
For transverse electric field (TE) 𝐸 𝑖𝑛 𝑝 𝑧 is always zero for the normal incidence but non zero for illumination at an arbitrary angle. Once the sample is illuminated by a polarized plane wave, the diffracted field passes through the high NA objective lens before it is measured by the polarization array sensor. The measured electric field, E 𝑚 , is then computed from the scattered field just before the objective by implementing the same procedure as above,

E 𝑚 = R 𝑇 𝑚 [𝐸 𝑠𝑥 𝐸 𝑠𝑦 𝐸 𝑠𝑧 ] 𝑇 . ( 19
)
Finally the information recorded by the polarization array sensor along various micro-analyzers axes can be computed as:

𝐸 𝑚𝑥 = E 𝑚 . â𝑥 and E my = E m .â y , (20)
where â𝑥 and â𝑦 are unit vectors along 0 • and 90 • micro-analyzer axis.

Hologram Demosaicking and Interpolation

The polarization array sensor we use consists in a CMOS sensor with polarization filters deposited on each photo-receptor, in groups of 4, oriented at 0, 45 or nearest-neighbor interpolation techniques [START_REF] Gao | Bilinear and bicubic interpolation methods for division of focal plane polarimeters[END_REF], we selected this approach since it provides smoother images with less interpolation artifacts [START_REF] Mihoubi | Survey of demosaicking methods for polarization filter array images[END_REF].

An important consideration for optimal image reconstruction is to fulfill Nyquist sampling requirement [START_REF] Lauer | New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope[END_REF]. For example, compared to our setup as described in [START_REF] Taddese | Optimizing sample illumination scanning in transmission tomographic diffractive microscopy[END_REF], sampling has to be tuned, when switching from a standard (polarization insensitive) CMOS sensor to a polarization sensitive camera, to satisfy Nyquist criterion for each sub-image. This is performed by adapting magnification in the detection arm of the interferometer (see also section 4). Doing so ensures that hologram demosaicking and interpolation, as previously described, delivers images containing maximal frequencies as allowed by the optical system.

Other Polarimetric Information

Unlike other 2-D polarimetric techniques that need insights about the thickness of the sample (since they only measure phase changes) [START_REF] Ferrand | Ptychography in anisotropic media[END_REF][START_REF] Aknoun | Quantitative retardance imaging of biological samples using quadriwave lateral shearing interferometry[END_REF], thanks to holographic imaging and numerical reconstruction based on first order Born approximation, PS-TDM can directly estimates the 3-D refractive index tensor 𝑛 of weakly scattering samples. This allows for estimating other polarimetric information by using further mathematical simplification. For instance, for a homogeneous anisotropic sample, it is possible to implement pixel-by-pixel diagonalization of the RI tensor to obtain two principal refractive indices, 𝑛 1 and 𝑛 2 . This can be given as:

𝑛(r) =       𝑛 𝑥 𝑥 𝑛 𝑥 𝑦 𝑛 𝑦 𝑥 𝑛 𝑦 𝑦       = R       𝑛 1 0 0 𝑛 2       R -1 =       𝑛 1 cos 2 (𝜃) + 𝑛 2 sin 2 (𝜃) (𝑛 2 -𝑛 1 ) cos(𝜃) sin(𝜃) (𝑛 2 -𝑛 1 ) cos(𝜃) sin(𝜃) 𝑛 1 sin 2 (𝜃) + 𝑛 2 cos 2 (𝜃)       , (21)
in which R is a rotation matrix, which is a 2-D form of R 𝑚 , given as:

R =       cos(𝜃) sin(𝜃) -sin(𝜃) cos(𝜃)       . ( 22
)
From the above expression, it is possible to obtain several 3-D polarimetric information, such as birefringence Δ𝑛, average RI 𝑛 𝑎𝑣𝑔 , and fast-axis orientation 𝜃:

Δ𝑛 = 𝑛 1 -𝑛 2 ( 23
)
𝑛 𝑎𝑣𝑔 = 𝑛 1 + 𝑛 2 2 (24) 𝜃 = 1 2 arctan 𝑛 𝑥 𝑦 + 𝑛 𝑦 𝑥 𝑛 𝑥 𝑥 -𝑛 𝑦𝑦 . (25)
Note that to get an orientation map, 𝜃 ∈ [-𝜋, 𝜋], one needs to use arctan 2 function and the computed map would be prone to phase wrapping problem. Under the Rytov approximation, the algorithm needs the phase as an input data, hence a phase unwrapping algorithm is implemented [START_REF] Volkov | Deterministic phase unwrapping in the presence of noise[END_REF].

As pointed out by one of the anonymous reviewers, a point of attention is the following: as light propagates through a birefringent sample, polarization rotates proportionally to the path-length.

This lead to a phase jump, forbidding detection to discern polarization states that are rotated by an angle 𝜙 from those rotated by an angle 𝜙+𝑛2𝜋 (n being an integer). To avoid this ambiguity, we here restrict investigations to small and low birefringent samples.

Image Reconstruction Simulations

Forward Model

Before experimental investigations, image simulations are conducted to validate the PS-TDM method by mimicking experimental conditions. Thus, Vectorial Beam Propagation Method (V-BPM) is implemented as a forward model to propagate the incident field through the considered phantoms. It provides a compromise between propagation accuracy and computation time, compared to finite-difference time-domain method (FDTD) and Born approximation [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF][START_REF] Shin | Tomographic measurement of dielectric tensors at optical frequency[END_REF][START_REF] Fan | Reconstructing complex refractive-index of multiply-scattering media by use of iterative optical diffraction tomography[END_REF].

Afterwards, the inversion algorithm that is used to reconstructed the refractive index tensor components from the vectorial scattered field is based on first order approximation. V-BPM necessitates dividing the 3-D sample into finitely small thickness 2-D slices, which allows to propagate the incident light beam through these slices.

The V-BPM model considering a slowly varying electromagnetic amplitude (SVEA) wave is

given as [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF]: The propagation is performed using the angular spectrum method.

e(𝑥, 𝑦, 𝑧 + 𝛿𝑧) = exp(𝑖2𝑘 𝑣 𝛿𝑛𝛿𝑧/cos 𝛼) refraction F -1 F e(𝑥, 𝑦, 𝑧) I exp -𝑖 𝑘 2 𝑥 + 𝑘 2 𝑦 𝑘 𝑧 𝛿𝑧 propagation , (26)

Simulation Parameters

The sample is assumed to be illuminated by a polarized plane wave of wavelength 𝜆 = 633 𝑛𝑚, and with unit amplitude. Simulations are performed in a 𝑥 × 𝑦 × 𝑧 = 28 × 28 × 28 𝜇m 3 volume.

The volume is sampled into 512 × 512 × 512 voxels of dimensions Δ𝑥 = Δ𝑦 = Δ𝑧 = 55 nm. In all simulations, the sample is assumed to be a homogeneous anisotropic sample, in which the RI tensor can be diagonalized into the principal refractive indices.

Using the V-BPM forward model, 100 holograms are generated for each polarization states using a rose scanning (8 overlapping petals, as given in [START_REF] Taddese | Optimizing sample illumination scanning in transmission tomographic diffractive microscopy[END_REF][START_REF] Taddese | Optimizing sample illumination scanning for reflection and 4pi tomographic diffractive microscopy[END_REF]) scheme within a cone of illumination whose polar angle is between 𝛼 𝑚𝑎𝑥 = [-67 • , 67 • ], corresponding to our experimental setup.

Such illumination scheme is chosen due to its simplified analytical expression, which helps to easily generate the incident beam at several directions. Once the diffracted field vectors are computed using V-BPM as forward model, the Jones matrix components are retrieved. Vectorial Rytov approximation is then implemented to calculate the synthetic object's 3-D RI tensor components, which are processed to give the principal refractive indexes.

Simulation Results

Image simulations have been performed using a home-developed Python code (Python 3.10.0). To do so, a spherical bead containing an ellipsoidal anisotropic structures is considered. The (central slices) of the computed refractive index distributions that are retrieved from the reconstructed RI tensor components are shown in Figure 3. One can observe blurring around those structures that are multiple scattering layers and around the edges, which results from the inversion algorithm based on Rytov approximation single-scattering assumption. Also, underestimation of the refractive index is observed, as can be confirmed from the RI profile plot along the 𝑦-axis. This is expected, since only 100 holograms are used. Even if we increase the number of holograms, employing vectorial Rytov approximation as inversion algorithm affects the image reconstruction efficiency. Such behavior is often observed when using simpler reconstruction methods such as Born or Rytov linear inversions [START_REF] Taddese | Optimizing sample illumination scanning in transmission tomographic diffractive microscopy[END_REF][START_REF] Taddese | Optimizing sample illumination scanning for reflection and 4pi tomographic diffractive microscopy[END_REF], and clearly constitutes a limit, but those methods have the advantage of being very fast, potentially allowing for real-time acquisition-reconstructiondisplay of the images, which would be interesting for high-throughput screening [START_REF] Bailleul | Tomographic diffractive microscopy: towards high-resolution 3-d real-time data acquisition, image reconstruction and display of unlabeled samples[END_REF]. One can also observe an elongation along the optical axis due to the known missing-cone problem, characteristic of all (non-confocal) transmission microscopes [START_REF] Streibl | Three-dimensional imaging by a microscope[END_REF]. Finally, note that due to intermixing of information, residual of the ellipsoidal structure from 𝛿𝑛 1 is observed in the 𝛿𝑛 2 image.

Experimental Setup

The experimental setup that we developed is an improvement of our previous systems [START_REF] Debailleul | Holographic microscopy and diffractive microtomography of transparent samples[END_REF][START_REF] Taddese | Optimizing sample illumination scanning in transmission tomographic diffractive microscopy[END_REF][START_REF] Bailleul | Tomographic diffractive microscopy: towards high-resolution 3-d real-time data acquisition, image reconstruction and display of unlabeled samples[END_REF],

based on the Mach-Zehnder interferometer, as depicted in Fig. 4. The laser beam (He-Ne, wavelength 𝜆 = 633 nm) is split into a reference beam and illumination beam, using a beam splitting optical fiber. In order to reduce polarization fluctuations of the reference beam, it is fed through a polarization maintaining optical fiber from Thorlabs (PM460-HP). Both reference and illumination beams are collimated. Then, a coherent plane wave exhibiting circular polarization is generated by using polarizing components. The polarization generating components, mounted in both illumination and reference arms, consist of a linear polarizer and a quarter waveplate. The circularly polarized illumination beam is then deflected by the tip/tilt mirror (Newport FSM300 Fast Steering Mirror), and angularly scans the sample, sandwiched between two microscope objectives (oil immersion 𝑛 𝑖𝑚𝑚 = 1.518, NA = 1.4 100× Olympus UPLFLN100XO2-2), the first serving as high-magnification condenser, the second as the detection objective.

Once the sample is illuminated by the plane waves with varying angular direction and polarization, the scattered field and reference beam are recombined by the recombining cube to make polarization sensitive holograms that are finally detected by the Polarized Array Sensor (PAS) from LUCID Vision Labs (PHX050S-PC). The PAS with micro-analyzers allows to record mosaicked holograms corresponding to 4 linearly polarization states, parallel, perpendicular, as well as at 45 • and 135 • polarization orientations, without the need for mechanical rotation of the analyzer. Lenses L 3 -L 4 are to be chosen so as to properly satisfy Nyquist sampling for each sub-image on the PAS.

Automation of the acquisitions have been realized using a homemade C++ acquisition software. This software uses the following libraries: eBUS SDK from Pleora Technologies [START_REF]eBUS SDK[END_REF] for managing hologram acquisitions via the PAS, and a Labjack (Labjack U3 DAC) [START_REF]Exodriver for labjack u3[END_REF] to control the FSM300 tip-tilt mirror, which delivers controlling signal to the FSM driver. This acquisition software also allows for quick visualization of the illumination scanning, as well as of the back focal plane of the detection objective, to ensure optimal scanning and detection. A GPU reconstruction code allows for a fast reconstruction and previsualization.

It should be noted acquisition multiplexing could also have been performed using a conventional sensor, considering two orthogonally polarized reference beams generated using a polarized beam splitter [START_REF] Verrier | Measuring enhanced optical correlations induced by transmission open channels in a slab geometry[END_REF] or a Wollaston prism [START_REF] Lee | Spatially multiplexed dielectric tensor tomography[END_REF]. This approach makes it possible to maintain the image field of view, with an imaging resolution limited by the microscope objective numerical aperture. To obtain the same results with a PAS, acquired images should be twice as big as the one obtained in a Fourier multiplexing framework [START_REF] Lee | Spatially multiplexed dielectric tensor tomography[END_REF]. Nevertheless, the choice of a PAS opens new outlooks such as Mueller imaging, or DIC imaging [START_REF] Verrier | 3d differential interference contrast microscopy using polarisation-sensitive tomographic diffraction microscopy[END_REF] with no modification of the experimental configuration, thus improving the versatility of the proposed arrangement. It also translates into a simpler experimental arrangement, as only the sensor is modified compared to a regular TDM experiment.

Experimental Results

Hologram pre-processing

Once mosaicked polarization sensitive holograms are recorded by the PSA, pixel-by-pixel demosaicking, and image interpolation are performed. Note that by making use of the rotation matrix, each interpolated hologram is then decomposed into its respective polarization. Same procedure applies for the holograms that are acquired without the sample (blank acquisitions) but only with the background medium, which are used as references to correct the holograms with the procedure presented in [START_REF] Liu | Tomographic diffractive microscopy and multiview profilometry with flexible aberration correction[END_REF]. Afterwards, each decomposed hologram is processed independently, similar to the scalar TDM reconstruction process.

Potato starch granules and silica beads

To validate our PS-TDM system, we have prepared a sample made of potato starch granules, which are well-investigated birefringent micro-objects [START_REF] Ortega-Ojeda | Interactions between amylose, native potato, hydrophobically modified potato and high amylopectin potato starches[END_REF], mixed with 5 𝜇m isotropic silica beads, embedded in an isotropic background medium (Eukitt from Sigma Aldrich) with refractive index 𝑛 0 = 1.49, and sandwiched between two cover glasses of thickness 170 𝜇m. To balance between image acquisition/reconstruction speed and image quality, the sample is illuminated with 400 angles using the optimized angular sample scanning scheme, 3D-UDHS from [START_REF] Taddese | Optimizing sample illumination scanning in transmission tomographic diffractive microscopy[END_REF][START_REF] Taddese | Optimizing sample illumination scanning for reflection and 4pi tomographic diffractive microscopy[END_REF]. Image demosaicking, then interpolations were conducted to retrieve the complex field corresponding to 𝑥 and 𝑦 polarization. Note that interpolation is critical, since a PAS introduces a phase shift and spatial-shear between two orthogonal polarized light. This effect can compromise PS-TDM image reconstruction, but it can be on the contrary beneficial in some cases, leading to a differential interference contrast (DIC) looking-like effect (see [START_REF] Verrier | 3d differential interference contrast microscopy using polarisation-sensitive tomographic diffraction microscopy[END_REF]).

Araneus diadematus silk fiber

So far, we illustrated benefits of our approach considering potato starch granules mixed with non birefringent beads, proving the ability to separate isotropic from anisotropic sample. In this subsection, we analyze the polarimetric information of a spider silk fiber. We selected Araneus diadematus silk fiber from dragline (radial line) of the web [START_REF] Römer | The elaborate structure of spider silk: structure and function of a natural high performance fiber[END_REF]. Such a sample is much smaller (1𝜇𝑚 or less) than potato starches, and known for exhibiting anisotropic structure at micrometer scale [START_REF] Little | Investigating the transverse optical structure of spider silk micro-fibers using quantitative optical microscopy[END_REF], which should be revealed using polarization sensitive TDM. The sample is immersed in Eukitt, and imaged with the same experimental conditions as that of Figs. 5, and 6

Obtained results are shown in Fig. 7. Such sample is thin and small enough to comply with both Born and weakly birefringent hypotheses. The reconstructed RI tensor image, shown in Fig. 7(a), emphasizes the anisotropic nature of the spider silk, which is proven by the presence of signal in the off-diagonal RI tensor components (𝑛 𝑥 𝑦 and 𝑛 𝑦 𝑥 ).

Further analysis is brought by computing the average RI 𝑛 𝑎𝑣𝑔 , the birefringence of the silk Δ𝑛, and the fast axis orientation 𝜃 that are given by Eqs. ( 24), [START_REF] Lauer | New approach to optical diffraction tomography yielding a vector equation of diffraction tomography and a novel tomographic microscope[END_REF], and [START_REF] Saba | Polarization-sensitive optical diffraction tomography[END_REF], respectively. Obtained results are illustrated in Figs. 7(b-d). Here, it should be noted that the rapid RI oscillations around the silk (from 1.45 to 1.5) are associated with a commonly known Rytov reconstruction artifact.

However, the average RI seems to be constant throughout the silk. Conversely, non-homogeneous structure is observed in the birefringence map of the silk fiber, as pointed out by a white arrow (see 7(c)).

We finally can mention the fast-axis orientation map (Fig. 7(d)) showing that orientation of the silk core is along the fiber axis, consistent with structural data [START_REF] Yazawa | Simultaneous effect of strain rate and humidity on the structure and mechanical behavior of spider silk[END_REF], but seems to differ from the orientation of its shell (see Fig. 2 in [START_REF] Strassburg | Functionalization of biopolymer fibers with magnetic nanoparticles[END_REF]). Interestingly, such orientation image results from computations, which eliminate reconstruction artifacts as seen on RI images. This could be of interest to obtain better sample discrimination from the background, but also possibly better resolved images.

Pinna nobilis oyster shell

To further demonstrate the performance of our PS-TDM, Pinna nobilis oyster shells are investigated.

Oyster shells are one example of a mesocrystalline material made out of calcium carbonate.

Structurally, the shell is made from single crystal calcite prisms, each prism being surrounded by organic walls (thin layer of amorphous calcium carbonate). Such structures are widely studied to hypothesize realistic biomineralization scenarios [START_REF] Ferrand | Quantitative imaging of anisotropic material properties with vectorial ptychography[END_REF][START_REF] Cuif | Biominerals and fossils through time[END_REF]. Such crystals are well suited for observation using 2-D microscopy (for example by 2-D vectorial ptychography [START_REF] Ferrand | Quantitative imaging of anisotropic material properties with vectorial ptychography[END_REF]), as they are of relatively adequate flatness and moderate thickness. Since they are also made of semi-transparent anisotropic material, they are well suited to be imaged by our PS-TDM.

The Pinna nobilis shell pieces were sampled in October 2013 from specimen living in the Mediterranean Sea (Giens, France). No animal was captured or killed. Small, millimeter-sized pieces were removed in situ from the spines covering the outer shell, in the vicinity of the shell growth margin, and preserved in a 70% ethanol/water solution (see [START_REF] Duboisset | Amorphous-to-crystal transition in the layer-by-layer growth of bivalve shell prisms[END_REF] for detailed explanations).

The reconstructed RI tensor components of the shell depicted in Fig. 8 These results are also confirmed by the computed average RI map, as shown in Fig. 8(c), as well as the fast-axis-orientation depicted Fig. 8(d).

Similar results have been obtained using optical vectorial ptychography (see Fig. 5 in [START_REF] Duboisset | Amorphous-to-crystal transition in the layer-by-layer growth of bivalve shell prisms[END_REF]), but, up-to-now, at lower magnification/resolution than with our system, which uses high-numerical aperture objectives, the drawback being a reduced field of view compared to [START_REF] Duboisset | Amorphous-to-crystal transition in the layer-by-layer growth of bivalve shell prisms[END_REF].

Conclusion

Experimental realization of quantitative 3-D polarimetric imaging approaches to image samples exhibiting low birefringence, such as biological samples, necessitates the development of highlysensitive optical systems. We have developed a simplified semi-automated high-resolution PS-TDM, which, thanks to the use of a PAS detector, multiplexes acquisitions, gaining a factor 2 in data collection. Full automation of the PS-TDM system could be done either by using motorized rotation stages or liquid crystal retarders in the polarization generation. Validation of our inversion algorithm based on Rytov approximation is done using synthetic data. The ability of our PS-TDM experimental setup to image birefringent specimens is validated using first a test samples made from birefringent and non-birefringent objects, then investigating a much smaller object such a spider silk fiber, and then tested onto biomineralized microcrystals.

Further analysis of our PS-TDM performances will require the use of calibrated samples.

Calibrated birefringent samples that are made for 3-D quantitative polarimetric imaging system, include liquid crystals [START_REF] Shin | Tomographic measurement of dielectric tensors at optical frequency[END_REF] and birefringent resolution targets such as for example those provided by Thorlabs (e.g. R2L2S1B). However, such samples are comprised of strong birefringence structures and may not be suitable for our application as they do not comply with our assumption of approximated Jones formalism approach. A highly-scattering 3-D sample with known RI deposited on silicon substrate is used in full polarization reflection TDM [START_REF] Zhang | Full-polarized tomographic diffraction microscopy achieves a resolution about one-fourth of the wavelength[END_REF][START_REF] Zhang | Far-field diffraction microscopy at 𝜆/10 resolution[END_REF]. In these works, the authors showed the advantage of full polarimetric measurements to improve the reconstructed image resolution, which has been done by combining scattered field information at multiple polarizations.

In general, fabricating such samples is not a simple task, as the phantoms would have to be weakly scattering, as well as weakly birefringent. For instance, 3-D printed phantoms that are fabricated using two-photon polymerization [START_REF] Krauze | 3d scattering microphantom sample to assess quantitative accuracy in tomographic phase microscopy techniques[END_REF][START_REF] Ziemczonok | 3d-printed biological cell phantom for testing 3d quantitative phase imaging systems[END_REF], used to characterize 3-D quantitative phase imaging techniques are presently limited to scalar system. It would be of great interest in the domain to have real 3-D birefringent test patterns to simultaneously be able to measure spatial resolution as well as RI sensitivity.

Our system could be of interest to study biological samples exhibiting birefringence at the cellular or sub-cellular level, such as spindle fibers, made of bundled microtubules, or actin distribution. Stress-induced birefringence at microscopic scale could also be investigated, such as in diatom or radiolarian shells, microcrystals, textile fibers, combustion residues. As possible extension, spectrally resolved, polarization-sensitive TDM should allow for a much-improved chemical sensitivity, which could be useful to identify micro-and nano-plastics, free-floating, and possibly even when internalized by cells.

circularly polarized beam, right and left circularly polarized light (LC and RC), described by their Jones vector,

𝐸 𝑖 = 1 √ 2       1 ±𝑖       . (27)
In the sample arm, the output field 𝑂 𝑖 𝑗 is expressed as cascaded Jones components:

𝑂 𝑖 𝑗 = 𝐽 𝜃 𝑎𝑛𝑙 𝐽 𝑠 𝐸 𝑖 , (28)
where 𝐽 𝜃 𝑎𝑛𝑙 is the Jones matrix of analyzer before the detector, and 𝐽 𝑠 is the anisotropic sample's complex Jones matrix:

𝐽 𝑠 =       𝐽 𝑥 𝑥 𝐽 𝑥 𝑦 𝐽 𝑦 𝑥 𝐽 𝑦𝑦       , 𝐽 𝜃 𝑎𝑛𝑙 =       0 0 cos 𝜃 sin 𝜃       (29)
Now using the 0 • and 90 • polarization, 𝑂 𝑖 𝑗 can be computed as:

𝑂 𝑖 𝑗 = 1 √ 2       0 0
𝐶 (𝐽 𝑥 𝑥 -𝑖𝐽 𝑥 𝑦 ) + 𝑆(𝐽 𝑦 𝑥 -𝑖𝐽 𝑦𝑦 ) 𝐶 (𝐽 𝑥 𝑥 + 𝑖𝐽 𝑥 𝑦 ) + 𝑆(𝐽 𝑦 𝑥 + 𝑖𝐽 𝑦𝑦 ),

      ( 30
)
where 𝐶 is cos 𝜃 and 𝑆 is sin 𝜃.

Similarly, the reference beam, 𝑅 𝑖 𝑗 , can be represented by 𝐽 𝜃 𝑎𝑛𝑙 and reference beam's Jones vector 𝐸 𝑟 𝑒 𝑓 :

𝑅 𝑖 𝑗 = 𝐽 𝜃 𝑎𝑛𝑙 I𝐸 𝑟 𝑒 𝑓 = 1 √ 2       0 0 𝐶 + 𝑖𝑆 𝐶 + 𝑖𝑆       (31)
It could also be possible to accounts for light source fluctuation or optical system inaccuracies in the reference arm. The best way to make sure that we are measuring only the sample contribution is to take a blank acquisition and subtract from both object and reference fields. Now, the polarization sensitive hologram that we record using the polarization sensitive camera is an interferogram between 𝑂 𝑖 𝑗 and 𝑅 𝑖 𝑗 . Recalling Eq. ( 13) and extracting the first diffraction order leads to: Extracting the +1 order:

𝑂 𝑖 𝑗 𝑅 * 𝑖 𝑗 = 1

  • polarization • 𝐸 𝑜 1𝑦 : right-circular illumination and detection at 90 • polarization • 𝐸 𝑜 2𝑥 : left-circular illumination and detection at 0 • polarization • 𝐸 𝑜 2𝑦 : left-circular illumination and detection at 90 • polarization

Fig. 1 .

 1 Fig. 1. Schematics representing the incident field's polarization for an oblique illumination of wave-vector k 𝑖 : 𝑧 -𝑎𝑥𝑖𝑠 is the normal incident beam direction. E 𝑥 𝑦 𝑖 is the incident field E 𝑖 projected at 𝑥𝑦 plane. 𝜃 and 𝜙 are respectively the polar and azimuthal angle.

Fig. 2 .

 2 Fig. 2. Schematics representation of PS-TDM demosaicking and interpolation process

Fig. 3 .

 3 Fig. 3. Simulation results of a spherical anisotropic bead, with two principal RIs (𝛿𝑛 1 and 𝛿𝑛 2 ). The images are the central slices of the 3-D object. The RI profile plot is taken from 𝛿𝑛 1 along the 𝑦-axis. Scale bar is 2 𝜇m.

Fig. 4 .

 4 Fig. 4. Experimental setup of high-resolution polarization sensitive TDM based on off-axis digital holographic microscope; CL: collimating lens, M 1,2 : Mirrors, PC 1,2 : Polarizing components, BS: beam spliter, PS-camera: polarization sensitive Camera, L 1,2,3,4 : Lenses, and RC: recombining component.

Fig. 5 .

 5 Fig. 5. Computed 2-D Jones tensor maps: (left) amplitude, and (right) phase values of potato starch granules (large spherical objects, visible in all components), and silica beads (not visible in the diagonal components). Scale bar is 10 𝜇m.

Figure 5

 5 Figure 5 shows the processed 2-D phase and amplitude maps of the Jones matrix elements of potato starch granules (anisotropic object) and silica beads (isotropic object) in the same field of view. Note that the isotropic silica beads are not visible in the off-diagonal Jones components. Unlike the silica beads, the potato starch granules structures are visible in all Jones components.

Fig. 6 .Figure 6

 66 Fig. 6. Reconstructed RI tensor components of potato starch granules, and silica beads in the same field of view; scale bar is 10𝜇𝑚. (See Visualization 1 for details)

Fig. 7 .

 7 Fig. 7. Reconstructed 3-D polarimetric information of Araneus diadematus spider silk. (a) RI tensor components, (b) birefringence, (c) average RI and (d) fast-axis orientation. Scale bar is 10𝜇𝑚 (a) and 2𝜇𝑚 (b-d).

6 .

 6 (a) clearly show the microstructure of the shell relative to the background medium. We observe an inhomogeneous structure of the prisms, and the RI tensor component of one prism varies from another. Note also the Region of Interest (ROI) depicted with white rectangles along the 𝑥 -𝑧 cross-section, where we observe bright colors (showing a positive RI contrast) relative to the background medium in both 𝑛 𝑥 𝑦 and 𝑛 𝑦 𝑥 RI tensor components. The same ROI shows dark structures in 𝑛 𝑥 𝑥 and 𝑛 𝑦𝑦 (showing a negative RI contrast). Such phenomena is also observed on potato starch granules, see Figure Other 3-D optical parameters are also computed to gain more insight about the polarimetric properties of such sample. The birefringence map in Fig. 8(b) reveals an optically inhomogeneous structure for the various prisms. We also observe smooth transition between adjacent prisms.

Fig. 8 .

 8 Fig. 8. Reconstructed 3-D polarimetric information of Pinna nobilis oyster shell. (a) RI tensor components, (b) birefringence, (c) average RI and (d) fast-axis orientation. Scale bar is 20 𝜇𝑚.

  2 [𝐶 (𝐽 𝑥 𝑥 -𝑖𝐽 𝑥 𝑦 ) + 𝑆(𝐽 𝑦 𝑥 -𝑖𝐽 𝑦𝑦 ) + 𝐶 (𝐽 𝑥 𝑥 + 𝑖𝐽 𝑥 𝑦 ) + 𝑆(𝐽 𝑦 𝑥 + 𝑖𝐽 𝑦𝑦 )] [𝐶 + 𝑖𝑆] (32)So 2 acquisitions are enough to compute all 4 Jones matrix components since the PS-Camera has mosaic images, containing both 0 and 90 degree polarization in a single acquisition. Thus, it provides a linear relationship between the measured complex scattered field and the birefringent sample's Jones tensor components.

	𝐽 𝑥 𝑥 = 𝐸 𝑜 11 + 𝐸 𝑜 12	(33)
	𝐽 𝑥 𝑦 = 𝑖[𝐸 𝑜 11 -𝐸 𝑜 12 ]	(34)

𝐽 𝑦 𝑥 = 𝑖[𝐸 𝑜 22 -𝐸 𝑜 21 ]

(35)

𝐽 𝑦𝑦 = 𝐸 𝑜 21 + 𝐸 𝑜

Acknowledgments

Julien Duboisset and Olivier Fauvarque are warmly acknowledged for their help during sample preparation of the Pinna nobilis sample. We also would like to thank Virginie Chamard and Patrick Ferrand for their fruitful discussions on the theoretical background on light polarization and polarization imaging.

Funding

This work is funded by Agence Nationale de la Recherche (ANR) under contracts ANR-18-CE45-0010 and ANR-19-CE42-0004 and by the Région Grand Est under contracts FRCR 18P-07855 and FRCR 19P-10656.

Data Availability Statement

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Disclosures

The authors declare no conflicts of interest.

Appendix

In Jones formalism, when light propagates through several optical elements, the variations of the electric field vector from beginning to the end can be expressed as the product of each associated Jones matrices or Jones vectors. Light in both reference and illumination arms is changed into