1. Introduction
Cavitation bubble clusters nucleate when the liquid pressure rapidly falls below a certain threshold. These clusters coherently oscillate and violently collapse to cause extreme energy concentration that leads to various critical consequences, and use in applications as diverse as injectors and pumps (e.g. Plesset & Ellis Reference Plesset and Ellis1955; Blander & Katz Reference Blander and Katz1975; Mørch Reference Mørch1980; Chen & Heister Reference Chen and Heister1994; Hashimoto et al. Reference Hashimoto, Yoshida, Watanabe, Kamijo and Tsujimoto1997; Prosperetti Reference Prosperetti2017), hydraulic machines (Arndt Reference Arndt1981), underwater propulsion and hydrofoils (e.g. Kubota, Kato & Yamaguchi Reference Kubota, Kato and Yamaguchi1992; Merkle, Feng & Buelow Reference Merkle, Feng and Buelow1998; Kunz et al. Reference Kunz, Boger, Stinebring, Chyczewski, Lindau, Gibeling, Venkateswaran and Govindan2000; Schnerr & Sauer Reference Schnerr and Sauer2001; Ganesh, Mäkiharju & Ceccio Reference Ganesh, Mäkiharju and Ceccio2016; Gnanaskandan & Mahesh Reference Gnanaskandan and Mahesh2016; Venning, Pearce & Brandner Reference Venning, Pearce and Brandner2022), medical ultrasound (e.g. Ikeda et al. Reference Ikeda, Yoshizawa, Tosaki, Allen, Takagi, Ohta, Kitamura and Matsumoto2006; Krimmel et al. Reference Krimmel, Colonius and Tanguay2010; Maxwell et al. Reference Maxwell, Wang, Cain, Fowlkes, Sapozhnikov, Bailey and Xu2011; Pishchalnikov, Williams & McAteer Reference Pishchalnikov, Williams and McAteer2011; Maeda et al. Reference Maeda, Kreider, Maxwell, Cunitz, Colonius and Bailey2015; Movahed et al. Reference Movahed, Kreider, Maxwell, Hutchens and Freund2016; Maeda Reference Maeda2020), surface cleaning (Verhaagen & Rivas Reference Verhaagen and Rivas2016; Yamashita & Ando Reference Yamashita and Ando2019), chemical synthesis (Suslick et al. Reference Suslick, Didenko, Fang, Hyeon, Kolbeck, McNamara, Mdleleni and Wong1999; Cairós & Mettin Reference Cairós and Mettin2017), and bio-inspired devices (Tang & Staack Reference Tang and Staack2019). Characterizing the dynamics is a challenge due to the complex interactions of bubbles involving disorders and stochasticity. Nuclei are typically micro-sized and polydisperse and randomly distributed. Their rapid, nonlinear oscillations are, except in controlled experiments (Bremond et al. Reference Bremond, Arora, Ohl and Lohse2006), practically not measurable. Molecular and hydrodynamics simulations can provide detailed insights into nucleation (Angélil et al. Reference Angélil, Diemand, Tanaka and Tanaka2014; Gallo, Magaletti & Casciola Reference Gallo, Magaletti and Casciola2021), while their time- and spatial scales have not reached those of practical cluster oscillations. Analyses have been made on the interaction dynamics in various regimes (Brennen Reference Brennen2014), yet no common knowledge has been established if the many-body coherence globally exists and if so scaling is possible, beyond the consensus that polydispersity induces strong disorders.
For the past decades, the Rayleigh–Plesset (R–P) equation and its variations have been actively explored to investigate the dynamics of single bubbles (Plesset Reference Plesset1949; Plesset & Prosperetti Reference Plesset and Prosperetti1977). Relatively few studies addressed the theory of clusters. By using the mean field approach to interacting bubbles modelled by the R–P equation, d'Agostino & Brennen (Reference d'Agostino and Brennen1989) derived a non-dimensional parameter that dictates the linear coherent oscillations of monodisperse clusters, the so-called ‘cloud interaction parameter’. Zeravcic, Lohse & Van Saarloos (Reference Zeravcic, Lohse and Van Saarloos2011) used the coupled R–P equations and identified disorders represented by the Anderson localization of acoustic energy in polydisperse, lattice-like clouds under weak excitation. We have recently extended the interaction parameter to the non-equilibrium, cavitating clusters under strong excitation by considering the effective interaction at excited states (Maeda et al. Reference Maeda, Maxwell, Colonius, Kreider and Bailey2018; Maeda & Colonius Reference Maeda and Colonius2019; Maeda & Maxwell Reference Maeda and Maxwell2021). To recall, we scale the mean kinetic energy of liquid induced by 
$N ({\gg }1)$ bubbles as (Maeda & Colonius Reference Maeda and Colonius2019)
\begin{equation} \langle K\rangle \sim \langle K_s \rangle (1+B_d),\end{equation}
where 
$K_s$ is the energy of a single bubble: 
$K_s=2{\rm \pi} \rho R^3_{b,c}\dot {R}^2_{b,c}$; 
$B_d$ is the parameter controlling the effective contribution of hydrodynamic inter-bubble interaction: 
$B_d = N{\langle \overline {R(t)} \rangle }/R_C$, and 
$R_{b,c}$ and 
$R_C$ denote the characteristic (reference) bubble radius and the cluster radius; 
$\rho _l$ is the liquid density; 
$\langle {\cdot } \rangle$ and 
$\overline {({\cdot })}$ denote time average during a period in which bubble dynamics are statistically stationary and the mean value about the bubbles (
$i=1,2,\ldots,N$), respectively. The scaling can be simply derived from the coupled R–P equation for the correlated (synchronized) limit of monodisperse bubbles (
$R_1=R_2\cdots =R_{b,c}$). Although realistic correlations are imperfect due to polydispersity and nonlienarity, 
$B_d$ was found to control well both the coherent dynamics of polydisperse cavitating clusters, and their acoustic emission in numerical simulations and experiments (Maeda & Colonius Reference Maeda and Colonius2019; Maeda & Maxwell Reference Maeda and Maxwell2021). Overall, previous studies indicate that the coherence can depend on both polydispersity and nonliearity in a non-separable manner, posing perplexing questions about the universality of scaling. The theoretical characterization of the nonlinear dynamics of disordered many-body systems is, in general, not a simple task. Meanwhile, greater computing power has enabled learning physics by analysing big data. Principal component analysis (PCA) is a powerful method for unsupervised learning which has seen recent success in characterizing the coherent physics of many-body and high-dimensional systems in fields ranging from quantum information to fluid dynamics (Milano & Koumoutsakos Reference Milano and Koumoutsakos2002; Holmes et al. Reference Holmes, Lumeley, Berkooz and Rowley2012; Lloyd, Mohseni & Rebentrost Reference Lloyd, Mohseni and Rebentrost2014; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017).
In this study, we introduce and demonstrate a method for unsupervised data decomposition to study the coherent bubble cluster dynamics by combining theory and PCA. PCA extracts dominant states and dynamical features, such as coherent quantum states and turbulent structures, and their amplitudes as the eigenfunction (feature) and the eigenvalue (variance) of the co-variance matrix of physical data. When applied to spatio-temporal data of dynamical systems, PCA is often denoted as proper orthogonal decomposition (POD). Those data can be properly weighted prior to PCA such that the variance becomes consistent with the norm induced by an energetic inner product of state variables (e.g. kinetic energy) (Lall, Marsden & Glavaški Reference Lall, Marsden and Glavaški1999; Rowley Reference Rowley2005). This weighting allows a physical interpretation that resulting features associated with a large variance are energetically dominant coherent structures. Proper weighting of Lagrangain bubble dynamics data is non-trivial since the linear variance of extracted features needs to account for the nonlinear interaction energy. For meaningful analysis, we introduce strategic pre-processing of the data prior to PCA such that the PC-variance becomes theoretically consistent with the energy modelled by the coupled R–P equation. Analysing simulation datasets of clusters, we show that the PCA can systematically extract not only coherent but also incoherent features whose magnitudes are respectively measured by the PC-variance and the entropy. We discover that the coherence is lost by disorders induced by polydispersity and nonlinearity, while under strong excitation, the underlying correlations are globally isolated in a single coherent feature whose variance (energy) is scaled by 
$B_d$, regardless of the disorders.
The remainder of this paper is as follows. In § 2, we describe the method. The first PC variance, spectral entropy and the coherence measure are introduced as quantifiable measures to characterize the coherent bubble dynamics from extracted features. In § 3, we verify and demonstrate the method using a numerical dataset of bubble clouds with different density and polydispersity parameters under various amplitudes of harmonic excitation. The deviation of the PC-variance from the physical energy is quantified for two weighting methods. The amplitude dependencies of the measures are quantified. The PC-spectra and their correlations with the coherent energy are analysed. Moreover, the extracted coherent dynamics is related to 
$B_d$ and its universality is discussed for cavitating clouds. In § 4, the physical significance of the method is discussed. In § 5, we state conclusions.
2. Methods
2.1. Principal component analysis of bubble dynamics data
For clusters modelled by the coupled R–P equations, observable dynamical variables are the bubbles’ radial velocities and radii. Consider a data matrix 
$\boldsymbol {Q}$ containing the 
$N_t$ snapshots of the radial velocities with a constant temporal interval,
\begin{equation} \boldsymbol{Q}=[\boldsymbol{q}_1,\boldsymbol{q}_2,\dots,\boldsymbol{q}_{Nt}], \end{equation}
where 
$\boldsymbol {q}_k$ denotes the vector containing the radial velocities of the 
$N$ bubbles at time 
$t_k$ (in the 
$k$th snapshot),
\begin{equation} \boldsymbol{q}_k=[\dot{R}_1(t_k),\dot{R}_2(t_k),\dots,\dot{R}_{N}(t_k)]^{\rm T}. \end{equation}
For later convenience, we also define the vector 
$\boldsymbol {r}_k$, containing the radii of the bubbles at the same instances:
\begin{equation} \boldsymbol{r}_k=[R_1(t_k),R_2(t_k),\dots,R_N(t_k)]^{\rm T}. \end{equation} PCA can be performed on 
$\boldsymbol {Q}$ by using the singular value decomposition (SVD) (e.g. Abdi & Williams Reference Abdi and Williams2010; Jolliffe & Cadima Reference Jolliffe and Cadima2016):
\begin{equation} \boldsymbol{Q}=\boldsymbol{U}\boldsymbol{\varSigma}\boldsymbol{X}^*. \end{equation}
The 
$i$th principal component (feature) is stored in the 
$i$th score matrix,
\begin{equation} \boldsymbol{\varPi}_{i}=\boldsymbol{U}\boldsymbol{\varSigma}_{i}, \end{equation}
where 
$\boldsymbol {\varSigma }_{i}$ contains the 
$i$th largest singular value and zeros elsewhere.
2.2. Weighted principal component analysis
Although this procedure for PCA is simple and straightforward, the physical meaning of the extracted features are obscure since PCA itself is not informed on the underlying dynamics of the system. A related issue of PCA for fluid flow data has been addressed in the context of POD. In POD of the snapshots of Eulerian fluid flow data, the state vectors consisting of the data matrix are often weighted such that the corresponding weighted inner product of the state becomes consistent with the kinetic energy of the original system. The dominant features (POD modes) can then be interpreted as energetically dominant coherent flow structures. Moreover, in our previous studies (Maeda & Colonius Reference Maeda and Colonius2019; Maeda & Maxwell Reference Maeda and Maxwell2021), the scaling of bubble cloud dynamics by 
$B$ was successfully demonstrated based on the total kinetic energy of liquid induced by interacting bubbles. These considerations motivate us to relate the energy to PCA for addressing the physics of bubble cloud.
Inspired by the POD, we consider weighting the present bubble dynamics data prior to PCA. To find an appropriate weight, we revisit the potential theory behind (1.1). The kinetic energy of the fluid induced by the spherical bubble oscillations is explicitly expressed as
\begin{equation} K =
2{\rm \pi}\rho_l\sum^N_{i=1} \Big[R_i^3\dot{R}_i^2+\sum^N_{j\ne
i}\frac{R_i^2R_j^2 \dot{R}_i
\dot{R}_j}{r_{ij}}\Big]+(H.O.T),\end{equation}
where 
$R_i$, 
$\dot {R}_i$ and 
$r_{ij}$ are the radius and the radial velocity of bubble 
$i$, and distance between the centres of bubble 
$i$ and 
$j$, respectively. The second term in the brackets represents the contribution of the long-range interactions. At each instant, 
$K$ can be expressed as a weighted inner product of 
$\boldsymbol {q}$:
\begin{equation} K=\boldsymbol{q}^T\boldsymbol{T}\boldsymbol{q}=(\boldsymbol{W}\boldsymbol{q})^{\rm T}(\boldsymbol{W}\boldsymbol{q}), \end{equation}where
\begin{equation} T_{ij}(\boldsymbol{r}_k) = \begin{cases} 2{\rm \pi}\rho R^3_i(t_k) & (i=j),\\[2pt] 2{\rm \pi}\rho \dfrac{R^2_i(t_k)R^2_j(t_k)}{r_{ij}} & (i\ne j). \end{cases} \end{equation}
The weight matrix, 
$\boldsymbol {W}$, can be obtained through the Cholesky factorization of 
$\boldsymbol {T}(\boldsymbol {r})$:
\begin{equation} \boldsymbol{T}(\boldsymbol{r})=\boldsymbol{W}\boldsymbol{W}^*. \end{equation}
This expression suggests the weighted data, 
$\boldsymbol {W}\boldsymbol {Q}$, is appropriate for use in PCA. However, 
$\boldsymbol {T}(\boldsymbol {r})$ and 
$\boldsymbol {W}$ are time dependent since 
${\boldsymbol {r}}$ can change in time during large-amplitude oscillations. The choice of 
$\boldsymbol {r}_k$ in the time series to define 
$\boldsymbol {W}$ is unclear. A straightforward choice is to use the temporal mean of the radius for each bubble, 
$\langle \boldsymbol {r}\rangle : \langle \boldsymbol {r}\rangle =\sum ^{Nt}_{k=1}{\boldsymbol {r}_k}/Nt$, but 
$T(\langle \boldsymbol {r}\rangle )$ is not a favourable approximation for 
$T(\boldsymbol {r})$ for bubbles with large-amplitude oscillations.
We address this obstacle by variable transformation as a means of pre-processing the data. The schematic is shown in figure 1. First, we transform the variables from 
$(\boldsymbol {q}, \boldsymbol {r})$ to 
$(\boldsymbol {\xi }, \boldsymbol {\eta })$, where 
$\boldsymbol {\xi }=[\xi _1, \xi _2,\ldots,\xi _N]^{\rm T}$ and 
$\boldsymbol {\eta }=[\eta _1, \eta _2,\ldots,\eta _N]^{\rm T}$, and 
$(\xi _i,\eta _i)=(q_ir_i, r_i)=(\dot {R}_iR_i,R_i)$ for 
$i\in [1,N]$. Here, 
$\xi _i$ is nothing but the velocity potential evaluated at the surface of bubble 
$i$. Supplemental discussions and justifications for this transformation are provided in Appendix A. The instantaneous energy of the system is expressed by a weighted inner product as
\begin{equation} K=\boldsymbol{\xi}^T\boldsymbol{P}(\boldsymbol{\eta})\boldsymbol{\xi}=(\boldsymbol{\varOmega}\boldsymbol{\eta})^{\rm T} (\boldsymbol{\varOmega}\boldsymbol{\eta}), \end{equation}
and 
$\boldsymbol {\varOmega }$ is the weight matrix satisfying
\begin{equation} \boldsymbol{P}(\boldsymbol{r})=\boldsymbol{\varOmega}\boldsymbol{\varOmega}^*, \end{equation}where
\begin{equation} P_{ij}(\eta) = \begin{cases} 2{\rm \pi}\rho \eta_i, & i=j,\\[2pt] 2{\rm \pi}\rho \dfrac{\eta_i \eta_j}{r_{ij}}, & i\ne j. \end{cases} \end{equation}
Second, we partially replace 
$\boldsymbol {\eta }$ with 
$\langle \boldsymbol {\eta }\rangle$ to approximate the system. Using 
$(\hat {\boldsymbol {\xi }},\hat {\boldsymbol {\eta }})$, the instantaneous energy of the approximate system is expressed as
\begin{equation} \hat{K}=\hat{\boldsymbol{\xi}}^T{\boldsymbol{P}}(\hat{\boldsymbol{\eta}})\hat{\boldsymbol{\xi}}.\end{equation}The temporal mean of the energy of the original system can then be approximated as
\begin{equation} \langle K\rangle \approx\langle\hat{K}\rangle =\langle\hat{\boldsymbol{\xi}}^T{\boldsymbol{P}(\hat{\boldsymbol{\eta}})}\hat{\boldsymbol{\xi}}\rangle \approx\langle\boldsymbol{\xi}^T\boldsymbol{P}(\langle\boldsymbol{\eta} \rangle)\boldsymbol{\xi}\rangle.\end{equation}
This is a critical result in the present context of PCA, since the energy is now related to the weighted inner product with the constant weight, 
$\boldsymbol {P}(\langle \boldsymbol {\eta } \rangle )$. Quantitative verification of this approximation is addressed through numerical experiments in the following section. Using the new set of variables, the PC decomposition is performed as
\begin{equation} \boldsymbol{\varOmega}\boldsymbol{\varXi}=\boldsymbol{U}_{\xi}\boldsymbol{\varSigma}_{\xi}\boldsymbol{X}_{\xi}^*.\end{equation}
The 
$i$th PC is stored in 
$\boldsymbol {\varPi }_{\xi,i}$: 
$\boldsymbol {\varPi }_{\xi,i}=\boldsymbol {U}_{\xi }\boldsymbol {\varSigma }_{\xi,i}$. We denote the 
$i$th largest singular value of 
$\boldsymbol {\varSigma }_{\xi }$ as 
$\sigma _i$. The degree of coherence can be measured by the normalized variance of the first PC:
\begin{equation} \hat{\sigma}_1^2=\sigma_1^2/\mathrm{tr}(\boldsymbol{\varSigma}^2), \end{equation}which represents the ratio of the energy occupied by the first PC to the total energy of the system. In the following sections, to distinguish the present method from the standard PCA, we denote (2.15) as the weighted principal component analysis of transformed data (WPCA-TD).
Figure 1. Schematic of the feature extraction by PCA from Lagrangian bubble dynamics data, 
$\boldsymbol {Q}$, after pre-processing. The variance of the resulting features is consistent with the interaction energy predicted by the coupled Rayleigh–Plesset equation.
2.3. Spectral entropy and coherence measure
We introduce two additional key measures to characterize the coherent physics of bubble clusters through WPCA-TD. First, to quantify the degree of incoherent bubble oscillations, we define the spectral entropy of the weighted co-variance matrix, 
$\boldsymbol {\varOmega }\boldsymbol {\varXi }$, as
\begin{equation} \hat{S}_{vN}=-\sum\frac{\hat{\sigma}^2_k(\ln\hat{\sigma}^2_k)}{\ln(N)}, \end{equation}
where 
$\ln _2(N)$ is the normalization factor. The spectral entropy of data is a discrete analogue of the Shannon entropy based on the spectrum of data and has been used to quantify the randomness of data in various applications (Kullback Reference Kullback1997; Alter, Brown & Botstein Reference Alter, Brown and Botstein2000; Hu et al. Reference Hu, Yin, Morris, Zhang and Wright2005; De Domenico & Biamonte Reference De Domenico and Biamonte2016). Aubry, Guyonnet & Lima (Reference Aubry, Guyonnet and Lima1991) introduced a similar definition of entropy for spatio-temporal signals of canonical fluid flows from their POD eigenvalues and used the entropy to characterize flow instabilities. In the present context, when bubbles are in perfect correlation, we expect to excite only the first PC capturing the entire energy of the system: 
$\hat {\sigma }^2_1=1$ and 
$\hat {\sigma }^2_k=0$ for 
$k: k\in [2,N]$, yielding 
$\hat {S}_{vN}=0$. In contrast, if the energy is equi-partitioned into all PCs, 
$\hat {\sigma }^2_k=1/N$ for all 
$k: k\in [1,N]$ and 
$\hat {S}_{vN}=1$. Second, we define the coherence measure 
$C$:
\begin{equation} C=\frac{B_1}{B_d}, \end{equation}where
\begin{equation} B_1=\frac{\sigma^2_1}{\sigma'^2_1}-1. \end{equation}
Here, 
$\sigma '^2_1$ is the first PC-variance excluding the contribution of interactions, obtained from the WPCA-TD of the same data using a diagonal weight matrix 
$\boldsymbol {\varOmega }'$ whose diagonal entries are those of 
$\boldsymbol {\varOmega }$. Further details of 
$\boldsymbol {\varOmega }'$ are provided in Appendix B. In the limit of perfect correlation (also see (1.1)),
\begin{equation} \langle K\rangle \approx \sigma^2_1 \approx (1+B_d)\sigma'^2_1\end{equation}and
\begin{equation} B_d\approx B_1. \end{equation}
Therefore, 
$C\approx 1$. This condition is typically realized in monodisperse clusters under weak (linear) oscillations (d'Agostino & Brennen Reference d'Agostino and Brennen1989). For real clusters under strong excitation, bubbles are not perfectly correlated and the energy can be distributed in broad features. In this regime, approximation (2.20) is not necessarily expected to hold and 
$C$ can take values far from unity. Phenomenologicaly speaking, 
$C$ quantifies from data the degree to which the most coherent mode of oscillation is represented by the cluster's mean-field interaction. A schematic of the process to obtain 
$C$ from data is provided in Appendix B.
3. Numerical experiments
3.1. Datasets
To verify and demonstrate the WPCA-TD, we use datasets of spherical clusters excited by 40 cycles of harmonic pressure excitation. Each cluster contains 
$O(10-10^3)$ bubbles with their initial radii following log-normal distributions with a reference radius of 
$R_{ref}=O(10) \mathrm {\mu }{\rm m}$; 
$\ln (R_0/R_{ref})=N(0,s_d)$, where 
$s_d$ is the lognormal standard deviation as the measure of polydispersity (Maeda & Colonius Reference Maeda and Colonius2019). We address 
$s_d=[0.1,0.3,0.5,0.7]$. In real bubble clouds, nuclei are expected to be polydisperse. Here, 
$s_d= 0.7$ may be a representative estimation based on previous studies (Katz Reference Katz1978; Ando, Liu & Ohl Reference Ando, Liu and Ohl2012; Maeda & Colonius Reference Maeda and Colonius2019) and is used unless noted. The smaller values of 
$s_d$ are used in some cases to quantify the effect of polydispersity. The bubbles are randomly distributed in the spherical region with a specified cluster radius. Similar parameters of clusters were previously simulated to compare with experiments (Maeda & Colonius Reference Maeda and Colonius2019; Maeda & Maxwell Reference Maeda and Maxwell2021). The density of bubbles is characterized by 
$B_0$, the value of 
$B_d$ at rest. The far-field pressure is given as 
$p_{\infty }(t)=p_0[1+A\sin (2{\rm \pi} ft)]$. The frequency of excitation is 
$f=\ $500 kHz unless noted, near the adiabatic resonant frequency of the reference bubble. The amplitude of excitation is defined by 
$A$, relative to the ambient pressure at 
$p_0=1.0$ atm. For each set of parameters, we compute an ensemble average by taking a mean of the results from 20 bubble clouds with distinct spatial placements of bubbles in the clouds. For data generation, we use mesh-free, coupled Keller–Miksis equations by modifying previous methods. Details of this method are provided in Appendix D.
3.2. Visualization of representative data
In figure 2, we show the evolution of representative quantities of a bubble cloud in the dataset with 
$(B_0,s_d,N)=(0.5,0.1,100)$ excited at 
$A=20$, during a stationary state. Figure 2(a) shows the radius of a representative bubble. The plot presents familiar features of cavitation bubbles forced by continuous, strong excitation, including fast events of collapse/rebound and slow growth/decay between them. During the slow phase, small-amplitude oscillations at excitation amplitude are evident. For the same bubble, figure 2(b) shows the evolution of the velocity potential obtained from the raw data, and those from the first and second dominant features extracted through the WPCA-TD. Compared with the radius, the raw data of the velocity potential looks much more symmetric about zero, even around the collapse events. The potential of the first feature presents a sinusoidal profile at the excitation frequency and with a constant amplitude which is close to the peak amplitude of the original data. The collapse events are not captured in this feature. The potential of the second feature is symmetric but has a much lower amplitude compared with the first feature. There is no clear similarity between the original data and the second feature, unlike that between the original data and the first-PC. Figure 2(c) shows the void fraction. The fraction oscillates with an amplitude of approximately 
$5.0\times 10^{-4}$ near 
$3.0\times 10^{-3}$ at the excitation frequency. Slow, small amplitude of variations are also observed. Figure 2(d) shows the kinetic energy of fluid induced by bubble oscillations. The energy oscillates approximately 
$0.75\ \mathrm {\mu }{\rm J}$ at doubled the excitation frequency. Although the frequency is as expected, the peak amplitude largely fluctuates in the window as well as the waveform is not symmetric, unlike the void fraction. The fluctuation and asymmetry indicated can be associated with the incoherent oscillations to the kinetic energy. For instance, if one bubble is expanding and another bubble is collapsing out of phase, their net contribution may not appear in the void fraction due to mutual cancellation but can appear in the energy. Overall, the averaged quantities, void fraction and the energy, are much smoother than the individual bubble dynamics. This can be trivially explained by the incoherence of violent collapse events among bubbles and coherence of the linear response against fundamental frequencies. Meanwhile, the quantitative nature of the coherent response is not predictable or easily analysable from these plots due to strong nonlinearity, especially under inter-bubble interactions with disordering factors including randomness of bubble position and polydispersity.
Figure 2. Evolution of representative quantities of a bubble cloud with 
$(B_0,s_d,N)=(0.5,0.1,100)$ excited at 
$A=20$, during a stationary state. (a) Radius of a representative bubble. (b) For the same bubble, 
$\xi$ obtained from the raw data, and those from the first and second dominant features extracted using WPCA-TD. (c) Void fraction. (d) Kinetic energy of the fluid induced by the bubble cloud.
Figure 3(a–c) shows the projected side-views of the three-dimensional bubble cloud of figure 2. The size of spheres can be interpreted as the mean energy of oscillations of bubbles at those locations in each feature. The overall distribution of the spheres of the first feature resembles that of the raw data, while the energy of the bubbles in the second feature is much smaller than that of the first feature. The plots therefore visually confirm that the first feature represents the most energetic mode of oscillations in the original data. This result also agrees with the observation of figure 3(b).
Figure 3. Projected side-views of the three-dimensional bubble cloud of figure 2. The size of the spheres denotes the root-mean-square amplitude of the velocity potential evaluated at the bubble surface at corresponding locations during the stationary state of oscillations, for (a) the raw data, (b) the first principal feature and (c) the second principal feature.
3.3. Error analysis
To show the effectiveness of the pre-processing, in figure 4, we plot the relative errors of the time-averaged kinetic energy of bubble clusters, one approximated using the original variable (
$\langle \boldsymbol {q}^T\boldsymbol {T}(\langle \boldsymbol {r}\rangle )\boldsymbol {q} \rangle$) and the other using the transformed variables (
$\langle \boldsymbol {\xi }^T\boldsymbol {P}(\langle \boldsymbol {\eta }\rangle )\boldsymbol {\xi } \rangle$), against the excitation amplitude. At 
$A<10^{-1}$, the error is nearly zero for both approximations. At 
$A>10^{-1}$, at which bubble dynamics become nonlinear, the error grows with 
$A$ for the former, while it remains small for the latter. This result confirms the improved approximation by the pre-processing.
Figure 4. Relative error of the mean kinetic energy of the fluid induced by clusters for the two approximations, 
$(\times ){:}\langle K\rangle \approx \langle \boldsymbol {q}^T\boldsymbol {T}(\langle \boldsymbol {q} \rangle )\boldsymbol {q}\rangle$ and 
$(\circ ){:}\langle K\rangle \approx \langle \boldsymbol {\xi }^T\boldsymbol {P}(\langle \boldsymbol {\eta }\rangle )\boldsymbol {\xi }\rangle$, against the excitation pressure amplitude with various polydispersities and values of 
$B_0$. Black, 
$(s_d,B_0)=(0.1,0.5)$; red, 
$(s_d,B_0)=(0.7,5.0)$.
3.4. Amplitude dependence of the key measures
The coherent dynamics critically depend on the excitation amplitude. Figure 5 shows the dependence of 
$\hat {\sigma }^2_1$, 
$\hat {S}_{vN}$ and 
$C$ against 
$A$ for various density and polidispersity parameters of clouds. This dependence is best highlighted in the result of a sparse, weakly polydisperse cluster (
$B_0=0.5, s_d=0.1$) in figure 5(a). The relative importance of the first PC, 
$\hat {\sigma }_1^2$, decays nearly monotonically from 0.9 to 0.2 through three distinct regimes. For 
$A\lessapprox 0.2$, 
$\hat {\sigma }^2_1\approx 1$, meaning the entire energy is captured by the first PC. Then 
$\hat {\sigma }^2_1$ rapidly decays to 0.4 and stays nearly constant up to 
$A\approx 3$. At 
$A>3$, 
$\hat {\sigma }^2_1$ decays again and stays nearly constant around 0.2. The decay indicates the decrease in the coherence with increasing 
$A$, and can be explained by the excitation of the nonlinear oscillations and cavitation triggered at 
$A\approx 1$ and above. Here, 
$\hat {S}_{vN}$ has a profile vertically mirrored to 
$\hat {\sigma }_1^2$; 
$\hat {S}_{vN}$ increases from approximately 0.1 to 0.5 through the three regimes, indicating more partitioning of the energy into multiple PCs and increase of incoherence, by increasing 
$A$. The mirrored profiles of 
$\hat {\sigma }_1^2$ and 
$\hat {S}_{vN}$ suggest that these parameters are complementary.
Figure 5. The first PC-variance (
$\hat {\sigma }_1^2$), the normalized von Neumann entropy (
$\hat {S}_{vN}$) and the coherence measure (
$C$), against the excitation amplitude with various initial density and polydispersity parameters, with 
$(B_0,s_d)= (a)$ 
$(0.5,0.1)$, (b) 
$(0.5,0.7)$, (c) 
$(5.0,0.1)$ and (d) 
$(5.0,0.7)$.
Remarkably, 
$C$ draws a square-well like profile, where 
$C\approx 1$ for both low (
$A<0.2$) and large (
$A>5$) amplitudes reaching a minimum value (
$C\approx 0.5$) at 
$A\approx 1$. This counter-intuitive result suggests that the energy of the first PC is scaled by the mean-field parameter 
$B_d$ regardless of the increase of incoherence for large 
$A$, and the scaling is lost only for the intermediate range at 
$1< A<5$. Comparisons with the result of dense, weakly polydisperse clusters 
$(B_0,s_d)=(5.0,0.1)$ shown in figure 5(b) highlights the effect of the density of bubbles. In figure 5(b), 
$\hat {\sigma }^2_1$ takes values near unity at small 
$A$, and decays at 
$A\approx 1.0$ to approximately 0.4 and stays nearly constant at 
$A>5.0$. At 
$A>10$, 
$\hat {\sigma }^2_1$ slightly grows against 
$A$. Here, 
$\hat {S}_{vN}$ has a profile vertically mirrored to 
$\hat {\sigma }_1^2$. The features of 
$\hat {\sigma }^2_1$ and 
$\hat {S}_{vN}$ are similar to those observed in figure 5(a), except that the transition occurs at a larger range of 
$A$. The sudden decay and growth of 
$\hat {\sigma }^2_1$ and 
$\hat {S}_{vN}$ can likewise be associated with the linear-to-nonlinear transition of bubble dynamics which results in incoherence. Meanwhile, the positive shift of the transition range of 
$A$ indicates that the dense bubble clouds tend to behave more coherently than sparse bubble clouds, agreeing with the previous theory (d'Agostino & Brennen Reference d'Agostino and Brennen1989; Maeda & Colonius Reference Maeda and Colonius2019). Here, 
$C$ is greater than 0.5 at almost all values of 
$A$, and is relatively more insensitive against 
$A$ compared with in figure 5(a), indicating that the cluster dynamics is moderately controlled by 
$B_d$. Figures 3(c) and 3(d) respectively correspond to sparse and dense clusters with strongly polydisperse nuclei; 
$(B_0,s_d)=(0.5,0.7)$ and 
$(5.0,0.7)$. These plots show clear differences from those of weekly polydisperse clusters in figures 5(a) and 5(b). For the sparse clusters (figure 5c), 
$\hat {\sigma }^2_1$ mildly decays from 0.8 to 0.7 at 
$0.1< A<0.2$, sharply decays to the minimum of 0.2 at 
$A\approx 0.4$ and then mildly grows to 0.3 at 
$A=1.2$. The sharp decay resembles that in figure 5(b), while the slope is milder. Additionally, 
$\hat {S}_{vN}$ has a profile mirrored to 
$\hat {\sigma }^2_1$, but it is less symmetric than in figures 5(a) and 5(b). Here, 
$\hat {S}_{vN}$ is nearly constant around 0.2 at 
$0.1< A<2.0$ and draws a concave curve with its maximum of approximately 0.5 at 
$A\approx 5.0$ followed by a smooth decay to 0.3 at 
$A=20$. The 
$C$ almost constantly increases from 0.2 to 
$O(1)$ throughout the plot, indicating that the coherent dynamics is controlled by 
$B_d$, only at large 
$A$, unlike the weakly polydisperse case of figure 5(a). The overall trend of the plots in figure 5(b) is similar to that in figure 5(a), although the changes of variables against 
$A$ are milder in figure 5(b).
Overall, the mutual trends of variables at 
$A>5$ are similar between figures 5(a) and 5(c), and between figures 5(b) and 5(d), indicating the decreasing influence of the initial polydispersity in the nonlinear regime. At 
$A< O(1)$, the polydisperse clouds tend to have smaller values of 
$\hat {\sigma }^2_1$, and larger values of 
$\hat {S}_{vN}$ and 
$C$. This can be explained by the enhancement of incoherent dynamics induced by polydispersity.
3.5. Principal component spectrum
To gain deeper insights to the meaning of 
$C$, in figures 6(a)–6(d), we show the PC-variances obtained at 
$A=2\times 10^{-2}$, 
$1.2$ and 
$20$, for the first 10 PCs, obtained from the sparse, weakly polydisperse clouds blueplotted in figure 5(a). The insets show the evolution of the square-root of the normalized total energy (
$\sqrt {K}$), and those of the first and the second PC-variances (
$\sqrt {K}_1$ and 
$\sqrt {K}_2$), during the four periods of excitation in statistically stationary states. As expected, in the linear regime (figure 6a), the first PC occupies nearly the entire energy and 
$\sqrt {K}_1$ evolves at the fundamental frequency. In the transition regime (figure 6b), the first PC occupies 40 % of the energy and the rest is partitioned into the sub-dominant PCs with a smooth decay. Here, 
$\sqrt {K}$ evolves more chaotic than the linear regime, as expected due to the nonlinear response of bubbles. Both 
$\sqrt {K}_1$ and 
$\sqrt {K}_2$ evolve with similar quasi-periodic profiles. We interpret that both PCs represent the coherent part of the energy.
Figure 6. (a–d) PC spectral profiles at 
$A=2\times 10^{-2}$, 1.2, 
$8.6$ and 20. Insets show the evolution of the square root of the normalized total energy (
$\sqrt {K}$, black) and those of the first (blue) and the second (red) PC-variances (
$\sqrt {K_1}$ and 
$\sqrt {K_2}$), during the four periods of harmonic excitation, with 
$\hat {t}$ being a non-dimensional time 
$\hat {t}=tf$. The 
$y$-axis of each inset is normalized by the maximum value of 
$\sqrt {K}$.
Interestingly, in the nonlinear regime (figure 6c), energy partitioning is non-smooth; the first PC occupies 40 % of the total energy similar to the transition regime, but the rest of the energy is broadly distributed into the other PCs with much smaller amplitudes. The evolution of 
$\sqrt {K}$ is non-periodic with noisy, fine structures of spikes. These spikes are expected due to the incoherent collapse events. The evolution of 
$\sqrt {K_1}$ is, in contrast, highly periodic and somewhat resembles that of figure 6(b). The evolution of 
$\sqrt {K_2}$ is more chaotic and less smooth than 
$\sqrt {K_1}$. The difference between the first-PC and the rest of the PCs suggest that, in this regime, only the first PC captures a major coherent feature and the rest of the PCs represent more the incoherent dynamics as broadband noise. Figure 6(d) shows the result of the nonlinear regime with a stronger excitation amplitude (
$A=20$). Overall, both the spectrum and the evolution of 
${K}$ look similar to that of figure 6(c), other than that the amplitude of the first PC-variance is increased to 60 %. The resemblance of figures 6(c) and 6(d) indicates that the dynamical features identified from these two plots are common in the nonlinear regimes.
The resemblance of the evolution of the first PCs in the linear and the nonlinear regimes can explain the recovery of 
$C$ at 
$A>5$ in figure 6(a). Although the overall dynamics are much more chaotic in the nonlinear regime, the contribution of the coherent interactions to the system's energy effectively appears only in the first PC in both regimes and therefore the relative contribution of the interaction to the first PC is commonly scaled by 
$B_d$. This result also implies that the underlying coherence in the nonlinear regime represents the perfect correlation (synchronized oscillations) like that of linear, monodisperse clouds.
3.6. Evaluation of the coherence measure
To assess the variation of the coherence measure dependent on the three regimes, in figures 7(a)–7(c), we plot 
$B_1$ against 
$B_d$ for 
$O(10^3)$ clusters with various values of 
$s_d$, 
$N$ and 
$R_C$, at the three distinct excitation amplitudes (
$A=2\times 10^{-2}$, 
$1.2$ and 
$20$). Appendix D summarizes the parameters used. With the weak excitation (figure 7a), bubble oscillations are in a linear regime and 
$B_d\approx B_0$. The data points are collapsed on the line of 
$C=1$ for 
$s_d=0.1$, while data points are scattered for the other values of 
$s_d$. This result is expected as 
$B_d$ was originally defined to scale the coherence of the monodisperse, perfectly correlated bubbles in the linear regime (d'Agostino & Brennen Reference d'Agostino and Brennen1989; Zeravcic et al. Reference Zeravcic, Lohse and Van Saarloos2011). With the intermediate excitation (figure 7b), the bubble oscillations are nonlinear. The data points are scattered from 
$C=1$, regardless of the value of 
$s_d$. In this regime, the result indicates that the coherence is lost due to the nonlinear dynamics, regardless of polydispersity. Surprisingly, with the strong excitation (figure 7c), the data points are collapsed on the line of 
$C=1$, meaning that the variance of the first PC is scaled by 
$B_d$, regardless of the parameters. This collapse is not observed for the second PC, regardless of the parameters (Appendix E). The results suggest that the scaling of 
$B_1$ is universal for cavitating clusters, which are typically excited at 
$O(1)$ MPa and above. Physically speaking, the scaling indicates that in cavitating clusters, the energy is partitioned into a single coherent mode of correlated (synchronized) oscillations and incoherent modes generalizing the aforementioned interpretation of figures 6(c) and 6(d). The coherent energy is controlled by the mean field originally derived for monodisperse, near-equilibrium bubbles. It is suggested that this partitioning and scaling are universal regardless of nuclei's polydispersity and the degree of nonlinearity. This finding explains the successful use of 
$B_d$ in characterizing and controlling seemingly disordered clusters in our previous studies (Maeda & Colonius Reference Maeda and Colonius2019; Maeda & Maxwell Reference Maeda and Maxwell2021). The isolation of coherence implies that the details of the microscopic scale as well as of the many-body interactions represented by higher-order PCs could be modelled as fast variables, which force the macroscopic (relevant) scales in the form of noise. A rigorous way to corroborate this hypothesis might be the Mori–Zwanzig projection operator formalism (Mori Reference Mori1965) by assuming that the fluctuations of fast and slow variables are uncorrelated. Although the choice of relevant slow variables may be unclear, the present approach is promising since the WPCA-TD can be seen as a prompter for more refined models of the evolution of slow variables.
Figure 7. 
$B_1$ against 
$B_d$ for 
$O(10^3)$ clusters with various values of 
$s_d$, 
$N$, 
$R_C$ and other physical parameters (Appendix D) in the (a) linear (
$A=2\times 10^{-2}$), (b) transition (
$A=1.2$) and (c) nonlinear (
$A=20$) regimes.
4. Significance of the WPCA-TD to inform on cavitation physics
In general, PCA can be considered as a mathematical technique for data decomposition and the physical meaning of extracted features can be left to one's interpretation. To the authors’ knowledge, studies of complex fluid flows employing data decomposition and feature extraction techniques often rely on one's intuition to discuss the physical meaning of features. In fact, in § 3.5, the analyses of the PC spectra required rigorous interpretations based on previous knowledge on cavitation physics. Meanwhile, the correlation between 
$B_1$ and 
$B_d$ identified in the last section shows that the present WPCA-TD can directly extract 
$B_d$ from data without additional signal processing like spectral filtering, which are often required for Fourier-based analysis, despite the presence of noisy and incoherent features. To generalize, WPCA-TD is not only a tool for data decomposition and feature extraction, but also can be used to directly inform on the non-dimensional number that controls the coherent physics of interests. This informative aspect of the WPCA-TD can provide a meaningful shortcut to address the physics of cloud cavitation.
To perform WPCA-TD, bubble dynamics data can be obtained from experimental measurements as well as other numerical approaches (e.g. Kameda & Matsumoto Reference Kameda and Matsumoto1996; Maeda & Colonius Reference Maeda and Colonius2018; Pishchalnikov et al. Reference Pishchalnikov, Behnke-Parks, Schmidmayer, Maeda, Colonius, Kenny and Laser2019). The cluster's shape can be arbitrary. Although a bubble's translation and deformation are neglected in the present numerical experiments, WPCA-TD can incorporate dynamical variables controlling these effects (e.g. Ilinskii, Hamilton & Zabolotskaya Reference Ilinskii, Hamilton and Zabolotskaya2007; Murakami, Gaudron & Johnsen Reference Murakami, Gaudron and Johnsen2020). Physically meaningful data decomposition requires the fine temporal resolution of individual bubble dynamics. In practice, such information would be difficult to obtain in experiments except for a small number of bubbles in a highly controlled environment. We thus emphasize that WPCA-TD would primarily be useful to process fine temporal resolution of numerical data.
5. Conclusion
In conclusion, to corroborate the coherent dynamics of bubble clusters, we introduced and demonstrated WPCA-TD, a method of PCA to comprehensively decompose the time series of Lagrangian bubble dynamics data into coherent dynamical features, in way similar to the modal decomposition of Eulerian flow field data. The data are pre-processed such that the PC-variance of the features becomes consistent with the hydrodynamic potential energy induced by bubble oscillations that is predicted by the R–P equation. By analysing simulation datasets of clusters under harmonic excitation, we demonstrated that the coherent energy and the degree of incoherence are respectively quantified by the variance and the spectral entropy. The coherence was lowered by disorders induced by the nuclei's polydispersity and nonlinear response of bubbles, as expected. Meanwhile, in cavitating regimes, underlying, correlated modes of oscillations were isolated in a single dominant feature. The variance of this feature was found to be controlled by the previously identified mean-field parameter, 
$B_d$, regardless of the disordering factors, indicating that the underlying coherent dynamics may be universal in cavitating clusters. These results suggest that the method can provide a simplified and comprehensive representation of complex bubble dynamics. Analogous to the use of POD for reduced-order modelling (ROM) of various flows, the dynamical features extracted by PCA may be used for ROM, without directly solving many-body interactions. Such a model may be promising for controlling cavitation and nucleation without tracking individual nuclei.
Acknowledgments
Both authors are appreciative of Professor Tim Colonius's mentoring during their days at Caltech. We are pleased to offer this paper as a small token of our appreciation on the occasion of his 60th birthday.
Funding
K.M. acknowledges support from SRB Co. Inc. and Purdue University. Some of the computation presented here used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF under grant TG-CTS190009, as well as Anvil at Purdue University through allocation PHY220130 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants #2138259, #2138286, #2138307, #2137603 and #2138296.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Variable transformation for a single bubble
In this section, we provide supplemental discussions and justification on the variable transformation used in the WPCA-TD, by considering the single bubble dynamics. By transformation 
$(\xi,\eta )=(R\dot {R}, R)$, with 
$\xi$ being nothing but the velocity potential evaluated at the bubble surface, the R–P equation can be expressed as
$$\begin{gather} \dot{\xi}=-\frac{1}{2}\left(\frac{\xi}{\eta}\right)^2+\frac{p_b-p_\infty(t)}{\rho}, \end{gather}$$
$$\begin{gather}\dot{\eta}=\frac{\xi}{\eta}, \end{gather}$$
where 
$p_b$ and 
$p_\infty$ are the pressure inside the bubble and that in infinity, respectively. The instantaneous energy of the system is expressed as
\begin{equation} K=2{\rm \pi}\rho \dot{R}^2R^3=2{\rm \pi}\rho \xi^2\eta. \end{equation}
Next, we approximate the system (A1). Given the time series data, the temporal average of 
$\eta$ can be computed as 
$\langle \eta \rangle$. Using 
$\langle \eta \rangle$, we approximate the system of (A1) on a space spanned by a set of new variables 
$(\hat {\xi },\hat {\eta })$:
$$\begin{gather}
\dot{\hat{\xi}}=-\frac{1}{2}\Bigg(\frac{\hat{\xi}}{\langle
\eta\rangle}\Bigg)^2+\frac{p_b-p_\infty(t)}{\rho},
\end{gather}$$
$$\begin{gather}\dot{\hat{\eta}}=\frac{\hat{\xi}}{\hat{\eta}}. \end{gather}$$
This system models well the dynamical features of the R–P equation including bifurcation, although 
$\hat {\xi }$ and 
$\hat {\eta }$ are now partially decoupled. Our explanation is the following. In the linear regime, the quadratic term of (A4) can be neglected and this system becomes identical with the original one. In the nonlinear, cavitating regime, the bubble size (
$\eta$) varies slowly near its peak after the explosive cavitation growth. Sporadic, fast collapse events do not influence much the temporal mean of the bubble size. Therefore, 
$\eta$ is close to 
$\langle \eta \rangle$ except at collapse events and, even if 
$\eta$ changes rapidly during the collapse and rebound, this change has a relatively small influence on the mean behaviour of the system at the time scale of (statistically) stationary bubble oscillations.
We can express the instantaneous energy of this approximate system as
\begin{equation} \hat{K}=2{\rm \pi}\rho \hat{\xi}^2\hat{\eta}. \end{equation}
With the same initial conditions (
$\xi =\hat {\xi }$ and 
$\eta =\hat {\eta }$ at 
$t=0$), the temporal mean of the energy of the original system is approximated as
\begin{equation} \langle K\rangle \approx\langle \hat{K}\rangle = 2{\rm \pi}\rho \langle\hat{\xi}^2 \hat{\eta}\rangle \approx 2{\rm \pi}\rho \langle\hat{\xi}^2\rangle\langle \eta\rangle. \end{equation}Note that Preston et al. (Reference Preston, Colonius and Brennen2007) used a POD-based analysis of single bubble dynamics for reduced-order modelling of heat and mass diffusion across the bubble interface. In the study, the temperature and concentration fields were obtained by solving partial differential equations and then represented by POD modes, following the POD/Galerkin framework. This is distinct from the present PCA (POD) of Lagrangian bubble dynamics data based on the direct transformation and projection of the Rayleigh–Plesset equation.
Appendix B. Schematic of the input-output procedure
Figure 8 shows the schematic of the present procedure to obtain the spectral entropy, PC-variance and coherence measure from the input data including the time series of the radius and radial velocity of Lagrangian bubbles.
Figure 8. Comprehensive schematic of the present procedure to obtain the spectral entropy, PC-variance and coherence measure from input data.
For computing 
$C$, we introduced the alternative weight, 
$\boldsymbol {\varOmega }'$. Here, 
$\boldsymbol {\varOmega }'$ is defined through 
$\boldsymbol {P}': \boldsymbol {P}'(\langle \boldsymbol {s}\rangle )=\boldsymbol {\varOmega }'\boldsymbol {\varOmega }'^*$, where 
$\boldsymbol {P}'$ is the diagonal matrix with its entries from 
$\boldsymbol {P}$:
\begin{equation} P'_{ij}(\eta) = \begin{cases} 2{\rm \pi}\rho \eta_i, & i=j,\\ 0, & i\ne j. \end{cases} \end{equation}
The corresponding inner product, 
$\boldsymbol {\xi }^T\boldsymbol {P}'\boldsymbol {\xi }$, represents the portion of the energy of the system excluding the contribution of inter-bubble interactions:
\begin{equation} K'= \boldsymbol{\xi}^T\boldsymbol{P}'\boldsymbol{\xi}, \end{equation}where
\begin{equation} K' = 2{\rm \pi}\rho_l\sum^N_{i=1} R_i^3\dot{R}_i^2+(H.O.T). \end{equation}Appendix C. Details of the K–M equation used for data preparation
In this section, we provide details of the formulation of the Keller–Miksis (K–M) equation for multiple bubbles, which is used to generate the Lagrangian bubble dynamics data. Various extensions of the Rayleigh–Plesset equation and its applications to a system of multiple bubbles are available (Takahira, Akamatsu & Fujikawa Reference Takahira, Akamatsu and Fujikawa1994; Doinikov Reference Doinikov2004; Ilinskii et al. Reference Ilinskii, Hamilton and Zabolotskaya2007; Yasui et al. Reference Yasui, Iida, Tuziuti, Kozuka and Towata2008). Our formulation can be derived from the K–M equation extended for multiple bubbles. First, we recall the general formulation for the oscillations of interacting spherical bubbles in weakly compressible liquid with arbitrary inter-bubble distances, whose derivation is, for instance, provided in Appendix 2.4 of Fuster & Colonius (Reference Fuster and Colonius2011). In the present study, we use a simplified version of this derivation. To recall, Fuster & Colonius (Reference Fuster and Colonius2011) describe the radial evolution of bubble 
$i$ as
\begin{equation} \ddot{R}_i\left(R_i\left(1-\frac{\dot{R}_i}{c}\right)\right)+\frac{3}{2}R^2_i\left(1-\frac{\dot{R}_i}{3c}\right)=F^*+I^*, \end{equation}
where 
$F^*$ and 
$I^*$ represent the forcing due to the external potential and the inter-bubble interaction, respectively, and are expressed as
\begin{equation} F^*=\frac{\partial\phi_\infty}{\partial t}\left(1-\frac{\dot{R}_i}{c}\right)+\frac{R_i}{c}\frac{\partial^2\phi_\infty}{\partial t^2}+H_i\left(1+\frac{\dot{R}_i}{c}\right)+\frac{R_i\dot{H}_i}{c} \end{equation}and
\begin{align} I^* &=\sum^N_{j\ne
i}\left[\left(1+\frac{\dot{R}_i}{c}\right)\frac{\partial\phi_j(R_i)}{\partial
t}\right] +\frac{R_i}{c}\sum^N_{j\ne i}
\Bigg[\Bigg(1+\frac{\dot{R}_j(t')}{c}\Bigg)\Bigg(1-\frac{\dot{R}_j}{c}\Bigg)\frac{\partial^2\phi_j(t'-R_j(t')/c)}{\partial
t'^2}\Bigg]
\end{align}
\begin{align} &\quad -\sum^N_{j\ne
i}\Bigg[\Bigg(1+\frac{\dot{R}_j(t')}{c}\Bigg)\frac{\partial\phi_j(t'-R_j(t')/c)}{\partial
t'}\frac{R_i}{R_j}\frac{\dot{R}_j}{c}\Bigg]\!.
\end{align}
Here, 
$\phi _\infty$ is the velocity potential of liquid at infinity, and 
$\phi _i(R_i)$ and 
$H_i$ are the potential and the enthalpy of liquid evaluated at the surface of bubble 
$i$; 
$t'$ is the retarded time defined as 
$t'=t-(d_{ij}-R_j)/c$, where 
$(d_{ij}-R_j)/c$ represents the travel time for the pressure wave to reach bubble 
$i$ from the surface of the bubble 
$j$ and 
$d_{ij}$ is the distance between the centres of those bubbles. In the sparse limit, the equation recovers the original Keller–Miksis equation (Keller & Miksis Reference Keller and Miksis1980).
To close the equations, several relations are considered. Using Bernoulli's equation, the potential for the bubble 
$i$ can be expressed as
\begin{equation}
\frac{\partial\phi_i(R_i)}{\partial
t}=-\Bigg(\frac{1}{2}\dot{R}^2_i+H_i+\sum^N_{j\ne
i}\frac{\partial\phi_j(R_i)}{\partial
t}+\frac{\partial\phi_\infty}{\partial t}\Bigg)\!.
\end{equation}
The velocity potentials at the bubble 
$i$ and bubble 
$j$ satisfy the following relation:
\begin{equation} \frac{\partial\phi_i(R_i)}{\partial t}=\frac{R_j(t')}{d_{ij}}\frac{\partial\phi_j(t'-R_j/c)}{\partial t'}. \end{equation}The enthalpy and the potential derivative are approximated as
$$\begin{gather} H_i\approx\frac{p_i-p_0}{\rho_{l,0}}, \end{gather}$$
$$\begin{gather}\frac{\partial\phi_\infty}{\partial t}\approx\frac{p_\infty-p_0}{\rho_{l,0}}. \end{gather}$$
Finally, 
$p_\infty$ is obtained using the information of the background Eulerian field computed on a mesh.
In the present study, we simplify 
$F^*$ and 
$I^*$ by invoking further approximations. First, consider that the bubble cluster size is smaller than the characteristic length-scale of the pressure wave in the field, 
$\lambda _c$. The inter-bubble distance in the cluster is naturally smaller than 
$\lambda _c$:
\begin{equation} d_{ij}-R_j\ll\lambda_c. \end{equation}
Dividing both sides with 
$c$,
\begin{equation} \frac{d_{ij}-R_j}{c}=t-t'\ll\frac{\lambda_c}{c}=\frac{1}{f_c}=T_c, \end{equation}
where 
$f_c$ and 
$T_c$ are the characteristic frequency and the period of the wave in the field. Therefore, the difference between 
$t$ and 
$t'$ is much smaller than the characteristic time scale of the dynamics. Given this knowledge, we approximate that 
$t'\approx t$. Second, we neglect the terms of the order of 
$(\dot {R}/c)^2$. Using these approximations, we can simplify 
$F^*$ and 
$I^*$ as
\begin{equation} F^*=\frac{p_i-p_\infty}{\rho_l}\left(1+\frac{\dot{R}_i}{c}\right)+\frac{R_i}{\rho_lc}\frac{\partial(\,p_i-p_\infty)}{\partial t} \end{equation}and
\begin{align} &I^* \end{align}
\begin{align} &\quad =-\sum^N_{j\ne
i}\frac{R_j}{d_{ij}}\Bigg(\frac{1}{2}\dot{R}^2_j+\ddot{R}_j\Bigg(R_j\Bigg(1-\frac{\dot{R}_j}{c}\Bigg)\Bigg)
+\frac{3}{2}\dot{R}^2_j\Bigg(1-\frac{\dot{R}_j}{3c}\Bigg)-\frac{p_j-p_\infty}{\rho_l}\frac{\dot{R}_j}{c}
\nonumber\\ &\qquad - \frac{R_j}{\rho_l
c}\frac{\partial(\,p_j-p_\infty)}{\partial t}\Bigg)\!.
\end{align}For harmonic excitation,
\begin{equation} p_\infty=p_a\sin(\omega t). \end{equation}We use the polytropic law to describe the pressure of the gas inside each bubble (Brennen Reference Brennen2014):
\begin{equation} p_i=p_{i,0}\left(\frac{R_{i,0}}{R_i}\right)^{3\gamma}, \end{equation}
where 
$\gamma$ is the constant polytropic exponent. Equations (C1), (C11) and (C13), together with (C14) and (C15), provide a complete system of ordinary differential equations for the radius of 
$N$ interacting bubbles, which can be readily solved with a given initial condition.
Appendix D. Summary of the parameters
Table 1 summarizes a set of parameters used to construct the database analysed the main manuscript. The number of bubbles 
$N$, cluster radius 
$R_c$, reference bubble radius 
$R_{ref}$, polydispersity measure 
$s_d$ and the forcing frequency were varied. The left table of parameters are used for the clusters in figure 4 (
$s_d=0.1,0.7$), figure 7(a) (all values of 
$s_d$) and figure 9(a) (all values of 
$s_d$). The right table of parameters are used for the clusters in figures 7(b), 7(c), 9(b) and 9(c).
Figure 9. 
$B_2$ against 
$B_d$ for 
$O(10^3)$ clusters with various values of 
$s_d$, 
$N$, 
$R_C$ and other physical parameters (Appendix D) in the (a) linear (
$A=2\times 10^{-2}$), (b) transition (
$A=1.2$) and (c) nonlinear (
$A=20$) regimes.
Table 1. Summary of the parameters used for the dataset of clusters analysed in the main manuscript.
Appendix E. Scaling of the second PC
Figure 9 shows the scaling of the interaction energy of the second dominant feature of bubble clouds in the dataset, 
$B_2=\sigma ^2_2/\sigma '^2_2 - 1$, for various values of excitation amplitude and nuclei polydispersity. For all parameters, data points are widely scattered and there is no clear correlation observed between 
$B_2$ and 
$B_d$, indicating that the second dominant feature does not represent coherent oscillations like those of the principal feature.
Open access
