Photoacoustic fluctuation imaging: theory and application to blood flow imaging

Photoacoustic fluctuation imaging, which exploits randomness in photoacoustic generation, provides enhanced images in terms of resolution and visibility, as compared to conventional photoacoustic images. While a few experimental demonstrations of photoacoustic fluctuation imaging have been reported, it has to date not been described theoretically. In the first part of this work, we propose a theory relevant to fluctuations induced either by random illumination patterns or by random distributions of absorbing particles. The theoretical predictions are validated by Monte Carlo finite-difference time-domain simulations of photoacoustic generation in random particle media. We provide a physical insight into why visibility artefacts are absent from second-order fluctuation images. In the second part, we demonstrate experimentally that harnessing randomness induced by the flow of red blood cells produce photoacoustic fluctuation images free of visibility artefacts. As a first proof of concept, we obtain two-dimensional images of blood vessel phantoms. Photoacoustic fluctuation imaging is finally applied in vivo to obtain 3D images of the vascularization in a chicken embryo.


INTRODUCTION
Photoacoustic (PA) imaging is a widely-spread biomedical imaging modality that makes use of ultrasound waves resulting from light absorption [1]. In acoustic-resolution PA imaging, arrays of ultrasound detectors are usually used to record PA signals generated by absorbing structures. PA reconstruction algorithms (such as delay-and-sum or backpropagation algorithms) are then used to provide images from a set of PA signals. PA images may, however, present artefacts that arise from coherence of PA waves and characteristics of the detection system such as geometry and frequency bandwidth. Limited-view artefacts occur when strongly anisotropic absorbing structures (such as blood vessels) coherently emit PA waves in some preferential directions, and these waves can not be detected due to the finite aperture and/or directivity of the ultrasound detectors. In such a case, some parts of the absorbing structure are not visible on the reconstructed PA image. Visibility artefacts also occur with resonant ultrasound detectors, which filter out low-frequency components of PA waves emitted by large absorbers (large as compared to the detection wavelength range). For instance, for large blood vessels, only the vessel boundaries may be visible in PA images.
Various solutions have been proposed to overcome the limited-view problem with finite-aperture detectors. One approach consists in enhancing the detection aperture by performing multiple PA acquisitions at different angles. This can be done by rotating the ultrasound probe around the object [2] or by spinning the object itself [3]. The detection aperture can also be augmented by placing acoustic reflectors at edges of the imaging zone [4,5]. Alternatively, the imaged object and the ultrasound probe can be placed inside a reverberating cavity [6][7][8]. This allows detecting PA waves emitted in all possible directions even when only a single-element probe is used [8].
All the mentioned techniques require that the whole range of imaging angles is accessible for detection. However, this is not necessarily feasible in real clinical environment. In addition, mechanical scanning of the detectors or sample increases the acquisition time and, therefore, degrades the temporal resolution. Methods that do not require detecting waves emitted in all directions have also been proposed. Sub-wavelength sparsely distributed absorbing particles act individually as isotropic PA sources. This can been exploited to remove limited-view artefacts: Dean-Ben et al. demonstrated that visibility in the limitedview scenario can be improved by means of the localization approach [9], or by using a nonlinear combination of tomographic reconstructions representing sparsely distributed moving particles [10]. Unless sparsely distributed contrast agents are used [11], these techniques remain limited in practice as blood consists of a concentrated suspension of red blood cells. One more approach consists in heating tissue locally with a focused ultrasound beam and thus generating artificial PA sources via the temperature dependence of the Gruneisen coefficient [12]. By scanning the focused ultrasound beam across the sample and accumulating the resulting PA images, the whole object is reconstructed. However, this approach is limited in practice by the safety thresholds of focused ultrasound. Similar to all other aforementioned approaches this method requires a very long acquisition time leading to a low temporal resolution. Recently, deep-learning based approaches have also been proposed as a way to palliate visibility artefacts from limited-view detection in photoacoustic imaging [13][14][15]. While powerful, deep-learning based imaging methods also have limitations. These include the availability of a training set and the performances inherently limited to unseen samples "close" enough to those of the training set, which is a severe practical limitation.
An approach based on multiple-speckle illumination was proposed by Gateau et al. [16]. In this experimental work, a random intensity distribution of speckle patterns, which changed from pulse to pulse, induced fluctuations in each pixel of the corresponding series of PA images. It was demonstrated experimentally with free-space-generated series of speckle patterns that a second-order fluctuation image provided a faithful representation of the absorbing distribution, free of limited-view and limited-bandwidth artefacts. However, exploiting optical speckle illumination for imaging tissue at depth turns out to be challenging since the fluctuation signal is expected to be very small when the speckle grain size is of the order of the optical wavelength. As a consequence, this technique has so far never been demonstrated in a close-to-clinical environment. Moreover, the work by Gateau et al. [16] provided no clear theoretical explanation for the visibility enhancement in fluctuation imaging with multiple speckle illumination. Since that work, there have been other studies based on PA signals fluctuations. In particular, such fluctuations were deployed to obtain super-resolution in PA imaging [17][18][19][20]. Most generally, PA fluctuation imaging utilizes some randomness in PA generation to provide enhanced images as compared to conventional PA imaging. While it has been successfully demonstrated experimentally for super-resolution imaging (induced either from multiple-speckle illumination [17][18][19] or from random distributions of absorbing particles [20]) and for visibility enhancement (with multiple-speckle illumination [16]), there is, however, no comprehensive theoretical description of fluctuation PA imaging to date.
In the first part of this work, we propose a unified theoretical framework relevant to both fluctuations induced by random illumination patterns and fluctuations induced by random distributions of absorbing particles. The theoretical considerations are validated in Monte-Carlo finite-difference time-domain (FDTD) simulations of PA generation in a random medium. In particular, we explain why visibility artefacts are absent in 2 nd -order fluctuation images. This property turns out to be completely independent of the origin of fluctuations. In addition, our theoretical results provide a quantitative comparison between fluctuations from multiple-speckle illumination and fluctuations from random distributions of absorbing particles. In the second part, we demonstrate experimentally that fluctuations induced by a blood flow can be exploited under physiological conditions to produce PA fluctuation images free of visibility artefacts, as opposed to conventional PA imaging. Two-dimensional images of vessel-mimicking phantoms flown with blood at physiological concentrations are obtained as a first proof-of-concept demonstration. Finally, the approach is validated in vivo by performing 3D imaging of a chicken embryo.

THEORY
In this section, we propose a unified theoretical framework for PA fluctuation imaging applicable to both fluctuations induced by random illumination patterns and fluctuations induced by random distributions of absorbing particles. We first provide theoretical expressions for the mean PA image and the fluctuation PA image. These expressions are then used to explain why the fluctuation images do not possess the visibility artefacts observed on the mean images (equivalent to conventional PA images). Finally, we discuss the dependence of the amplitude of the fluctuation image on characteristic statistical properties (mean, variance, characteristic size) of the random process inducing fluctuations (illumination patterns or distribution of absorbing particles). In particular, we compare the case of randomly-generated multiple-speckle illumination to that of randomly-distributed red blood cells, and discuss quantitatively why flowing red blood cells lead to PA fluctuations much larger than those produced by multiple speckle illumination deep in tissue.

A.1. Theoretical framework
PA imaging provides images of the absorbed energy density α that we describe in this work as: In the expression above, F 0 has the dimension of a light fluence, µ 0 has the dimension of an absorption coefficient; f (r) is a dimensionless binary function that describes the structure containing optical absorbers (for instance blood vasculature with red blood cells) : f (r) = 1 inside the imaged structure, f (r) = 0) outside; g k (r) is the k th realization of the random process inducing fluctuations: it may either describe a spatial distribution of light intensity or a spatial distribution of absorbers. For random illumination patterns, such as speckle patterns [16], F 0 × g k (r) (with g k being a continuous real-valued function) is the associated random fluence corresponding to each laser shot k, and the absorption coefficient µ 0 is considered homogenous in the absorbing structure f (r). For random distributions of absorbing particles, such as red blood cells, µ 0 × g k (r) is the associated heterogeneous absorption coefficient for each laser shot k while g k (r) is a binary function equal to 1 inside the particles, and 0 elsewhere. From the most general perspective, g k (r) × f (r) thus describes the spatial distribution of absorbed energy. We assume that g k (r) is a spatially stationary and isotropic random process. Its relevant statistical properties in the context of our work are the mean η =< g k (r) > k and the variance σ 2 g =< (g k (r) − η) 2 > k as first-order statistical properties and its autocovariance as a second-order statistical property. We note that C(0) = σ 2 g . We further assume that the normalized autocovariance C( r 1 − r 2 )/C(0) is a function with a finite volume V g and characteristic width D g (with D n g = V g , n=2 or n=3 depending on the relevant dimensionality): For multiple-speckle illumination, V g (resp. D g ) is of the order of the characteristic volume (resp. size) of a speckle grain. For random distributions of mono-dispersed absorbing particles, V g (resp. D g ) is of the order of the particle characteristic volume (resp. size).
We assume that each PA reconstruction A k (r) corresponding to laser shot k can be expressed as the convolution between α k (r) and the point spread function (PSF) h(r) of the imaging system: where Γ is the Grüneisen parameter. The convolution in Eq. 3 implies a space-invariant PSF, which is a valid assumption within the small fields of view that we deal with in this work. We consider the most general framework where the PSF h(r) is either a real-valued function or a complex-valued function. A real-valued PSF corresponds to reconstruction based on realvalued RF signals, whereas a complex-valued PSF corresponds to reconstruction based on complex-valued signals. Complex-valued reconstruction provides a straightforward way to remove radio-frequency oscillations by using the modulus of the complex-valued reconstruction as the final image, and is a widely used approach in medical ultrasound imaging. In the context of photoacoustic imaging, an illustration of the differences between real-valued only images and complexmodulus images may be found in a previous publication from our group [20].

A.2. Mean photoacoustic image
We consider the mean PA image defined by E[A](r) =< A k (r) > k , where < . > k denotes an ensemble average (average over laser shots in practice). As η =< g k (r) > k , E[A](r) can be straightforwardly expressed as: In the case of random distributions of absorbing particles, η =< g k (r) > k corresponds to the volume fraction of absorbers, and ηµ 0 is the average absorption coefficient of the imaged structure. In the case of random illumination patterns and a homogeneously absorbing structure, < g k (r)F 0 > k = ηF 0 corresponds to the average light fluence inside the imaged structure, and µ 0 is the absorption coefficient of the homogeneously absorbing imaged structure. Most generally, Eq. (4) corresponds to conventional PA reconstruction of a homogeneously absorbing structure described by f (r). Eq. (4) is independent of the real-valued or complex-valued nature of h(r). For real-valued reconstruction, E[A](r) directly represents the mean PA image, whereas for complex-valued reconstruction, the mean PA image is rather represented by the modulus |E[A](r)|.

A.3. Photoacoustic fluctuation image
In the context of this work, we define the PA fluctuation image as a second-order fluctuation image, defined as the square root of the variance image (r)] * > k , for both real-valued and complex-valued reconstructions. Higher-order PA fluctuation images may also be defined [20], but are beyond the scope of this work. The variance image σ 2 [A](r) can be expressed by use of the autocovariance function C(r 1 , r 2 ): As the fundamental and major assumption of our work, we now assume that the characteristic size D g of the random distribution is much smaller than the shortest characteristic sizes of both h(r) and f (r). In the context of random distributions of red blood cells, D g is of the order of 10 µm. In the context of speckle illumination, D g is the typical size of a speckle grain and is therefore not larger than 1 µm for visible or near-infrared light. Moreover, the shortest characteristic length scale of h(r) for resonant transducers is the wavelength λ at the central frequency.
Our assumption is therefore fully verified for vessels at least several tens of µm in diameter and central frequencies up to several tens of MHz (acoustic wavelengths being not smaller than several tens of µm), which corresponds to the most of practical situations of interest. Under this assumption, C(r 1 , r 2 ) may be replaced in the integral by This leads to the following final expression for the PA fluctuation image, valid for both real-valued and complex-valued PSF: While it has already been shown by our group that the PA fluctuation image could be written as a convolution involving the square of the PSF [17,20], Eq. (5) for the first time explicitly links the amplitude of the fluctuation image to statistical properties of the random process at stake (through its mean η, its standard deviation σ 2 g and the volume of its normalized autocovariance V g ). In particular, Eq. (5) remains valid when the characteristic size of f (r) is smaller than that of h(r), which is of interest for super-resolution imaging: as opposed to our previous works on fluctuation-based super-resolution imaging [17,20], we here provide a quantitative prediction of the second-order superresolution photoacoustic fluctuation image, which is relevant for quantitative analysis of such images. In this work, we however focus on the consequences of Eq. (5) on the visibility problem.

B. Consequences on visibility artefacts
Eq. (5 ) indicates that the fluctuation image involves a convolution with the square of the PSF. This is a well-known result from the Super-resolution Optical Fluctuation Imaging (SOFI) approach initially proposed for super-resolution imaging relying on blinking fluorophores. The SOFI approach permits obtaining super-resolved images based on the fact that the PSF to the 2 nd power is sharper than the PSF itself [21]. Our group adapted the SOFI approach to super-resolution PA imaging, with fluctuations induced from either multiple-speckle illumination [17] or moving absorbers [20]. Although this convolution with the squared PSF resulting in the 2 nd −order fluctuation image has been known both for optical and PA imaging, its effect was analyzed only from the resolution perspective [17,20]. However, as we will now discuss, it also has a direct impact in the context of both the limited-view and resonant-transducer problems. To illustrate our discussion, we first use two-dimensional (2D) numerical simulations showing how fluctuation imaging removes visibility artefacts. The numerical simulations were carried out with the 2D version of a finite-difference time-domain (FDTD) freely-available code developed by our group [22]. The simulations output consisted of a set of RF signals detected by a linear transducer array similar to that used in the experiments (central frequency 15 MHz, bandwidth 80%, N = 128 elements, pitch 110 µm). The input of each simulation consisted of a source image defined on a cartesian grid with a spatial step of 2.5 µm.

B.1. Illustration with simulation results
Three different structures f (r) were used to illustrate how PA fluctuation imaging removes artefacts of conventional PA imaging, as illustrated in Fig. 1 (top row). Similar structures used in the phantoms experiments are represented in Fig.4. The simulation input for each realization of the random process consisted of a random source image g k (r) × f (r). The distribution g k (r) was represented by an image of super-pixels (group of q × q pixels on the simulation grid), that was randomly assigned binary values 0 or 1 following a Bernouilli distribution with mean value η. The simulations permitted producing the full range of η =< g k > = 0 -100%. In this model, η represents the surface fraction of the absorbed light and can be varied from 0 to 100 %. Fig.1 shows some example results obtained in our Monte-Carlo simulations. To obtain these results, the surface fraction of absorbed energy was set to η = 50 % and the size of each individual light absorbing patch was 20 µm × 20 µm (16 × 16 superpixel). Fig.1 shows the modulus images obtained from complexvalued signals. Fig.1 (column a) illustrates the limited-view problem: conventional PA imaging (middle row) only shows features parallel to the probe, while PA fluctuation imaging (bottom row) renders the full structure correctly. Fig.1 (column b) illustrates the limited-bandwidth problem: conventional PA imaging (middle row) only shows the upper and lower boundaries (high spatial frequencies) of the imaged object and does not reveal its inner content (low spatial frequencies), while PA fluctuation imaging (bottom row) correctly renders the full structure. Fig.1 (column c) illustrates both problems and confirms the outperformance of fluctuation imaging (bottom row). Such outperformance have first been demonstrated experimentally in our group several years ago in the specific context of multiplespeckle illumination [16], but without any clear theoretical explanation.
To illustrate the principle of the proposed approached, the simulations results presented in Fig.1 were obtained without any additional sources of fluctuations (such as detection noise or laser pulse energy fluctuations). Additional simulation results presented in the Supplementary Information (section 4) confirms that our approach also works very efficiently under noisy environment: pulse laser energy fluctuations may be compensated for via SVD filtering; the fluctuation of interest may be detected on top of a constant noise background provided that a sufficiently large number N of measurements are used to estimate the fluctuations image. Our simulation results with noise show that N may range in practice from a few units to a few hundred (see Fig. S3 of the Supplementary Information), depending on the ratio between the fluctuation of interest and the noise level (FNR, see section 3 of the Supplementary Information.

B.2. Mechanism of visibility enhancement in fluctuation imaging
Based on the expression of the fluctuation image involving a convolution with the squared PSF, we now explain why PA fluctuation imaging removes visibility artefacts. The PSF corresponding to the images shown in Fig.1 was computed in the same simulation setting by placing a point source in the center of the field of view (FOV) at the distance of z = 15 mm from the probe. We note h r (r) the real-valued PSF and h c (r) its complex-valued counterpart for the considered 1D linear transducer array (aligned along the horizontal direction as in Fig.4). Fig.2 shows the PSF h(r) and its square modulus |h| 2 (r) along and h 2 (r) for limited-view photoacoustic imaging with resonant transducer, for both the real-valued PSF (h r (r)) and the complex-valued PSF (h c (r)). Top row: spatial domain. Bottom row: Fourier domain. All functions are normalized to one. The PSF was derived from FDTD-computed signals, with the same 1D linear array as that used for Fig. 1, with a point source 15 mm away from the transducer.
with the corresponding spatial Fourier transforms.
The bipolar nature of h r (r) and h c (r) for resonant transducers leads to low spatial frequency components of the objects being filtered out, which can be explained with showing the PSF in the Fourier space. Moreover, the anisotropy of h r (r) and h c (r) leads to a selective detection of only spatial frequencies within the numerical aperture of the transducer. Previous analyses of |h| 2 (r) in the context of superresolution PA imaging used its extended support in the Fourier space as compared to that of h(r), to explain the obtained superresolution (see Supplementary Information in [17,20]). However, it also follows from the support of h 2 r (r) and |h c | 2 (r) in the Fourier space (see Figs.2.(b 2 ) and Fig.2.(d 2 )) that imaging with the squared modulus of the PSF preserves low spatial frequencies in all directions, as opposed to h(r). As a consequence, the fluctuation image does not suffer from the limited-view and resonant-bandwidth artefacts observed in conventional imaging.  1 )). Interestingly, the conclusions above drawn from the analysis in the Fourier space may also be obtained in the physical space by considering the difference between coherent and incoherent summation. Since h(r) is a bipolar function, its convolution with absorbing objects which are always positive may lead to destructive interference, which in practice occurs when dealing with large objects or those that emit PA waves away from the detectors. As opposed to h(r), h 2 (r) is positive only and can therefore not lead to destructive interference. In other words, conventional imaging involves coherent summation of bipolar PA waves from positive sources, that may lead to destructive interference, whereas fluctuation imaging involves summation of positive fluctuations that can never interfere destructively. One can draw a parallel with optics where amplitudes sum up coherently while intensities sum up incoherently.

C. Amplitude of the fluctuation image
The mechanism discussed and illustrated above with Monte-Carlo simulations and a simple random model is independent of the nature of fluctuations. These fluctuations may arise either from random heterogeneous illumination patterns (such as speckle patterns) or from random distributions of optical absorbers (such as red blood cells). However, the amplitude of the fluctuation image fundamentally depends on the statistical properties of the relevant random process (mean η, standard deviation σ 2 g , typical length scale D g ), as indicated by Eq.5. To illustrate our discussion, we consider the fluctuation amplitude expected from structures much larger than the PSF. In this case, f 2 (r) * |h| 2 (r) r |h| 2 (r )dr = ||h|| 2 , and the fluctuation amplitude is given by: , as a consequence of the resonant bandwidth artefact (bipolar nature of h(r)). In practice, PA fluctuation amplitude must be larger than other fluctuations such as the the detection noise to be measurable. The fluctuation amplitude is in part determined by the sensitivity of the system (quantified through ||h|| 2 ), and by the statistical properties of the random distribution of absorbed energy (η, σ g and V g ).
We first demonstrate by use of 2D numerical simulations that the fluctuation amplitude depends on η and V g = D 2 g as predicted by our theory. To this end, we shall use the simple random model already introduced in section B.1. Then we provide more specific expressions for the fluctuation amplitude in the case of multiple speckle illumination and in the case of random distributions of absorbing particles. We finally discuss why fluctuations from absorbing particles such as red blood cells are expected to be much larger than those from multiple speckle illumination deep in tissue. We recall that equations 5 and 6 are valid only for characteristic lengths D g much smaller than the shortest dimension of the PSF (about 50 µm as half a wavelength in the context of our work, see Fig. 2).

C.1. Validation with simulations
To validate our theory with simulations, we consider a simple model of a random distribution of absorbed energy for which V g (and thus D g ) is independent of η. To do so, we use the model already introduced in section B.1, for which g k is defined on a cartesian grid with pixel/voxel size D. Each pixel of the grid may acquire a random binary value 0 or 1 from a Bernoulli distribution with mean value of η. We recall that η defined as such corresponds to the average surface/volume fraction η =< g k (r) > k of the corresponding random medium. The variance for the Bernoulli distribution is given by σ 2 g = η × (1 − η). Moreover, it can be shown straightforwardly for this random medium that C(r , , with n being the space dimension (n=2 or n=3) and Λ D being the 1D unit triangular function with support 2 × D. From Eq.2, one obtains that V g = D n (i.e. D g = D), independently of η. The expression for the fluctuation amplitude for this random distribution of absorbed energy becomes : To validate the predictions of Eq.7, we ran 2D simulations (n = 2) for the disk structure (much larger than the PSF) with  Fig.3, the dependence of the fluctuation amplitude on η (for fixed D = 10 µm) and on the pixel size D (for fixed η = 50 %) measured from the FDTD simulations (N = 100 random realizations for each set of parameters) exactly matches the one following from Eq.7. In particular, the dependence on η follows η × (1 − η) and η = 50 % is the value that maximizes the fluctuation amplitude. As expected, the size dependence diverges from the theory when D approaches half a wavelength (50 µm).

C.2. Multiple speckle illumination
We consider here a homogeneously absorbing structure f (r) (constant absorption coefficient µ 0 ), illuminated by random multiple speckle illumination described by the fluence F k (r) = F 0 × g k (r). We further assume that the speckle is fully developed, so the statistical properties of the fluence distribution are well-known: the intensity/fluence follows an exponential distribution, with σ F =< F >, where < F >= ηF 0 is the average intensity/fluence. D g = D s is a measure of the speckle grain size, and is independent of the average fluence ηF 0 . V s = D n s is the speckle grain area (i.e. n = 2) or volume (i.e. n = 3) depending on the relevant imaging configuration. The fluctuation amplitude in this case is thus given by: While it has long been acknowledged that detecting fluctuations from speckle illumination required large enough speckle grains, Eq.8 provides for the first time a quantitative theoretical prediction of the expected fluctuation amplitude as a function of all relevant parameters, including the geometry of the imaged object and the PSF of the imaging system. For a fixed average fluence, Eq.8 shows that the fluctuation amplitude is proportional to the square root of the area/volume of the speckle grain (depending on the relevant dimensionality).

C.3. Random distribution of absorbing particles
We now consider randomness that arises from the sample rather than from the illumination. In this case, the fluence is constant (F 0 ), and the absorption distribution is described by µ k (r) = µ 0 × g k (r). More specifically, we consider randomness induced by random distributions of absorbing particles of a given characteristic volume V p (characteristic size D p = V 1/n p ) inside the structure to image f (r), such as red blood cells flowing in blood vessels. In this case, g k (r) is a binary function (g k (r) = 1 inside absorbing particles, g k (r) = 0 outside) and η =< g k (r) > k simply represents the average volume fraction of absorbers. For blood, µ 0 is the absorption coefficient of pure hemoglobin (i.e. of a single red blood cell), and ηµ 0 is the average absorption coefficient of the whole blood as a suspension of red blood cells. For g k (r) which is a binary random variable with η =< g k (r) > k (Bernoulli distribution), the variance is , regardless of the spatial structure of the random medium. In this case, the expression for the fluctuation image becomes: The characteristic volume V g (η), defined by Eq.2, is of the order of V p , but its exact value generally depends on the volume fraction η, especially at high volume fraction where the positions of neighbors may become correlated. The dependence of σ[A](r) to the volume of the absorber (for a fixed η) is exactly the same as the dependence to the size of the speckle grain: the larger the particles, the larger the amplitude of the fluctuation image, within the limit that the particle size should be small compared to the PSF. For a given particle size and shape, the value of η that maximizes the fluctuation amplitude is given by the maximum of η(1 − η)V g (η), whose exact form depends on V g (η). In the simplest case where V g is independent of η, our theory predicts that the fluctuation amplitude is proportional to η × (1 − η), which is maximized for η = 50 %. It is out of the scope of this work to investigate the dependance of V g (η) for various types of media. Nonetheless, we demonstrate in the Supplementary information that V g (η) may be expressed as a function of V p through the so-called packing factor: The packing factor W(η) has been used in the field of ultrasound imaging to study the non linear dependence of echogeneicity as a function of η in particles media [23][24][25][26], including blood [24,[26][27][28]. By use of the expression above, the PA fluctuation may thus also be expressed as which interestingly has the same dependence on η × W(η) as that reported for the backscattering coefficient [26], and agrees with recent experimental results [29] (see Supplementary Inf.)

C.4. Multiple speckle illumination versus particles fluctuations
As pointed out earlier, it is the ratio between the amplitude of the fluctuations of interest and the amplitude of other fluctuations (laser pulse energy fluctuations, thermal noise) that determines the ability to detect the relevant fluctuations experimentally. The theoretical expressions of the fluctuation images obtained for the case of multiple speckle illumination and random distributions of absorbers provide a quantitative comparison between the two situations. From Eqs. 8 and 9, one gets: where it is assumed that the average fluence and the absorption coefficient are identical in both cases. Applied to red blood cells (V p ∼ 100 µm 3 , η ∼ 50 %, W(50%) ∼ 16 % [26] and Supp. Inf) and and speckle grains in the multiple scattering regime deep in tissue (D s ∼ λ 2 ∼ 0.25 µm in the visible range), the relation above (with n = 3) indicates that fluctuations arising from red blood cells as absorbing particles are at least one order of magnitude larger than fluctuations expected from multiple speckle illumination.

D. Conclusions from the theory
In this part, we provided a unified theoretical framework to predict the amplitude of PA fluctuation images, whether fluctuations arise from random illuminations or from random distributions of absorbers. As a first result of our theoretical investigation, we explained why fluctuation PA imaging palliates the visibility artefacts emerging in conventional limited-view and resonant-bandwidth PA imaging. As a second and major result, we provided with Eq.5 a general expression of the fluctuation amplitude as a function of statistical properties of the random process inducing fluctuations. We consequently demonstrated that fluctuations from red blood cells randomly distributed in blood vessels are expected to be one or two orders or magnitude larger than those expected from multiple speckle illumination of blood vessels at depth in tissue. We recall that this conclusion holds because optical speckle grains at depth are much smaller (half the optical wavelength) than both the acoustic PSF (h(r)) and the vessel size ( f (r)), as the main assumption of our theory. These theoretical results quantitatively explain why it has recently been possible to experimentally detect fluctuations from flowing red blood cells (in the context of super-resolution imaging [20]), whereas the use of speckle-induced fluctuations imaging remain limited to proof-of-concept experiments where speckle grains could be made large enough via free-space propagation. As a general conclusion, our theoretical results suggest that fluctuations from flowing red blood cells at physiological concentrations (around 50 %) may be exploited to remove the visibility artefacts of conventional PA imaging. This has never been demonstrated experimentally so far, and is the objective of the second part of our work.

A. Phantoms experiments
As a first experimental demonstration of the proposed approach, we obtained cross-sectional images of the three different types of structure introduced in the simulation part. In each experiment, PA signal fluctuations were produced by a controlled flow of blood at physiological concentration.

A.1. Experimental setup
The three types of experimental configurations are shown in Fig. 4. The first experiment consisted in imaging a C-shaped structure formed by a polycarbonate (PC) capillary (inner diameter D = 100 µm, wall thickness w = 20 µm, Paradigm Optics Inc., Vancouver, USA ) while the second and the third experiments consisted in imaging a cross-section of a glass tube (inner diameter D = 1 mm, wall thickness w = 100 µm, Capillary Tube Supplies Ltd, UK) aligned within the imaging plane (2 nd experiment) and perpendicularly to the imaging plane (3 rd experiment). Both samples and the probe were held in a water tank. For each experiment, blood flow was induced through the tubes and capillary via a syringe pump (KDS Legato 100, KD For each PA acquisition, the raw data was available as a radio-frequency (RF) frame containing the signals received by each transducer element of the CMUT array measured with a sampling frequency f s = 62.5 MHz. The total number of acquired RF frames was M 1 = 1,000 (C-shaped capillary), M2 = 1,000 (glass tube, parallel orientation) and M3 = 10,000 (glass tube, perpendicular orientation). Given the 100 Hz PRF, the total acquisition times for each experiment was T 1 = T 2 = 10 s and T 3 = 100 s. A longer acquisition time was chosen for the 3rd experiment to ensure that the number of independent red blood cells (RBC) configurations induced by the blood flow was similar in all the three experiments. In contrast to the 1st and the 2nd experiments, the blood flow in the 3rd experiment was oriented perpendicular to the imaging plane. However, the characteristic PSF size in the direction perpendicular to the imaging (xz) plane was about ten times larger (of the order of 1 mm) than the dimension of the PSF within the imaging plane (of the order of 100 µm). With the chosen blood velocity of 1 cm/s, the RBC distribution over 100 µm (respectively 1 mm) distance is typically renewed every 10 ms (respectively 100 ms). Consequently, each of T 1 , T 2 and T 3 corresponds to 1000 independent configurations during the measurement time.

A.2. Data processing
For each of the three configurations, N = 1000 independent frames over the total acquisition time were used for data processing. For each experiment, N = 1000 complex-valued RF frames were computed from the N = 1000 real-valued frames via a Hilbert transform along the time axis. Complexvalued PA reconstructions A k (x, z) (with k=1...N) were computed by applying standard delay-and-sum beamforming reconstruction. The mean reconstruction was computed as E[A](x, z) = 1 N ∑ M k=1 A k (x, z), and the mean image was defined as |E[A](x, z)|. While the fluctuations of interest were those induced by the blood flow, other sources of fluctuations were present in the experiments, including laser pulse energy fluctuations and detection noise. Because laser pulse energy fluctua-tions affect identically all PA sources, they may be compensated by pulse-to-pulse energy monitoring. However, this requires additional hardware and signal processing.
Here, we used spatio-temporal singular value decomposition (SVD) as a key processing step to filter out laser pulse energy fluctuations (PEF). The SVD approach was first introduced in ultrasound imaging to discriminate tissue and blood motion [30] and was used by our group [20] later on to extract relevant fluctuations in fluctuation-based super-resolution PA imaging. In brief, SVD decomposes the initial data into a basis of spatiotemporal singular vectors. By choosing the singular vectors corresponding to relevant fluctuations, one can suppress signals with different spatiotemporal behaviour such as tissue motion, laser and electronics noise, etc. In our experiments, the first singular vectors (with the highest singular values) corresponded to laser PEF. In the Supplementary Information, we illustrate on both simulation and experimental results how filtering out the first singular values indeed provides a very efficient way to remove the effect of parasitic laser PEF (see Supplementary information, section 4.A, Fig. S2). We also illustrate the effect of SVD filtering as a function of the range of selected singular values (see Supplementary information, section 4C, Fig. S5). For all experimental images, the captions indicate the range of singular values that were selected (noted SVD cutoffs from now on).
For each set of N = 1000 SVD-filtered PA images, A SVD k (x, z), the standard deviation image was computed through its variance as : We note that the expression above is similar to the one used in ultrasound power Doppler (in which particular case the mean value of the SVD filtered images is zero) [30].

A.3. Results
The experimental results are shown in Fig. 5. Row 1 shows the conventional photoacoustic image. Row 2 shows the fluctuation images that are obtained without the SVD step. Row 3 shows the fluctuation images obtained after SVD filtering. These results demonstrate that fluctuations induced by a blood flow at physiological concentration can be exploited to suppress visibility artefacts of conventional PA imaging. Both limited-view and resonant-bandwidth artefacts are indeed well suppressed, in agreement with the simulation results shown in Fig.1. The residual signal (clutter) observed at the bottom of the images b 2 and c 2 is likely due to acoustic reverberation induced by the tube made of glass. This artefact is thus related to the sample rather than the proposed method. The slight inhomogeneity of the fluctuation amplitude within the C-shaped capillary (image a 2 ) may be caused by aberration induced by the capillary and/or an imperfect alignment of the structure in the imaging plane. As confirmed by simulations (see Fig. S2 of the supplementary material), the laser pulse energy fluctuations (PEF, measured with a powermeter around 3% relative rms value) completely overwhelms the fluctuation image if not compensated for. Our results illustrate that removing the first singular values efficiently filters out the effect of PEF, and avoid demanding hardware-based monitoring and compensation of those fluctuations. Fig. S5 of the supplementary material further illustrates that once the low singular values are removed to filter out the  (Fig. S3, bottom row), which show that depending on the desired image quality, N could be reduced from a a few hundred to a few tens of images, increasing the temporal resolution by the same factor. In addition, the 1D profiles through the fluctuation image presented in Fig. S3 also clearly illustrate that the fluctuation of interest appears on top of a noise fluctuation background caused by detection noise. On the principle, provided that N is sufficiently large, any fluctuation of interest can be detected on top of the noise background, thanks to the additivity of the variance of independent noise sources. Importantly, we note that quantitative measurements (not considered here) requires that the noise background can be estimated and substracted from the variance image. In all the results presented in this paper, the colormap was set such that black corresponded to the background noise on the fluctuation image.
The dependence on N illustrate the trade-off between acquisition time/temporal resolution and quality of the fluctuation image, which in practice will depend on the level of fluctuation of interest as compared to noise. For the phantom experiments shown here, the peak to noise ratio in the signal domain was around 25-30, but the fluctuation of interest (from blood motion) were about half the noise rms value, which corresponded to a quite realistic situation. Thanks to both the beamforming step (which averages out RF noise by summation over the 128 channels) and the statistical estimation of the fluctuation along N realisations, the fluctuation of interest can in the end be detected in the image domain on top of the noise background even if weaker.

B. In vivo experiments
The proposed approach is further illustrated with a preliminary demonstration of 3D imaging of a live chicken embryo. The chorioallantoic membrane (CAM) of the chicken embryo is composed of small blood vessels whose role is to capture oxygen through the shell and to supply it to the developing embryo. Vessels are spread on a large surface and can have 3D orientation in the vicinity of the embryo, and, therefore, provide a relevant model to assess the feasibility and performance of the proposed approach. Fig. 6. Schematic of the experimental setup for 3D in vivo imaging of chicken embryo.

B.1. Experimental setup and methods
The experimental configuration is illustrated on Fig.6. A custom 256-channel spherical transducer array (8 MHz, bandwidth 80%, f = 35 mm, f /D = 0.7) was designed by our group and fabricated by Imasonics (Voray-sur-L'Oignon, France). This transducer was connected to the same Verasonics acquisition system used in our phantom experiments. An in-house coupling cone, filled with water and closed with a latex membrane, was used to couple the transducer to the sample. Fertilized eggs were obtained from a local farm and placed in an incubator for 6 days.
After making a small hole through the shell using a scalpel, 1.5 mL of egg white was removed with a syringe and the top part of the shell was cut with scissors. Warm Phosphate-Buffered Saline (PBS) solution was poured into the shell to ensure acoustic contact with the latex membrane. The sample was held by a hammock made of food wrap. The hammock was attached to a bowl containing water maintained at 36 • C using a hot plate. The sample was illuminated with λ = 730nm light via a custom fiber bundle (Ceramoptec, Germany) going through the transducer, resulting in a fluence at the top surface of the sample of approximately 3 mJ/cm 2 . M = 768 frames were acquired at 100 Hz and saved to the computer for off-line data processing. The associated average power density, 300 mW/cm 2 , was of the order of the ANSI limit of 230 mW/cm 2 . Offline processing was performed exactly as in the 2D phantom experiments, with delay-and-sum beamforming, SVD filtering and statistical computations.

B.2. Results
The top row of Fig.7 shows the maximum intensity projections (MIP) of the mean image in the three principal directions. These images present strong deterministic background features usually called clutter. In our experiment, this clutter is likely due to the limited number of transducers (N = 256) sparsely distributed over the ultrasound probe. Despite the relatively large detection aperture ( f /D = 0.7), several structures oriented predominantly along the z direction (upward and downward vessels) are missing or very weak in this conventional reconstruction, because of the limited-view problem. In contrast, MIPs of the 3D fluctuation images (bottom row of Fig.7) reveal vessels with all possible orientations. This preliminary experiment also suggests that, in addition to removing the visibility artefacts, extracting the flow-induced fluctuations also allows removing the clutter that appears presumably due to the sparse nature of our probe (inducing a 3D PSF with grating lobes). For the results presented on Fig.7, all the N = 768 frames (measurement time 7.68 seconds) were used to compute the fluctuation images, with SVD cutoffs = . In the supplementary information, we present results obtained from different values of N (Fig. S3) and also illustrate the effect of SVD filtering as a function of singular values cutoffs (Fig. 4C). As observed for the phantoms experiments, filtering out the first singular values is essential to isolate efficiently the fluctuations from blood from those from laser PEF. Fig. S3 shows that N = 128 (measurement time 1.28 seconds) still provides a relatively good image quality, and illustrate the tradeoff between acquisition time/temporal resolution and quality of the fluctuation image.
The purpose of this preliminary experiment was to illustrate the performance of the method in a biologically-relevant situation, with realistic blood flow. A detailed description and investigation of the performance of our 3D imaging setup are being carried out and was out of the scope of the work presented here.

DISCUSSION AND CONCLUSION
In conclusion, we demonstrated both theoretically and experimentally that fluctuations induced by blood flow can be exploited to palliate visibility artefacts of conventional limitedview and resonant-bandwidth PA imaging. We first propose a theoretical framework that provides a quantitative expression of the fluctuation image as a function of the system PSF and statistical properties of the random medium at stake. Our theoretical framework applies to fluctuations induced either by fluctuating illumination patterns or by random distributions of absorbing particles (including flowing red blood cells). Our theoretical predictions suggested in particular that visibility problems could be overcome through exploiting the flow of blood at physiological concentration (overcoming visibility problems had been demonstrated previously only experimentally and those experiments relied on multiple speckle illuminations under unrealistic conditions).
We first provided two-dimensional images obtained with vessel-mimicking structures flown with blood at physiological concentration. Finally, we presented a preliminary experiment demonstrating 3D fluctuation imaging in vivo of the vasculature of a chicken embryo CAM model. We emphasize that the proposed approach does not require additional hardware and can be implemented on conventional PA equipment since it relies on simple statistical processing.
The proposed technique has two main potential limitations inherent to its very principle, which are fundamentally related as one may be overcome at the expense of the other : temporal resolution and sensitivity. As discussed in section 3.A.3, and illustrated on Fig. S3 of the Supplementary Information, any fluctuation of interest can in principle be detected on top of any noise background (sensitivity), provided that the number N of acquisitions is made large enough, at the cost of the temporal resolution.
What we illustrate in the supplementary material (see sections 3 and 4) is that the relevant number to quantify the level of fluctuation is the fluctuation to noise ratio (FNR), rather than the conventional signal to noise ratio on the PA signals. As a rule of thumb, as illustrated by our phantom experiments and additional simulation results provided in Section 4 of the Supplementary Information (Fig. S2), for 128 channels, a FNR on the order of 1 requires N ∼ a few hundreds, while a FNR ∼ 10 requires N ∼ a few tens. Importantly, this estimate is valid for N uncorrelated frames, which was the case for our simulations and phantoms experiments. For laser pulse repetition rate high enough, or equivalently for blood flow slow enough, there might be correlations between RF frames, in which case the robustness in the statistical estimation is affected. In particular, it should be kept in mind that samples with different dynamics (such as vascular networks with different vessel sizes and different blood velocities) may be reconstructed with a variable robustness. For instance, an robust reconstruction of blood microcapillaries may require a longer acquisition time that for larger vessels, if the laser repetition period is short enough to resolve the dynamics in the microcapillaries.
In the case of our proof-of-concept in vivo experiment, the chicken embryo was mostly transparent, which is a relatively favorable situation as compared to imaging deep into tissue. The SNR (SNR = peak to noise ratio in the signals domain) in this experiment was on the order of 50, about twice higher than for the phantom experiments. However, the fluence value was approximately 3 mJ/cm 2 , whereas this could be increased by a factor ten for deep tissue experiments. In section 6 of the Supplementary material, we predict that the amplitude of the photoacoustic fluctuation image is typically 40 times (for our 15 MHz linear array) to 100 times (for our 8 MHz spherical 2D array) lower than the amplitude of the conventional photoacoustic image, when imaging blood vessels 25 µm (linear array) and 50 µm (spherical 2D array). Converted into the signal domain, this pre-dicts that a FNR about 1 requires a SNR of about 20 for our linear array and a SNR of about 50 for our 2D spherical array, which corresponds to the cases of both our phantom experiments and in vivo experiments).
Although drawing firm conclusions relative to imaging deep in tissue must involve further in vivo investigations in small animal or humans, our results on the chicken embryo suggest that the technique has the potential to be applied for deep tissue imaging, provided that a sufficient fluctuation to noise ratio (FNR) can be reached, i.e. that the number N of acquisitions is sufficiently large. Our estimates above indicate that a SNR around a few tens should provide blood flow a FNR about 1, which should provide very good images for N about a few hundreds.
When compared to deep-learning based approaches to solve the visibility problem, which can be implemented in a singleshot mode, our approach has a lower temporal resolution. It however has the merit of extreme simplicity, and can be implemented straightforwardly, as opposed to deep-learning approaches which requires the availability of a training set, and some quite specific know-how.
The proposed theoretical framework could also be potentially extended to investigate spatio-temporal cross-correlations of fluctuation images, which carry quantitative information on the flow of absorbing particles, and thus may contribute to the field of PA Doppler measurements. Finally, although we focused here on the application to the visibility problem, our theoretical prediction of the second-order fluctuation image will be also be useful to quantitatively exploit fluctuation-based super-resolution imaging.

FUNDING INFORMATION
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 681514-COHERENCE).

DISCLOSURES
The authors declare no conflicts of interest.
By integrating the equation above, and by use of A n * u p dr = A n (r)dr × u p (r)dr = A n (r)V 2 p dr, one gets the following expression for V g : The pair correlation function B(r) widely used in condensed matter physics [1], also introduced in the context of ultrasound physics [2][3][4][5][6][7], is related to the random microscopic density by the following expression [7]: where m = η V p is the average concentration of particles. By combining Eqs. S5 and S6, and noting that σ 2 g = η × (1 − η), one finally gets the following expression for V g : corresponds to the value for q = 0 of the so-called structure factor and is known in ultrasound backscattering theory as the packing factor [2][3][4][5][6][7]. Finally, one gets the following simple relationship between V g , V p , η and W(η):

COMPARISON TO ULTRASONIC BACKSCATTERING
From the above expression of V g , the amplitude of the photoacoustic variance image can be expressed as We note that the dependence of the photoacoustic variance as a function of η and W(η) is identical to that of the ultrasound backscattering coefficient BSC reported in the field of ultrasound imaging of blood as [6]: where µ s = σ s η V p is the ultrasound scattering coefficent (σ s being the scattering cross-section of a single particle).
It was out of the scope of our experimental work here to confront our theory with respect to its dependency on the volume fraction. However, it was recently reported from photoacoustic experiments on blood flow that the dependence of the amplitude of photoacoustic fluctuations (expressed as variance) as a function of blood concentration agreed with the dependence predicted from the backscattering theory as ηW(η) [8], which thus support our theoretical results. In these experiments, the best agreement was obtained with the expression of the packing factor of hard spheres (m=3), given by [6]: With η ∼ 50 % (typical volume fraction for blood), the formula above gives W ∼ 16 %, which is the value that was used in the main text to compare the fluctuations from blood to that from multiple speckle experiments. A detailed comparison of Eqs. S9 and S10 is out of purpose at this stage, given the quite different conditions for which these expressions have been obtained: Eq. S10 for the backscattering coefficient is obtained under the assumption of harmonic ultrasound, and do not take into account any specific detection geometry or resolution, while Eq. S9 is a theoretical prediction of a variance image, including both the properties of the pulsed imaging system (through its PSF) and the statistical properties of the medium.

FLUCTUATIONS OF INTEREST vs PARASITIC FLUCTUATIONS
In this second part of the Supplementary information, we investigate and discuss the influence of parasitic sources of fluctuation on the estimation of the fluctuation of interest. In photoacoustic fluctuation imaging, the two main sources of parasitic fluctuations are the electronic noise present on the measured RF signals, and the shot-to-shot energy fluctuation of the laser source. Temporal jitter between the ultrasound acquisition and laser pulses may also be an additional source of fluctuation, depending on the level of jitter as compared to the ultrasound time scale. In our experiments, the temporal jitter was on the order of 1 ns, and its effect was negligible as compared to that of the shot-to-shot energy fluctuations given our ultrasound frequency range, as verified through simulations. We thus limit our discussion here to the effects of the electronic noise and shot-to-shot laser energy fluctuations.
We first introduce the relevant quantities used to measure the various levels of fluctuations, in both the signal and image domains. We then investigate the effects of the electronic noise and shot-to-shot laser energy fluctuations on the performances of photoacoustic fluctuation imaging. In particular, we investigate the influence of the number N of images to acquire on the quality of the fluctuation image, on both simulations and experimental results. We also explain and illustrate in some details the effect of the SVD processing step, including a sensitivity analysis to the choice of SVD filtering cutoffs. Finally, based on the theoretical expressions provided in the main manuscript, we predict and compare quantitatively the amplitude of the conventional photoacoustic signals/images to that of the fluctuation signals/images, in the context of imaging small blood vessels (with both our linear and 2D arrays). In Fig. S1, we define the relevant noise, fluctuation and peak amplitude, in both the signals domain and the image domain. We used the fact that the sum of independent random variables has a total variance given by the sum of each variance. In conventional ultrasound signal processing, the SNR would be defined as P/n. In the context of fluctuation imaging, the useful "signal" is rather a fluctuation signal, noted F. To avoid confusion, we defined the peak-to-noise ration PNR = P/n, the fluctuation to noise ratio FNR = F/n. As illustrated in the image domain, the relevant "SNR" quantity for fluctuation imaging is given by

DEFINITIONS OF RELEVANT NOISE-RELATED QUANTITIES
While the value of FNR BF in the image domain is similar to the value of FNR in the signals domain, they are strictly speaking not identical as the noise is usually uncorrelated between channels, whereas there may be some spatial coherence in the fluctuation from real sources. As independent fluctuations are summed incoherently, whereas PA signals are summed coherently during the beamforming step, the peak-to-fluctuation ratio PFR in the signals domain is usually extremely different from PFR BF = P BF F BF defined from the comparison between the mean image and the fluctuation image. The possibility to estimate robustly a fluctuation image depends fundamentally on the FNR (i.e. on the detection noise, laser pulse energy, transducer sensitivity, etc.) while the PFR only depends on the absorbers distribution and the system PSF.

FLUCTUATION OF INTEREST VS PARASITIC FLUCTUATIONS
A. removing pulse energy fluctuations with SVD filtering. In Fig. S2, left panel, we provide simulations images that were obtained in extremely noisy conditions: SNR = 5, FNR = 0.5 (RF noise twice larger than the fluctuations of interest), and relative pulse energy fluctuations of 50%. The second raw illustrate the results of computing the fluctuation image without prior SVD filtering. In this condition, the fluctuation is completely dictated by the PEF, and what is seen is the average photoacoustic source distribution that fluctuates with the PEF. When just the first singular value is removed, the effect of the PEF is completely filtered out, and is image is visually identical to that observed without PEF in the simulation. The same effect is observed for the phantom experiment on all three types of structures (as already shown in Fig.5 of the main manuscript). The only difference between simulations and experiments is that it is necessary to remove more than one singular value to filter out the PEF in the experimental situations. We have no clear explanation for this at that stage. The sensitivity to the choice of SVD filtering cutoffs is illustrated further below (Fig.S5). Fig.S3 illustrates both the effects of the electronic noise and of the number N of acquisition. As a key conclusion, the 1D profiles through the fluctuation images shows that the object appears on top of the noise background. The level of noise is actually only dictated by the electronic noise (value 1 in the image space, as our choice of unit). The ability to detect the object on top of the noise is actually depending on the statistical power with which both the noise and the fluctuation of interest are estimated: for large FNR values, the noise is negligible, and the only fluctuation comes from the fact that N is finite. When noise is significant (FNR on the order of 1 or lower), it is always possible to detect the additional fluctuation from the object, but at the cost of a large N. For the phantom results shown at the bottom row (corresponding to the experiment used for Fig. 5), the SNR was about 25, and the FNR about 0.8. Fig.S3 clearly demonstrates that the choice of N is dictated by the FNR value, and a trade-off between temporal resolution and robust estimation of the object. Fig.S4 illustrates this trade-off for the chicken embryo experiments.   The main message from this sensitivity plot is that removing the first singular values is essential to filter out the effect of laser PEF. While removing just one the first singular value allowed to perfectly removed the effect of the PEF in simulations, it is necessary to remove more for the experiments. We have no clear explanation for this at that stage. One other important conclusion is that filtering out high singular values has little effect on the final image. We observe on quantitative analysis of the fluctuation images that both the noise and the object are decreased when high singular values are filtered out. We observed that depending on the image quality criterion, filtering some high singular values can slightly improve the image quality. This explain why we present images with both low and high singular values cutoff. Fig.S5D, obtained for the phantom experiment (B), illustrates how a CNR metrics (defined as CNR = Fluctuation in object -background fluctuation level spatialvarianceo f ensemblevarianceinbackground+spatialvarianceo f theensemblevarianceinobject ) is affected by the choice of SVD cutoffs. This figure illustrate that from a CNR point of view, the key choice is that of the low SVD cutoffs: too low would leave some PEF, while too high degrades the object of interest. This metrics is very marginally affected by the choice of the high SVD cutoff.

PFR BF PREDICTION FOR VESSEL-LIKE STRUCTURE
In this section, we use equations (4) and (5) of the main manuscript to numerically compute the ration R(r) between the mean photoacoustic image and the fluctuation image, for vessel like structures with various orientation, in 3D.
This ratio allows a quantitative comparison between the two modalities in given situations defined by a 3D object f (r) and the 3D PSF, and corresponds to a local value (computed pixel wise) of the PFR BF as defined in section 3. A blood flow is considered with η = 50%, W(50%) = 16% and V p = 100µm 3 . The signals to reconstruct the 3D PSF were computed for our linear array (L22-8v) using Field-2 [9][10] and measured for our custom 2D array. We computed (S12) for the case of a horizontal and a vertical thin vessel (λ/4 diameter, λ/4 = 24µm for the linear array @15MHz and λ/4 = 47µm for the 2D array @8 MHz) (see Fig.S6a-d). With this object geometry, conventional mean PA image suffers from visibility artefacts only for the vertical vessel as shown in Fig.S6b-e. The fluctuation images show no visibility artefact (Fig.S6c-f). Images are normalized by the maximum of the conventional mean images. R values are evaluated by taking the ratio of the images at the center of the object r obj for the different types of vessel, horizonal and vertical. Table S1 indicates the computed values of R. As an indication, we also computed R for the case of a transverse (Transv.) vessel for the 2D imaging case, ie. going across the imaging plane. The ratio between the fluctuation at any configuation r obj and the value of E[A](r hori ) in the absence of visibility artefact is also computed at the bottom line. The fluctuation is about 40x smaller than the mean PA image in the in-plane 2D imaging case (vessel diameter λ/4, 15 MHz). In the transverse flow case, the anisotropy of the PSF induces a different value around 100. In 3D imaging, the PSF is isotropic and the values for both the horizontal and vertical vessels are close to 100. It should be kept in mind that these values vary with the frequency/geometry of the probe and the considered vessel size. However, the quantitative evaluation of R(r obj ) offers a way to predict the fundamental limits in the detection of PA fluctuations at depths in the image space, as this ratio is independent of the fluence.