• We analyze the ground state of a Bose-Einstein condensate in the presence of higher-order interaction (HOI), modeled by a modified Gross-Pitaevskii equation (MGPE). In fact, due to the appearance of HOI, the ground state structures become very rich and complicated. We establish the existence and non-existence results under different parameter regimes, and obtain their limiting behaviors and/or structures with different combinations of HOI and contact interactions. Both the whole space case and the bounded domain case are considered, where different structures of ground states are identified.
  • We propose a new fourth-order compact time-splitting ($S_\text{4c}$) Fourier pseudospectral method for the Dirac equation by splitting the Dirac equation into two parts together with using the double commutator between them to integrate the Dirac equation at each time interval. The method is explicit, fourth-order in time and spectral order in space. It is unconditional stable and conserves the total density in the discretized level. It is called a compact time-splitting method since, at each time step, the number of sub-steps in $S_\text{4c}$ is much less than those of the standard fourth-order splitting method and the fourth-order partitioned Runge-Kutta splitting method. Comparison among $S_\text{4c}$ and many other existing time-splitting methods for the Dirac equation are carried out in terms of accuracy and efficiency as well as long time behavior. Numerical results demonstrate the advantage in terms of efficiency and accuracy of the proposed $S_\text{4c}$. Finally we report the spatial/temporal resolutions of $S_\text{4c}$ for the Dirac equation in different parameter regimes including the nonrelativistic limit regime, the semiclassical limit regime, and the simultaneously nonrelativisic and massless limit regime.
  • In this paper, we systematically review mathematical models, theories and numerical methods for ground states and dynamics of spinor Bose-Einstein condensates (BECs) based on the coupled Gross-Pitaevskii equations (GPEs). We start with a pseudo spin-1/2 BEC system with/without an internal atomic Josephson junction and spin-orbit coupling including (i) existence and uniqueness as well as non-existence of ground states under different parameter regimes, (ii) ground state structures under different limiting parameter regimes, (iii) dynamical properties, and (iv) efficient and accurate numerical methods for computing ground states and dynamics. Then we extend these results to spin-1 BEC and spin-2 BEC. Finally, extensions to dipolar spinor systems and/or general spin-F (F>=3) BEC are discussed.
  • We present a regularized finite difference method for the logarithmic Schr\"odinger equation (LogSE) and establish its error bound. Due to the blow-up of the logarithmic nonlinearity, i.e. $\ln \rho\to -\infty$ when $\rho\rightarrow 0^+$ with $\rho=|u|^2$ being the density and $u$ being the complex-valued wave function or order parameter, there are significant difficulties in designing numerical methods and establishing their error bounds for the LogSE. In order to suppress the round-off error and to avoid blow-up, a regularized logarithmic Schr\"odinger equation (RLogSE) is proposed with a small regularization parameter $0<\varepsilon\ll 1$ and linear convergence is established between the solutions of RLogSE and LogSE in term of $\varepsilon$. Then a semi-implicit finite difference method is presented for discretizing the RLogSE and error estimates are established in terms of the mesh size $h$ and time step $\tau$ as well as the small regularization parameter $\varepsilon$. Finally numerical results are reported to confirm our error bounds.
  • Based on the two-dimensional mean-field equations for pancake-shaped dipolar Bose-Einstein condensates in a rotating frame with both attractive and repulsive dipole-dipole interaction (DDI) as well as arbitrary polarization angle, we study the profiles of the single vortex state and show how the critical rotational frequency change with the s-wave contact interaction strengths, DDI strengths and the polarization angles. In addition, we find numerically that at the `magic angle' $\vartheta=\arccos(\sqrt{3}/3)$, the critical rotational frequency is almost independent of the DDI strength. By numerically solving the dipolar GPE at high rotational speed, we identify different patterns of vortex lattices which strongly depend on the polarization direction. As a result, we undergo a study of vortex lattice structures for the whole regime of polarization direction and find evidence that the vortex lattice orientation tends to be aligned with the direction of the dipoles.
  • We study asymptotically and numerically the fundamental gap -- the difference between the first two smallest (and distinct) eigenvalues -- of the fractional Schr\"{o}dinger operator (FSO) and formulate a gap conjecture on the fundamental gap of the FSO. We begin with an introduction of the FSO on bounded domains with homogeneous Dirichlet boundary conditions, while the fractional Laplacian operator defined either via the local fractional Laplacian (i.e. via the eigenfunctions decomposition of the Laplacian operator) or via the classical fractional Laplacian (i.e. zero extension of the eigenfunctions outside the bounded domains and then via the Fourier transform). For the FSO on bounded domains with either the local fractional Laplacian or the classical fractional Laplacian, we obtain the fundamental gap of the FSO analytically on simple geometry without potential and numerically on complicated geometries and/or with different convex potentials. Based on the asymptotic and extensive numerical results, a gap conjecture on the fundamental gap of the FSO is formulated. Surprisingly, for two and higher dimensions, the lower bound of the fundamental gap depends not only on the diameter of the domain, but also the diameter of the largest inscribed ball of the domain, which is completely different from the case of the Schr\"{o}dinger operator. Extensions of these results for the FSO in the whole space and on bounded domains with periodic boundary conditions are presented.
  • We establish uniform error bounds of a finite difference method for the Klein-Gordon-Zakharov system (KGZ) with a dimensionless parameter $\varepsilon \in (0,1]$, which is inversely proportional to the acoustic speed. In the subsonic limit regime, i.e. $0<\varepsilon \ll 1$, the solution propagates highly oscillatory waves in time and/or rapid outgoing initial layers in space due to the singular perturbation in the Zakharov equation and/or the incompatibility of the initial data. Specifically, the solution propagates waves with $O(\varepsilon)$-wavelength in time and $O(1)$-wavelength in space as well as outgoing initial layers in space at speed $O(1/\varepsilon)$. This high oscillation in time and rapid outgoing waves in space of the solution cause significant burdens in designing numerical methods and establishing error estimates for KGZ. By adapting an asymptotic consistent formulation, we propose a uniformly accurate finite difference method and rigorously establish two independent error bounds at $O(h^2+\tau^2/\varepsilon)$ and $O(h^2+\tau+\varepsilon)$ with $h$ mesh size and $\tau$ time step. Thus we obtain a uniform error bound at $O(h^2+\tau)$ for $0<\varepsilon\le 1$. The main techniques in the analysis include the energy method, cut-off of the nonlinearity to bound the numerical solution, the integral approximation of the oscillatory term, and $\varepsilon$-dependent error bounds between the solutions of KGZ and its limiting model when $\varepsilon\to0^+$. Finally, numerical results are reported to confirm our error bounds.
  • We present a new approach for predicting stable equilibrium shapes of crystalline islands on flat substrates, as commonly occur through solid-state dewetting of thin films. The new theory is a generalization of the widely used Winterbottom construction i.e., an extension of the Wulff construction for particles on substrates). This approach is equally applicable to cases where the crystal surface energy is isotropic, weakly anisotropic, strongly anisotropic and "cusped". We demonstrate that, unlike in the classical Winterbottom approach, multiple equilibrium island shapes may be possible when the surface energy is strongly anisotropic. We analyze these shapes through perturbation analysis, by calculating the first and second variations of the total free energy functional with respect to contact locations and island shape. Based on this analysis, we find the necessary conditions for stable equilibria and exploit this through a generalization of the Winterbottom construction to identify all possible stable equilibrium shapes. Finally, we propose a dynamical evolution method based on surface diffusion mass transport to determine whether all of the stable equilibrium shapes are dynamically accessible. Applying this approach, we demonstrate that islands with different initial shapes may evolve into different stationary shapes and show that these dynamically-determined stationary states correspond to the predicted stable equilibrium shapes, as obtained from the generalized Winterbottom construction.
  • We study asymptotically and numerically the fundamental gaps (i.e. the difference between the first excited state and the ground state) in energy and chemical potential of the Gross-Pitaevskii equation (GPE) -- nonlinear Schrodinger equation with cubic nonlinearity -- with repulsive interaction under different trapping potentials including box potential and harmonic potential. Based on our asymptotic and numerical results, we formulate a gap conjecture on the fundamental gaps in energy and chemical potential of the GPE on bounded domains with the homogeneous Dirichlet boundary condition, and in the whole space with a convex trapping potential growing at least quadratically in the far field. We then extend these results to the GPE on bounded domains with either the homogeneous Neumann boundary condition or periodic boundary condition.
  • Experiments, theory and atomistic simulations show that finite triple junction mobility results in non-equilibrium triple junction angles in evolving polycrystalline systems. These angles have been predicted and verified for cases where grain boundary migration is steady-state. Yet, steady-state never occurs during the evolution of polycrystalline microstructures as a result of changing grain size and topological events (e.g., grain face/edge switching - "$T_1$" process, or grain disappearance "$T_2$" or "$T_3$" processes). We examine the non-steady evolution of the triple junction angle in the vicinity of topological events and show that large deviations from equilibrium and/or steady-state angles occur. We analyze the characteristic relaxation time of triple junction angles $\tau$ by consideration of a pair of topological events, beginning from steady-state migration. Using numerical results and theoretical analysis we predict how the triple junction angle varies with time and how $\tau$ varies with triple junction mobility. We argue that it is precisely those cases where grain boundaries are moving quickly (e.g., topological process in nanocrystalline materials), that the classical steady-state prediction of the finite triple junction mobility triple junction angle is inapplicable and may only be applied qualitatively.
  • We study analytically and numerically stability and interaction patterns of quantized vortex lattices governed by the reduced dynamical law -- a system of ordinary differential equations (ODEs) -- in superconductivity. By deriving several non-autonomous first integrals of the ODEs, we obtain qualitatively dynamical properties of a cluster of quantized vortices, including global existence, finite time collision, equilibrium solution and invariant solution manifolds. For a vortex lattice with 3 vortices, we establish orbital stability when they have the same winding number and find different collision patterns when they have different winding numbers. In addition, under several special initial setups, we can obtain analytical solutions for the nonlinear ODEs.
  • We propose an efficient and accurate parametric finite element method (PFEM) for solving sharp-interface continuum models for solid-state dewetting of thin films with anisotropic surface energies. The governing equations of the sharp-interface models belong to a new type of high-order (4th- or 6th-order) geometric evolution partial differential equations about open curve/surface interface tracking problems which include anisotropic surface diffusion flow and contact line migration. Compared to the traditional methods (e.g., marker-particle methods), the proposed PFEM not only has very good accuracy, but also poses very mild restrictions on the numerical stability, and thus it has significant advantages for solving this type of open curve evolution problems with applications in the simulation of solid-state dewetting. Extensive numerical results are reported to demonstrate the accuracy and high efficiency of the proposed PFEM.
  • We present several numerical methods and establish their error estimates for the discretization of the nonlinear Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter $0<\varepsilon\ll 1$ which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. We begin with the conservative Crank-Nicolson finite difference (CNFD) method and establish rigorously its error estimate which depends explicitly on the mesh size $h$ and time step $\tau$ as well as the small parameter $0<\varepsilon\le 1$. Based on the error bound, in order to obtain `correct' numerical solutions in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the CNFD method requests the $\varepsilon$-scalability: $\tau=O(\varepsilon^3)$ and $h=O(\sqrt{\varepsilon})$. Then we propose and analyze two numerical methods for the discretization of the nonlinear Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their $\varepsilon$-scalability is improved to $\tau=O(\varepsilon^2)$ and $h=O(1)$ when $0<\varepsilon\ll 1$ compared with the CNFD method. Extensive numerical results are reported to confirm our error estimates.
  • We present a uniformly accurate finite difference method and establish rigorously its uniform error bounds for the Zakharov system (ZS) with a dimensionless parameter $0<\varepsilon\le 1$, which is inversely proportional to the speed of sound. In the subsonic limit regime, i.e., $0<\varepsilon\ll 1$, the solution propagates highly oscillatory waves and/or rapid outgoing initial layers due to the perturbation of the wave operator in ZS and/or the incompatibility of the initial data which is characterized by two nonnegative parameters $\alpha$ and $\beta$. Specifically, the solution propagates waves with $O(\varepsilon)$- and $O(1)$-wavelength in time and space, respectively, and amplitude at $O(\varepsilon^{\min\{2,\alpha,1+\beta\}})$ and $O(\varepsilon^\alpha)$ for well-prepared ($\alpha\ge1$) and ill-prepared ($0\le \alpha<1$) initial data, respectively. This high oscillation of the solution in time brings significant difficulties in designing numerical methods and establishing their error bounds, especially in the subsonic limit regime. A uniformly accurate finite difference method is proposed by reformulating ZS into an asymptotic consistent formulation and adopting an integral approximation of the oscillatory term. By adapting the energy method and using the limiting equation via a nonlinear Schr\"{o}dinger equation with an oscillatory potential, we rigorously establish two independent error bounds and obtain error bounds at $O(h^2+\tau^{4/3})$ and $O(h^2+\tau^{1+\frac{\alpha}{2+\alpha}})$ for well-prepared and ill-prepared initial data, respectively, which are uniform in both space and time for $0<\varepsilon\le 1$ and optimal at the second order in space. Numerical results are reported to demonstrate that our error bounds are sharp.
  • We analyze rigorously error estimates and compare numerically spatial/temporal resolution of various numerical methods for the discretization of the Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter $0<\varepsilon\ll 1$ which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. We begin with several frequently used finite difference time domain (FDTD) methods and obtain rigorously their error estimates in the nonrelativistic limit regime by paying particular attention to how error bounds depend explicitly on mesh size $h$ and time step $\tau$ as well as the small parameter $\varepsilon$. Based on the error bounds, in order to obtain `correct' numerical solutions in the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the FDTD methods share the same $\varepsilon$-scalability on time step: $\tau=O(\varepsilon^3)$. Then we propose and analyze two numerical methods for the discretization of the Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their $\varepsilon$-scalability on time step is improved to $\tau=O(\varepsilon^2)$ when $0<\varepsilon\ll 1$. Extensive numerical results are reported to support our error estimates.
  • We propose and rigourously analyze a multiscale time integrator Fourier pseudospectral (MTI-FP) method for the Dirac equation with a dimensionless parameter $\varepsilon\in(0,1]$ which is inversely proportional to the speed of light. In the nonrelativistic limit regime, i.e. $0<\varepsilon\ll 1$, the solution exhibits highly oscillatory propagating waves with wavelength $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively. Due to the rapid temporal oscillation, it is quite challenging in designing and analyzing numerical methods with uniform error bounds in $\varepsilon\in(0,1]$. We present the MTI-FP method based on properly adopting a multiscale decomposition of the solution of the Dirac equation and applying the exponential wave integrator with appropriate numerical quadratures. By a careful study of the error propagation and using the energy method, we establish two independent error estimates via two different mathematical approaches as $h^{m_0}+\frac{\tau^2}{\varepsilon^2}$ and $h^{m_0}+\tau^2+\varepsilon^2$, where $h$ is the mesh size, $\tau$ is the time step and $m_0$ depends on the regularity of the solution. These two error bounds immediately imply that the MTI-FP method converges uniformly and optimally in space with exponential convergence rate if the solution is smooth, and uniformly in time with linear convergence rate at $O(\tau)$ for all $\varepsilon\in(0,1]$ and optimally with quadratic convergence rate at $O(\tau^2)$ in the regimes when either $\varepsilon=O(1)$ or $0<\varepsilon\lesssim \tau$. Numerical results are reported to demonstrate that our error estimates are optimal and sharp. Finally, the MTI-FP method is applied to study numerically the convergence rates of the solution of the Dirac equation to those of its limiting models when $\varepsilon\to0^+$.
  • We propose a sharp-interface continuum model based on a thermodynamic variational approach to investigate the strong anisotropic effect on solid-state dewetting including contact line dynamics. For sufficiently strong surface energy anisotropy, we show that multiple equilibrium shapes may appear that can not be described by the widely employed Winterbottom construction, i.e., the modified Wulff construction for an island on a substrate. We repair the Winterbottom construction to include multiple equilibrium shapes and employ our evolution model to demonstrate that all such shapes are dynamically accessible.
  • We derive rigorously one- and two-dimensional mean-field equations for cigar- and pancake-shaped Bose-Einstein condensates (BEC) with higher order interactions (HOI). We show how the higher order interaction modifies the contact interaction of the strongly confined particles. Surprisingly, we find that the usual Gaussian profile assumption for the strongly confining direction is inappropriate for the cigar-shaped BEC case, and a Thomas-Fermi type profile should be adopted instead. Based on the derived mean field equations, the Thomas-Fermi densities are analyzed in presence of the contact interaction and HOI. For both box and harmonic traps in one, two and three dimensions, we identify the analytical Thomas-Fermi densities, which depend on the competition between the contact interaction and the HOI.
  • We introduce and analyze a novel mean-field model for polariton condensates with velocity dependence of the effective polariton mass due the photon and exciton components. The effective mass depends on the in-plane wave vector k, which at the inflection point of the lower polariton energy branch becomes infinite and above negative. The polariton condensate modes of the new mean-field theory are now sensitive to mass variations and for certain points of the energy dispersion the polariton condensate mode represents fractional quantum mechanics. The impact of the generalized kinetic energy term is elucidated by numerical studies in 1D and 2D showing significant differences for large velocities. Analytical expressions for plane wave solutions as well as a linear waves analysis show the significance of this new model.
  • A multiscale time integrator Fourier pseudospectral (MTI-FP) method is proposed and analyzed for solving the Klein-Gordon-Schr\"{o}dinger (KGS) equations in the nonrelativistic limit regime with a dimensionless parameter $0<\varepsilon\le1$ which is inversely proportional to the speed of light. In fact, the solution to the KGS equations propagates waves with wavelength at $O(\varepsilon^2)$ and $O(1)$ in time and space, respectively, when $0<\varepsilon\ll 1$, which brings significantly numerical burdens in practical computation. The MTI-FP method is designed by adapting a multiscale decomposition by frequency to the solution at each time step and applying the Fourier pseudospectral discretization and exponential wave integrators for spatial and temporal derivatives, respectively. We rigorously establish two independent error bounds for the MTI-FP at $O(\tau^2/\varepsilon^2+h^{m_0})$ and $O(\varepsilon^2+h^{m_0})$ for $\varepsilon\in(0,1]$ with $\tau$ time step size, $h$ mesh size and $m_0\ge 4$ an integer depending on the regularity of the solution, which imply that the MTI-FP converges uniformly and optimally in space with exponential convergence rate if the solution is smooth, and uniformly in time with linear convergence rate at $O(\tau)$ for $\varepsilon\in(0,1]$ and optimally with quadratic convergence rate at $O(\tau^2)$ in the regime when either $\varepsilon=O(1)$ or $0<\varepsilon\le \tau$. Thus the meshing strategy requirement (or $\varepsilon$-scalability) of the MTI-FP is $\tau=O(1)$ and $h=O(1)$ for $0<\varepsilon\ll 1$, which is significantly better than classical methods. Numerical results demonstrate that our error bounds are optimal and sharp. Finally, the MTI-FP method is applied to study numerically convergence rates of the KGS equations to the limiting models in the nonrelativistic limit regime.
  • In this paper, we propose a regularized Newton method for computing ground states of Bose-Einstein condensates (BECs), which can be formulated as an energy minimization problem with a spherical constraint. The energy functional and constraint are discretized by either the finite difference, or sine or Fourier pseudospectral discretization schemes and thus the original infinite dimensional nonconvex minimization problem is approximated by a finite dimensional constrained nonconvex minimization problem. Then an initial solution is first constructed by using a feasible gradient type method, which is an explicit scheme and maintains the spherical constraint automatically. To accelerate the convergence of the gradient type method, we approximate the energy functional by its second-order Taylor expansion with a regularized term at each Newton iteration and adopt a cascadic multigrid technique for selecting initial data. It leads to a standard trust-region subproblem and we solve it again by the feasible gradient type method. The convergence of the regularized Newton method is established by adjusting the regularization parameter as the standard trust-region strategy. Extensive numerical experiments on challenging examples, including a BEC in three dimensions with an optical lattice potential and rotating BECs in two dimensions with rapid rotation and strongly repulsive interaction, show that our method is efficient, accurate and robust.
  • In this paper, we propose efficient and accurate numerical methods for computing the ground state and dynamics of the dipolar Bose-Einstein condensates utilising a newly developed dipole-dipole interaction (DDI) solver that is implemented with the non-uniform fast Fourier transform (NUFFT) algorithm. We begin with the three-dimensional (3D) Gross-Pitaevskii equation (GPE) with a DDI term and present the corresponding two-dimensional (2D) model under a strongly anisotropic confining potential. Different from existing methods, the NUFFT based DDI solver removes the singularity by adopting the spherical/polar coordinates in Fourier space in 3D/2D, respectively, thus it can achieve spectral accuracy in space and simultaneously maintain high efficiency by making full use of FFT and NUFFT whenever it is necessary and/or needed. Then, we incorporate this solver into existing successful methods for computing the ground state and dynamics of GPE with a DDI for dipolar BEC. Extensive numerical comparisons with existing methods are carried out for computing the DDI, ground states and dynamics of the dipolar BEC. Numerical results show that our new methods outperform existing methods in terms of both accuracy and efficiency.
  • We study dimension reduction for the three-dimensional Gross-Pitaevskii equation with a long-range and anisotropic dipole-dipole interaction modeling dipolar Bose-Einstein condensation in a strong interaction regime. The cases of disk shaped condensates (confinement from dimension three to dimension two) and cigar shaped condensates (confinement to dimension one) are analyzed. In both cases, the analysis combines averaging tools and semiclassical techniques. Asymptotic models are derived, with rates of convergence in terms of two small dimensionless parameters characterizing the strength of the confinement and the strength of the interaction between atoms.
  • We present efficient and accurate numerical methods for computing the ground state and dynamics of the nonlinear Schr\"odinger equation (NLSE) with nonlocal interactions based on a fast and accurate evaluation of the long-range interactions via the nonuniform fast Fourier transform (NUFFT). We begin with a review of the fast and accurate NUFFT based method in \cite{JGB} for nonlocal interactions where the singularity of the Fourier symbol of the interaction kernel at the origin can be canceled by switching to spherical or polar coordinates. We then extend the method to compute other nonlocal interactions whose Fourier symbols have stronger singularity at the origin that cannot be canceled by the coordinate transform. Many of these interactions do not decay at infinity in the physical space, which adds another layer of complexity since it is more difficult to impose the correct artificial boundary conditions for the truncated bounded computational domain. The performance of our method against other existing methods is illustrated numerically, with particular attention on the effect of the size of the computational domain in the physical space. Finally, to study the ground state and dynamics of the NLSE, we propose efficient and accurate numerical methods by combining the NUFFT method for potential evaluation with the normalized gradient flow using backward Euler Fourier pseudospectral discretization and time-splitting Fourier pseudospectral method, respectively. Extensive numerical comparisons are carried out between these methods and other existing methods for computing the ground state and dynamics of the NLSE with various nonlocal interactions. Numerical results show that our scheme performs much better than those existing methods in terms of both accuracy and efficiency.
  • We propose a sharp interface model for simulating solid-state dewetting where the surface energy is (weakly) anisotropic. The morphology evolution of thin films is governed by surface diffusion and contact line migration. The mathematical model is based on an energy variational approach. Anisotropic surface energies lead to multiple solutions of the contact angle equation at contact points. Introduction of a finite contact point mobility is both physically based and leads to robust, unambiguous determination of the contact angles. We implement the mathematical model in an explicit finite difference scheme with cubic spline interpolation for evolving marker points. Following validation of the mathematical and numerical approaches, we simulate the evolution of thin film islands, semi-infinite films, and films with holes as a function of film dimensions, Young's angle $\theta_i$, anisotropy strength and crystal symmetry, and film crystal orientation relative to the substrate normal. We find that the contact point retraction rate can be well described by a power-law, $l \sim t^n$. Our results demonstrate that the exponent $n$ is not universal -- it is sensitive to the Young's angle $\theta_i$ (and insensitive to anisotropy). In addition to classical wetting (where holes in a film heal) and dewetting (where holes in a film grow), we observe cases where a hole through the film heals but leave a finite size hole/bubble between the continuous film and substrate or where the hole heals leaving a continuous film that is not bonded to the substrate. Surface energy anisotropy (i) increases the instability that leads to island break-up into multiple islands, (ii) enhances hole healing, and (iii) leads to finite island size even under some conditions where the Young's angle $\theta_i$ suggests that the film wets the substrate.