# Spontaneous Formation of Star-Shaped Surface Patterns
in a Driven Bose-Einstein Condensate

K. Kwon  Department of Physics, Korea Advanced Institute of Science and
Technology, Daejeon 34141, Korea  K. Mukherjee  Department of Physics, Indian
Institute of Technology Kharagpur, Kharagpur-721302, India  S. Huh  Department
of Physics, Korea Advanced Institute of Science and Technology, Daejeon 34141,
Korea  K. Kim  Department of Physics, Korea Advanced Institute of Science and
Technology, Daejeon 34141, Korea  S. I. Mistakidis  Center for Optical Quantum
Technologies, Department of Physics,University of Hamburg, Luruper Chaussee
149, 22761 Hamburg, Germany  D. K. Maity  Department of Physics, Indian
Institute of Technology Kharagpur, Kharagpur-721302, India  P. G. Kevrekidis
Department of Mathematics and Statistics, University of Massachusetts,
Amherst, Massachusetts 01003-4515, USA  S. Majumder  Department of Physics,
Indian Institute of Technology Kharagpur, Kharagpur-721302, India  P.
Schmelcher  Center for Optical Quantum Technologies, Department of
Physics,University of Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany
The Hamburg Centre for Ultrafast Imaging, University of Hamburg,Luruper
Chaussee 149, 22761 Hamburg, Germany  J.-y. Choi<EMAIL_ADDRESS>Department of Physics, Korea Advanced Institute of Science and Technology,
Daejeon 34141, Korea

###### Abstract

We observe experimentally the spontaneous formation of star-shaped surface
patterns in driven Bose-Einstein condensates. Two-dimensional star-shaped
patterns with $l$-fold symmetry, ranging from quadrupole ($l=2$) to heptagon
modes ($l=7$), are parametrically excited by modulating the scattering length
near the Feshbach resonance. An effective Mathieu equation and Floquet
analysis are utilized, relating the instability conditions to the dispersion
of the surface modes in a trapped superfluid. Identifying the resonant
frequencies of the patterns, we precisely measure the dispersion relation of
the collective excitations. The oscillation amplitude of the surface
excitations increases exponentially during the modulation. We find that only
the $l=6$ mode is unstable due to its emergent coupling with the dipole motion
of the cloud. Our experimental results are in excellent agreement with the
mean-field framework. Our work opens a new pathway for generating higher-lying
collective excitations with applications, such as the probing of exotic
properties of quantum fluids and providing a generation mechanism of quantum
turbulence.

††preprint: APS/123-QED

Introduction.– Spontaneous pattern formation is frequently encountered in
various research fields, including chemistry, biology, nonlinear optics, and
cosmology Cross and Hohenberg (1993); Maini  _et al._ (1997); Petrov  _et al._
(1997). Faraday waves constitute one the earliest and most celebrated examples
thereof that can be observed when a fluid in a vessel is subject to a vertical
periodic modulation Faraday (1831). The underlying mechanism of these
phenomena is the existence of instabilities, manifested in the related
nonlinear hydrodynamic equations. The instabilities characterize a dominant
wavelength that breaks the spatial and temporal symmetries of the system Cross
and Hohenberg (1993). This has applications in measuring the intrinsic
properties of fluids like density and surface tension Noblin  _et al._ (2005);
Shen  _et al._ (2010). Moreover, recent experiments with tracer particles in
the Faraday waves have revealed the emergence of two-dimensional (2D)
turbulence Xia  _et al._ (2011); von Kameke  _et al._ (2011); Francois  _et
al._ (2014) and self-organization of flows with various patterns Alarcón  _et
al._ (2020), thus extending the research scope of parametrically driven
systems.

Bose-Einstein condensates (BECs) of atomic gases offer a fertile platform for
transferring the relevant knowledge of nonlinear dynamics in classical
settings to the realm of quantum many-body systems Engels  _et al._ (2007);
Nguyen  _et al._ (2019); Zhang  _et al._ (2020); Tsubota  _et al._ (2013);
Kevrekidis  _et al._ (2015); Madeira  _et al._ (2020). One-dimensional Faraday
waves have been indeed observed in BECs under the periodic modulation of the
transverse trap frequency of an elongated condensate Engels  _et al._ (2007)
or of the $s$-wave scattering length near the magnetic Feshbach resonance
Nguyen  _et al._ (2019). Increasing the modulation strength with low driving
frequency, irregular patterns (granulation) are generated characterized by
fairly sizable quantum fluctuations Nguyen  _et al._ (2019), and bearing
features of quantum turbulence Lode  _et al._ (2021). However, the majority of
experimental efforts has been performed in one-dimension Engels  _et al._
(2007); Nguyen  _et al._ (2019). As such, various correlation patterns
emerging from non-linear wave mixing Zhang  _et al._ (2020), surface
excitations either from parametrically-driven multi-component systems Maity
_et al._ (2020) or from quantum fluctuations, such as quantum capillary waves
in optical lattices Rath  _et al._ (2011), and the relation of Faraday waves
to turbulent behavior in higher dimensions Seman  _et al._ (2011); Madeira
_et al._ (2020) are yet a largely unexplored territory.

Figure 1: (Color online) Experimental observation of star-shaped condensates.
(a) A BEC of 7Li atoms is prepared in a pancake shaped trap consisting of a
red-detuned optical trap for the axial confinement and a magnetic trap,
induced by Feshbach magnetic field curvature, for the radial confinement. The
scattering length is modulated by oscillating the magnetic field near the
Feshbach resonance. The mean modulation amplitude $\bar{a}_{m}$ is defined as
the average of upper ($a_{m}^{+}$) and lower ($a_{m}^{-}$) modulation peaks,
$\bar{a}_{m}=(a_{m}^{+}+a_{m}^{-})/2$. (b-g) Representative images (single-
shots) of the ensuing regular polygons with $D_{l}$ symmetry triggered by the
periodic modulation. The modulation frequencies are 84 Hz ($D_{2}$), 104 Hz
($D_{3}$), 119 Hz ($D_{4}$), 132 Hz ($D_{5}$), 147 Hz ($D_{6}$), and 161 Hz
($D_{7}$), respectively. The mean modulation amplitude $\bar{a}_{m}=19a_{B}$
($a_{m}^{+}=21a_{B}$ and $a_{m}^{-}=17a_{B}$). The scale bar in the $l=2$ mode
represents 100 $\mu{}m$. (h) The orientation angle $\phi^{*}$ for pentagon-
shaped BECs, defined in (e), is measured with 400 consecutive experimental
runs. The histogram displays the corresponding occurrence of the angle with
bin size of 3∘. (i) Hexagon and (j) heptagon shaped patterns obtained by
solving the 3D GPE (see main text).

In this Letter, we report the controlled generation of surface modes of
different wavenumbers in atomic quantum fluids. Specifically, 2D regular
polygons exhibiting an $l$-fold symmetry (from $l=2$ to $l=7$) develop from
radially symmetric condensates by modulating the atomic interactions near the
Feshbach resonance. The observed surface patterns have no preferred
orientation and oscillate sinusoidally according to the modulation frequency.
The associated spatial and temporal symmetry breaking phenomenon can be
understood in terms of the hydrodynamic (parametric) instability, where an
effective Mathieu equation describes the stability boundaries of the
individual patterns. Probing the instability regions, we accurately measure
the dispersion relation of the surface modes of the harmonically trapped
superfluid. The experimental results show an excellent agreement with the
predictions of the three-dimensional (3D) mean-field Gross-Pitaevskii equation
(GPE).

Experimental and Theoretical Setup.– The experiment is initiated by producing
a BEC of 7Li atoms in the $\ket{F=1,m_{F}=1}$ state near the Feshbach
resonance Kim  _et al._ (2019). The scattering length is set to
$a_{bg}=138(6)a_{B}$ ($a_{B}$ is the Bohr radius). The condensate resides in a
highly anisotropic harmonic trap consisting of a tight confining optical trap
in the axial direction Huh  _et al._ (2020) and a weak magnetic trap
constituting the radially symmetric confinement [Fig. 1(a)]. The trap
frequencies are measured to be
$[\omega_{r},\omega_{z}]=2\pi\times[29.4(2),725(5)$] Hz. Then, we apply a
sinusoidally oscillating magnetic field to the pancake shaped condensate,
which modulates the scattering length $a_{s}(t)$ of the atoms [Fig. 1(a)].
Following a modulation time $t$, we take in-situ absorption images under the
Feshbach magnetic field and measure the atomic density distribution.

After 1 s of modulation, the condensate boundary is strongly deformed,
displaying 2D regular polygon patterns along the $xy$-plane [Fig. 1(b)-(g)].
The generation of the surface modes is mostly driven by the scattering length
modulation as the oscillation amplitude of the radial trap frequency is non-
sizable, i.e., of about 0.05%. Moreover, the regular polygons show no
preferred orientation in the horizontal plane, manifesting the spatio-temporal
symmetry breaking phenomenon [Fig. 1(h)].

These surface modes are equally reproduced within the full 3D GPE framework
[Fig.1(i)-(j)] which reads

$\begin{split}i\hbar\frac{\partial}{\partial
t}\Psi&(x,y,z,t)=\bigg{(}-\frac{\hbar^{2}}{2m}\nabla_{\bf{r}}^{2}+\frac{1}{2}m\omega_{r}^{2}(x^{2}+y^{2}\\\
&+\lambda^{2}z^{2})+\frac{4\pi\hbar^{2}a_{s}(t)}{m}\absolutevalue{\Psi}^{2}\bigg{)}\Psi(x,y,z,t),\end{split}$
(1)

where $\int{\bf{dr}}{|\Psi|}^{2}=N$ and
$\nabla_{\bf{r}}^{2}\equiv\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z}^{2}$.
Also, $m$ and $\lambda=\omega_{z}/\omega_{r}$ represent the atomic mass and
the aspect ratio of the trap.

To emulate the thermal fraction in the experiments (less than 10$\%$), we
consider a weak amplitude perturbation to the ground state, $\Psi_{G}(x,y,z)$,
of the BEC. The initial wave function
$\Psi_{{\rm{in}}}(x,y,z)=\Psi_{G}(x,y,z)[1+\epsilon\delta(x,y,z)]$ Zheng  _et
al._ (2013a); Proukakis and Jackson (2008). Here, $\delta(x,y,z)$ is a
Gaussian random distribution having zero mean and variance unity produced by
using the so-called Box-Mueller algorithm Box and Muller (1958); Zheng  _et
al._ (2013b) and $\epsilon\ll 1$ mimics the thermal fraction being, herein, of
the order of $\epsilon\sim 0.1$, see Supplemental Material sup .

Subsequently, we let the system [described by Eq. (1)] evolve upon considering
a periodic modulation of the scattering length of the form
$a_{s}(t)=a_{bg}+\bar{a}_{m}\cos(\omega_{m}t)$, where $\bar{a}_{m}$ and
$\omega_{m}$ are the mean modulation amplitude and frequency, dictated by the
experiment. Figures 1(i) and (j) show characteristic density profiles,
$n(x,y,t)=\int dz\absolutevalue{\Psi(x,y,z,t)}^{2}$, of the $D_{6}$ and
$D_{7}$ star patterns, respectively, obtained within the GPE framework. In the
simulations, the surface modes are robust to the experimental imperfections,
such as anisotropy in the radial confinement ($\sim$ 3%), oscillations of the
radial curvature ($\sim$0.05%), and high harmonics in the scattering length
modulation protocol.

Results and Discussion.– To understand the underlying mechanism of the surface
deformation, we study the frequency dependence of the surface modes at a fixed
mean modulation amplitude $\bar{a}_{m}=19a_{B}$. We characterize the polygon
shaped BECs with $D_{l}$ symmetry, shown in Fig. 1(b)-(g), by the Fourier
amplitude $\mathcal{F}_{l}$ of the condensate radius over the azimuthal angle.
The $\mathcal{F}_{l}$ quantifies the displacement of the condensate boundary
with $l$-fold symmetry sup . The $D_{l}$ shaped-BECs varying from ellipses
($l=2$) to regular heptagons ($l=7$) can be identified. Figure 2(a) displays
the spectral peak $\mathcal{F}_{l}$ of each mode under various modulation
frequencies. The surface modes ($l=2-7$) are only excited at certain driving
frequency intervals, in a way strongly reminiscent of the tongues in the
Mathieu equation Kovacic  _et al._ (2018). The resonance curves are
asymmetric, resembling the response of a Duffing oscillator, which is well
represented for $l\leq 4$ modes. When we reduce the modulation amplitude, the
resonance spectra become more symmetric and acquire a narrower width sup ,
highlighting the role of non-linear interactions during the surface
deformation.

Figure 2: (Color online) Hydrodynamic instability of the surface excitations.
(a) Spectral peak of various $l$-fold star patterns, created after 1 s of
modulation, as a function of the modulation frequency. The spectral peak for
the hexagon ($l=6$ mode) is taken after $t=0.3$ s due to the involved dipole
instability (see main text). The filled circles designate the experimental
results and solid lines refer to the GPE predictions; notice the very good
agreement between the two. Each data point is averaged over 5-10 independent
experimental realizations, and the error bars denote the standard deviation of
the mean. (b) Floquet stability tongues for different modulation strengths
$\bar{a}_{m}/a_{bg}$ and frequencies $\omega_{m}/(2\pi)$ characterizing
distinct $l$-fold patterns.

The onset of resonance behavior of the surface excitations can be unveiled by
a Mathieu equation analysis for the amplitude of the density deformation, see
details in Ref. sup . Since the observed surface modes have no radial nodes,
we assume a density disturbance of the form $\delta
n=\zeta_{l}(t)r^{l}e^{il\phi}$. At short modulation times the deformation is
small, and we arrive at the Mathieu equation for $\zeta_{l}$ after linearizing
the hydrodynamic equations of the superfluid,

$\ddot{\zeta}_{l}(t)+\omega^{2}_{l}\big{[}1+\frac{\bar{a}_{m}}{a_{bg}}\cos(\omega_{m}t)\big{]}\zeta_{l}(t)=0,$
(2)

where $\omega_{l}=\sqrt{l}\omega_{r}$. This equation represents a
parametrically driven oscillator with a natural frequency $\omega_{l}$, having
a series of resonances at $\omega_{m}=2\omega_{l}/n$, where $n$ is an integer.

Within Floquet theory, a solution
$\zeta_{l}(t)=e^{(s+i\alpha\omega_{m})t}\sum_{k=-\infty}^{\infty}\zeta^{(k)}_{l}e^{ik\omega_{m}t}$
is sought, where $s$ is its growth rate and $\alpha$ is the Floquet exponent
Eckardt (2017). For $s>0$ the system is dynamically unstable and pattern
formation takes place at the surface. Setting $s=0$, we provide the marginal
stability boundaries of the $D_{l}$ symmetric patterns in Fig. 2(b). The
stability diagram is composed of a series of resonant tongues, where the
system exhibits star-shape patterns if $\bar{a}_{m}$ and $\omega_{m}$ reside
inside or at the boundaries of a specific tongue. Otherwise, the BEC cloud
solely performs a collective breathing motion. The spectrum in Fig. 2(a) can
be interpreted as the intersection of the instability tongues at modulation
strength $\Gamma=\bar{a}_{m}/a_{bg}\simeq 0.14$. Notice the not only
qualitative but also quantitative match of the instability tongues between
theoretical analysis, numerical findings and experimental results and the weak
deviation of the latter two when the nonlinear effects become more prominent
as discussed above. Including a dissipative term $\gamma\dot{\zeta}_{l}$ to
Eq. (2) lifts the tongues, suppressing pattern formation under a threshold
amplitude. The dissipation rate $\gamma=2\pi\times 1.8~{}\textrm{Hz}$, best
matching the experimentally measured threshold amplitudes, is used. The
temporal dynamics from the 3D GPE presents sub-harmonic oscillations
$\omega_{m}/2$ of the surface modes, leading to the Floquet exponent
$\alpha=1/2$.

Figure 3: (Color online) Dispersion relation of the surface modes for
superfluids in a harmonic trap. The spectral peak value for each $l$ mode is
obtained by fitting a Gaussian function to each resonance spectrum with small
driving amplitude. Solid line represents a guide to the eye for the dispersion
law of the surface mode in trapped superfluids. The inset shows the ratio of
the measured resonance frequency and the hydrodynamic predictions. The error
bars indicate a 95% confidence interval of the fit. Uncertainty of the radial
trap frequency is marked by the shaded region.

The theoretical investigations provide a deeper insight into the spontaneous
pattern formation. The natural frequency ($\omega_{l}=\sqrt{l}\omega_{r}$) of
the Mathieu equation is the dispersion of the surface excitation modes of
superfluids in a harmonic potential Stringari (1996). It indicates that the
observed star-shaped BECs subject to driving originate from the parametric
excitation of the surface mode with high multipolarity. In other words, one
can infer the dispersion laws by studying the resonance spectrum of each mode.
To measure the resonance frequencies of the surface excitations, we
investigate the surface mode spectrum with marginal mean modulation amplitude
$(\bar{a}_{m}=8-9a_{B})$ and sufficient condensate atom number
$(N=4.1(2)\times 10^{6})$. Figure 3 depicts the measured resonance frequencies
$\omega_{\text{res}}(l)$ of different modes up to $l=7$, which show a
remarkable agreement with the predicted square-root scaling dispersion. Within
the parameter regime that we operate, other effects from beyond mean-field
Pitaevskii and Stringari (1998), dipolar interactions Pollack  _et al._
(2009), and finite particle number Stringari (1996); Pitaevskii and Stringari
(1998); sup  are negligible. The observed small deviations ($\sim 2\%$) might
be attributed to the impact of the modulation on the resonance spectrum sup
and trap imperfections such as anharmonicity of the optical dipole trap Holten
_et al._ (2018).

Figure 4: (Color online) Time-evolution of the triangular surface mode. (a)
Dynamics of the growth rate of the triangular mode under a 104 Hz modulation.
The fluctuations of the spectral peak ($l=3$) increase exponentially in the
course of time. Each data point is a single experimental realization. Inset:
The averaged spectral peak $\langle{\mathcal{F}}_{3}\rangle$ over a single
oscillation time interval, $[t,t+9.5]$ ms. Dashed line designates an
exponential fit to the data (dark blue circle), and the solid line is obtained
from the GPE simulation. (b) Zooming in at 1 s reveals the oscillation of the
surface mode with 104 Hz driving frequency. Each data point corresponds to the
mean over four independent experimental runs, and the error bars indicate the
standard deviation of the mean. Inset: Absorption images during the first
oscillation period. The oscillation frequency is sub-harmonic, a result that
is confirmed by the GPE calculations, see also Ref. sup .

Another characteristic feature of the parametric excitations is the
exponential growth of the associated unstable modes as described by the
solution of Eq. (2). Focusing on the $l=3$ triangular mode, as a
representative example of this phenomenology, we investigate the time-
evolution of the spectral peak $\mathcal{F}_{3}(t)$ at resonance driving
$\omega_{m}=2\pi\times 104$ Hz. Initially, small amplitude fluctuations of the
condensate radius build up with no clear patterns, having
$\mathcal{F}_{3}\simeq 0$. After 300 ms of modulation, the azimuthal angular
symmetry of the condensate boundary breaks, rapidly forming a regular triangle
with sharp edges at $t=600$ ms. Averaging the spectral peak within one period
of oscillation $\langle\mathcal{F}_{3}\rangle$, we observe a clear
manifestation of the parametric instability via the exponential growth
dynamics [Fig. 4(a) inset], where the characteristic growth rate increases for
higher $l$ symmetry modes and driving amplitude sup .

As the evolution settles into a periodic pattern, the triangular surface mode
undergoes a regular oscillation characterized by the external driving
frequency [Fig. 4(b)]. The actual dynamics is sub-harmonic, a fact that is
also confirmed within the GPE calculations. The peaks oscillate 90 degree out-
of-phase with respect to the driving field, reflecting the dynamics of the
Mathieu equation under resonant frequency driving. When turning off the
periodic modulation after the development of the surface mode (not shown), the
latter experiences a relaxation towards a symmetric shape, satisfying
$\langle\mathcal{F}_{3}(t)\rangle\sim e^{-t/\tau_{3}}$, with a characteristic
time scale $\tau_{3}\approx 780$ ms. The decay time decreases with the mode
number such that $\tau_{7}\approx 130$ ms for the $l=7$ mode. It was also
found to decrease with higher temperature Fedichev  _et al._ (1998); Onofrio
_et al._ (2000).

Lastly, we would like to comment on the dynamics of the hexagon mode ($D_{6}$
pattern). Unlike the other surface excitations, this mode is found to be
unstable due to its emergent coupling with the dipole motion of the cloud.
This behavior is also found by the mean-field simulations. In the present
study, we focus on the dispersion law of the surface modes, such that the
resonance spectra in Fig. 2(a) and Fig. 3 are obtained for relatively short
evolution times ($t\approx 300$ ms) and in particular before the dipole motion
destabilizes the hexagonal pattern. Further details of the long-time dynamics
of $D_{6}$ patterns are provided in Ref. sup . This observation motivates
further efforts to unveil possible signatures of turbulent properties of the
surface modes.

Conclusions.– We observe experimentally and analyze theoretically the
generation of star-shaped surface patterns in a BEC due to the Faraday wave
instability induced by the periodic modulation of the scattering length.
Quantitative monitoring of the patterns enables to identify the growth rate of
the parametric instability and to measure the dispersion relation of the
surface excitations of superfluids in a harmonic trap, in very good agreement
with theoretical predictions and numerical computations. Since our
experimental method does not require special engineering to shape the
condensates, it can be applied to various quantum fluids in a broad context,
such as Fermi gases Altmeyer  _et al._ (2007); Holten  _et al._ (2018) and
dipolar quantum fluids Tanzi  _et al._ (2019), as well as exciton-polariton
BECs Whittaker  _et al._ (2017). Moreover, relevant ideas extend to binary
mixtures or quantum droplets, where the resonance spectrum can be utilized to
extract the interfacial tension of the superfluid boundary Maity  _et al._
(2020). By increasing the modulation strength, a transition to granulation and
turbulent behavior, of interest in its own right Lode  _et al._ (2021), can be
also studied in a two-dimensional condensate.

###### Acknowledgements.

We acknowledge discussions with Jongchul mun and thank to Junhyeok Hur and
Haejun Jung for critical reading of the manuscript. K. Kwon, S. Hur, and J.-y.
Choi is supported by the Samsung Science and Technology Foundation BA1702-06
and National Research Foundation of Korea (NRF) Grant under Projects No.
2020R1C1C1010863. K. Kim is supported by KAIST UP program. K. M. and D.K.M.
acknowledge MHRD, Govt. of India for the research fellowship. S. I. M.
gratefully acknowledges financial support in the framework of the Lenz-Ising
Award of the University of Hamburg. This material is based upon work supported
by the National Science Foundation under Grant No. PHY-1602994 and under Grant
No. DMS-1809074 (P.G.K.).

## References

  * Cross and Hohenberg (1993) M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
  * Maini  _et al._ (1997) P. Maini, K. Painter, and H. P. Chau, J. Chem. Soc. Faraday Trans. 93, 3601 (1997).
  * Petrov  _et al._ (1997) V. Petrov, Q. Ouyang, and H. L. Swinney, Nature (London) 338, 655 (1997).
  * Faraday (1831) M. Faraday, Phil. Trans. R. Soc. Lond. 121, 299 (1831).
  * Noblin  _et al._ (2005) X. Noblin, A. Buguin, and F. Brochard-Wyart, Phys. Rev. Lett. 94, 166102 (2005).
  * Shen  _et al._ (2010) C. L. Shen, W. J. Xie, and B. Wei, Phys. Rev. E 81, 046305 (2010).
  * Xia  _et al._ (2011) H. Xia, D. Byrne, G. Falkovich, and M. Shats, Nat. Phys. 7, 321 (2011).
  * von Kameke  _et al._ (2011) A. von Kameke, F. Huhn, G. Fernández-García, A. P. Muñuzuri, and V. Pérez-Muñuzuri, Phys. Rev. Lett. 107, 074502 (2011).
  * Francois  _et al._ (2014) N. Francois, H. Xia, H. Punzmann, S. Ramsden, and M. Shats, Phys. Rev. X 4, 021021 (2014).
  * Alarcón  _et al._ (2020) H. Alarcón, M. Herrera-Muñoz, N. Périnet, N. Mujica, P. Gutiérrez, and L. Gordillo, Phys. Rev. Lett. 125, 254505 (2020).
  * Engels  _et al._ (2007) P. Engels, C. Atherton, and M. A. Hoefer, Phys. Rev. Lett. 98, 095301 (2007).
  * Nguyen  _et al._ (2019) J. H. V. Nguyen, M. C. Tsatsos, D. Luo, A. U. J. Lode, G. D. Telles, V. S. Bagnato, and R. G. Hulet, Phys. Rev. X 9, 011052 (2019).
  * Zhang  _et al._ (2020) Z. Zhang, K.-X. Yao, L. Feng, J. Hu, and C. Chin, Nature Phys. 16, 652 (2020).
  * Tsubota  _et al._ (2013) M. Tsubota, M. Kobayashi, and H. Takeuchi, Phys. Rep. 522, 191 (2013).
  * Kevrekidis  _et al._ (2015) P. Kevrekidis, D. Frantzeskakis, and R. Carretero-González,  _The Defocusing Nonlinear Schrodinger Equation_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015).
  * Madeira  _et al._ (2020) L. Madeira, M. Caracanhas, F. dos Santos, and V. Bagnato, Annual Review of Condensed Matter Physics 11, 37 (2020).
  * Lode  _et al._ (2021) A. U. J. Lode, M. C. Tsatsos, P. G. Kevrekidis, G. D. Telles, D. Luo, R. G. Hulet, and V. S. Bagnato, arXiv:2103.07479 (2021).
  * Maity  _et al._ (2020) D. K. Maity, K. Mukherjee, S. I. Mistakidis, S. Das, P. G. Kevrekidis, S. Majumder, and P. Schmelcher, Phys. Rev. A 102, 033320 (2020).
  * Rath  _et al._ (2011) S. P. Rath, B. Spivak, and W. Zwerger, Phys. Rev. Lett. 107, 155703 (2011).
  * Seman  _et al._ (2011) J. Seman, E. Henn, R. Shiozaki, G. Roati, F. Poveda-Cuevas, K. Magalhães, V. Yukalov, M. Tsubota, M. Kobayashi, K. Kasamatsu, and V. Bagnato, Laser Phys. Lett. 8, 691 (2011).
  * Kim  _et al._ (2019) K. Kim, S. Huh, K. Kwon, and J.-y. Choi, Phys. Rev. A 99, 053604 (2019).
  * Huh  _et al._ (2020) S. Huh, K. Kim, K. Kwon, and J.-y. Choi, Phys. Rev. Research 2, 033471 (2020).
  * Zheng  _et al._ (2013a) H. Zheng, Y. Hao, and Q. Gu, J. Phys. B: At. Mol. and Opt. Phys. 46, 065301 (2013a).
  * Proukakis and Jackson (2008) N. P. Proukakis and B. Jackson, J. Phys. B: At. Mol. and Opt. Phys. 41, 203002 (2008).
  * Box and Muller (1958) G. E. P. Box and M. E. Muller, The Annals of Mathematical Statistics 29, 610 (1958).
  * Zheng  _et al._ (2013b) H. Zheng, Y. Hao, and Q. Gu, Journal of Physics B: Atomic, Molecular and Optical Physics 46, 065301 (2013b).
  * (27) See Supplemental Material at [URL] for the description of 1) Data analysis; 2) Dynamics of the patterns; 3) Details on the mean-field simulations; 4) Instability growth rate; 5) Mathieu equation and Floquet theory; 6) Finite-size effect; 7) Effect of the modulation amplitude; and 8) Unstable dynamics of the $D_{6}$ pattern.
  * Kovacic  _et al._ (2018) I. Kovacic, R. Rand, and S. Mohamed Sah, Applied Mechanics Reviews 70 (2018).
  * Eckardt (2017) A. Eckardt, Rev. Mod. Phys. 89, 011004 (2017).
  * Stringari (1996) S. Stringari, Phys. Rev. Lett. 77, 2360 (1996).
  * Pitaevskii and Stringari (1998) L. Pitaevskii and S. Stringari, Phys. Rev. Lett. 81, 4541 (1998).
  * Pollack  _et al._ (2009) S. E. Pollack, D. Dries, M. Junker, Y. P. Chen, T. A. Corcovilos, and R. G. Hulet, Phys. Rev. Lett. 102, 090402 (2009).
  * Holten  _et al._ (2018) M. Holten, L. Bayha, A. C. Klein, P. A. Murthy, P. M. Preiss, and S. Jochim, Phys. Rev. Lett. 121, 120401 (2018).
  * Fedichev  _et al._ (1998) P. O. Fedichev, G. V. Shlyapnikov, and J. T. M. Walraven, Phys. Rev. Lett. 80, 2269 (1998).
  * Onofrio  _et al._ (2000) R. Onofrio, D. S. Durfee, C. Raman, M. Köhl, C. E. Kuklewicz, and W. Ketterle, Phys. Rev. Lett. 84, 810 (2000).
  * Altmeyer  _et al._ (2007) A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R. Geursen, M. Bartenstein, C. Chin, J. H. Denschlag, and R. Grimm, Phys. Rev. Lett. 98, 040401 (2007).
  * Tanzi  _et al._ (2019) L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature (London) 574, 382 (2019).
  * Whittaker  _et al._ (2017) C. E. Whittaker, B. Dzurnak, O. A. Egorov, G. Buonaiuto, P. M. Walker, E. Cancellieri, D. M. Whittaker, E. Clarke, S. S. Gavrilov, M. S. Skolnick, and D. N. Krizhanovskii, Phys. Rev. X 7, 031033 (2017).
  * Crank and Nicolson (1947) J. Crank and P. Nicolson, Math. Proc. Camb. Philos. Soc 43, 50–67 (1947).
  * Antoine  _et al._ (2013) X. Antoine, W. Bao, and C. Besse, Comput. Phys. Commun 184, 2621 (2013).
  * Hill (2002) D. F. Hill, Physics of Fluids 14, 158 (2002).
  * Madelung (1927) E. Madelung, Z. Phys. 40, 322 (1927).
  * Ames (2008) W. F. Ames,  _Numerical methods for partial differential equations_, 2nd ed. (Boston Academic Press, 2008).
  * Dalfovo  _et al._ (1997) F. Dalfovo, S. Giorgini, M. Guilleumas, L. Pitaevskii, and S. Stringari, Phys. Rev. A 56, 3840 (1997).
  * Staliunas  _et al._ (2002) K. Staliunas, S. Longhi, and G. J. de Valcárcel, Phys. Rev. Lett. 89, 210406 (2002).
  * Hill and Eaves (2010) R. J. A. Hill and L. Eaves, Phys. Rev. E 81, 056312 (2010).
  * Goldman and Dalibard (2014) N. Goldman and J. Dalibard, Phys. Rev. X 4, 031027 (2014).
  * Barone  _et al._ (1977) S. R. Barone, M. A. Narcowich, and F. J. Narcowich, Phys. Rev. A 15, 1109 (1977).
  * Benjamin  _et al._ (1954) T. B. Benjamin, F. J. Ursell, and G. I. Taylor, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 225, 505 (1954).
  * Kumar (1996) K. Kumar, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 452, 1113 (1996).
  * Kevrekidis  _et al._ (2017) P. G. Kevrekidis, W. Wang, R. Carretero-González, and D. J. Frantzeskakis, Phys. Rev. Lett. 118, 244101 (2017).
  * Katsimiga  _et al._ (2020) G. C. Katsimiga, S. I. Mistakidis, T. M. Bersano, M. K. H. Ome, S. M. Mossman, K. Mukherjee, P. Schmelcher, P. Engels, and P. G. Kevrekidis, Phys. Rev. A 102, 023301 (2020).
  * Katsimiga  _et al._ (2021) G. C. Katsimiga, S. I. Mistakidis, P. Schmelcher, and P. G. Kevrekidis, New Journal of Physics 23, 013015 (2021).
  * Note (1) Note that we have performed a BdG analysis Katsimiga  _et al._ (2020) for a radial 2D symmetric trap confirming this argument.

Supplemental Material: “Spontaneous Formation of Star-Shaped Surface Patterns
in a Driven Bose-Einstein Condensate”

##  I Data analysis

Figure 5: (Color online) Characterizing surface modes with $D_{3}$ symmetry.
(a) in-situ absorption image for the $l=3$ surface mode. For a given angle
$\phi$ in the $xy$-plane, we obtain the surface boundary $R$ by bilinear fit.
(b) The boundary $R$ as a function of the angle $\phi$ for the triangular-
shaped condensate. (c) Angular spectrum of the surface boundary (note the
logarithmic scale).

The analysis process for characterizing the star-shaped surface patterns is
displayed in Fig. 5. Using a bilinear fit function, we obtain the radial
boundary of the condensate from its center-of-mass and estimate the radius $R$
as a function of the polar angle $\phi$ [Fig. 5(a)]. After taking a Fast
Fourier Transformation (FFT) to the condensate boundary, we normalize the
result to the number of data points, $N_{d}$, in the $R(\phi_{i})$, and obtain
the normalized Fourier spectrum of the angular mode,

$\mathcal{F}_{l}=\frac{1}{N_{d}}\absolutevalue{\sum_{i=0}^{N_{d}-1}R(\phi_{i})e^{-i2\pi
li/N_{d}}}.$ (3)

For the triangular symmetric surface pattern ($D_{3}$), the angular spectrum
displays a clear primary peak at $l=3$ and higher order peaks at integer
multiples ($l=3n$, $n$ is an integer) [Fig. 5(c)]. In the experiment, we use
the spectral peak $\mathcal{F}_{l}$ at the first peak position ($n=1$) of the
angular spectrum.

##  II Dynamics of the patterns

As explicated in the main text, the subharmonic nature of the star patterns
can be directly inferred from the Mathieu equation and the Floquet analysis.
Here, we also demonstrate this subharmonicity of the patterns within the full
3D Gross-Pitaevskii equation (GPE) by monitoring the density evolution of the
7Li BEC along the $xy$-plane, namely $n(x,y,t)=\int
dz\absolutevalue{\Psi(x,y,z,t)}^{2}$. Particularly, we consider a 7Li BEC
consisting of $N=4\times 10^{6}$ atoms which interact with scattering length
$a_{bg}=138.6a_{B}$ and are confined in the experimentally dictated 3D
harmonic trap. The system is initialized in its ground state configuration
[Fig. 6(a)]. The dynamics is induced upon considering a weak perturbation to
the ground state in order to emulate the experimentally observed thermal
fraction (see next section) and following a time-periodic driving of the
scattering length $a_{s}$ as explained in the main text.

We remark that the overall phenomenology, namely the pattern formation, is
exclusively caused by the periodic modulation of $a_{s}$ and it is not altered
in the absence of the thermal fraction. The latter accelerates the formation
of the patterns for a specific driving amplitude within a fixed evolution
time. Moreover, it is worth mentioning that due to its relatively strong
confinement the transversal direction plays only a minor role for the
discussed phenomena and thus we can safely conclude that the dynamics along
the $z$-direction is essentially “frozen”. Indeed, probing the density
$n(z,t)=\int dxdy\absolutevalue{\Psi(x,y,z,t)}^{2}$ we can deduce that it
undergoes a weak amplitude breathing motion (not shown for brevity). The above
facts have been confirmed both experimentally and theoretically within the 3D
GPE framework.

Figure 6: (Color online) (a)-(e) Snapshots of the density profiles $n(x,y,t)$
integrated over the $z$-direction at few select time-instants (see the
legends) in the course of an evolution simulation showcasing the generation of
the $D_{3}$ symmetric pattern. The 7Li BEC comprises of $N=4\times 10^{6}$
particles. The dynamics is triggered by applying a periodic modulation of the
scattering length $a_{s}$ characterized by a frequency $\omega_{m}=2\pi\times
104$ Hz and amplitude $\bar{a}_{m}=19a_{B}$. The colormap represents the
height of the bosonic density in units of $\mu m^{-3}$.

Below, we showcase the dynamical generation of the three-fold ($D_{3}$)
pattern realized for a driving frequency $\omega_{m}=2\pi\times 104$ Hz and
amplitude $\bar{a}_{m}=19a_{B}$, see Fig. 6(b)-(f). Closely inspecting the
time-evolution of the $D_{3}$ configuration we can clearly discern its
subharmonic nature. For instance, its lobes at $t=1067.5$ ms [Fig. 6(d)] and
at $t=1077$ ms [Fig. 6(e)] are rotated by an angle $\phi=\pi/6$ and $\phi=0$
respectively as compared to the one at $t=1058$ ms [Fig. 6(c)]. This dynamical
behavior and rotation of the patterns occurs for even longer timescales, see
for instance $n(x,y,t=1200~{}{\rm ms})$ depicted in Fig. 6(f). As a result, we
can deduce that the natural angular frequency of $D_{3}$ is $\omega_{l}\approx
2\pi\times 52$ Hz, which corresponds to half of the driving frequency, i.e.,
$\omega_{l}=\omega_{m}/2$ and thus confirms the subharmonicity of the $D_{3}$
structure. Finally, let us note that a similar dynamical response takes place
for all the star-shape patterns (not shown for brevity) realized at the
appropriate driving frequencies presented in Fig. 2 of the main text.

##  III Details of the 3D mean-field simulations

To simulate the periodically driven non-equilibrium dynamics and analyze the
spontaneous star-shape pattern formation of the quasi-2D 7Li BEC, we
numerically solve the full 3D time-dependent GPE of Eq. (1) in the main text
using a split-time Crank-Nicolson method Crank and Nicolson (1947); Antoine
_et al._ (2013). In particular, we reduce Eq. (1) of the main text into a
dimensionless form by performing a re-scaling of the spatial coordinates
$x^{\prime}\rightarrow x/a_{\rm osc}$, $y^{\prime}\rightarrow y/a_{\rm osc}$
and $z^{\prime}\rightarrow z/a_{\rm osc}$, the time
$t^{\prime}\rightarrow\omega_{x}t$ and the wave function
$\Psi^{\prime}(x^{\prime},y^{\prime},z^{\prime},t^{\prime})\rightarrow\sqrt{(N/a^{3}_{\rm
osc})\Psi^{(}x,y,z,t)}$. Here, $a_{\rm osc}=\sqrt{\hbar/(m\omega_{x})}$ is the
harmonic oscillator length along the $x$-direction. Recall that
$\omega_{x}=\omega_{y}\equiv\omega_{r}$. The ground state, $\Psi_{\rm
G}(x,y,z)$, of the BEC is obtained by numerically solving Eq. (1) in imaginary
time, $t=i\tau$. The normalization of the wave function is ensured at each
iteration of the imaginary time propagation by applying the transformation
$\Psi(x,y,z,t)\rightarrow\Psi(x,y,z,t)/\norm{\Psi(x,y,z,t)}$.

As discussed in the main text, experimentally a small amount of thermal atoms
(less than $10\%$) relative to the condensed atoms is present. To take into
account this thermal fraction we consider a weak perturbation modeled by a
noise term associated with a random perturbation on top of the ground state
$\Psi_{\rm G}(x,y,z)$ Zheng  _et al._ (2013a); Proukakis and Jackson (2008).
This normally distributed random perturbation emulates the experimentally
observed deformations to the condensate configuration due to thermal effects.
Therefore the initial wave function, $\Psi_{\rm in}(x,y,z)$, employed for the
propagation reads

$\displaystyle\Psi_{\rm in}(x,y,z)=\Psi_{\rm
G}(x,y,z)[1+\epsilon\delta(x,y,z)].$ (4)

In this expression, $\delta(x,y,z)$ represents a normally distributed random
perturbation with zero mean and variance unity which is generated by the so-
called Box-Mueller algorithm Box and Muller (1958). Also, $\epsilon\ll 1$ is
the amplitude of the perturbation which is taken to be $\epsilon=0.1$. Having
at hand the initial state of Eq. (4), we then propagate the GPE in real time,
under the influence of the scattering length modulation introduced in the main
text, up to $1500$ ms. To certify that our results are not affected by the
choice of the random distribution we utilize an average of the time-evolved
wave function over a sample of twenty different realizations for each specific
dynamical scenario. We have also checked that our results do not change upon
enlarging the sample of the noise term e.g. from twenty to forty. Our
numerical simulations are carried-out in a three-dimensional grid with hard-
wall boundary conditions. The hard-walls are located at $x_{\pm}=\pm 30a_{\rm
osc}$, $y_{\pm}=\pm 30a_{\rm osc}$, and $z_{\pm}=\pm 6a_{\rm osc}$. Their
location, of course, does not impact our results since for all of our
numerical simulations the BEC density resides within the spatial interval
$-16a_{\rm osc}<x,y<16a_{\rm osc}$ and $-0.8a_{\rm osc}<z<0.8a_{\rm osc}$.

##  IV Instability growth rate

To shed light into the growth rate of the Faraday-type instability (resulting
in the star-shape density patterns) and also its dependence on the driving
amplitude we monitor the spectral strength associated to the emerging
azimuthal defomation (see also main text) of the 7Li BEC. The dynamics of
$\mathcal{F}_{l}(t)$ is presented in Fig. 7(a), (b) regarding the $D_{3}$ and
$D_{7}$ patterns being realized for $\omega_{m}/(2\pi)=104$ Hz and $161$ Hz
respectively. It is evident that $\mathcal{F}_{l}(t)$ undergoes a vigorous and
rapidly fluctuating behavior, while its overall trend is reminiscent of a
typical instability growth rate Hill (2002).

Indeed, by employing the “moving average” of $\mathcal{F}_{l}(t)$ in terms of
time it is possible to smoothen its highly fluctuating behavior and
subsequently capture the growth rate of the involved instability. This “moving
average” being denoted in the following as $\mathcal{F}^{\rm avg}_{l}(t_{i})$
corresponds to the mean value of $\mathcal{F}_{l}(t)$ within the time-interval
$\delta t_{i}=2\pi/\omega_{m}$ containing $n$ time-instants. It reads

$\displaystyle\mathcal{F}^{\rm
avg}_{l}(t_{i})=\frac{\mathcal{F}_{l}(t_{i})+\mathcal{F}_{l}(t_{i-1})+\dots+\mathcal{F}_{l}(t_{i-(n-1)})}{n}.$
(5)

The respective $\mathcal{F}^{\rm avg}_{l}(t)$ is also provided in Fig. 7(a),
(b) for the $D_{3}$ and $D_{7}$ patterns. Apparently it features an
exponential growth until $t\approx 1000$ ms ($t\approx 600$ ms) for the
$D_{3}$ ($D_{7}$) configuration thus capturing the build-up of the Faraday
instability characterized by the azimuthal deformation of the BEC.
Consequently, $\mathcal{F}^{\rm avg}_{l}(t)$ exhibits a relatively low
amplitude oscillatory behavior around a certain value indicating that the
effect of the instability becomes less drastic. The underlying instability
growth rate can be determined via performing a fit of the $\mathcal{F}^{\rm
avg}_{l}(t)$ (until it acquires its maximum value) to the function
$Ae^{t/\tau}$, with $1/\tau$ quantifying the aforementioned instability growth
rate.

Figure 7: (Color online) Time-evolution of the spectral strength
$\mathcal{F}_{l}(t)$ (blue line) and its “moving average” over a single
driving period (orange line) for the (a) $D_{3}$ and (b) $D_{7}$ patterns. (c)
The growth rate $1/\tau$ and (d) the maximum spectral strength,
$\mathcal{F}^{\rm max}_{l}(t)$, as a function of the driving amplitude
$\bar{a}_{m}$ for different symmetric patterns (see legends). The individual
$D_{l}$ patterns are realized at $\omega_{m}=2\pi\times(84,104,119,147,158)$
Hz for $l=2,3,4,6,7$ while $\bar{a}_{m}=19a_{B}$ fixed. The 7Li BEC consists
of $N=4\times 10^{6}$ atoms. The dynamics is initiated by applying a periodic
modulation of the scattering length $a_{s}$ (see main text).

The instability growth rate $1/\tau$ for different $D_{l}$ patterns with
respect to $\bar{a}_{m}$ is illustrated in Fig. 7(c). We observe that $1/\tau$
increases linearly for larger $\bar{a}_{m}$ and also for higher-fold symmetric
patterns. This behavior clearly suggests that lower-fold configurations
require longer time to appear and therefore their associated instability is
weaker. To expose the impact of the driving amplitude $\bar{a}_{m}$ on the
development of the Faraday instability we show in Fig. 7(d) the maximum value
$\mathcal{F}^{\rm max}_{l}$ of the spectral strength $\mathcal{F}_{l}(t)$ for
different $\bar{a}_{m}$ and also specific values of $\omega_{m}$ corresponding
to a particular $l$-fold pattern ranging from $D_{2}$ to $D_{7}$. We remark
that $\mathcal{F}^{\rm max}_{l}$ for fixed $\omega_{m}$ (or equivalently
$l$-fold pattern) essentially provides a measure of the instability strength;
if $\mathcal{F}^{\rm max}_{l}\neq 0$ implies that the instability develops
otherwise $\bar{a}_{m}$ is not able to seed it (at least in the experimentally
relevant time window). Indeed, depending on the $l$-fold symmetric pattern we
find that a specific threshold driving amplitude is required in order to
initiate the underlying pattern formation. This threshold driving amplitude is
larger for lower-fold symmetric patterns and it lies very close to the one
predicted by the Floquet analysis (see the resonant tongues depicted in Fig.
2(b) of the main text) and also identified experimentally. Furthermore, for
lower-fold symmetric patterns the global maximum of $\mathcal{F}_{l}$, i.e.
$\mathcal{F}^{\rm max}_{l}\neq 0$, is larger [Fig. 7(c)], see also the
stability diagram displayed in Fig. 2(a) of the main text. This means that the
overall azimuthal deformation is more dramatic for lower $D_{l}$
configurations.

##  V Mathieu equation and Floquet theory

In the following, we elaborate on the derivation of the Mathieu equation [Eq.
(2) in the main text] which provides an effective single-particle description
of the deformation amplitude at the surface of the BEC. Then, we discuss the
application of the Floquet analysis to the Mathieu equation used to determine
the emergent instability tongues. Since the dynamical pattern formation takes
place in the weakly confined $xy$-plane and it is decoupled from the
$z$-direction, we are able to factorize the 3D wave function as
$\Psi(x,y,z,t)=\psi(x,y,t)\tilde{\phi}(z)$. Indeed, as it has been already
argued the dynamics along the transverse $z$-direction, besides undergoing a
very weak amplitude breathing motion, it essentially remains “frozen”. This
allows us to consider a Gaussian (ground-state) ansatz for $\tilde{\phi}(z)$,
namely
$\tilde{\phi}(z)=(m\omega_{z}/(\pi\hbar))^{1/4}e^{-m\omega_{z}z^{2}/(2\hbar)}$.
Substituting these ansatze into the 3D GPE [Eq. (1) in the main text] and
performing the integration over $z$, we can straightforwardly obtain the
following 2D GPE subjected to a periodically driven scattering length. Namely

$\displaystyle\begin{split}i\hbar\frac{\partial\psi}{\partial
t}=&-\frac{\hbar^{2}}{2m}(\partial_{x}^{2}+\partial_{y}^{2})\psi+\frac{1}{2}m\omega_{r}^{2}(x^{2}+y^{2})\psi\\\
&+g_{2D}\left(1+\frac{\bar{a}_{m}}{a_{bg}}\cos(\omega_{m}t)\right)\absolutevalue{\psi}^{2}\psi,\end{split}$
(6)

where $g_{2D}=g\sqrt{m\omega_{z}/(2\pi\hbar)}$ is the effective 2D interaction
strength and $g=4\pi\hbar^{2}a_{bg}/m$.

Next, we define the so-called Madelung transformation Madelung (1927)
$\psi=ne^{i\theta}$, where $n(x,y,t)$ and $\theta(x,y,t)$ are the 2D density
and phase (or velocity potential) of the condensate respectively and also
employ the superfluid velocity $\mathbf{v}=\frac{\hbar}{m}\nabla\theta$ Ames
(2008). Inserting these expressions into Eq. (6), it is possible to arrive at
the hydrodynamic form of the 2D GPE described by the following coupled set of
equations

$\frac{\partial n}{\partial
t}=-\frac{\hbar}{m}\nabla\cdot\left(n\nabla\theta\right),$ (7)

and

$\begin{split}m\frac{\partial\mathbf{v}}{\partial
t}&=-\nabla\bigg{[}\frac{1}{2}mv^{2}+\frac{1}{2}m\omega_{r}^{2}(x^{2}+y^{2})\\\
&+ng_{2D}\left(1+\frac{\bar{a}_{m}}{a_{bg}}\cos(\omega_{m}t)\right)-\frac{\hbar^{2}}{2m\sqrt{n}}\nabla^{2}\sqrt{n}\bigg{]}.\end{split}$
(8)

The wave function of the initial (ground) state is real while the phase is
spatially uniform and therefore at $t=0$ it holds that
$\mathbf{v}=\nabla\theta=0$. In particular, $\theta$ and $n(x,y,t)$ retain
their equilibrium values if the system is let to evolve in time in the absence
of any external perturbation. However, under the influence of the time
periodic modulation the equilibrium density is weakly modified as
$n\rightarrow n+\delta n$ while $\mathbf{v}\neq 0$ since $\theta$ becomes
spatially dependent. Note also that $\delta n\ll n$ represents a weak
amplitude distortion of the BEC surface, an assumption that is valid for small
modulation amplitudes ($\bar{a}_{m}\ll a_{bg}$) and short timescales.
Linearizing Eqs. (7)-(8) in terms of $\delta n\ll n$ and $\mathbf{v}\ll 1$
they can be cast into the following compact form

$m\frac{\partial^{2}\delta n}{\partial
t^{2}}=g_{2D}\left(1+\frac{\bar{a}_{m}}{a_{bg}}\cos(\omega_{m}t)\right)\nabla\cdot(n\nabla\delta
n).$ (9)

This describes the time-evolution of the density deformation due to the time
periodic perturbation. Since the pattern formation occurs in a circular
geometry it can be more conveniently addressed by adopting polar-coordinates,
i.e., $r=\sqrt{x^{2}+y^{2}}$ and $\phi=\tan^{-1}(y/x)$. In this sense, we
assume that $\delta n(t)=\zeta_{l}(t)r^{l}e^{il\phi}$, with $\zeta_{l}(t)$
being the deformation amplitude of the $l$-fold pattern. Furthermore,
operating at the large particle number limit (this is the case also for the
present experiment) it is reasonable to utilize the Thomas-Fermi (TF)
approximation. Within the latter, the equilibrium density of the BEC is
$n=m\omega_{r}^{2}(R_{TF}^{2}-r^{2})/(2g_{2D})$, with $R_{TF}$ denoting the TF
radius of the cloud Dalfovo  _et al._ (1997). Using these expressions for
$\delta n$ and $n$ into Eq. (9) we obtain the following Mathieu equation for
the amplitude of the density deformation

$\displaystyle\ddot{\zeta}_{l}(t)+\omega^{2}_{l}\big{[}1+\frac{\bar{a}_{m}}{a_{bg}}\cos(\omega_{m}t)\big{]}\zeta_{l}(t)=0.$
(10)

Here, $\omega_{l}=\sqrt{l}\omega_{r}$ is the natural angular frequency of the
$l$-mode pattern. Importantly, this Mathieu equation is a reduction of the
higher dimensional nonlinear GPE by using linear stability analysis into an
equation involving a single degree-of-freedom Staliunas  _et al._ (2002). It
addresses the dynamics of the amplitude of the deformation at the BEC surface
due to the periodic modulation of the scattering length Maity  _et al._ (2020)
and the subsequent formation of standing waves characterizing the Faraday-like
instability Shen  _et al._ (2010); Hill and Eaves (2010); Faraday (1831).

Figure 8: (Color online) (a) The peak of the spectral strength as identified
by the FFT amplitude of $R(\phi)$ [see also Eq. (3)] for varying driving
frequency and two distinct particle numbers (see legend). (b) Breathing mode
frequency $\omega_{\rm br}/\omega_{r}$ as a function of the particle number
$N$. The percentage of the $\omega_{\rm br}$ deviation from the hydrodynamic
limit (here $N=10^{7}$) is indicated by $\mathcal{Z}$. The inset illustrates
the time-evolution of $\sigma_{R}(t)$ and the peak spectral strength
manifesting the exponential growth due to the instability (prior to an
eventual saturation) for $\omega_{m}=2\pi\times 119$ Hz and $N=4\times
10^{6}$.

Note that according to Eq. (10) the pattern formation occurs for any
$\bar{a}_{m}>0$. However, as we have experimentally observed and confirmed
within the GPE simulations [Fig. 7(d)] there exists a finite threshold value
of $\bar{a}_{m}$ which is required to generate the patterns. This threshold
amplitude is larger for lower symmetric patterns, see Fig. 7(d). It depends on
the considered total evolution time and the thermal fraction $\epsilon$.
Notably, increasing one of these two parameters enhances the possibility of
observing a specific pattern, thereby reducing the required threshold driving
amplitude. Accordingly, it is desirable to construct an extended form of Eq.
(10) whose solution captures the pattern specific threshold driving amplitude.
For this reason we include a phenomenological dissipative term
$\gamma\dot{\zeta}_{l}$ in Eq. (10). Naturally, the value of $\gamma$
decreases for a smaller threshold driving amplitude. As we have demonstrated
in the main text [Fig. 2(b)], owing to this dissipative term the minima of the
resonant tongues representing the pattern specific threshold driving
amplitudes are shifted to larger values when compared to the $\gamma=0$ case.
In fact, one can deduce the value of $\gamma=2\pi\times 1.8$ Hz by comparing
these shifts to those obtained from the experiment. According to Floquet
theory, Eq. (10) admits a solution of the form

$\zeta_{l}(t)=e^{(s+i\alpha\omega_{m})t}\sum_{k=-\infty}^{\infty}\zeta^{(k)}_{l}e^{ik\omega_{m}t},$
(11)

where $s$ is the growth rate of the deformation amplitude and $\alpha$ is the
Floquet exponent Goldman and Dalibard (2014); Barone  _et al._ (1977); Eckardt
(2017). Inserting Eq. (11) into Eq. (10) we arrive at the eigenvalue problem

$\displaystyle\begin{split}\bigg{(}\frac{2(k+\alpha)^{2}\omega^{2}_{m}}{\omega^{2}_{l}}-&\frac{2i\gamma(k+\alpha)\omega_{m}}{\omega^{2}_{l}}-2\bigg{)}\zeta^{(k)}_{l}\\\
&=\frac{\bar{a}_{m}}{a_{bg}}\big{(}\zeta^{(k+1)}_{l}+\zeta^{(k-1)}_{l}\big{)},\end{split}$
(12)

whose solutions determine the Floquet stability tongues in the parameter space
($\bar{a}_{m}$,$\omega_{m}$) as depicted in Fig. 2(b) of the main text. These
are also known as Arnold tongues Benjamin  _et al._ (1954); Kumar (1996) in
nonlinear dynamical systems characterizing their stability. Namely, for a
given parameter space ($\mathcal{X}$, $\mathcal{Y}$) when residing outside a
specific tongue the system is stable otherwise it is unstable. Interestingly,
the term $2i\gamma(k+\alpha)\omega_{m}/\omega^{2}_{l}$ in Eq. (12) determines
the corresponding threshold amplitude of the $l$-fold symmetric pattern. It is
apparent that for a fixed value of $\gamma$ this term decreases with
increasing $l$ thus confirming our experimental observation that a smaller
threshold driving amplitude is required for the generation of higher-fold
symmetric patterns.

##  VI Finite-size effect

The star-shape patterns appear due to the activation of the surface modes
which are collective modes of the equilibrium state Dalfovo  _et al._ (1997);
Ames (2008); Stringari (1996). Recall that in the TF limit the spectrum of the
collective excitations of the ground state is $\omega_{\rm
exc}=\omega_{r}[l+2n_{r}(1+l)+2n_{r}^{2}]^{1/2}$ Stringari (1996), where
$n_{r}$, $l$ are the radial and azimuthal wavenumbers respectively. For
$n_{r}=0$ this spectrum captures only the surface modes. This prediction turns
out to be very accurate for low-lying modes, but a deviation from this limit
can be noticeable regarding higher $l$-modes and small atom numbers Stringari
(1996); Dalfovo  _et al._ (1997). That is, the nature of excitations undergo a
crossover from collective to single-particle character. In this section, we
present the numerical and experimental studies of the finite-size effect on
the surface mode spectrum.

Figure 9: (Color online) Finite-size effect on the $D_{7}$ symmetric pattern.
(a) The resonance spectrum is measured after 1.2 s of modulation with
$\bar{a}_{m}=16a_{B}$ for two different atom numbers (see legend). Solid lines
represent the Gaussian fit curve to the data. (b) Peak frequency for the
heptagon ($l=7$) at distinct condensate atom numbers $N$. The shaded area
depicts the GPE predictions (right axis) including systematic uncertainties
from the experimental determination of the radial trap frequency
($\omega_{r}=2\pi\times 29.4(2)$ Hz).

Fig. 8(a) displays the resonance spectrum of $l$-fold symmetric patterns for
two different atom numbers, $N=4\times 10^{6}$ and $N=3\times 10^{4}$. It
becomes apparent that the resonant driving frequencies of each $l$-fold
pattern increase with decreasing particle number $N$. Interestingly, the
aforementioned spectral shift is more pronounced for higher-fold symmetric
patterns which can therefore serve as better indicators for probing the
crossover from finite particle number to the TF regime. The shift of
eigenfrequencies can be demonstrated by determining the corresponding
Bogoliubov-de-Gennes linearization spectrum of the GPE Kevrekidis  _et al._
(2017); Katsimiga  _et al._ (2020, 2021). However, for a 3D setup (as
considered herein) this is a computationally challenging task 111Note that we
have performed a BdG analysis Katsimiga  _et al._ (2020) for a radial 2D
symmetric trap confirming this argument..

Note that, irrespectively of the particle number, at the early stages of the
driven dynamics the BEC undergoes a radial expansion and contraction i.e. a
breathing mode. Subsequently, in the long-time dynamics, the Faraday
instability induces the spatial symmetry breaking phenomenon Maity  _et al._
(2020). The latter is related to the activation of the surface modes manifest
in the form of star-shape patterns. Thus the behavior of the breathing mode is
inherently related to the general inset illustration of the star patterns,
namely the radial oscillation gradually transforms into an azimuthal
deformation capturing the onset of the Faraday instability, see inset of Fig.
8(b). The aforementioned mechanism motivates the investigation of the
breathing mode behavior with respect to the particle number in order to
understand the impact of $N$ on the natural angular frequencies of the star-
shape patterns.

To monitor the breathing motion of the cloud we resort to the second moment of
the center-of-mass coordinate which allows us to determine the instantaneous
spatial extent of the BEC

$\displaystyle\sigma^{2}_{R}(t)=\int
dxdydz(x^{2}+y^{2}+z^{2})\absolutevalue{\Psi(x,y,z,t)}^{2}.$ (13)

The time-evolution of this quantity for $N=4\times 10^{6}$ and
$\omega_{m}=2\pi\times 104$ Hz is presented in the inset of Fig. 8(b). As it
can be seen, it features an oscillatory behavior of constant amplitude until
$t\approx 600$ ms. Afterwards its fluctuation amplitude decreases due to the
appearance of the star-shape patterns and finally saturates around $t=1000$ ms
where the patterns become stabilized. The underlying breathing frequency
$\omega_{\rm br}$ corresponds to the one with larger Fourier amplitude of
$\sigma^{2}_{R}(t)$ which is shown in Fig. 8(b) for different particle
numbers. Evidently, $\omega_{\rm br}$ decreases with increasing $N$ and in
particular exhibits a saturation tendency after the TF region is reached, see
Fig. 8(b). Thus, we can conclude that the shift of the natural angular
frequencies of the patterns stems from the increase of the underlying
collective excitation frequencies for smaller particle number. This
observation is consistent with the analytical predictions of previous studies
Pitaevskii and Stringari (1998); Stringari (1996). It also yields that the
finite-size effects are suppressed in our system since $N>(a_{\rm
osc}/a_{bg})^{2}\log(R_{TF}/a_{\rm osc})$ holds when $N>10^{6}$, a result that
is also confirmed by our GPE calculations [Fig. 8(b)]. To explicitly expose
deviations of the breathing frequency, $\omega_{\rm br}$, for a specific $N$
from the hydrodynamic limit ($N\rightarrow\infty)$ prediction, we employ
$\mathcal{Z}=[\omega_{\rm br}(N)-\omega_{\rm
br}(N\rightarrow\infty)/\omega_{\rm br}(N\rightarrow\infty)]\%$. As a
representative example of the hydrodynamic limit we take $N=10^{7}$, and
inspect $\mathcal{Z}$ for varying $N$. As can be seen, the deviation from the
hydrodynamic limit for $N>10^{6}$ is less than $0.2\%$, thus further showing
that finite-size effects are vanishingly small.

We also experimentally study the finite-size effect on the spectral peak shift
of the $l=7$ mode for different atom number in the condensate [Fig. 9]. In the
experiment, we set the mean modulation amplitude $\bar{a}_{m}=16a_{B}$ to
observe a $D_{7}$ symmetry surface mode for low particle numbers. A gradual
blue-shift is observed when we have less than $2\times 10^{6}$ atoms in the
condensate. Such a trend is also captured by the GPE calculations. The
absolute resonance frequencies in the experiment exhibit an offset from the
GPE predictions simply due to the blue-shift arising for strong modulations.
For our experimental parameters, the mode number for the crossover from the
finite to the TF atom limit Dalfovo  _et al._ (1997) is
$l_{c}=2^{1/3}(R_{TF}/a_{\rm{osc}})^{4/3}\simeq 40$ for $N=10^{6}$.
Apparently, it is much larger than the $l=7$ of the theory prediction. As such
the observed small shift in the resonance frequencies is reasonable. According
to these investigations, we conclude that finite-size effects in the
dispersion spectrum of the surface excitations [Fig. 3 in the main text] are
very small and in particular less than $0.4\%$ for $N\geq 2\times 10^{6}$
atoms.

Figure 10: (Color online) (a) Resonance spectrum of the square mode ($l=4$)
with increasing modulation amplitude (from light blue to dark blue). The
vertical offsets have been added to the modulation spectrum for clarity. (b)
The lower and upper bound frequencies of $D_{4}$ for varying mean modulation
amplitude $\bar{a}_{m}$. Inside each specific frequency interval the $D_{4}$
surface mode is well developed (unstable); otherwise the star pattern
formation is suppressed (stable region). The dashed lines are linear fits to
the boundary frequencies.  Figure 11: (Color online) Time-evolution of the
surface mode with $D_{6}$ symmetry. (a) The spectral peak of the $l=6$ mode
gradually disappears in the course of the modulation time. (b) Simultaneously,
the amplitude of the dipole motion increases after $t\sim 400$ ms. The solid
lines correspond to sinusoidal fits. Insets show snapshot images of the
$D_{6}$ mode at (i) $t=398$ and (ii) $t=597$ ms, respectively.

Finally, let us comment on the emergence of beyond mean-field corrections
manifesting in the breathing frequency of the BEC cloud. As we argue below
these effects, that we theoretically determine herein, are adequately small
and they pose a challenge to experimental detection. Recall that the breathing
mode frequency in the hydrodynamic limit is given by $\omega^{2}_{\rm
br}=\omega_{r}^{2}[2+(3/2)\lambda^{2}-(1/2)\sqrt{(9\lambda^{4}-16\lambda^{2}+16)}]$
with $\lambda=\omega_{z}/\omega_{r}$. For $\omega_{r}=2\pi\times 29.4$ Hz,
$\lambda=725/29.4$, we thus have $\omega_{\rm br}=2\pi\times 53.67$ Hz. On the
other hand, the beyond-mean-field correction Pitaevskii and Stringari (1998)
reads $\delta\omega_{\rm br}/\omega_{\rm
br}=[63\sqrt{\pi}/128]\sqrt{a_{bg}^{3}n(0)}$, where
$a_{bg}^{3}n(0)=[15^{2/5}/(8\pi)](N^{1/6}a_{bg}/a_{\rm osc})^{12/5}$. In this
sense, for our setup with $N=4\times 10^{6}$, $\omega_{r}=2\pi\times 29.4$ Hz,
$a_{bg}=138.6a_{B}$, and $a_{\rm{osc}}=6.96~{}\mu m$, the anticipated relative
deviation of the breathing frequency due to the presence of quantum
fluctuations is $\delta\omega_{br}\approx 2\pi\times 0.092$ Hz. Evidently,
from our simulations we have $\omega_{\rm br}=2\pi\times 53.80$ Hz, see Fig.
8(b). Therefore, the small deviation $\sim 0.09$ Hz of the breathing mode
frequency observed in our simulations from the hydrodynamic limit is
attributed to the presence of beyond mean-field contributions.

##  VII Effect of the modulation amplitude

The resonance spectrum of the surface modes shows a strong dependence on the
mean modulation amplitude, $\bar{a}_{m}$. As an example of this behavior,
Figure 10(a) displays the surface mode spectrum with $l=4$ at different
$\bar{a}_{m}$. With higher modulation amplitude, the spectrum is broadened, as
explicated in the instability diagram of the main text [Fig. 2(b)], and its
shape is deformed to highly asymmetric with respect to the peak modulation
frequency. To identify the resonance frequency $\omega_{\text{res}}$ of the
surface mode, we reduce the mean modulation amplitude so that the resonance
spectrum becomes narrower and symmetric. For example, the resonance frequency
of the $D_{4}$ symmetry mode is extracted from the Gaussian fit to the
spectrum with $\bar{a}_{m}=9a_{B}$, which gives
$\omega_{\text{res}}=2\pi\times 119.17(5)$ Hz. Alternatively, we could
estimate the resonance frequency from the experimentally obtained stability
diagram [Fig. 10(b)]. Displaying the upper and lower bound frequencies of the
surface mode, we are able to estimate the resonance frequency ($2\pi\times
118.64(5)$ Hz) at vanishingly small modulation amplitude. The resonance
frequency obtained from the stability diagram shows a better agreement with
the hydrodynamic prediction ($l=2,3$ modes exhibit a similar behavior). This
might indicate that the finite $\bar{a}_{m}$ still affects the resonance
spectrum, which could lead to the minor deviations of the measured resonance
frequency with respect to the hydrodynamic theory prediction.

##  VIII  Unstable dynamics of the $D_{6}$ pattern

In contrast to the other surface modes, the hexagonal-shaped BEC is found to
be dynamically unstable after long-time modulation. Fig. 11(a) displays the
time-evolution of the spectral peak associated with $D_{6}$ symmetry,
$\mathcal{F}_{6}$, at a driving frequency of 146 Hz. The star-shaped
oscillation gradually disappears, and simultaneously the shape of the atomic
cloud is distorted, such that the center-of-mass of the condensate is
displaced from the origin. We characterize this displacement by utilizing the
third moment of the density distribution $\mu_{3}=\int r^{3}n~{}d^{2}r$. The
dynamics of this quantity is presented in Fig. 11(b).

At the early stages of the pattern formation ($t=300$ ms), the $l=6$ surface
mode possesses a well-defined six-fold rotational symmetry at its center-of-
mass. However, after $t=400$ ms, the dipole motion of the cloud is enhanced
and thus acts against the involved star-shaped oscillation. Eventually, only
the dipole motion is left without noticeable surface mode after 600 ms of
modulation. This manifests the emergence of a mode coupling between the
ensuing $l=6$ surface mode and the dipole mode. Similar results are also
demonstrated in the numerical simulations within the 3D GPE. This implies that
the mode-coupling phenomenon might not be solely attributed to experimental
imperfections. Although a deeper understanding of its physical origin is still
lacking, we believe that it deserves further investigation to unveil possible
signatures of turbulent properties of the surface modes, which are of interest
in their own right, e.g., as pathways towards turbulence generation.