2.1 Field solutions

The continuity relation in equation (2.1a,b ) implies that the pressure $p$ is a harmonic field. Similarly, the momentum equation shows that $\boldsymbol{s}$ is a gradient of harmonic functions. Also, the presence of two interfaces indicates that two different classes of spherical basis solutions for the Laplace equation should be considered to describe $p$ and $\boldsymbol{s}$ . Firstly, as the fluid domain is truncated by the outer surface, the fields should include regular solutions which remain finite in a finite domain but become infinite far away from the origin. Moreover, the presence of the bubble allows us to take into account singular solutions with an infinite value at the centre of the cavity, which is precluded from the domain of interest by the inner interface. Thus, we draw inspiration from Stokesian dynamics (Sangani & Mo Reference Sangani and Mo1996; Bhattacharya, Mishra & Bhattacharya Reference Bhattacharya, Mishra and Bhattacharya2010; Navardi & Bhattacharya Reference Navardi and Bhattacharya2013) to propose a general solution for equation (2.1a,b )

Here, $\unicode[STIX]{x1D719}_{lm}^{i+}$ and $\unicode[STIX]{x1D719}_{lm}^{i-}$ are the regular and singular basis functions in spherical coordinates with origin at the centre of the $i$ th interface. Accordingly, these are defined by

with $Y_{lm}^{(i)}$ being the spherical harmonics on the surface of the $i$ th sphere where the angular variations are defined by two usual quantum numbers $l=1,2,3,\ldots$ and $m=0,\pm 1,\pm 2,\ldots ,\pm l$ . The corresponding time-dependent amplitudes are $\unicode[STIX]{x1D6FC}_{lm}^{i\pm }$ , and overdots in equation (2.3) represent second-order temporal derivatives.

We use only $\unicode[STIX]{x1D719}_{lm}^{0+}$ and $\unicode[STIX]{x1D719}_{lm}^{1-}$ in (2.3) and (2.4). This is because the regular solution appearing due to the truncation of the domain by the external surface should be naturally expanded around its centre. Similarly, the singular function introduced due to the exclusion of the bubble volume should be expressed in the corresponding spherical coordinates. Still, the solutions of $p$ and $\boldsymbol{s}$ can have alternative representations than the one presented in (2.3) and (2.4). Such flexibility exists due to two transformation relations.

The regular basis functions can be expanded in Taylor series equating these to a linear combination of spherical solutions centred around another point. Hence,

where $T_{r}(R|l^{\prime }m^{\prime },lm)$ are the transformation coefficients which represent the Taylor series expansion of $\unicode[STIX]{x1D719}_{lm}^{0+}$ at the centre of the cavity. At this point, we choose the $z$ axis to be aligned with the line joining the centres of the drop and the bubble. Under such a construction, $T_{r}$ uncouples $m$ and $m^{\prime }$ , and assumes the following simplified form

Here, $\tilde{a}=a_{1}/a_{0}$ , $\tilde{R}=R/a_{0}$ and $\unicode[STIX]{x1D6FF}_{mm^{\prime }}$ is the Kronecker delta implying the block diagonalization of different azimuthal modes defined by $m$ .

Similar transformation between the singular solutions can also be derived using Taylor series expansion, although the presence of the singularity limits its validity by the convergence radius. Hence, for points outside the sphere of radius $R$ around the centre of the drop:

The new transformation coefficients are obtained as

where the earlier choice of the $z$ -axis is taken for granted.

The complete description of the dynamics requires evaluation of the amplitudes $\unicode[STIX]{x1D6FC}_{lm}^{i\pm }(t)$ . This is achieved by satisfying the boundary conditions in (2.2), where we use the completeness of the solutions in (2.3) and (2.4) along with the flexibility given by (2.6)–(2.9).

2.2 Matrix formulation

We combine the boundary condition (2.2), the general solutions (2.3) and (2.4) and transformations (2.6)–(2.9) to derive two sets of ordinary differential equations describing $\unicode[STIX]{x1D6FC}_{lm}^{0+}$ and $\unicode[STIX]{x1D6FC}_{lm}^{1-}$ simultaneously. The first set of these relations comes from the surface tension dictated stress discontinuity at the outer surface of the drop

The second group represents the same condition at the interface of the embedded bubble:

These equations should be considered for all possible $l$ and $m$ .

We cast equations (2.10) and (2.11) into a matrix form which has a block-diagonal structure, uncoupling mutual effects between two modes with different $m$ . Accordingly, one finds

Here, components in column $|\!\unicode[STIX]{x1D6FC}_{m}^{i\pm }\!\rangle$ are amplitudes for specific $m$ but different $l$ . Elements of square matrices $[T^{f}]$ and $[T^{r}]$ are the transformation coefficients $T_{f}(l,l^{\prime }|R,m)$ and $T_{r}(l,l^{\prime }|R,m)$ . Also, $[\tilde{T}^{f}]$ and $[\tilde{T}^{r}]$ are related to $[T^{f}]$ and $[T^{r}]$

whereas $[D^{0}]$ and $[D^{1}]$ have only non-zero components in diagonal positions so that

The identity operator in equation (2.12) is shown as $[I]$ .

We combine $|\!\unicode[STIX]{x1D6FC}_{m}^{0+}\!\rangle$ and $|\!\unicode[STIX]{x1D6FC}_{m}^{1-}\!\rangle$ into a single column notation $|\!\unicode[STIX]{x1D6FC}_{m}\!\rangle$ which is the transpose of $\{\langle \!\unicode[STIX]{x1D6FC}_{m}^{0+}\! |\langle \!\unicode[STIX]{x1D6FC}_{m}^{1-}\!|\}$ so that equation (2.12) can be written in a more concise form:

where

We refer to $[M(R,m)]$ and $[S(R,m)]$ as the mass and spring matrices, respectively. The spectra of the interfacial pulsation of the system are directly dependent on $[M(R,m)]$ and $[S(R,m)]$ . For a given $m$ , temporal variation in any $k$ th mode of vibration is expressed as a sinusoidal behaviour:

where $\unicode[STIX]{x1D714}_{km}$ is the corresponding natural frequency. When this variation is inserted into (2.15), $|\!\unicode[STIX]{x1D6FC}_{km}(0)\!\rangle$ and $\unicode[STIX]{x1D714}_{km}^{2}$ can be identified as the eigen vectors and eigen values of the matrix $[P]=[M]^{-1}[S]$ . In our computation, we construct $[M]$ , $[S]$ and $[P]$ for arbitrary bubble position $R$ and quantum number $m$ . Then, the eigen values of $[P]$ are calculated to reveal the prominent frequencies of the system.

It is to be noted here that the described matrix formulation reveals an inherent coupling between the dynamics of the two interfaces. This interdependence is quantified by the off-diagonal elements of $[M]$ and $[S]$ matrices represented by the blocks $[T^{f}]$ , $[T^{r}]$ , $[\tilde{T}^{f}]$ , $[\tilde{T}^{r}]$ . As the components of the off-diagonal blocks vary as power functions of $\tilde{a}$ , the effect of one surface on another becomes negligible for small void fractions. This is manifested in our prior work (Bhattacharya Reference Bhattacharya2016) where it is seen that deformability of the bubble is insignificant for $\tilde{a}<0.3$ . Such behaviour is corroborated by experimental studies also (Obreschkow et al. Reference Obreschkow, Kobel, Dorsaz, Bosset, Nicollier and Farhat2006).

2.3 Test for spectral convergence

For a specific $m$ value, $l$ can be any positive integer as long as $l\geqslant |m|$ . Therefore, $[M]$ , $[S]$ and $[P]$ are infinite dimensional matrices in reality. However, for practical simulations, their dimensionality has to be truncated by enforcing a maximum value for $l$ given by $l_{mx}+|m|-1$ . This spectral truncation makes the dimensionality of all relevant matrices $2l_{mx}$ due to the $2\times 2$ block shown in (2.16). As a result, we expect $2l_{mx}$ number of eigen values for $[P]$ leading to the same number of natural frequencies.

For a concentric geometry, these are classified into $(+)$ and $(-)$ configurations in an earlier paper (Bhattacharya Reference Bhattacharya2016) accounting for a twofold increase in the number of possible modes compared to the pure drop analysed by Rayleigh. For $(+)$ states, the outer surface deforms relatively less, and maintains perfect alignment with the bubble. In contrast, the $(-)$ states show stronger waves at the drop periphery pulsating in the perpendicular orientation with respect to the cavity. The two types of interfacial waves are referred to as the bubble and the sloshing modes, respectively.

For an eccentrically placed system also, one can identify these modes, each of which reveals $l_{mx}$ number of frequencies after the spectral truncation. These values in each set for given $m$ are arranged in ascending order numbered by a subscript $k$ , and referred to as $\unicode[STIX]{x1D714}_{km}^{\pm }$ according to their origin from either the $(+)$ or $(-)$ group.

We are interested in the first few of these values for $k=2,3$ , where it is to be noted that $\unicode[STIX]{x1D714}_{1m}^{\pm }$ are identically $0$ for translational symmetry of the system (Bhattacharya Reference Bhattacharya2016). The spectral truncation, however, introduces errors in the desired outputs. Fortunately, $\unicode[STIX]{x1D714}_{km}^{\pm }$ for any given $m$ should have negligible error due to the truncation in matrix dimension as long as $l_{mx}$ is adequately greater than $k$ . Evidence of such spectral convergence indicates the correctness in both mathematical analysis and numerical simulations.

Accordingly, we define a non-dimensional frequency $\tilde{\unicode[STIX]{x1D714}}_{km}^{\pm }$ as

In figure 2, computed values of $\tilde{\unicode[STIX]{x1D714}}_{km}^{\pm }$ are plotted as functions of $l_{mx}$ for a non-axisymmetric perturbation represented by $m=1$ . We also choose two geometric configurations with different size and position of the bubble. In both, we try to ensure close proximity between two interfaces, as this potentially causes slower spectral convergence. Even then, the curves show a fast saturation of computed results for all cases when $l_{mx}\geqslant 7$ . The figure also reveals slower saturation of $(-)$ modes for larger bubbles, as solutions with higher $l$ are needed to properly represent the misaligned waves at two closely located but eccentrically placed interfaces. Conversely, $(+)$ modes converge relatively later for a smaller cavity because its deflation normal to the surface makes the peripheral variations more difficult capture by basis functions.

Figure 2.

Computed values of normalized natural frequencies for $k=2,3$ and $m=1$ are plotted as functions of spectral truncation order $l_{mx}$ for two geometric configurations with (1) $R=0.4a_{0}$ and $a_{1}=0.5a_{0}$ (thick line), as well as (2) $R=0.25a_{0}$ and $a_{1}=0.7a_{0}$ (thin line).

2.4 Parametric study of natural frequency revealing fine-structure split

For a concentric position of the bubble with $R=0$ , matrices $[M]$ , $[S]$ and $[P]$ do not have any dependence on the quantum number $m$ . This happens because the dynamic equation for a spherically symmetric unperturbed geometry should not have any systemic azimuthal variation. Consequently, a number of vibrational modes with different $m$ exhibit the exact same natural frequency when the centres of the drop and the cavity coincide. We refer to this phenomenon as azimuthal degeneracy.

For an eccentrically placed bubble, however, the azimuthal degeneracy disappears, because $[M]$ , $[S]$ and $[P]$ depend on $m$ when $R\neq 0$ . As a result, one can see small variations in natural frequencies $\unicode[STIX]{x1D714}_{km}^{\pm }$ for the modes with different $m$ but same $k$ . This is the fine-structure split in the wave spectra similar to the one seen in atomic physics, where small changes in energy level occur for different states with different quantum numbers.

In figures 3 and 4, we illustrate the aforementioned fine-structure split in the frequency spectra for two different size ratios $\tilde{a}$ . For all computations, $l_{mx}$ is considered as $15$ , so that convergence in connection to spectral truncation is not an issue. In each panel, a specific $k$ is chosen, and different $\unicode[STIX]{x1D714}_{km}^{+}$ or $\unicode[STIX]{x1D714}_{km}^{-}$ are plotted as functions of normalized $R$ for different $m$ values. For concentric position, $k$ matches exactly with the quantum number $l$ . Hence, for each $k$ , we expect $(2k+1)$ number of modes with different $m$ as it can take the values $0,\pm 1,\pm 2,\ldots ,\pm k$ . Among these possible $m$ values, however, modes with $+m$ and $-m$ have the exact same frequency and exact same contour of the deformed interface – these are only separated by a temporal phase. Hence, for each $k$ , one can find $k+1$ number of fine deviations in frequencies for different $m$ values. This dictates the number of curves for respective $k$ values in figures 3 and 4.

Figure 3.

Normalized natural frequencies for $l=2$ and $l=3$ as well as plus and minus configurational modes are plotted as functions of position of the bubble normalized by its maximum possible value when the size ratio $\tilde{a}=0.5$ and $m=0$ (solid line), $m=1$ (dashed line), $m=2$ (dash-dotted line) and $m=3$ (dotted line).

Figure 4.

Same as figure 3 with size ratio $\tilde{a}$ is $0.7$ instead of $0.5$ .

To highlight the fine-structure split, $\tilde{\unicode[STIX]{x1D714}}_{km}^{\pm }$ are divided by their values $\tilde{\unicode[STIX]{x1D714}}_{km}^{0\pm }$ for concentric position at $R=0$ . The normalizing factor $\tilde{\unicode[STIX]{x1D714}}_{km}^{0\pm }$ is well described in our recently published paper (Bhattacharya Reference Bhattacharya2016) and corroborated in an earlier paper (Tsamopoulos & Brown Reference Tsamopoulos and Brown1987). For example, we find $\tilde{\unicode[STIX]{x1D714}}_{2m}^{0-}$ to be $2.74$ and $2.29$ for $\tilde{a}=0.5$ and $\tilde{a}=0.7$ . Similarly, our calculations show $\tilde{\unicode[STIX]{x1D714}}_{2m}^{0+}$ to be $10.11$ and $7.32$ for similar size ratios. These values match exactly with the ones reported in Tsamopoulos & Brown (Reference Tsamopoulos and Brown1987). In figures 3 and 4, we use these quantities to renormalize the presented curves which start from $1$ . Then, a gradual departure for different $m$ values with increasing $R$ can be seen in the plots. For small departure from the spherically symmetric geometry, the frequencies should vary as parabolic functions of $R$ . This fact is evident in figures 3 and 4. We analyse this behaviour in more detail in a perturbation analysis presented in § 3, where the general numerical simulation is independently corroborated. Also, we consult the results from the aforementioned earlier work focusing on a concentric geometry to check whether the new findings match exactly with the known values when $R\rightarrow 0$ .

When $R$ approaches $a_{0}-a_{1}$ , two interfaces nearly touch each other at the location of minimum separation. It should be noted that touching surfaces lead to the rupture of the film (Patzer & Homsy Reference Patzer and Homsy1975), and the entire analysis is no longer valid from physical considerations. Still, we want to investigate this situation from a purely mathematical perspective. For closely situated interfaces, it is well known that methods based on basis function expansions typically have issues with spectral convergence (Bhattacharya et al. Reference Bhattacharya, Mishra and Bhattacharya2010; Navardi & Bhattacharya Reference Navardi and Bhattacharya2013). In Stokesian dynamics (Bhattacharya Reference Bhattacharya2008; Navardi & Bhattacharya Reference Navardi and Bhattacharya2010a ,Reference Navardi and Bhattacharya b ,Reference Navardi and Bhattacharya c ), this problem is addressed by lubrication analysis, where the predominance of the contributions to hydrodynamic stress from the contact region is presumed. In the present problem, however, the frequency of oscillation is not dictated by the lubrication flow due to its negligible inertia. As a result, one cannot provide a quantitative explanation for spectral behaviour in such a situation. Despite this difficulty, the perturbation analysis in § 3 can qualitatively show why the frequency should vary linearly from a fixed value when $R\rightarrow (a_{0}-a_{1})$ .

The frequencies for different $m$ but the same $k$ diverge considerably with increasing $R$ . Still, their differences are significantly smaller than the difference between two modes with unequal $k$ or different configurational modes $(+)$ or $(-)$ . Thus, the spectral plot for interfacial waves would show band structures in frequency space with values for the same $k$ and configurational modes grouping together in close proximity. The width and location of these bands change with $R$ and $\tilde{a}$ . Such a realization can lead to an efficient characterization of multiphase drops. In § 4, we outline an algorithm for this purpose.

3.1 Expansion of matrices representing small change in bubble position

Our perturbation analysis analytically describes the change in the frequency values for a specific mode due to a small change in the bubble position $R$ . This requires an expansion of the relevant eigen value equation involving the matrix $[P]$ defined in § 2.2.

Let us assume that the $i$ th eigen value $\unicode[STIX]{x1D706}_{0}^{(i)}$ for the eigen vector $|\!e_{0}^{(i)}\!\rangle$ of $[P]$ is known for a specific position of the bubble at $R=R_{0}$ . The quantities associated with this specific location are denoted by the subscript $0$ , so that $[P]$ at $R_{0}$ is referred to as $[P_{0}]$ which satisfies

The question is what the change in the eigen value would be if the position of the bubble were perturbed by a short distance $\unicode[STIX]{x1D6FF}$ .

In the new position also, one can write the similar eigen value equation

where the subscript $0$ is dropped to indicate the new position. We consider a perturbative expansion of equation (3.2) treating $\unicode[STIX]{x1D6FF}$ as a small parameter

where the subscript $k$ represents the $k$ th-order term in the respective perturbations.

We insert the expansion given in (3.3) into (3.2) to derive a hierarchy of equations corresponding to different perturbative orders. The leading zeroth-order relation is already given by (3.1). The next order for $i=1$ is governed by the following

The row of eigen vectors $\langle f_{0}^{(j)} |$ for the matrix $P_{0}$ is also the reciprocal basis for the column vectors. Hence, $\langle f_{0}^{(j)} |$ satisfies two conditions simultaneously

We contract equation (3.4) from the left with $\langle f_{0}^{(j)}\! |$ . When the properties of $\langle f_{0}^{(j)}\! |$ in (3.5) are utilized, the predominant correction in the eigen values can be derived:

This is the slope by which the eigen value $\unicode[STIX]{x1D706}^{(i)}$ varies with the position of the bubble.

The frequency of the $i$ th mode, renamed here as $\unicode[STIX]{x1D714}^{(i)}$ instead of $\unicode[STIX]{x1D714}_{km}^{\pm }$ , is the square root of $\unicode[STIX]{x1D706}^{(i)}$ . Let us now assume $R_{0}\rightarrow (a_{0}-a_{1})$ with the two interfaces nearly touching each other. Then, one expects a linear dependence of $\unicode[STIX]{x1D714}^{(i)}$ on the thickness $h$ of the so called lubrication region around minimum separation. Thus, the variation would be

where the constant $\unicode[STIX]{x1D6FA}_{1}^{(i)}$ is given by

Unfortunately, the leading-order eigen values under touching conditions with $h\rightarrow 0$ is not readily available from a simplified mathematical analysis. One has to solve the field equations in the tangent sphere coordinate system so that the corresponding matrix formulation ultimately yields these quantities. However, such a result would not have any physical meaning, as the bubble will simply rupture in such a configuration. This is why we did not pursue a detailed calculation when the two interfaces nearly touch each other. Still, from a purely mathematical perspective, the perturbation analysis provides a qualitative agreement, as all curves in figures 3 and 4 obey the functional form in equation (3.7). This means that the change $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}^{(i)}$ in frequency should scale as

instead of the scaling $\sqrt{\unicode[STIX]{x1D6FE}/(\unicode[STIX]{x1D70C}a_{0}^{3})}$ seen in (2.18).

3.2 Asymptotic results under the limit of nearly concentric geometry

The first-order correction in equation (3.6) should not be enough to explain the leading-order variation of frequencies with $R$ for a nearly concentric geometry. The reason is based on the reflectional symmetry in the $x$ – $y$ plane making $\unicode[STIX]{x1D706}_{1}^{(i)}$ equal to zero when $R_{0}=0$ . In other words, if the bubble is displaced upward or downward by the same distance with respect to the centre of the drop, the resulting two systems will be essentially identical apart from the configurational orientation. Thus, in both cases, the wave frequencies should be exactly the same, making the eigen value an even function of $R$ . This is why we have seen a parabolic dependence of all $\unicode[STIX]{x1D714}_{km}^{\pm }$ in figures 3 and 4 for $R\rightarrow 0$ . Quantitative description of this variation requires an analysis of the second-order perturbation.

The second-order relation for the eigen value equation is

We contract equation (3.10) from the left with $\langle f_{0}^{(i)}\! |$ , and use its properties to find

From (3.10), one can conclude

Here, the superscript $BI$ denotes partial or block inverse of singular matrix $[P_{0}-\unicode[STIX]{x1D706}_{0}^{(i)}I]$ with its non-null block being identified and inverted without changing its null space

Combining, equations (3.11) and (3.12), the second-order correction in the eigen value is obtained

where $\unicode[STIX]{x1D706}_{1}^{(i)}=0$ is presumed for $R_{0}=0$ .

Fortunately, for $R_{0}=0$ , all matrices $[P_{0}]$ , $[P_{1}]$ , $[P_{2}]$ have simple forms. Our analysis first constructs $[P_{0}]$ , $[P_{1}]$ , $[P_{2}]$ in terms of the perturbation terms in the mass matrix $[M]$ and spring matrix $[S]$ . Accordingly, the expansions for $[M]$ and $[S]$ are defined

and the $i$ th-order contributions $[M_{i}]$ and $[S_{i}]$ are explicitly determined using the $R$ -dependent expression of the respective matrices. The explicit forms of the $i$ th-order terms are given below

and for $i>0$

Here, the components of the off-diagonal blocks of $[M_{i}]$ are only non-zero when $l=l^{\prime }+i$ for which

Similarly, only non-zero elements in the off-diagonal blocks of $[S_{i}]$ can be evaluated by relating them to $[M_{i}]$

In contrast, the diagonal blocks of the matrices are either the identity matrix $[I]$ or the null matrix $[0]$ or the matrices defined in (2.14).

For our problem, $[P]$ is the product of the inverse of $[M]$ and $[S]$ . This definition yields:

For $R_{0}=0$ , we checked that the specific forms of $\langle f_{0}^{(i)}\! |$ , $[P_{1}]$ and $|\!e_{0}^{(i)}\!\rangle$ make $\unicode[STIX]{x1D706}_{1}^{(i)}=0$ according to (3.6). Consequently, $\unicode[STIX]{x1D706}_{2}^{(i)}$ contributes as the most prominent correction. Then, we can express the non-dimensional frequency in the following form

where the dimensionless parameter $\tilde{\unicode[STIX]{x1D6FA}}_{km}^{\pm }$ represents the quadratic behaviour in $R$

Here, the superscript $(i)$ is interchangeably used in place of quantum numbers.

Figure 5.

The parameter $\tilde{\unicode[STIX]{x1D6FA}}_{lm}$ computed from (3.24) according to perturbation theory is plotted as a function of $\tilde{a}$ for plus and minus spectral modes as well as $l=2$ and $l=3$ , when $m=0$ (solid line), $m=1$ (dotted line), $m=2$ (dashed line) and $m=3$ (dash-dotted line). The results from the general simulation presented in § 2 are shown as large dotted points for a few selected size ratios.

In figure 5, we present the calculated $\tilde{\unicode[STIX]{x1D6FA}}_{km}^{\pm }$ according to the perturbation analysis as a function of the size ratio. The results from the asymptotic theory are presented as the line curves for different quantum numbers. In contrast, the corresponding values from the general simulation, discussed in § 2, are shown as isolated data points for a few specific $\tilde{a}$ . The latter match exactly with the former, establishing an important corroboration between the two independent approaches, validating both.

4.1 Outline of the detection and characterization algorithm

The detection method would require three major steps. Firstly, the interfacial pulsation has to be recorded in an experiment (Saffren, Elleman & Rhim Reference Saffren, Elleman and Rhim1982). Secondly, the relevant spectral quantities have to be evaluated from the stored data. Finally, the convenient features of experimentally obtained values should be properly matched with the theoretical prediction of a possible system by an efficient algorithm.

The first two steps are conceptually straightforward. One can use a high-speed camera to view the pulsation of the drop. As a result, its radial deformation $s^{ex}$ at the outer surface can be recorded as a series of two-dimensional images taken at subsequent time intervals. The observed displacement from the average contour has to be stored as a function of time and polar coordinate $\unicode[STIX]{x1D6FD}$ on the viewing plane. Then, the frequency spectra of the interfacial dynamics can be experimentally determined by a dual Fourier transform $C_{\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})$

Complete spectral description of the system requires $C_{\unicode[STIX]{x1D708}}$ to be evaluated for $\unicode[STIX]{x1D708}=0,1,2,\ldots$ , as natural frequencies for all $l$ values would be revealed by $C_{\unicode[STIX]{x1D708}}(\unicode[STIX]{x1D714})$ . Consequently, spectral plots of $C_{\unicode[STIX]{x1D708}}$ in the frequency domain will show several spikes grouped together in bands according to our theory of fine-structure split.

This spectral manifestation can be used to predict the volume fraction $\unicode[STIX]{x1D719}$ and position $R$ of the cavity inside the drop by matching the observation with the theoretical predictions. For this purpose, an efficient algorithm can be built after identifying useful features of the variation of frequency values for different modes in parametric space. The strategy would be to recognize quantities which only vary with one system parameter while showing negligible dependence on others. From the spectral results in § 2, we identify such signature features which are discussed in the subsequent subsections.

4.2 Convenient features of fine spectral structure indicating the volume fraction

Figure 6.

Computed $\tilde{\unicode[STIX]{x1D714}}_{kk}^{-}$ are plotted as functions of volume fraction with $R/a_{0}=0.2$ (thin line) and $R/a_{0}=0.4,0.6,0.8$ (thick line) for $k=2,3$ .

The efficient characterization algorithm requires recognition of spectral features which only varies with either of position $R$ or volume fraction $\unicode[STIX]{x1D719}=\tilde{a}^{3}$ . Such quantities are readily available in the plots presented in § 2.4. One can see that the frequencies for the $(-)$ modes corresponding to $k=|m|$ have negligible variation with $R$ . This is especially true for higher values of $l$ . The slow dependence in the fields along $z$ -direction for $k=m$ is the reason behind the observed behaviour. Thus, one can identify $\unicode[STIX]{x1D714}_{kk}^{-}$ to be a signature variable which should be conveniently used to find the volume fraction $\unicode[STIX]{x1D719}$ .

The approximate invariance of $\unicode[STIX]{x1D714}_{kk}^{-}$ with $R$ is further revealed in figure 6. We choose four different values of $R$ and present the non-dimensional $\tilde{\unicode[STIX]{x1D714}}_{22}^{-}$ and $\tilde{\unicode[STIX]{x1D714}}_{33}^{-}$ as functions of volume fraction of the cavity. The curves in the figures are of different lengths because allowable volume fractions are limited for an eccentrically situated bubble. For $k=m=2$ , three of these curves merge with each other altogether. For $k=m=3$ , the invariance with $R$ is even more prominent – there, all four curves collapse onto each other showing no change whatsoever for variation in $R$ .

This implies that $\tilde{\unicode[STIX]{x1D714}}_{33}^{-}$ can be treated as a signature parameter which solely depends on $\unicode[STIX]{x1D719}$ . So if this variation is tabulated properly, one can predict the volume fraction from the observed value of the $\tilde{\unicode[STIX]{x1D714}}_{33}^{-}$ .

However, in order to characterize the bubble, $\tilde{\unicode[STIX]{x1D714}}_{33}^{-}$ has to be located in the experimentally observed spectral data. The second band in the spectral plot corresponds to the $k=3$ and $(-)$ configurational mode. The last peak in that cluster represents the frequency value for $m=k=3$ . Thus, if one locates this quantity in the experimental data acquired following the procedure in § 4.1, the volume fraction for the cavity can be detected.

4.3 Detection of bubble position from bandwidth

Once the volume fraction is known, it is easy to find the position of the bubble from various other quantities. Among all others, the most reliable and convenient parameter could be the width $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}_{k}^{\pm }$ of the first few prominent bands. This quantity is the interval in the spectral space in which the fine-structure split ensures a number of modes for the same specific $l$ to be clustered together. One can evaluate the width of this band from the difference between the maximum and minimum of closely situated frequencies for different $m$ but the same $l$ .

In figure 7, we present bandwidths for four most relevant modes as functions of the normalized position of the bubble for two different volume fractions. The plots show non-trivial variation in $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}_{k}^{\pm }$ with $R$ .

Figure 7.

Bandwidths of non-dimensional frequency are plotted as functions of normalized position for $(-)$ configurational modes with $k=2$ (solid line) and $k=3$ (dash-dot line) as well as $(+)$ configurational modes with $k=2$ (dashed line) and $k=3$ (dotted line), where the size ratio is $\tilde{a}=0.5$ (left) and $\tilde{a}=0.7$ (right).

The bands for specific $k$ can be identified in the experimentally acquired spectral data by the number of fine structures in that group. This number should be exactly the same as $k+1$ . Also, there would be two bands with same number of fine structures. Between these two, the one with the lower value represents the $(-)$ configurational mode, whereas the $(+)$ mode corresponds to the larger number. In this way it is possible to locate $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D714}_{k}^{\pm }$ in the observed frequency spectra.

These quantities can reveal the position of the bubble in the following way. Firstly, the theory outlined in § 2 can be used to compute the bandwidth as a function of $R$ for the known value of the already detected volume fraction. Then, an error $\unicode[STIX]{x1D716}$ can be defined as

where superscript ‘ $ex$ ’ describes the experimentally observed bandwidth. One can choose a few representative $k$ so that $\unicode[STIX]{x1D716}$ relates to a mean square departure in the simulated result from the experimental values. In the end, an error versus position plot can indicate the location of the cavity by presuming that it corresponds to the minimum  $\unicode[STIX]{x1D716}$ .

It is to be noted that the position of the bubble may change in time due to various reasons. Firstly, it can move due to the presence of an external potential such gravity. In that case, a low Bond number limit would ensure a very slow temporal change in $R$ . Secondly, the system itself might be unstable, causing departure due to nonlinear shape oscillation (Shaw Reference Shaw2006). However, the time required for such translation is always much larger than the temporal scale dictated by the interfacial frequencies (Tsamopoulos & Brown Reference Tsamopoulos and Brown1987). Thus, for either case, one can apply a multiple time scale analysis to describe the complex dynamics. In such a formulation, the leading-order spectra for nominal fluctuations would be similar to our findings. Then, additionally, the natural frequencies would also alter slowly over a longer period where this slow evolution can indicate the translational dynamics inside the liquid domain.

Table 1.

The comparison between the values of the system parameters used as inputs for the dynamic evolution and the estimations of the same quantities obtained from the detection algorithm.

We perform a numerical experiment to confirm the viability of the described detection technique. Accordingly, a bubble–drop system is considered where an input configuration with specific $\tilde{a}$ and $R$ values is assigned. Then, it is assumed that $126$ different modes corresponding to $k$ values varying from $2$ to $10$ are involved in the interfacial pulsation. These include waves with different $m=0,1,\ldots ,k$ and configurational orientations $+$ and $-$ for each $k$ . The non-dimensional amplitudes for each of the modes are randomly generated numbers between $0.01$ and $0.05$ . Once these random values are fixed, the deformation at the outer surface is calculated using the field solutions equation (2.4) given in § 2.1. Consequently, the projection of the instantaneous image of the deforming drop is computed taking the $y$ – $z$ surface to be the two-dimensional viewing plane. The radial deviation $s^{ex}$ of the image is recorded as a function of time and the polar angle on the viewing plane. Next, the spectral response $C_{\unicode[STIX]{x1D708}}$ is calculated using the dual Fourier transform in (4.1), revealing the natural frequencies of the system as the resulting peaks in the spectral space. Finally, $\tilde{a}$ and $R$ are estimated from these data by exploiting known features of $\tilde{\unicode[STIX]{x1D714}}_{33}^{-}$ as well as minimizing the errors according to (4.2) in the observed and computed bandwidths.

The estimated and the input values of $\tilde{a}$ and $R$ are presented in table 1. We select $5$ different systems for our case study so that it covers the possible parametric space properly. Accordingly, the first simulation has moderate $R$ and $\tilde{a}$ identified as $0.25$ and $0.45$ , respectively. Then, we choose two extreme values for $R$ as $0.1$ and $0.4$ representing a slightly eccentric bubble and one nearly touching the external surface. Similarly, two additional $\tilde{a}$ values are also accounted for to evaluate the effectiveness of our scheme for both small and large cavities. The extreme values of one parameter are taken into account keeping the other parameter moderate so that ultimately $5$ different configurations can include all possible situations.

The results for the numerical experiment show the correctness of the detection algorithm. A majority of the case studies show less than $1\,\%$ relative error between the actual input and the estimated prediction. The error is lower for larger bubbles and for more eccentric locations. This is understandable because such configurations enhance the interactions between the interfaces, leading to a more accurate detection due to more prominent surface dynamics. The comparative exploration is a necessary consistency test inferring the outlined diagnostic technique to be robust and implementable as long as the caveats of our analysis are satisfied and experiments are accurate enough.