Microtomographic particle image velocimetry measurements of viscoelastic instabilities in a three-dimensional microcontraction

Abstract Viscoelastic flow through an abrupt planar contraction geometry above a certain Weissenberg number ($Wi$) is well known to become unstable upstream of the contraction plane via a central jet separating from the walls and forming vortices in the salient corners. Here, for the first time, we consider three-dimensional (3-D) viscoelastic contraction flows in a microfabricated glass square–square contraction geometry. We employ state-of-the-art microtomographic particle image velocimetry to produce time-resolved and volumetric quantification of the 3-D viscoelastic instabilities arising in a dilute polymer solution driven through the geometry over a wide range of $Wi$ but at negligible Reynolds number. Based on our observations, we describe new insights into the growth, propagation and transient dynamics of an elastic vortex formed upstream of the 3-D microcontraction due to flow jetting towards the contraction. At low $Wi$ we observe vortex growth for increasing $Wi$, followed by a previously unreported vortex growth plateau region. In the plateau region, the vortex circulates around the jet with a period that decreases with $Wi$ but an amplitude that is independent of $Wi$. In addition, we report new out-of-plane asymmetric jetting behaviour with a phase-wise dependence on $Wi$. Finally, we resolve the rate-of-strain tensor $\boldsymbol{\mathsf{D}}$ and ascribe local gradients in $\boldsymbol{\mathsf{D}}$ as the underlying driver of circulation via strain hardening of the fluid in the wake of the jet.

developing computational models capable of studying highly elastic flows (Afonso et al. 2011;Pimenta & Alves 2017;Alves, Oliveira & Pinho 2021). Under negligible inertia (i.e. Reynolds numbers Re 1), for Weissenberg numbers (Wi = λγ , where λ is the fluid relaxation time andγ the shear rate) beyond a critical value Wi c ≈ 0.5, pipe flow moving towards a contraction becomes sufficiently elastic that it separates from the upstream walls, forming a central 'jet' that enters the constriction and vortices around the mouth of the constriction (McKinley et al. 1991;Rothstein & McKinley 1999). Initially the corner vortices are static in their placement as Wi is increased, but they eventually grow in size until a Hopf bifurcation characterized by a periodic fluctuation of the vortex separation point occurs. For certain contraction ratios β (β = w 0 /w c , the length scale ratio of uniform channel w 0 to the contraction width w c ; see figure 1a) a lip vortex may form for Wi < Wi c (Giesekus 1968), as summarized by Rothstein & McKinley (2001), developing into a corner vortex with increasing Wi. For Wi > Wi p Wi c , the corner vortices become increasingly unsteady for increasing Wi (McKinley et al. 1991;Rothstein & McKinley 1999, and may lead to chaos through period-or frequency-doubling routes (McKinley et al. 1991;Afonso et al. 2011). Rothstein & McKinley (2001) studied the flow of a polymer solution through an axisymmetric contraction-expansion geometry with a contraction ratio β = 4. For Wi = 3.5, they noted a jet of flow separating from the corner vortex and moving through the contraction as the vortex itself precessed in the azimuthal direction. Based on laser Doppler velocimetry (LDV) measurements, both the jet and the vortex precession shared the same fundamental frequency. From visual observations, they postulated that the jet was sustained along the upstream length of the vortex, such that the region of flow formed a helical path towards the contraction. Alves, Pinho & Oliveira (2005, 2008 performed flow visualization to support their numerical results with a square-square contraction (β = 4) channel to quantify a nonlinear vortex growth with Wi, showing a local decrease in vortex size near Wi = 10 but subsequent monotonic growth for higher Wi. Both works described an unsteady azimuthal flow instability for Wi > 52 with a fundamental frequency that scaled linearly with Wi. Interestingly, they reported a flow inversion whereby at low Wi particles entered the vortex from the corner-normal plane (z = ±y, with x downstream) and exited into the contraction at the midplane (z = 0), but reversed this pattern at higher Wi. Through simulations they showed that the flow inversion Wi was model-dependent (in some cases Wi ≤ 5.8) and that the phenomenon could be attributed to strong extensional effects, as it was not captured in purely shear-thinning simulations. Finally, they noted that a fully three-dimensional (3-D) study was required.
Experiments by Sousa et al. (2009Sousa et al. ( , 2011 for a square-square contraction with various β described similar azimuthal vortex precession above Wi = 62 for β = 4, and they characterized the fundamental frequency via flow visualization. They observed that for β > 4, the fundamental frequency was independent of Wi, while it grew with log(Wi) for β = 4. Both studies saw similar flow inversion behaviour as Alves et al. (2008), inverting at Wi = 2 for β = 4. Alves et al. (2008) and Sousa et al. (2011) compared experiments against various numerical approaches, but restricted their studies to steady-state solutions at high Wi due to the numerical burden of solving 3-D transient viscoelastic flows (see the review article of Alves et al. (2021)). Afonso et al. (2011) were able to solve a β = 4 contraction flow in transience up to Wi = 100 in two dimensions, and in steady state up to Wi = 5000 in three dimensions, via the log conformation approach (Fattal & Kupferman 2005) using the Oldroyd-B (Oldroyd 1950) and Phan-Thien-Tanner (PTT; Thien & Tanner 1977) constitutive models. They noted the expected upstream vortex growth with Wi, and used the two-dimensional (2-D) transient solutions to show frequency-doubling behaviour leading to a chaotic-like regime. In three dimensions, they showed that the contraction flow included regions of strong extension and calculated a steady-state flow inversion, supporting the experimental findings of Alves et al. (2008) and Sousa et al. (2009).
The onset and subsequent dynamics of this elastic flow instability are highly sensitive to the contraction geometry and fluid rheology, and it has been shown by Alves & Poole (2007) that the extensional rheology in particular dictates the instability rather than shear-thinning effects. For planar contraction geometries, the numerical simulations of Alves & Poole (2007) showed that Wi p 2 ∼ β, as predicted by the Pakdel-McKinley criterion (McKinley, Pakdel & Öztekin 1996). This indicates that the onset of instability is driven by elastic stress growth on the curved streamlines near the contraction throat.
More recently, viscoelastic contraction flow has received attention at the microscale, where inertia can often be neglected and elastic effects at high Wi become dominant even for low-viscosity fluids (see the thorough review provided by Rodd et al. (2007)). In the negligible-Re regime, there is an apparent void of knowledge regarding the transient nature of 3-D flow at moderate to high Wi. Although several studies have provided a qualitative visual description of intriguing 3-D viscoelastic flow phenomena in 3-D contractions, due to the difficulty of resolving such flows at the microscale and the immense computational burden of solving such flows numerically, this flow type remains to be quantitatively described.
923 R6-3 To date, elastic contraction flow has been studied primarily via localized or planar measurements such as particle image velocimetry (PIV), LDV or streak imagery. Global pressure measurements have also been employed to provide insights into drag reduction. In recent years, tomographic particle image velocimetry (TPIV) has received increasing attention as a method whereby 3-D flow volumes can be resolved via the reconstruction of a particle-laden flow from overlapping lines of sight, followed by cross-correlation between subsequent particle volumes (Elsinga et al. 2006a). This method can also be applied at the microscale (µ-TPIV), with multiple lines of sight provided by stereomicroscopy. Holographic particle image velocimetry (HPIV) has also shown success in taking volumetric viscoelastic flow measurements in microscale geometries (Qin et al. 2019(Qin et al. , 2020, reporting a bistable negative wake ahead of a cylinder and out-of-plane instability modes along the flow separatrix of a cross-channel. However, HPIV is quite limited in terms of volume depth and reconstruction resolution compared to µ-TPIV (Schäfer & Schröder 2011). Nonetheless, the novel results of Qin et al. (2019Qin et al. ( , 2020 suggest that, despite decades of research on fundamental viscoelastic flows, deep insights are still to be elucidated once out-of-plane dynamics are captured. Here, using a dilute solution of a high-molecular-weight polymer, we report the first investigations of 3-D viscoelastic contraction flow at the microscale using the µ-TPIV method. We focus on the range of Wi encompassing the transition from the vortex growth regime (which is accompanied by the growth of a steady central jet) to the onset of periodic vortex precession (which is accompanied by the circulation of the jet). We demonstrate that the circulation of the jet has a phase-wise asymmetry that depends on the nominal Wi. By fully resolving the 3-D velocity field, we can assess the true velocity-gradient tensor and thus the rate-of-strain tensor. We show that the precession of the corner vortex is driven by the central jet continuously retreating from regions of increased rate of strain, and hypothesize that the underlying driving mechanism is localized strain hardening of the polymer solution.

Flow cell and viscoelastic fluid
The experiments were conducted in a square-sectioned contraction-expansion flow cell (figure 1a) fabricated from fused silica glass via selective laser-induced etching (Gottmann, Hermans & Ortmann 2012) using a commercial LightFab 3-D printer (LightFab GmbH). This process can resolve features on the micrometre scale, with a surface root-mean-square (r.m.s.) value of approximately 1 µm (Pimenta et al. 2020). We measured the channel width and height from an X-ray microtomography scan (figure 1a) as w 0 = 860 ± 10 µm outside the contraction and w c = 255 ± 5 µm inside the contraction, yielding a contraction ratio of β = 3.4. Figure 1(a) displays our dimensionless coordinate system, where each component is reduced by w c /2 (e.g. x * = 2x/w c ). Flow was stepped to each velocity by two syringe pumps in a push-pull configuration.
The viscoelastic test fluid was a polymeric solution composed of 107 parts per million (ppm) partially hydrolysed polyacrylamide (HPAA; M W = 18 MDa, Polysciences Inc., USA) in a solvent of 85 wt% glycerol and 15 wt% deionized water. The refractive index of the fluid is closely matched to that of the fused silica flow cell. An Anton-Paar MCR 502 stress-controlled rheometer was used with a cone-and-plate geometry (50 mm diameter, 1 • angle) to characterize the shear viscosity of the fluid under steady shear. Figure 1(b) shows that the fluid is weakly shear thinning and has a zero-shear viscosity of η 0 = 184 mPa s. 923 R6-4 The relaxation time of the fluid was measured as λ = 0.65 s using a capillary breakup extensional rheometer (Haake CaBER, Thermo Fisher Scientific; see Anna & McKinley 2001) fitted with endplates of diameter d 0 = 6 mm. We plot the ratio of the apparent extensional viscosity η E,app to the zero-shear viscosity η 0 against the accumulated Hencky strain H = 2 ln(d 0 /d(t)) in figure 1(c). The fluid exhibits strong strain hardening with η E,app ≈ 400η 0 at high strains. The zero-shear viscosity η 0 and the characteristic length scale w c /2 are used to calculate Re = ρU c w c /2η 0 and Wi = λγ = 2λU c /w c (where U c is the average flow velocity in the contraction). In the present work, Re 10 −2 and as such is considered negligible.
The flow was recorded as double-frame images captured at 12 Hz, with a flow-rate-dependent time interval between laser pulses t such that no particle moved more than 8 pixels. Frames were preprocessed with local background subtraction and Gaussian smoothing at 3 × 3 pixels. A 3-D calibration was performed by capturing reference images of a microgrid at the planes z = ±450 µm and z = 0 µm, fully encompassing the depth of the flow cell (w 0 ), and a coordinate system was interpolated between these planes using a third-order polynomial. Particle positions in three dimensions were reconstructed from the images using four iterations of the Fast MART (Multiplicative Algebraic Reconstruction Technique) algorithm implemented in the commercial PIV software DaVis 10.1.2 (Lavision GmbH). Fast MART initializes the particle volume using the multiplicative line-of-sight routine (MLOS; Worth & Nickels 2008;Atkinson & Soria 2009), followed by iterations of Sequential MART (SMART; Atkinson & Soria 2009). We concluded the algorithm with five iterations of the motion tracking enhancement (MTE) method (Novara, Batenburg & Scarano 2010;Lynch & Scarano 2015) to reduce spurious 'ghost' particles which arise from randomly overlapping lines of sight (Elsinga, Van Oudheusden & Scarano 2006b) and thus do not correlate in time.
Volume self-calibration (Wieneke 2008) was employed to improve the accuracy of reconstruction. Particle displacements between particle volumes were obtained using a multigrid iterative cross-correlation technique, with the final pass at 32 × 32 × 32 voxels with 75 % overlap for a vector grid of 31 µm. To reduce measurement noise, the vector field was spatiotemporally filtered with a second-order polynomial regression across neighbourhoods of 5 3 vectors in space and extended through five increments in time for a total kernel size of 5 4 points. This polynomial regression visibly reduced measurement noise while predominantly preserving space-time resolution owing to the small kernel size: the filtering wavelength is substantially smaller than the flow dynamics reported in this work. This is a common approach to denoising TPIV data (Scarano & Poelma 2009;Elsinga et al. 2010;Schneiders, Scarano & Elsinga 2017). Ultimately we resolved 1600 flow volumes per recording: a duration of 133 s (200λ).
Uncertainty quantification for TPIV is a topic of ongoing research (Atkinson et al. 2011;Sciacchitano 2019), but a priori comparisons of experimental velocity fields to direct numerical solutions or synthetic reconstructions have yielded a TPIV uncertainty of the order of 0.1-0.3 pixels (Atkinson et al. 2011). As an a posteriori assessment of our measurement quality, we validated conservation of mass for time-averaged flow volumes (Zhang, Tao & Katz 1997). Flow divergence ∇ · u was calculated at each experimentally determined velocity vector u, and the error assessed relative to an assumption of incompressible flow (relative error ζ = (∇ · u) 2 /tr(∇u · ∇u)). The value ζ averaged 0.25 for 3 ≤ Wi ≤ 87. Divergence error relative to the magnitude of vorticity was 0.14 at Wi = 87, in good agreement with the value of 0.2 obtained from a three-camera TPIV experiment by Kempaiah et al. (2020).

Results and discussion
Measurements were taken over a range of flow velocities encompassing 1.5 ≤ Wi ≤ 87 for a region of interest upstream of the contraction. Note that the downstream (expansion) side was imaged in a separate series of experiments, but flow remained steady across the Wi range investigated. Throughout the discussion of the flow field kinematics, we non-dimensionalize lengths by w c /2 and velocities by U c . Deformation rates are hence reduced by 2U c /w c . Times are non-dimensionalized either by the fluid relaxation time λ or by the fundamental period of circulation in the case of periodic flows. All non-dimensional quantities are indicated by a superscript '*'.

Steady flow at low Wi
As flow approaches the contraction at Wi ≤ 5, a steady separation point forms where the flow separates from the walls of the channel and passes as a central jet through the constriction. Evidently, the minimum Wi achieved in our experiments precludes possible observation of lip vortices, which occur at lower Wi than the onset of jetting behaviour and the appearance of corner vortices (e.g. Rothstein & McKinley 2001). Additionally, although we do not track particle positions in our measurements, it was visually observed that flow entered the corner vortex from the midplane, i.e. the inverted state found by Alves et al. (2008) and Sousa et al. (2009), indicating strong extensional effects. Figure 2 presents (a) average isosurfaces and (b) a midplane x * -y * slice of the streamwise velocity u * x for Wi = 5. The pink isosurface shows u * x = 1/3 (i.e. the central jet), while the grey surface marks u * x = 0 (i.e. the edges of the recirculating regions). Flow separation follows the upstream down zero-crossing of u * x ; the central jet is characterized by positive u * x , while the corner vortices drive negative u * x backflow along the walls towards the separation point.

The onset of periodic instability
For increasing Wi, the corner vortex propagates upstream but remains steady, until, for Wi > Wi p , the flow transitions to a periodic instability characterized by axial fluctuation of the upstream separation point and a circumferential precession of the central jet (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2021.620 for animations of the jet for 5 ≤ Wi ≤ 44). Note that a similar helical flow pattern was previously postulated by Rothstein & McKinley (2001) for an axisymmetric contraction. They observed it as azimuthal motion of a local jet separating from the lip near the corner vortex, and suggested that the jet lined the boundary of the corner vortex. Our results validate their hypothesis: a jet starts from the upstream separation point, propagates through the  Our results expand on prior experiments in axisymmetric abrupt contractions which saw the vortex size increase monotonically for increasing Wi, although at a lower Wi range than probed here due to larger length scales involved (e.g. Wi < 5 for McKinley et al. (1991) and Wi < 8 for Rothstein & McKinley (1999)). In planar microfluidic contraction flow experiments, Rodd et al. (2007) reported a decrease in L * v for a single data point at the maximum Wi = 24 achieved, but were unable to extrapolate the trend. Numerical work by Comminal et al. (2016) revealed an L * v plateau accompanied by periodic vortex annihilation for Wi > 14 in a 2-D contraction, which they attributed to an accumulation of elastic strain upstream of the contraction. However, they acknowledged that their use of the Oldroyd-B constitutive model (Oldroyd 1950) lacks physical mechanisms (such as finite extensibility) which would otherwise limit elastic stress.

Out-of-plane jet dynamics
As shown in the x * -t * space-time diagrams in figure 3(a-e), the flow instability is strongly periodic for Wi ≥ 11. Thus, to collapse the dataset and further reduce noise, we deploy time-synchronous averaging (TSA) to reduce the time dimension to a single average cycle. This method is further discussed in Bechhoefer & Kingsley (2009). We isolated the time series of velocity magnitude at a single point in the volume, and then used the local maxima in the time series at that point to segregate each period at all points in the volume. The cycles are mapped to a normalized phase time φ * ranging from 0 to 1 (i.e. one cycle), and averaged into f s T bins, where f s is the sampling frequency (12 Hz) and T = T * /λ is the average period (T * shown in figure 3f ). Henceforth, phase-averaged quantities will be marked by .
We used the phase-binned data to investigate the 3-D trajectory of the central jet passing through the toroidal corner vortex. We found the maximum flow velocity in each y * -z * plane along the x * axis for −L * v ≤ x * ≤ 0, and took the y * -z * coordinates of the peak velocity as the location of the core of the jet in x * . Thus, we extract x * -y * -z * trajectories of the circulating inner jet as it approaches the contraction. To quantify fluctuations in the jet velocity throughout phase time, we calculate the normalized fluctuating velocity along In this way Û is not affected by the accelerating flow velocity as the jet approaches the contraction, and is normalized to compare between different Wi values.
Data of Û are presented in figure 4(a) for Wi ≥ 5, along with snapshots of the location of the jet and planar projections of the swept area of the circulating jet core in figure 4(b,c) for Wi = 44 and Wi = 87 (trajectories in figure 4(b,c) are coloured by their φ * value in figure 4(a) for comparison). Note that we observe a phase-wise asymmetry of the jet in both the normalized fluctuating component of velocity and the jet location. For each Wi, the jet starts at the same place in the channel at φ * = 0; however, the directionality of the precession about x * is seemingly random between experiments. Maxima and minima in Û are evident for Wi = 11 near φ * = 0.18, 0.88 and 0.56, and are apparent but slightly degraded by Wi = 22. The fact that the velocity fluctuation is small (less than 5 %) and almost symmetric about φ * = 0.5 indicates that the cyclical fluctuation in the 11 ≤ Wi ≤ 22 regime is likely to be geometric in origin as the channel has an approximately 5 µm (or approximately 0.02w c ) variation between the width and height. At higher Wi (e.g. Wi = 44 and Wi = 87), the jet had the same directionality as for Wi = 22 and was mapped to the  same locations in phase time, but strikingly the dual Û peaks are eliminated. Instead, we observe a single minimum and maximum for Û and the values occur at different phase times from those seen at lower velocities. Thus it is inferred that the peak fluctuation is independent of the y * -z * coordinate between different Wi.
To relate velocity fluctuations within the jet to the spatial distribution, we present in figure 4(b,c) trajectories of the jet core position in x * -y * -z * over phase time. Here we see that the minimum Û for Wi = 44 coincides with a skewed curvature of the jet towards +z * near φ * = 0.42, and the jet is comparatively tighter towards z * = 0, which aligns in time with the maximum Û near φ * = 0.78: the distribution of the jet varies inversely with the fluctuating velocity such that the flow rate is preserved. A similar effect is seen for Wi = 87, but with a shifted phase timing of Û , placing the asymmetry more in the x * -y * plane. Reconfiguration to single-sided asymmetry in the jet at Wi ≥ 44 coincides with a phenomenon previously discussed: the growth plateau in L * v from 11 ≤ Wi ≤ 22 despite a decreasing core period of the circulating jet (figure 3f ). Whereas the symmetrically swirling jet appears to prohibit vortex growth, a break to strong asymmetry above Wi = 22 proves a preferable route and L * v again increases for increasing Wi.
3.4. The role of the rate-of-strain tensor Dilute solutions of high-molecular-weight polymers are known to shear-thin under shear flow and to strain-harden under extensional deformations (Tirtaatmadja & Sridhar 1993;Solomon & Muller 1996) of the channel. We can gain qualitative insight into the relevance of strain hardening to the dynamics of our microcontraction flow by considering the local rate-of-strain tensor D = (∇ u + ∇ u T )/2. Here we rely solely on the µ-TPIV measurements to obtain the velocity vector field before calculating D. We take the magnitude of D as γ = √ 2(D : D) (reduced as γ * = γ w c /(2U c )) to highlight regions of high rate of strain where strain hardening of the HPAA is more likely.
We extracted y * -z * slices at an arbitrary x * = −8 for Wi = 87 to show a simplified perspective on the relationship between 0.5 γ * and the location of the jet (coloured contours of 0.5 U * , where U * = U /U c ) in figure 5(a-d). The individual panels progress along phase time φ * , with each plane including the 0.5 U * contour from the following phase time. Red arrows connect circular markers for the locations of maximum velocity on the plane between time steps with the tail and head markers showing the present and future positions, respectively. Two trends are noted: First, while the jet (the solid coloured contour) stays centred about γ * , the jet location forward in time is always towards the exterior of the γ contour. Second, a phase-wise asymmetry manifests in the distribution of γ * : φ * = 0.12 has low asymmetry in Û (figure 5b) as well as in the distribution of γ ; at φ * = 0.62, Û is highly deviated and this aligns with a strong mismatch in γ * about the circumference of the jet. Therefore, it appears that a mismatch in rate of strain about the jet can strongly influence the phase-wise progression of the jet as it circulates, with positive and negative fluctuations in Û accompanied by an imbalanced distribution of γ * .
In figure 6, we present isosurfaces of 0.5 γ * and 0.5 U * , the maximum rate of strain and flow velocity, for Wi = 87. Figure 6(a-d) show four time steps throughout a circulation, where the volume of high rate of strain forms a band about the circumference of the core of the jet (the pink isosurface). Moreover, the γ * isosurface extends further upstream on one side of the jet for all time steps, i.e. the rate of strain is greater along one side of the jet. An animated loop of the jet circulating with the rate-of-strain volume is shown in supplementary movie 2. We directly compare the rate of strain forward in time in figure 6(e) as projections of γ * from (a-d) from the −x * direction. Moving clockwise from the φ * = 0.12 isosurface, each surface forward in time presents a decrease in the rate of strain in the clockwise direction or, in other words, a perpetual retreat of the central jet from regions of increased rate of strain. The dynamics of elastic contraction flow has been reported to be sensitive to the extensional rheology (Rothstein & McKinley 2001; Alves & Poole 2007), and since our HPAA solution strain hardens (figure 1c), we can infer from the flow kinematics that the central jet moves about the contraction driven by gradients of strain-hardening HPAA, which will tend to follow the regions of increased rate of strain. This sheds new light onto the likely significance of extensional rheology for local dynamics in viscoelastic contraction flow, as experimentally describing flow topology via directly resolving the rate-of-strain tensor has not been achieved hitherto for elastic flow instability.

Conclusions
Using µ-TPIV, we have experimentally resolved for the first time the highly 3-D dynamics of viscoelastic flow through a square-square microcontraction at low Re and high Wi. We captured steady and periodic flow instability for a central jet of fluid passing through a toroidal vortex pinned about the contraction entrance, and observed a new vortex growth plateau where the period of instability decreases with Wi but the enveloped vortex volume stagnates. For low 11 ≤ Wi ≤ 22, the jet circulates with a symmetric cyclical velocity fluctuation likely originating from geometric imperfections. This region coincides with the vortex growth plateau. At higher Wi ≥ 44, a strong asymmetry in the jet forms, which corresponds to exiting the growth plateau. This would indicate that the asymmetric mode provides a preferable route to vortex growth. We determined the first experimental mapping of the full rate-of-strain tensor to transient dynamics for viscoelastic flow instabilities. We relate regions of increased rate of strain and the extensional rheology of the fluid to the directionality of the circulating jet. Gradients of strain hardening in the fluid provide a likely circulation mechanism as the jet retreats from locally strain-hardened regions of the flowing viscoelastic polymer solution, allowing us to gain new insight into the significance of extensional rheology to local dynamics of viscoelastic contraction flow.
Funding. We gratefully acknowledge the support of the Okinawa Institute of Science and Technology Graduate University (OIST) with subsidy funding from the Cabinet Office, Government of Japan. We also acknowledge financial support from the Japanese Society for the Promotion of Science (JSPS, Grant Nos. 21K14080, 18K03958, 18H01135 and 21K03884).