**2. Methods**

As shown in Figure 1, an imaging target is placed in the region of interest (ROI), which is surrounded by a skull-like scattering layer. Under the illumination of the pulse laser, the optical absorber will radiate ultrasound waves to the surrounding medium as a result of the PA e ffect. After the acoustic waves propagate through the random scattering layer, a passive ultrasound annular array with 576 elements detects the ultrasound signal. The annular array is then divided into *N*S sectors for data processing and imaging. Each sector has *N* = 45 elements and can be approximated as a linear array so that the filtered imaging method proposed in our previous studies can be applied for the image reconstruction [26,27]. For each sector array, the recorded *N*-channel signals are arranged in a vector **H**(*t*) = [*H*1(*t*),*H*2(*t*), ... ,*Hn*(*t*), ... ,*HN*(*t*)] where *n* = 1,2, ... ,*N*.

**Figure 1.** The schematic of the scenario considered in this study. For any target placed in the region of interest (ROI) within a scattering layer, its propagation distances to the *n*th transducers of a sector array, and the equivalent linear array are *l*1 and *l*2, respectively. The corresponding delay is then given by *dn* = (*l*1 − *l*2)*F*s/*<sup>c</sup>*, where *c* is the sound speed and *F*s is the sampling frequency.

The sector array can be approximated to a linear array by calculating the corresponding delay. For any imaging target, its propagation distances to the *n*th transducers of the two arrays are *l*1 and *l*2, respectively. The corresponding delay is *dn* = (*l*1 − *l*2)*F*s/*<sup>c</sup>*, where *c* is the sound speed and *F*s is the sampling frequency. The time domain signals for the linear array are then **J**(*t*) = [*H*1(*t* − *d*1),*H*2(*t* − *d*2), ... ,*Hn*(*<sup>t</sup>* − *dn*), ... ,*HN*(*<sup>t</sup>* − *dN*)]. The signals can be written in their frequency domain by applying a Fourier transform to **<sup>P</sup>**(*T*,*f*) = [*P*1(*T*,*f*),*P*2(*T*,*f*), ... ,*Pn*(*T*,*f*), ... ,*P <sup>N</sup>*(*T*,*f*)], where *Pn*(*T*,*f*) is the short-time Fourier transform of the *n*th channel signal in the time window [ *T* − Δ*t*/2, *T* + Δ*t*/2], *T* is the time of flight, and Δ*t* is the width of the window. **<sup>P</sup>**(*T*,*f*) contains both the direct wave **<sup>P</sup>**D(*<sup>T</sup>*,*f*) and the scattering wave **<sup>P</sup>**S(*<sup>T</sup>*,*f*). When the concentration of the scatterers is su fficiently high, the scattering wave occupies a dominant position and the imaging quality is greatly degraded due to the randomness of the scattering process.

To reduce the influence of the randomly scattering, a matrix filtering method can be applied to extract the direct wave from the received signals of a linear array [26]. Firstly, a response matrix can be constructed as:

$$\mathbf{K} = \mathbf{P}\mathbf{P}^{\mathrm{T}} = \underbrace{\mathbf{P}^{\mathrm{D}}(\mathbf{P}^{\mathrm{D}})^{\mathrm{T}}}\_{\mathrm{Coherence,K}^{\mathrm{C}}} + \underbrace{\left[\mathbf{P}^{\mathrm{D}}(\mathbf{P}^{\mathrm{S}})^{\mathrm{T}} + \mathbf{P}^{\mathrm{S}}(\mathbf{P}^{\mathrm{D}})^{\mathrm{T}} + \mathbf{P}^{\mathrm{S}}(\mathbf{P}^{\mathrm{S}})^{\mathrm{T}}\right]}\_{\mathrm{Random,K}^{\mathrm{R}}} \tag{1}$$

Depending on whether or not contain the term **P**S, the matrix **K** can be divided into two parts: **K**<sup>R</sup> and **K**C, where **K**<sup>R</sup> is a random term due to the random nature of the scatterer distribution and the scattering paths and **K**<sup>C</sup> is independent of the scattering paths and shows a deterministic relation of the phase between its elements along the antidiagonal:

$$\beta\_{\emptyset} = \frac{k\_{m-q,m+q}^{\mathbb{C}}}{k\_{m,m}^{\mathbb{C}}} = \exp\left[jk(qw)^2/Z\right] \tag{2}$$

where *k* is the wave number, *w* is the pitch size, and *Z* is the depth of the time window. The phase relation along the antidiagonal only depends on the channel position *qw*, indicating that the direct waves from the optical absorbers manifest themselves as a particular coherence on the antidiagonals of the matrix.

Using the coherent property, we can extract the direct wave by the following steps. Firstly, the original matrix **K** is rotated 45◦ counterclockwise so that the coherent property occurs on the columns. The elements are divided into two parts based on their symmetry characteristic and are expanded to two new matrices: **A**1 and **A**2 [26]. The matrices **A**1 and **A**2 are called indi fferently **A** as they are filtered in the same way. A filtering matrix **F** = **CC**† is then applied to extract the coherent part from matrix **A** as: **A**<sup>F</sup> = **FA** = **FA**<sup>C</sup> + **FA**R, where **C** = [exp(*jky*<sup>2</sup> 1/2*Z*), ... , exp(*jky*<sup>2</sup> *u*/2*Z*), ... , exp(*jky*<sup>2</sup> *N*A /2*Z*)] with *u* = 1, 2, ... , *N* A*,* and *yu* = [*xm* + *xn*+(*N* A − 1) *w*] is the characteristic spacing of direct waves. The coherent part remains unchanged during the filtering process, as **FA**<sup>C</sup> = **A**<sup>C</sup> and the random part decreases by a factor (N + 1)/2. Finally, a filtered matrix **K**<sup>F</sup> with the original coordinate can be obtained by applying the inverse process of the first step to matrices **A**1 F and **A**2 F.

Using a time reversal operator, we can recover the image of the optical absorbers from the matrix **K**<sup>F</sup> corresponding to the single sector area [28,29]. The singular value decomposition is applied to the matrix **K**<sup>F</sup> = **<sup>U</sup>**F**Λ**F**V**F†, and each optical absorber is mainly associated with one non-zero singular λi F contained in diagonal matrix **Λ [30]**. Therefore, by numerically back-propagating the singular vector, the image of the *i*th singular valued at a depth of *Z* = *cT* can be achieved by the time reversal operator **I***i* = λi <sup>F</sup>|**G** × **V**|, where **G** with the component *gm*= exp(*jkrm*)/ √*rm* is the propagating operator between the array and the focal plane in a homogeneous media, and *rm* is the distance between the *m*th element in the array and the point in the focal plane. By superimposing the imaging results of multiple sectors, we can further improve the imaging quality and reduce the influence of the inhomogeneous aberrating effect of the scattering layer. For the wanted imaging point A, we calculated the image value *Im*<sup>A</sup> for all the linear arrays with *m* = 1, 2, ... , *N*s. Taking the first sector area as an example, the image value

for points corresponding to the array can be obtained by applying the time reversal operator. If point A is within the imaging area, we can find two imaging points, B and C, closest to point A in the focal plane and then perform a weighted evaluation based on the distance from point A to points B and C:

$$I\_1^{\Lambda} = \frac{I\_1^{\text{B}} d\_{\text{AC}} + I\_1^{\text{C}} d\_{\text{AB}}}{d\_{\text{AC}} + d\_{\text{AB}}} \tag{3}$$

where *d*AC and *d*AB represent the distances between A, B and A, C, respectively. If point A is outside the imaging area, the imaging value is then recorded as 0. When we calculate the imaging value of point A for all sector areas, its final value is the result of taking the summation and average to all the non-zero values: *I*<sup>A</sup> = *Ns <sup>m</sup>*=1 *<sup>I</sup>*A*m*/(*Ns* − *Nz*), where *N*Z is the number of zero values. As such, we can fully use the valid information received by all the transducers of the annular array for image reconstruction of each point in the ROI.
