{'fraction_non_alphanumeric': 0.036606756997229074, 'fraction_numerical': 0.01110913389429799, 'mean_word_length': 4.684682080924856, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structure and stability of two-dimensional Bose-Einstein condensates under both harmonic and lattice confinement

In this work, we study pancake-shaped Bose-Einstein condensates confined by both a cylindrically symmetric harmonic potential and an optical lattice with equal periodicity in two orthogonal directions. We first identify the spectrum of the underlying two-dimensional linear problem through multiple-scale techniques. Then, we use the results obtained in the linear limit as a starting point for a nonlinear existence and stability analysis of the lowest energy states, emanating from the linear ones, in the nonlinear problem. Two-parameter continuations of these states are performed for increasing nonlinearity and optical lattice strengths, and their instabilities and temporal evolution are investigated. It is found that the ground state as well as one of the excited states are either stable or weakly unstable for both attractive and repulsive interatomic interactions.


Introduction
The last decade has witnessed a tremendous amount of research efforts in the physics of atomic Bose-Einstein condensates (BECs) [1,2]. The study of BECs has yielded a wide array of interesting phenomena, not only because of very precise experimental control that exists over the relevant experimental procedures [3,4], but also because of the intimate connections of the description of dilute-gas BECs with other areas of physics, such as superfluidity, superconductivity, lasers and coherent optics, nonlinear optics, and nonlinear wave theory. Of particular emphasis in much experimental and theoretical work is the setting of a BEC trapped in periodic potentials, usually combined with an additional harmonic trapping potential. From the standpoint of nonlinear interactions, mathematical descriptions of BECs held in purely harmonic traps are now well known. Nevertheless, apart from studies focusing on the transition between superfluidity and an insulator state [5], relatively little attention has been given to an understanding of the varieties of many-body states that may possibly exist with intermediate lattice strengths, where phase coherence is maintained across the sample. A more complete understanding of BEC behavior in such lattice potentials is relevant and important to current work with BECs, and to an even broader array of topics, in particular discrete nonlinear optics and nonlinear wave theories. Such regimes of BEC physics are experimentally and theoretically accessible, and comparisons between theoretical and experimental results are certainly possible. Here, we present a theoretical stability examination of BECs with either attractive or repulsive interatomic interactions in a combined harmonic and periodic potential.
Many of the common elements between BECs and other areas of physics, and in particular optics, originate in the existence of macroscopic coherence in the many-body state of the system. Mathematically, BEC dynamics are therefore often accurately described by a mean-field model, namely a partial differential equation of the nonlinear Schrödinger (NLS) type, the so-called Gross-Pitaevskii (GP) equation [1,2,3]. The GP equation is particularly successful in drawing connections between BEC physics, nonlinear optics and nonlinear wave theories, with vortices and solitary waves examples of common elements between these areas. The GP equation is a classical nonlinear evolution equation (with the nonlinearity originating from the interatomic interactions) and, as such, it permits the study of a variety of interesting nonlinear phenomena. These phenomena have primarily been studied by treating the condensate as a purely nonlinear coherent matter-wave, i.e., from the viewpoint of the nonlinear dynamics of solitary waves. Relevant studies have already been summarized in various books (see, e.g., Ref. [4]) and reviews (see, e.g., Refs. [6] for bright matter-wave solitons, [7,8] for vortices in BECs, [9] for dynamical instabilities in BECs, [10,11] for nonlinear dynamics of BECs in optical lattices).
On the other hand, many static and dynamic properties of BECs confined in various types of external potentials can be studied by starting from the non-interacting limit, where the nonlinearity is considered to be negligible. The basic idea of such an approach is that in the absence of interactions the GP equation is reduced to a linear Schrödinger equation for a confined single-particle state; in this limit, and in the case of, e.g., a harmonic external potential, the linear problem becomes the equation for the quantum harmonic oscillator characterized by discrete energies and corresponding eigenmodes [12]. Exploiting this simple physical picture, one may then use analytical and/or numerical techniques for the continuation of these linear eigenmodes supported by the particular type of the external trapping potential into nonlinear states as the interactions become stronger. This idea has been explored at the level of one-dimensional (1D) [13] and higher-dimensional states [14] in the case of a harmonic trapping potential, where nonlinear stationary modes were found from a continuation of the (linear) states of the quantum harmonic oscillator. The same problem has been studied in the framework of the so-called Feshbach resonance management technique in Ref. [15], where a linear temporal variation of the nonlinearity was considered. The continuation of the linear states to their nonlinear counterparts has also been explored from the point of view of bifurcation and stability theory [16]. Finally, in the same spirit but in the twodimensional (2D) setting, radially symmetric nonlinear states of harmonically trapped "pancake-shaped" condensates were recently investigated in Ref. [17].
Importantly, all of the above studies provide a clear physical picture of how genuinely nonlinear states of harmonically confined BECs (such as dark and bright matter-wave solitons in 1D or ring solitons and vortices in 2D), are connected to and emanate from the eigenmodes of the quantum harmonic oscillator. Similar considerations also hold for BECs confined in optical lattices. In this case, pertinent nonlinear stationary states (such as spatially extended nonlinear Bloch waves, truncated nonlinear Bloch waves, matter-wave gap solitons in 1D, and gap vortices in 2D and 3D), can be understood by the structure of the band-gap spectrum of the linear Bloch waves supported in the non-interacting limit (see, e.g., Ref. [18] and references therein). However, there are only few studies for condensates confined in both harmonic and optical lattice potentials, and these are basically devoted to the dynamics of particular nonlinear structures (such as dark [19] and bright [20] solitons in 1D, and vortices in 2D [21]). Thus, the structure of condensates confined in such superpositions of harmonic and periodic potentials remains, to the best of our knowledge, largely unexplored.
Nevertheless, such a study is particularly relevant to current work with BECs, and even suggests new avenues for exploration. In particular, one might ask whether the addition of a weak optical lattice might increase the stability of excited states that are known to be unstable in harmonic traps. Stability, if it is found, may add new realistic options for the engineering of new quantum states of BECs. It may also be interesting to investigate the Mott insulator transition by starting from a stable excited state of a weak optical lattice. Also, the transport of excited states (which is not discussed in this paper) through a lattice structure may have application in future BEC interferometry experiments. Finally, the advances of far-off-resonant optical trapping techniques allow for the creation of strongly pancake-shaped condensates that may be confined by harmonic and spatially periodic components, and we expect that the theoretical considerations described here may be directly explored using current experimental techniques.
Our aim in the present work is to contribute to this direction and study the structure and the stability of a pancake-shaped condensate confined by the combination of a harmonic trap and a periodic potential, with periodicity in two orthogonal directions. We will adopt the above mentioned approach of continuation of linear states to nonlinear ones, thus providing a host of interesting solutions that have not been explored previously and yet should be tractable within presently available experimental settings. In particular, our analysis starts by first considering the non-interacting limit. In this regime, we employ a multiscale perturbation method (which uses the harmonic trap strength as a formal small parameter) to find the discrete energies and the corresponding eigenmodes of the pertinent single-particle Schrödinger equation with the combined harmonic and periodic potential. We then use this linear limit as a starting point for initializing a 2D "nonlinear solver" that identifies the relevant stationary nonlinear eigenstates as a function of the chemical potential (i.e., the nonlinearity strength) and of the optical lattice depth.
Once the basic structure of the condensate is found, we subsequently perform a linear stability analysis of the nonlinear modes that can be initiated by the noninteracting ground and first few excited states. When nonlinear states are found to be unstable, we use direct numerical simulations to study their dynamics and monitor the evolution of the relevant instability. Essential results that will be presented below are the following: for a fixed harmonic trap strength, there exist certain regions in the parameter plane defined by the chemical potential and the optical lattice depth, where not only the ground state, but also excited states are stable or only weakly unstable. Particularly, an excited state with a shape resembling an out-of-phase matter-wave soliton pair (for attractive interactions) is found to persist for long times, being stable (weakly unstable) for attractive (repulsive) interatomic interactions. Thus, the ground state and the aforementioned excited state have a good chance to be observed in a real experiment with either attractive or repulsive pancake BECs. Similar conclusions can be drawn for more complex states, such as a quadrupolar one which may also be stable in the attractive case, however, higher excited states are typically more prone to instabilities, as is shown in our detailed numerics below.
The paper is organized as follows: In Section II we present the model and study analytically the non-interacting regime. The continuation of the linear states to the nonlinear ones, as well as the stability properties of the nonlinear states are presented in Section III. Finally, in Section IV, we summarize our findings and present our conclusions.

The model and its analytical consideration
At sufficiently low temperatures, and in the framework of the mean-field approach, the condensate dynamics can be described by the order parameter Ψ(r, t). We assume that the condensate is kept in a highly anisotropic trap, with the transverse (x, y) and longitudinal (z) trapping frequencies chosen so that ω x = ω y ≡ ω ⊥ ≪ ω z . In such a case, the condensate has a nearly planar, so-called "pancake" shape (see, e.g., Refs. [22] for relevant experimental realizations), which allows us to assume a separable wave function, Ψ = Φ(z)ψ(x, y), where Φ(z) is the ground state of the respective quantum harmonic oscillator. Then, averaging of the underlying three-dimensional (3D) GP equation along the longitudinal direction z [23] leads to the following 2D GP equation for the transverse component of the wave function (see also Refs. [8,9,4]): Here, ∆ ≡ ∂ 2 x + ∂ 2 y is the 2D Laplacian, m is the atomic mass, and g 2D = g 3D / √ 2πa z is an effective 2D coupling constant, where g 3D = 4πh 2 a/m (a being the scattering length) and a z = h/mω z is the longitudinal harmonic oscillator length. Finally, the potential V ext (x, y) in the GP Eq. (1) is assumed to consist of a harmonic component and a square 2D optical lattice (OL) created by two pairs of interfering laser beams of wavelength λ: In the above expression, r 2 ≡ x 2 + y 2 , while the optical lattice is characterized by two parameters, namely its depth V 0 and its periodicity d = π/k = (λ/2)/ sin(θ/2), where θ is the angle between the two beams that create the x-direction lattice, and between the two beams that create the y-direction lattice.
Measuring length in units of a L = d/π, time in units of ω −1 L =h/E L , and energy in units of E L = 2E rec =h 2 /ma 2 L (where E rec is the lattice recoil energy), the GP Eq. (1) can be put into the following dimensionless form: In the normalized GP Eq. (3), the wavefunction is rescaled as , the parameter s is given by s = sign(g 2D ) = ±1 (with s = +1 or s = −1 corresponding, respectively, to repulsive or attractive interatomic interactions), while the normalized trapping potential V (x, y) is now given by: In the above equation, the normalized lattice depth V 0 is measured in units of 4E rec , while the normalized harmonic trap strength is given by where a ⊥ = h/mω ⊥ is the transverse harmonic oscillator length. Note that using realistic parameter values (see, e.g., Ref. [24]), namely, a lattice periodicity 0.3 µm, a recoil energy E rec /h ∼ 6KHz (assuming an atomic mass corresponding to 87 Rb), and using a transverse trap frequency ω ⊥ = 2π × 5Hz, the parameter Ω is of order of 10 −4 ; thus, it is a natural small parameter of the problem. Our analysis starts by considering the non-interacting limit s → 0, in which the GP equation becomes a linear Schrödinger equation. Then, seeking stationary localized solutions of the form ψ(x, y, t) = exp(−iE m,n t)u m,n (x, y) [where E m,n are discrete energies and u m,n (x, y) are the corresponding linear eigenmodes], and rescaling spatial variables by √ Ω, we obtain the following equation: The next step is to separate variables through u(x, y) = u m (x)u n (y) and split the energy into E m,n = E m + E n , to obtain two 1D eigenvalue problems of the same type as above, namely, and a similar one for y (with x replaced by y and the subscript m replaced by n).
We now restrict ourselves to the physically relevant regime of 0 < Ω ≪ 1 as discussed above. In this case, we may use ν ≡ √ Ω as a formal small parameter and develop methods of multiple scales and homogenization techniques [16] in order to obtain analytical predictions for the linear spectrum. In particular, introducing the fast (i.e., rapidly varying) variable X = x/ν and rescaling the energy as ε m = E m /Ω, the eigenvalue problem of Eq. (7) is expressed as follows: where (note that V 0 /Ω was treated as an O(1) parameter). Additionally, we consider a formal series expansion (in ν) for u m and ε m [16], namely, To this end, substitution of this expansion in the eigenvalue problem of Eq.
while the corresponding eigenfunction is given by: where Combining the results in the two orthogonal directions x and y yields a total energy eigenvalue and a corresponding eigenfunction which up to normalization factors can be written as: The particularly appealing feature of this expression is that it allows us to combine various (ground or excited) states in the x direction with different ones along the y direction. In this work we will focus on the simplest possible combinations of m, n ∈ {0, 1, 2} and examine the various states generated by the combinations of these, which we hereafter denote |m, n . In particular, below we will focus on the ground state |0, 0 and the excited states |1, 0 , |1, 0 + |0, 1 , |1, 1 , and |2, 0 . Our aim is to investigate which of these states persist in the nonlinear regime in both cases of attractive and repulsive interatomic interactions, and study the stability of these states in detail. Notice that we will only illustrate (by an appropriate curve in the numerical results that follow) the linear limit, E m,n of Eq. (15) above, for the case of attractive interactions.

The non-interacting limit
We start by examining the validity of the above analytical predictions concerning the linear limit of the problem, namely, Eqs. (15) and (16). The results are summarized in Fig. 1. Panel (a) shows the 1D harmonic oscillator energy spectrum (for Ω = 0.1) and compares it with the energy spectrum obtained from numerical and approximate theoretical [see Eq. (13)] solutions of the combined harmonic and optical lattice potential for V 0 = 0.3. Panel (b) offers a similar comparison but for a larger lattice depth, namely V 0 = 0.5. One can clearly see that the theoretical calculation approximates very accurately the numerical results for the first few states (i.e., n = 0, 1, 2), while deviations become more significant for higher-order excited states. Panels (c), (d) and (e) show the zeroth, first and second eigenfunction of the purely harmonic potential (thick solid line), as well as of the harmonic trap and optical lattice potential as found numerically (thin solid line) and analytically [given by Eq. (14)]. The green dash-dotted line represents the form of the combined potential. We once again note the good agreement of our analytical results above in comparison with the full numerical computation.

The approach for the interacting case
We now consider the full nonlinear problem of Eq. (3). In the following, we will monitor the two-dimensional, V 0 × µ, parameter space (where µ is the chemical potential of the relevant modes) for nonlinear excitations that stem from the linear spectrum of the problem. We perform the relevant analysis first in the case of attractive interatomic interactions and then in the case of repulsive ones. Notice that the parameter Ω will be fixed to a relatively large value, namely Ω = 0.1; this is done for convenience in our numerical simulations (such "large" values of Ω correspond to smaller condensates that can be analyzed numerically with relatively coarser spatial grids), but we have checked that our results remain qualitatively similar for smaller values of Ω (results not shown here). Nevertheless, it should be noted that even such a value of Ω, together with the considered range of values of the other normalized parameters (chemical potential, number of atoms, lattice depth, etc -see below) is still physically relevant. For example, our choice may realistically correspond to a 87 Rb condensate containing ∼ 15, 000 atoms, confined by a harmonic potential with frequencies ω z = 20ω ⊥ = 2π × 240Hz (so that ω ⊥ = 2π × 12 Hz) and an optical lattice potential with a periodicity d ≈ 3µm. The recoil energy in this case is E rec /h = 60 Hz (so that ω L = 2π × 120 Hz, giving Ω = 0.1), and a lattice depth of V 0 = 0.3 corresponds to ∼1.2 E rec . To further set the scale for the simulations described below, such a BEC with repulsive interactions in the purely harmonic trap (where the optical lattice is not applied) would have a chemical potential of µ ∼ 0.5 in our dimensionless units. In the following sections, stationary solutions of the full nonlinear problem are sought in the form where µ m,n (which is the nonlinear analog of the energy E m,n found in the non-interacting limit) represents the chemical potential. Note that we will henceforth avoid using subscripts m and n when the meaning is clear, in the interest of avoiding notational clutter.

Existence and Stability
We begin by looking at the case of attractive interatomic interactions. The most fundamental solutions are those belonging to the |0, 0 branch, which represents the ground state of the system and is shown in Fig. 2. The top left panel of this figure shows the diagnostic that we will typically use to follow the V 0 × µ surface, namely the rescaled number of particles N(V 0 , µ) = |u m,n | 2 dxdy as a function of the chemical potential µ introduced above, and the optical lattice depth V 0 . Essentially, the grayscale values in this plot correspond to the number of atoms needed to obtain a particular chemical potential with a particular lattice depth; lighter values correspond to more atoms. As N becomes smaller, through the appropriate variation of µ, we approach the linear limit so one expects the solution to degenerate to the corresponding linear eigenmode (for µ tending to the corresponding eigenvalue of the linear problem). Figure 2 shows our observations for this fundamental branch, which seems to disappear for µ m,n (V 0 ) ≈ E m,n (V 0 ) = − 1 2 V 2 0 + 1 − 1 4 V 2 0 Ω (n + m + 1). Naturally the surface degenerates to its linear limit for µ m,n (V 0 ) → E m,n (V 0 ) and the number of particles is a decreasing function of µ, contrary to what is the case in the repulsive nonlinearity (see below). This is a well-known difference between the two cases that has been documented elsewhere (see, e.g., Refs. [15,16]). It is relevant to indicate here that although there is a direct correspondence between the atom number N and the chemical potential µ, we opt to illustrate our results as a function of µ (and V 0 ), since the latter is the relevant parameter entering the mathematical setup of the problem and it is the one for which we developed an analytical prediction in the linear limit. Nevertheless, we also give N(V 0 , µ), so as to associate in each case the relevant chemical potential (and lattice strength) for a given configuration with the corresponding physical quantity, i.e., the atom number.
It is important to highlight here that that the numerical computations have been performed in a domain of size 201 × 201, with ∆x = ∆y = 0.15. The size of the grid weakly affects the value of the respective eigenvalues. In particular, the energy eigenvalues in the non-interacting limit of the one-dimensional decoupled eigenvalue problem (with V 0 = 0.5) for the n = 0 and n = 1 mode that we report below are E 0 ≈ −0.01453 and E 1 ≈ 0.07639 respectively for this coarser domain, while for the a domain of size 3001 with ∆x = 0.01 they become −0.01408 and 0.07692 respectively (computations show the discrepancy is uniformly smaller for smaller values of V 0 ). This feature will weakly affect the quantitative aspects of the results that follow (in essence, providing an error bar in the estimates below of ∆µ ≈ 5 × 10 −4 and similar for the eigenvalues λ introduced below), but is essentially necessary, given the limitations of standard eigenvalue solvers for such big Jacobian eigenvalue problems.
The linear stability of the solutions is analyzed by using the following standard ansatz for the perturbation, where λ = λ r + iλ i is the generally complex eigenvalue (subscripts r and i denote, respectively, the real and imaginary parts of λ) of the ensuing Bogoliubov-de Gennes equations [1,2,3,4], and (a, b) T is the corresponding eigenvector. The real part λ r of the eigenvalue then determines the growth rate of the potential instability of the solution, since for positive real values the perturbation will grow exponentially in time. The imaginary part λ i denotes the oscillatory part (frequency) of the relevant eigenmode. The top right of Fig. 2 depicts the stability domain, denoted by S(V 0 , µ) = max λ (λ r ), in terms of the maximum real part of all eigenvalues as a function of the lattice depth V 0 and the chemical potential µ. This quantity S is a particularly important one from a physical point of view since if a perturbation to the system has initially a projection p(0) onto the most unstable eigenmode of the linearization, then this perturbation will grow in time according to ||p(t)|| = ||p(0)|| exp(St) while the solution is sufficiently close in space to the original profile. Hence, given the initial condition profile [which determines p(0)] and S, we can determine for an unstable configuration the time t required for the instability to manifest itself, i.e., for p(t) to become of the order of the solution amplitude.
It is important to note, in connection to our numerical linear stability results, that the |0, 0 branch can become unstable for µ < µ cr (V 0 ) (see top right panel of Fig. 2) due to the appearance of a real pair of eigenvalues. This instability for large N is something that may be expected in the case of attractive interactions under consideration, as the corresponding 2D GP equation for an a homogeneous BEC (i.e., without any external potential) is well-known to be subject to collapse [25]. However, it should also be  expected that very close to the linear limit the growth rate of the instability is essentially zero (cf. with the top left panel of Fig. 2 of Ref. [17] for V 0 = 0 which is not shown here). Essentially, the potential appears to stabilize the solitary wave against dispersion in this regime (i.e., close to the linear limit), but cannot stabilize it against the catastrophic collapse-type instability. Furthermore, in the presence of the optical lattice we can observe that there is always a range of chemical potentials for which the condensate is stabilized, in accordance with what was originally suggested in Ref. [26]. Furthermore, even in the 3D case it is in principle possible to arrest collapse by appropriate choices of the parameters [27].
Next, we consider real-valued solutions with m + n = 1. The |1, 0 state (again in the case of attractive interactions) is shown in Fig. 3. This branch is always unstable, due to up to two real eigenvalue pairs and one complex quartet. A typical example of the branch in the bottom panels of the figure reveals this instability.
The |1, 0 + |0, 1 configuration for the attractive interactions case is shown in     4. This configuration turns out to be unstable in a large fraction of the regime of parameters considered due to a quartet of complex eigenvalues. However, remarkably, as V 0 and µ are increased and decreased respectively, it is possible to actually trap this state in a linearly stable form (eliminating the relevant oscillatory instability). This indicates that it would be of particular interest to try to identify such a state (which resembles an out-of-phase soliton pair) in a real experiment. Also, as expected, the solution degenerates to its linear counterpart as µ → E. Images of a typical unstable solution and its complex quartet are shown along with a stable solution from the top right-hand region of the two-parameter diagram. We should also note in passing that states in the form of |1, 0 +i|0, 1 would produce a vortex waveform; however since such states have been studied in some detail earlier in Ref. [21] in a similar setting (i.e., in the presence of an external potential containing both harmonic and lattice components), we do not examine them in more detail here. We now turn to solutions featuring m + n = 2. First, we consider the |1, 1 branch for the attractive case in Fig. 5. In this case the solution may possess between one and three complex eigenvalue quartets in its linearization (the middle panel of the figure shows a particular unstable case where there are two such quartets). However, once again, there exists a region in the right side of the relevant parameter space [i.e., for appropriate (V 0 , µ)] where the solution is found to be linearly stable and all potential oscillatory instabilities are suppressed. The bottom panel of Fig. 5 shows such a linearly stable case of the quadrupolar configuration, which, again, should be experimentally accessible.
Finally, we consider the state |2, 0 , as depicted in Fig. 6. This configuration is highly unstable throughout our parameter space, with up to four real pairs and one complex quartet of eigenvalues. A typical example of the unstable configuration and its spectral plane of eigenvalues is shown in the bottom panel of Fig. 6.

Dynamics
Now we corroborate our existence and stability results (for the attractive interactions case) with an investigation of the actual dynamics of typical unstable solutions selected from the above families. For each case, the particular solution presented in the corresponding figures is perturbed with random noise distributed uniformly between −0.05 and 0.05 and integrated over time. It is important to note that although a random perturbation is used here to "emulate" the experimental noise, the system is deterministic and the sole relevant feature of any (generic) random perturbation is its projection onto the most unstable eigendirection(s) of the perturbed solution profile. These projections, as indicated above, will grow (determinstically) according to the corresponding growth rate. For the time propagation, we implement a standard 4th-order Runge-Kutta integrator scheme where we have numerical consistency and stability for the conservative time step of ∆t = 10 −3 . The results are compiled in Fig. 7. Panel (a) in Fig. 7 depicts the catastrophic instability of the ground state, |0, 0 , which is subject to collapse, occurring at t ≈ 25. Panel (b) shows similar behavior for the |1, 0 state, in which the two lobes appear to self-focus independently, although one eventually prevails and collapses for t ≈ 35. It is very interesting to note that while the |1, 0 state collapses, its more stable superposition with the |0, 1 state survives for longer times, as expected, and also eventually merges into a ground-state-like (single pulse) configuration, which was found to have a number of atoms just on the unstable side of the boundary of stability for such a structure. The resulting state actually survives for very long times, oscillating within one of the wells where it originally collected itself [see the panel (c)], apparently stabilized by the ensuing oscillations. Panels (d) and (e) show, respectively, the relatively rapid break up and subsequent collapse of the |1, 1 and |2, 0 states.

Repulsive interatomic interactions
3.4.1. Existence and Stability Now we will investigate the results pertaining to repulsive interatomic interactions for the same linear states examined above. We once again start with the ground state |0, 0 branch shown in Fig. 8. The top left panel of Fig. 8 shows the number of atoms N(V 0 , µ) as in the previous section. However, the linear stability S(V 0 , µ) for this case is omitted because, as may be expected, this branch is stable throughout the parameter space, in contrast to its attractive counterpart (which is subject to collapse).
We now turn to excited states with m + n = 1. Figure 9 shows features similar to the previous one, but now for the state |1, 0 . This branch is found to always be unstable due to the appearance of up to three real eigenvalue pairs. The top panels depict the dependence of the number of atoms N(V 0 , µ) (left), and instability growth rate S(V 0 , µ) (right) on the lattice depth V 0 and the chemical potential µ. A sample profile and its eigenvalue spectrum are given in the bottom panels, indicating the presence, in this case, of three real eigenvalues pairs. The next state we consider is the |1, 0 + |0, 1 state which is presented in Fig. 10. This state always possesses a quartet of complex eigenvalues, and up to two additional pairs of real eigenvalues, and is thus unstable for all µ. It it worth noticing, however, that the instability weakens to relatively benign small magnitude complex quartets for intermediate values of the chemical potential, roughly µ ∈ (0.4, 0.9), and large lattice depths, V 0 > 0.3. This suggests that such a configuration should be long-lived enough that it could be observable in experiments with repulsive condensates.
Next, we consider the states with n + m = 2 (again for repulsive interatomic interactions), starting with the |1, 1 branch in Fig. 11. The branch is also always It is interesting to note that such states are reminiscent of the domain-walls presented in Ref. [28] (here the domain-wall is imposed by the difference in phase), which however were found as potentially stable structures in multi-component condensates.
Next, the case of the |2, 0 state is shown in Fig. 12. Here, there are up to three complex quartets along with three real pairs of eigenvalues. It is notable that for higher values of the lattice depth, these states are deformed as the lattice "squeezes" the central maximum separating the two minima (see middle left panel of the figure). The contour plot shown in the bottom right panel suggests that further increase of lattice depth may lead to a new configuration altogether when the two local maxima eventually pinch off of  each other. This deformation is a direct consequence of the presence of the (repulsive) nonlinearity, which results in drastically different configurations as compared to the linear limit of the structure.

Dynamics
We performed numerical simulations to investigate the evolution of typical unstable states in the case of repulsive interatomic interactions, using similar time-stepping schemes as discussed above in the case of attractive interactions. Apart from the ground state, all excited states presented in the previous section were predicted to be unstable, and this is confirmed in this section. In the particular case of the state |1, 0 + |0, 1 which was found to be weakly unstable (see bottom row of Fig. 10), the instability takes a considerable time to manifest itself. The evolution of this state is depicted in panel (b) of Fig. 13; it is clearly seen that the state persists up to t ≈ 300. On the other hand, panel (a) shows the evolution of the state |1, 0 which persists up to t ≈ 40, while the bottom row of panels shows the dynamics of the (c) |1, 1 state and of the (d) |2, 0 , both persisting also up to t ≈ 40. All of these excited states degenerate  into ground-state-like configurations. Notice that transient vortex-like structures seem to appear during this process but they do not persist in the eventual dynamics and are hence not further discussed here.

Conclusions and discussion
In summary, we have studied the structure and the stability of a pancake-shaped condensate (with either attractive or repulsive interatomic interactions) confined in a potential with both a harmonic and an optical lattice component. Starting from the non-interacting limit, and exploiting the smallness of the harmonic trap strength, we have employed a multiscale perturbation method to find the discrete energies and the corresponding eigenmodes of the pertinent 2D linear Schrödinger equation. Then, we used the results found in this linear (non-interacting) limit in order to identify states persisting in the nonlinear (interacting) regime as well. This investigation revealed that the most fundamental states (emanating from combinations of the ground state and the first few excited states in the two orthogonal directions of the optical lattice) can indeed be continued in the nonlinear regime. To demonstrate this continuation, we used two-parameter diagrams involving the effective strength of the nonlinearity (through the chemical potential) and the optical lattice depth. Excited states were typically found to be unstable. The instability was found to result in either wavefunction collapse or a robust single-lobed structure in the case of attractive interactions; on the other hand, in the case of repulsive interactions, the instability was always found to lead to the ground state of the system. Nevertheless, noteworthy exceptions of stable or very weakly unstable states were also revealed. These include the |1, 0 + |0, 1 and the |1, 1 states in the case of attractive interatomic interactions. Moreover, in the case of repulsive interactions, the same state, |1, 0 +|0, 1 , was found (in certain parameter regimes) to be only very weakly unstable. Direct numerical simulations confirmed that the instability of this state is indeed weak, and it manifests itself at large times, an order of magnitude larger than the ones pertaining to the manifestation of instabilities of other excited states. Thus, it is clear that the obtained results suggest that the state |1, 0 + |0, 1 has a good chance to be observed in experiments with either an attractive or a repulsive pancake condensate. It is especially important to highlight that these states are stabilized (or quasi-stabilized) only in the presence of a sufficiently strong optical lattice; hence the latter potential plays a critical role in determining the stability of the states presented herein.
As described in Section III B, the parameters used in our analysis have been chosen in order to facilitate convenience of the numerical computations, while also within range of experimentally achievable limits of atom number, chemical potential, harmonic oscillator frequencies, and optical lattice depth and periodicity. Furthermore, we believe that the stability and spatial structure of the states examined here can be examined experimentally. For example, we imagine utilizing a BEC held in a pancakeshaped harmonic trap, created by an optical field. By using an optical trap rather than a magnetic trap, the scattering length of a BEC may be adjusted using a Feshbach resonance. We envision that an optical lattice potential is ramped on and superimposed on a BEC with an interatomic scattering length tuned to be near zero. Once the lattice has reached the desired depth, the scattering length can be further adjusted with a magnetic field to be either positive or negative (the latter option would need to be within a region of stability that does not result in collapse of the BEC). Finally, phase imprinting techniques [29] can be used to generate the desired phase profile of the BEC. By optically examining the state of the BEC at various points in time after phase profile imprinting, the stability of the generated states can then be examined and compared with our numerical results and stability analysis. For example, with an optical lattice frequency of ω L = 2π × 120 Hz (as in the example of Section III B), the time unit of our dynamical evolution plots is 1.3 ms. This implies that for the cases we have examined, signatures of instability would be typically visible on the experimentally feasible 10 to 100 ms timescale. We therefore believe that our predictions could be examined with current experimental techniques.
There are various directions along which one can extend the present considerations. A natural one is to extend the analysis to fully 3D condensates and examine the persistence and stability of higher-dimensional variants of the presented states. A perhaps more subtle direction is to consider a different basis of linear eigenfunctions in the 2D problem, namely instead of the Hermite-Gauss basis used here, to focus on the Laguerre-Gauss basis of the underlying linear problem with the parabolic potential. Under such a choice, it would be interesting to examine how solutions of that type, in-cluding one-node and multi-node ring-like structures (see, e.g., Ref. [17] and references therein), are deformed in the presence of the lattice and how their stability is correspondingly affected. Finally, as discussed above, it appears that the setting considered herein should be directly accessible to present experiments with pancake-shaped BECs. In view of that, it would be particularly relevant to examine which ones among the structures presented in this work can survive for evolution times that are of interest within the time scales of an experiment.