*2.2. Converter*

Due to intermittent problems with generating solar photovoltaic power, it is possible that the DC voltage produced is not large enough to support further research. It needs a DC-DC boost converter attached to it in order to increase to the appropriate level. The boost converter is modeled using typical calculations based on the system rating. An inverter witha2L voltage source performs the DC-AC conversion. The inverter's output is represented as follows:

$$P\_{inverter,out} = \eta\_{inverter} \times P\_{DC} \tag{6}$$

*ηinverter* = efficacy of the inverter;

*PDC* = power outcome from DC-DC boost converter;

*Pinverter*,*output* = Power delivered from the inverter.


*4.1.1 Reliability analysis*

In terms of reliability testing, the scale modified in this study adopted Cronbach's α to confirm the internal consistency of the scale. According to Sharma [41], Cronbach's α suggests that the construct has good reliability criteria if it is greater than .7. The Cronbach's α values for each construct of the modified scale in this study are shown in **Table 1**. The α values of self-efficacy, learning interest, learning anxiety, Learner's Perceived Value, and student engagement were all in the range of .8–.95, which meets the recommended criteria. Therefore, the measuring tools in this study had good internal consistency.


**4. Description of the numerical method**

Velocity, pressure, and temperature of the steam, as well as the temperature of the tube wall with the complex cross section, will be computed using the ANSYS/CFX v13.0 software [25].

The direct and inverse problems will be solved. In the inverse problem, the steam temperature at the superheater inlet and the heat transfer coefficient on the flue-gas side will be determined using measured steam temperatures at selected locations in the superheater. If the flue-gas side heat transfer coefficient is determined based on the measured steam temperature at the superheater outlet, the secant method is used. In this case, it is assumed that the inlet steam temperature is also measured. The Levenberg–Marquardt method is used to solve the least squares problem when the inlet steam temperature and the flue-gas side heat transfer coeffi‐ cient are determined based on the measured steam temperature at the selected points along the steam flow path. At every iteration step, a direct problem was solved using the ANSYS/CFX v13.0 software. Time-averaged three-dimensional equations of conservation of mass, momentum and energy are solved using an element-based finite volume method. The heat transfer in the flowing steam was modeled using equations of mass, momentum, and energy conservation which can be written as follows [25]:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{U}) = 0 \tag{19}$$

$$\frac{\partial(\rho \mathbf{U})}{\partial t} + \nabla \cdot (\rho \mathbf{U} \otimes \mathbf{U}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \mathbf{S}\_{\mathbf{M}} \tag{20}$$

$$\frac{\partial(\rho h)}{\partial t} - \frac{\partial p}{\partial t} + \nabla \cdot (\rho \,\mathrm{U}h) = \nabla \cdot (k \nabla T) + \mathbf{U} \cdot \nabla p + \tau : \nabla \mathbf{U} + \mathbf{S}\_{\mathrm{E}} \tag{21}$$

where the stress tensor *τ* is related to the strain rate by

$$
\tau = \mu (\nabla \mathbf{U} + (\nabla \mathbf{U})^T - \frac{2}{3} \delta \nabla \cdot \mathbf{U}) \tag{22}
$$

Temperature distribution in double omega tube walls (**Figure 18**) is governed by the heat conduction equation

$$\mathcal{L}\_w \rho\_w \frac{\partial T\_w}{\partial t} = \nabla \cdot \left[ k\_w(T\_w) \nabla T\_w \right] \tag{23}$$

Two equation turbulence model SST (Shear Stress Transport), which was developed based on the Wilcox *k*–*ω* and *k*–*ε* turbulence models [53], was used to model the flow turbulence. The baseline *k*–*ω* model is a two equations turbulence model that offers a good compromise between numerical effectiveness and computational accuracy. This model is a combination of the *k*–*ω* and *k*–*ε* models.

The Wilcox *k*–*ω* model is given by two equations:

**4. Description of the numerical method**

240 Numerical Simulation - From Brain Imaging to Turbulent Flows

energy conservation which can be written as follows [25]:

*t* r

where the stress tensor *τ* is related to the strain rate by

t m

*t t* r

¶ ¶

conduction equation

*t* r

( ) ( ) *p*

( ) ( ) () : *h p h kT p*

¶ ¶ - +Ñ× =Ñ× Ñ + ×Ñ + Ñ +

2 ( () ) 3 *T*

Temperature distribution in double omega tube walls (**Figure 18**) is governed by the heat

[() ] *w w w ww w T c kT T*

¶ =Ñ× Ñ

*t* r

d

= Ñ + Ñ - Ñ× **UU U**

r

r

¶

Velocity, pressure, and temperature of the steam, as well as the temperature of the tube wall with the complex cross section, will be computed using the ANSYS/CFX v13.0 software [25].

The direct and inverse problems will be solved. In the inverse problem, the steam temperature at the superheater inlet and the heat transfer coefficient on the flue-gas side will be determined using measured steam temperatures at selected locations in the superheater. If the flue-gas side heat transfer coefficient is determined based on the measured steam temperature at the superheater outlet, the secant method is used. In this case, it is assumed that the inlet steam temperature is also measured. The Levenberg–Marquardt method is used to solve the least squares problem when the inlet steam temperature and the flue-gas side heat transfer coeffi‐ cient are determined based on the measured steam temperature at the selected points along the steam flow path. At every iteration step, a direct problem was solved using the ANSYS/CFX v13.0 software. Time-averaged three-dimensional equations of conservation of mass, momentum and energy are solved using an element-based finite volume method. The heat transfer in the flowing steam was modeled using equations of mass, momentum, and

( )0

 t

**U U U S** (20)

 t

**U U USE** (21)

¶ (23)

**U** (19)

(22)

r

¶ +Ñ× =

¶ + Ñ × Ä = -Ñ + Ñ × + ¶ **M**

$$\frac{\partial \left(\rho k\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} \left(\rho U\_j k\right) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left(\mu + \frac{\mu\_t}{\sigma\_{k1}}\right) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + P\_k - \beta^\* \rho k o o \tag{24}$$

$$\frac{\partial \left(\rho oo\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{j}} \left(\rho U\_{j}\rho o\right) = \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \left(\mu + \frac{\mu\_{r}}{\sigma\_{a1}}\right) \frac{\partial oo}{\partial \mathbf{x}\_{j}} \right] + \alpha\_{1} \frac{oo}{k} P\_{k} - \beta\_{1} {}^{\*}\rho oo^{2} \tag{25}$$

To develop the baseline *k*–*ω* turbulence model, the transformed *k*–*ε* model is also used:

$$\frac{\partial \left(\rho k\right)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_{j}} \left(\rho U\_{j} k\right) = \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \left(\mu + \frac{\mu\_{i}}{\sigma\_{k2}}\right) \frac{\partial k}{\partial \mathbf{x}\_{j}} \right] + P\_{i} - \beta^{\*} \rho k o o \tag{26}$$

$$\frac{\partial(\rho o)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j}(\rho U\_j o) = $$

$$= \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_r}{\sigma\_{a2}} \right) \frac{\partial o}{\partial \mathbf{x}\_j} \right] + 2\rho \frac{1}{\sigma\_{a2}} \frac{\partial k}{\partial \mathbf{x}\_j} \frac{\partial o}{\partial \mathbf{x}\_j} + \alpha\_2 \frac{o}{k} P\_k - \beta\_2 \rho o^2 \tag{27}$$

Multiplying Eqs. (24) and (25) of the Wilcox model by weighting function *F*1, the transformed *k*–*ε* Eqs. (26) and (27) by a weighting function (1 – *F*1) and adding the transformed equations gives the baseline *k*–*ω* model:

$$\frac{\partial(\rho k)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\rho U\_j k) = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_r}{\sigma\_{\mathbb{k}3}} \right) \frac{\partial k}{\partial \mathbf{x}\_j} \right] + P\_k - \beta^\* \rho k o o + P\_{\mathbb{k}b} \tag{28}$$

$$\begin{aligned} \frac{\partial(\rho o)}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_j} (\rho U\_j o) &= \\ = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \mu + \frac{\mu\_t}{\sigma\_{o3}} \right) \frac{\partial o}{\partial \mathbf{x}\_j} \right] + \left( 1 - F\_i \right) 2 \rho \frac{1}{\sigma\_{o2}} \frac{\partial k}{\partial \mathbf{x}\_j} \frac{\partial o}{\partial \mathbf{x}\_j} + \alpha\_3 \frac{o}{k} P\_k - \beta\_3 \rho o^2 + P\_{ob} \end{aligned} \tag{29}$$

The baseline *k*–*ω* model has the advantages of the Wilcox *k*–*ω*, and the *k*–*ε* but is not able to predict onset and amount of flow separation from smooth surfaces [53]. The main reason for this drawback is neglecting of the transport of the turbulent shear stress. This causes an overprediction of the eddy-viscosity. In order to eliminate this deficiency, the shear stress transport (SST) turbulence model was developed [53], in which the turbulent viscosity *νt* is limited by the following equation:

$$\nu\_i = \frac{\alpha\_i k}{\max(\alpha\_i o, SF\_1)}\tag{30}$$

$$\begin{array}{c} \frac{\partial \ T\_w}{\partial n} = 0 \\\\ \frac{\partial T\_w}{\partial n} = 0 \\\\ \frac{\partial T\_w}{\partial n} \end{array}$$

$$\sigma\_w^{\text{h}\_\mathbf{g}}$$

**Figure 19.** One-fourth of a tube cross section.

Because of the symmetry, only one-fourth of a tube cross section was analyzed (**Figure 19**). Heat is transferred by convection from the combustion chamber to the tube through the surface E–F. The surfaces A–F and B–C–D–E are thermally insulated. The tube inner surface A–B is cooled by the flowing live steam. The following boundary conditions were assumed for the tube wall (**Figure 19**):

Solving Inverse Heat Transfer Problems When Using CFD Modeling http://dx.doi.org/10.5772/63807 243

$$\left.T\right|\_{A-B} = T\_w \left|\_{A-B}, \quad k \right. \frac{\left.\hat{\mathcal{O}}T}{\left.\hat{\mathcal{O}}\,\mathbf{n}}\right|\_{A-B} = -k\_w \left.\frac{\left.\hat{\mathcal{O}}T\_w}{\left.\hat{\mathcal{O}}\,\mathbf{n}}\right|\_{A-B}\right|\right. \tag{31}$$

$$k\_w \frac{\partial T\_w}{\partial n}|\_{B-C-D-E} = 0\tag{32}$$

$$k\_w \frac{\partial T\_w}{\partial n}|\_{E-F} = h\_g \left(T\_g - T\_w|\_{E-F}\right) \tag{33}$$

$$k\_{\boldsymbol{w}} \frac{\partial \boldsymbol{T}\_{\boldsymbol{w}}}{\partial \boldsymbol{n}} \mid\_{\mathcal{A}-\boldsymbol{F}} = \mathbf{0} \tag{34}$$

Also, the fluid temperature and the steam mass flow rate or steam velocity are known at the inlet of the tube. The steam pressure is given at the outlet of the tube. The direct conjugate heat transfer problem described by Eqs. (19) and (29) with appropriate boundary conditions was solved using the ANSYS/CFX software.


**4.2 Microenvironment immunomodulation**

MSCs have proven to regulate the immune response through cell-to-cell contact and by the secretion of soluble mediators including cytokines, prostaglandins, enzymes, and proapoptotic and antiapoptotic molecules [27, 37–39]. Different studies involving MSC transplantation in exacerbated immune response models such as peritonitis and ulcerative colitis ameliorate inflammation by reducing the expression levels of pro-inflammatory cytokines such as interleukin-1 beta, interleukin-12, interleukin-6, and tumor necrosis factor-α (TNFα). In addition, these cells exert a decrease of the classical phenotype M1 marker and an increase of the alternative phenotype M2, as well as a marked macrophage reprogramming from M1 to M2 [17, 40–42]. Moreover, MSCs can suppress T cell activation and proliferation by downregulating the expression of costimulatory molecules on the surface of dendritic cells [43], interleukin-10, transforming growth factor-B (TGFβ), nitric oxide, and indoleamine 2,3-dioxygenase enzyme production in response to inflammation as well as interleukin-2 absorption [37, 44, 45]. Similarly, B cell activation can be disturbed, and the regulatory B cell phenotype can be promoted [46, 47].

#### **4.3 Tissue repairing and regeneration**

MSCs participate in repairing many tissues, mostly by the secretion of TGFβ and vascular endothelial growth factor (VEGF) to promote angiogenesis [48], extracellular matrix remodeling, and reduction of the scar formation in chronic wounds [49, 50]. Similarly, in pathologies where gliosis, demyelination, and neuroinflammation occur, MSCs have shown neuroprotective activities such as vascular stabilization and angiogenesis by tight junction protein expression [51], neuronal suppression apoptosis [52, 53], glia hypertrophy prevention [54], and promotion of myelinization by the activation of oligodendrocyte precursor cells [53, 55]. Additionally, MSCs promote synaptic transmission [56], neurite outgrowth, and axonal sprouting mostly by excretion of trophic factors including brain-derived neurotrophic factor (BDNF) and nerve growth factor (NGF) [12, 57]. SCI models motor skills are increased [58] inclusively; bladder and erectile dysfunction improvement have been reported [59, 60].

#### **4.4 Multipotential capacity**

Originally it was believed that MSCs could exclusively differentiate into cells from the mesodermal lineage [23, 24]; however, in the last 20 years many authors have proven that a proper microenvironment can promote greater plasticity and that MSCs from different sources can differentiate into dermal, neural, or glial cells in vitro and in vivo when they are exposed to neurotrophic factors and specific cytokines [61–63]. Smooth muscle and endothelial cells derived from MSCs can be detected and improve heart functions in ischemic myocardium models [64]. Also, skin, articular cartilage, and bone regeneration have been reported, but mostly when MSCs are combined with natural and artificial scaffolds or when genetically modified [65–67]. In order to achieve CNS regeneration, different sources of MSCs and culture methods have been tested in murine SCI models; however, transplants have demonstrated to improve functional recovery by differentiation into neurons, astrocytes, but mostly oligodendrocytes [54, 55, 58].

### **4.5 Exosomes as mechanisms for cell-to-cell communication and drug delivery vehicles**

Exosomes are extracellular vesicles released by many cells, including MSCs. Their length is between 30 and 100nm, and they can bind to cells through receptor and ligand interaction or by fusion with the target cell membrane to deliver high amounts of cytokines, growth factors, microRNAs, and mRNAs capable of modifying peptide and protein synthesis [60]. Thus, the derived MSC exosomes could be the most attractive therapy in SCI models since neuronal differentiation from MSCs remains poor [58, 68] and several studies show that the regenerative and anti-inflammatory mechanisms are mostly mediated by paracrine factors [27, 36, 69–71]. Furthermore, MSC exosomes are natural drug delivery vehicles that can be modified and produced in high quantities [28, 72, 73]. MSCs have proven to be safe in many different preclinical studies; however, clinical trials involving exosomes in SCI therapy are not yet recruited due to the fact that optimal MSC culture conditions and protocols for exosome isolation is still to be established (www.clinicaltrials.gov).


**Modulational Instability, Inter-Component Asymmetry, and Formation of Quantum Droplets in One-Dimensional Binary Bose Gases**

#### **Thudiyangal Mithun 1,2,\* ,†, Aleksandra Maluckov 3,\* ,†, Kenichi Kasamatsu 4,† , Boris A. Malomed 5,† and Avinash Khare 6,†**

Received: 9 December 2019; Accepted: 14 January 2020; Published: 18 January 2020

**Abstract:** Quantum droplets are ultradilute liquid states that emerge from the competitive interplay of two Hamiltonian terms, the mean-field energy and beyond-mean-field correction, in a weakly interacting binary Bose gas. We relate the formation of droplets in symmetric and asymmetric two-component one-dimensional boson systems to the modulational instability of a spatially uniform state driven by the beyond-mean-field term. Asymmetry between the components may be caused by their unequal populations or unequal intra-component interaction strengths. Stability of both symmetric and asymmetric droplets is investigated. Robustness of the symmetric solutions against symmetry-breaking perturbations is confirmed.

**Keywords:** quantum droplet; binary Bose–Einstein condensate; modulational instability

#### **1. Introduction**

The mean-field (MF) theory of weakly interacting dilute atomic gases rules out formation of a liquid state [1,2]. However, it has been recently shown that a liquid phase arises if one takes into account beyond-MF effects originating from quantum fluctuations around the MF ground state of weakly interacting binary (two-component) Bose gases [3]. A fundamental property that allows one to interpret this phase as a fluid is incompressibility: It maintains a limit density which cannot be made larger (see details below), hence adding more atoms leads to spatial expansion of the state. Another fundamental feature of this quantum-fluid phase is that it facilitates self-trapping of *quantum droplets* (QDs), which are stabilized by the interplay between the contact MF interaction and the beyond-MF Lee–Huang–Yang (LHY) correction [4]. Binary Bose–Einstein condensates (BECs) with competing intra- and inter-component MF interactions of opposite signs offer a remarkable possibility for the generation of QDs, as proposed by Petrov [3]. This possibility was further elaborated in various settings, including different effective dimensions [5–20]. In particular, the dynamics of QDs with the flat-top (FT) or Gaussian shape, which correspond to large or relatively small numbers of particles, respectively, was addressed in the framework of the one-dimensional (1D) reduction of the model [20]. The theoretical prediction was followed by experimental creation of QDs in mixtures of two different atomic states of 39K, with quasi-2D [21,22] and fully 3D [23,24] shapes (see also recent reviews [25,26]). Very recently, the creation of especially long-lived QDs was reported in a heteronuclear 41K-87Rb

system [27]. Another theoretically predicted and experimentally realized option for the creation of QDs makes use of the single-component condensate with dipole–dipole interactions [28–35]. It is relevant to mention that the formation of multiple droplets was also predicted and experimentally observed as an MF effect in strongly nonequilibrium (turbulent) states of BECs [36].

Collective modes of QDs are a subject of special interest, as they reveal internal dynamics of the droplets [20,24,32,37,38]. In particular, the stable existence of the QDs is secured if the particle-emission threshold lies below all excitation modes, hence a perturbation in the form of such modes will not cause decay of the droplet.

We aim to address issues that are related to the creation of QDs in the 1D setting and were not addressed in previous works. First, we consider modulational instability (MI) of spatially uniform plane-wave (PW) states, in the framework of the coupled system of Gross–Pitaevskii (GP) equations with the LHY corrections, for the two-component MF wave function of the binary condensate. This is the system which was originally derived in [5]. Recently, MI has been experimentally demonstrated in BECs with attractive interactions [39–41]. Other examples of the MI are provided by the binary BEC with the linear Rabi coupling or the spin-orbit coupling [41,42], and by a system combining the MF and LHY terms [43]. The linear-stability analysis, followed by direct simulations of the corresponding GP equations, shows that the lower branch of the PW states exhibits MI, the instability splitting the PW into a chain of localized droplet-like structures. Secondly, we address properties of the QDs in the binary condensate in the framework of the two-component GP system, without assuming effective inter-component symmetry, which reduces the system to a single-component GP equation. The asymmetry implies different MF self-repulsion coefficients in the two components, and/or unequal norms in them. Although properties of QDs have been studied by using the two-component GP system in some papers [6,11,14,15,17,18], the explicit asymmetry of the system parameters has not been addressed, except for [14] in which the situation for 39K-39K and 23Na-87Rb atomic mixtures have been considered. We conclude that the population difference between the components does not significantly affect density profiles of QDs in the system with equal MF self-repulsion strengths in the two components. On the other hand, we find that profiles of the QD solutions are essentially asymmetric when the self-repulsion coefficients are different in the components. Generally, the numerical findings corroborate stability of the known symmetric states against symmetry-breaking perturbations. We also address the MI of the two-component system, and demonstrate that chains of asymmetric QDs can be generated by the MI-induced nonlinear evolution.

The paper is organized as follows. In Section 2 we introduce the model and discuss conditions necessary for the formation of the droplets. Section 3.1 deals with the single-component version of the symmetric system. We consider various solutions admitted by it (PW, FT, periodic, etc.), and apply the linear-stability analysis of the PW solution to assess the MI, in a combination with direct simulations. In Section 3.2, we address the stability of asymmetric droplets, as well as the formation of droplets in the two-component asymmetric system via the MI. The paper is concluded by Section 4. Additional symmetric and asymmetric exact and approximate analytical solutions are presented in Appendices.

#### **2. Model and Methods**

We consider the 1D model of the two-component condensate with coefficients of the intra-component repulsion, *g*1 > 0 and *g*2 > 0, and inter-component attraction, *g*12 < 0. In the weak-interaction limit, the corresponding energy density, which includes the MF terms and LHY correction, was derived in [5]:

$$\mathcal{E}\_{\rm ID} = \frac{\left(\sqrt{\mathcal{g}\_1}\rho\_1 - \sqrt{\mathcal{g}\_2}\rho\_2\right)^2}{2} + \frac{\mathcal{g}\delta\mathfrak{g}\left(\sqrt{\mathcal{g}\_2}\rho\_1 + \sqrt{\mathcal{g}\_1}\rho\_2\right)^2}{\left(\mathcal{g}\_1 + \mathcal{g}\_2\right)^2} - \frac{2\sqrt{\mathcal{m}}\left(\mathcal{g}\_1\rho\_1 + \mathcal{g}\_2\rho\_2\right)^{3/2}}{3\pi\hbar},\tag{1}$$

where *m* is the atomic mass (the same for both components), *ρj* = |Ψ*j* | 2 (*j* = 1, 2) is the density of the *j*-th component, represented by the MF wave function Ψ*j* , and

$$\mathbf{g} \equiv \sqrt{\mathbf{g}\_1 \mathbf{g}\_2} \qquad \delta \mathbf{g} \equiv \mathbf{g}\_{12} + \mathbf{g}. \tag{2}$$

The last term in Equation (1) represents the LHY correction. Derivation of Equation (1) assumes that the binary BEC is close to the point of the MF repulsion-attraction balance, with |*δg*| *g* . In experiments, *δg* may be tuned to be both positive and negative [21–23].

Equation (1) is valid in the case of tight confinement applied in the transverse dimensions, which makes the setting effectively one-dimensional. In the 3D case, the LHY term ∼ −*ρ* 3/2 (for *ρ*1 = *ρ*2 ≡ *ρ*) is replaced by one ∼ +*ρ* 5/2. A detailed consideration of the crossover from 3D to 1D [12,44,45] in the two-component system is a problem which may be a subject of a separate work. Here, it is relevant to compare the symmetric version of Equation (1) for the energy density with that recently presented in [12]. It demonstrates that an accurately derived LHY contribution to the energy density of the 1D system contains, in addition to the *ρ* 3/2 term which was derived in [5], a term ∼ *ρ* 2 , which can be absorbed into the mean-field energy density, and a higher-order term ∼ *ρ* 3 , which was omitted in the analysis reported in [12]. A conclusion formulated in that work is that the energy density originally derived in [5] is literally valid if the ratio of the mean-field energy to that of the transverse confinement takes values ≤ 0.03. For typical experimental parameters, this implies that the difference between absolute values of scattering lengths of the mean-field intra-component repulsion and inter-component attraction should be ≤ 1 nm, which may be achieved in the experiment. The 1D QDs originate from the balance of the second term in Equation (1), corresponding to the weakly repulsive MF interaction, with *δg* > 0, and the LHY term, which introduces effective attraction in the 1D setting, on the contrary to the repulsion in the 3D setting [5,20].

The energy functional, R +∞ −∞ E1D*dZ*, gives rise to the system of GP equations, which include the LHY correction,

$$\begin{split} i\hbar \frac{\partial \Psi\_{1}}{\partial T} &= -\frac{\hbar^{2}}{2m} \frac{\partial^{2} \Psi\_{1}}{\partial Z^{2}} + (\mathcal{g}\_{1} + \mathbb{G}\mathcal{g}\_{2})|\Psi\_{1}|^{2}\Psi\_{1} - (1 - \mathbb{G})\mathcal{g}|\Psi\_{2}|^{2}\Psi\_{1} - \frac{\mathcal{g}\_{1}\sqrt{m}}{\pi\hbar} \sqrt{\mathcal{g}\_{1}|\Psi\_{1}|^{2} + \mathcal{g}\_{2}|\Psi\_{2}|^{2}\Psi\_{1}} \\ i\hbar \frac{\partial \Psi\_{2}}{\partial T} &= -\frac{\hbar^{2}}{2m} \frac{\partial^{2} \Psi\_{2}}{\partial Z^{2}} + (\mathcal{g}\_{2} + \mathbb{G}\mathcal{g}\_{1})|\Psi\_{2}|^{2}\Psi\_{2} - (1 - \mathbb{G})\mathcal{g}|\Psi\_{1}|^{2}\Psi\_{2} - \frac{\mathcal{g}\_{2}\sqrt{m}}{\pi\hbar} \sqrt{\mathcal{g}\_{1}|\Psi\_{1}|^{2} + \mathcal{g}\_{2}|\Psi\_{2}|^{2}\Psi\_{2}} \end{split} (3)$$

where *T* and *Z* are the time and coordinate measured in physical units, and parameter

$$\mathbf{G} = \frac{\mathbf{2g\delta g}}{(\mathbf{g\_1} + \mathbf{g\_2})^2} \tag{4}$$

measures the deviation from the MF repulsion–attraction balance point, see Equation (2). The normalization of the components of the wave function is determined by numbers of bosons in each component:

$$N\_{\hat{j}} = \int\_{-\infty}^{+\infty} |\Psi\_{\hat{j}}|^2 dZ. \tag{5}$$

Further, rescaling

$$\left(\frac{\text{mg}^2}{\hbar^3}\right)T \equiv t, \quad \left(\frac{\text{mg}}{\hbar^2}\right)Z \equiv z, \quad \left(\frac{\hbar}{\sqrt{\text{mg}}}\right)\Psi\_{1,2} \equiv \psi\_{1,2} \tag{6}$$

casts Equation (3) in the normalized form,

$$\begin{aligned} i\frac{\partial\psi\_1}{\partial t} &= -\frac{1}{2}\frac{\partial^2\psi\_1}{\partial z^2} + (P + GP^{-1})|\psi\_1|^2\psi\_1 - (1 - G)|\psi\_2|^2\psi\_1 - \frac{P}{\pi}\sqrt{P|\psi\_1|^2 + P^{-1}|\psi\_2|^2}\psi\_1, \\ i\frac{\partial\psi\_2}{\partial t} &= -\frac{1}{2}\frac{\partial^2\psi\_2}{\partial z^2} + (P^{-1} + GP)|\psi\_2|^2\psi\_2 - (1 - G)|\psi\_1|^2\psi\_2 - \frac{1}{\pi P}\sqrt{P^{-1}|\psi\_2|^2 + P|\psi\_1|^2}\psi\_2. \end{aligned} \tag{7}$$

where parameter

$$P \equiv \sqrt{\frac{\mathcal{G}\_1}{\mathcal{G}\_2}} = \frac{\mathcal{G}\_1}{\mathcal{g}} \tag{8}$$

determines the asymmetry of the system, in the case of *P* 6= 1. Note that, as concerns stationary solutions with chemical potentials *µ*1,2, sought for as

$$
\psi\_{1,2}(z,t) = \exp(-i\mu\_{1,2}t)\phi\_{1,2}(z),
\tag{9}
$$

states with mutually proportional components, *φ*1(*z*) = *Kφ*2(*z*), are only possible in the fully symmetric case with *P* = 1, *µ*1 = *µ*2, and *K* = 1. In previous works [5,20], 1D solutions for QDs were considered only in the framework of the single GP equation which corresponds to symmetric system Equation (7) with *P* = 1 and *ψ*1 = *ψ*2.

#### **3. Modulation Instability Versus QDs**

In this section, we address MI of PWs in both symmetric and asymmetric GP systems, and relate it to formation of the QDs in the binary bosonic gas. To the best of our knowledge, this is the first work aiming to associate the MI with the formation of the 1D droplets in the system with unequal components. We first consider MI in the framework of the single-component reduction of the symmetric version of Equation (7), after briefly reviewing stationary solutions of the GP equation. Next, we extend the analysis for the two-component GP system, which makes it possible to produce asymmetric QDs, starting from the MI of asymmetric PW states.

#### *3.1. The Single-Component GP Model*

Under the single-component reduction of the binary system, with *g*1 = *g*2 ≡ *g* and *ψ*1 = *ψ*2 ≡ *ψ*, Equation (1) simplifies to [5]

$$
\varepsilon\_{\rm{1D}} \equiv \frac{\hbar^4}{m^2 \mathcal{g}^3} \mathcal{E}\_{\rm{1D}} = \frac{\delta \mathcal{g}}{\mathcal{g}} n^2 - \frac{2^{5/2}}{3\pi} n^{3/2},\tag{10}
$$

with the single dimensionless density, *n* = |*ψ*| 2 ≡ *h*¯ 2/*mg ρ*. Assuming a spatially uniform state, the equilibrium density and the corresponding chemical potential are given by

$$m\_0 = \frac{8}{9\pi^2} \left(\frac{\mathcal{g}}{\delta \mathcal{g}}\right)^2, \; \mu\_0 = -\frac{4}{9\pi^2} \frac{\mathcal{g}}{\delta \mathcal{g}}.\tag{11}$$

Density *n*0 corresponds to the minimum of the energy per particle, *∂n* - *n* −1 *ε*1D(*n*) = 0, and *µ*0 is negative for *δg*/*g* > 0. The corresponding single GP equation is

$$i\frac{\partial\psi}{\partial t} = -\frac{1}{2}\frac{\partial^2\psi}{\partial z^2} + \frac{\delta g}{g}|\psi|^2\psi - \frac{\sqrt{2}}{\pi}|\psi|\psi\tag{12}$$

with normalization condition R +∞ −∞ |*ψ*(*z*)| 2*dz* = *N*, where *N* ≡ *N*1 = *N*2 is the number of atoms in each component.

Although coefficient *δg*/*g* can be scaled out in Equation (12), as done in [20], we keep it here as a free parameter. This option is convenient for the subsequent consideration of the MI, treating *δg*/*g* and density *n* as independent constants, which may be matched to experimentally relevant parameters.

Below, we address two stationary solutions of Equation (12). One is the QD bound state of a finite size, which was studied in detail in [5] and [20]. The other solution is the PW with uniform density. Here, we briefly recapitulated basic properties of these solutions for the completeness of the presentation. In Section 3.1.3 we address the MI of the PWs and associate it with the spontaneous generation of chains of localized modes. Additional families of exact analytical solutions of Equation (12) are given in Appendix A.

#### 3.1.1. The Droplet Solution

As shown in [5,20,46], at *δg*/*g* > 0 Equation (12) gives rise to an exact soliton-like solution representing a QD, maintained by the balance between the effective cubic self-repulsion and quadratic attraction:

$$\psi(z,t) = \frac{Ae^{-i\mu t}}{1 + B\cosh(\sqrt{-2\mu}z)}, \quad A = \sqrt{n\_0}\frac{\mu}{\mu\_0}, \quad B = \sqrt{1 - \frac{\mu}{\mu\_0}}.\tag{13}$$

This solution exists in a finite range of negative values of the chemical potential *µ*0 < *µ* < 0, featuring the FT shape at 0 < *µ* − *µ*0 |*µ*0|, with size *L* ≈ (−2*µ*0) −1/2 ln h (1 − *µ*/*µ*0) −1 i [5,20]. A typical density profile of the FT solution is displayed in the inset of Figure 1. At *µ* = *µ*0, the size of the droplet diverges, and the solution carries over into the delocalized PW with uniform density, *n* = *n*0. The fact that the density of the condensate filling the FT state cannot exceed the largest value, *n*0, implies its incompressibility. For this reason, the condensate may be considered as a fluid, as mentioned above. With the increase of *µ* from *µ*0 towards *µ* = 0, the maximum density of the localized mode,

$$n\_{\text{max}} \equiv n(z=0) = n\_0 \left(\frac{\mu}{\mu\_0}\right)^2 \left(1 + \sqrt{1 - \frac{\mu}{\mu\_0}}\right)^{-2},\tag{14}$$

monotonously decreases from *n*0 to 0. The QD's FWHM size, defined by condition *n* (*z* = *L*FWHM/2) = *n* (*z* = 0) /2, also shrinks at first with increasing *µ*, attaining a minimum value (*L*FWHM)min ≈ 2.36/ √ −*µ*0 at *µ*/*µ*0 ≈ 0.776. Further increase of *µ* towards *µ* = 0 makes the QD broader, its width diverging as *L*FWHM ≈ 1.71/√ −*µ* at *µ* → −0.

**Figure 1.** The maximum density *n*max ≡ *n*(*z* = 0) in the FT (flat-top) state, as per Equation (14), and the PW (plane-wave)/density are displayed as functions of *µ* by the red-solid and blue-dashed curves, respectively, for *δg*/*g* = 0.05. In this case, Equation (11) yields *n*0 = 36.025 and *µ*0 = −0.900633. The PW solution includes upper and lower branches corresponding to *n* ±, as given by Equation (20), the lower one (marked by circles) being subject to the MI (modulational instability). The spinodal point is one with coordinates (*µc*, *nc*). For other values of *δg*/*g*, the plot can be generated from the present one by rescaling. The inset shows the density profile of the FT solution for *δg*/*g* = 0.05 and *µ* = *µ*0 + 0.00001, very close to the delocalization limit (the transition to PW).

The norm of the exact QD solutions given by Equation (13) is

$$N(\mu) = n\_0 \sqrt{-\frac{2}{\mu\_0}} \left[ \ln \left( \frac{1 + \sqrt{\mu/\mu\_0}}{\sqrt{1 - \mu/\mu\_0}} \right) - \sqrt{\frac{\mu}{\mu\_0}} \right]. \tag{15}$$

It satisfies the well-known Vakhitov–Kolokolov (VK) necessary stability criterion [47],

$$\frac{dN(\mu)}{d\mu} = -\frac{n\_0}{\mu\_0^2} \sqrt{-\frac{\mu}{2}} \frac{1}{1 - \mu/\mu\_0} < 0,\tag{16}$$

due to *µ*0 < 0 and

$$0 < \mu/\mu\_0 < 1.\tag{17}$$

Full stability of the QD family has been verified by direct simulations of the evolution of perturbed QDs in the framework of Equation (12).

It is relevant to mention that the exact solution of Equation (13) is valid too at *δg*/*g* < 0, when the cubic term in Equation (12) is self-attractive, like the quadratic one. In that case, *µ*0 is positive, as per Equation (11), while the chemical potential of the self-trapped state remains negative, as the solution of Equation (13) may exist only at *µ* < 0. Then, it follows from Equation (13) that the soliton-like mode exists for all values of *µ* < 0 (unlike the finite interval Equation (17), in which the solution exists for *δg*/*g* > 0), and it does not feature the FT shape. Rather, with the increase of −*µ*, it demonstrates a crossover between the KdV-soliton shape ∼ sech2 p −*µ*/2*z* and the nonlinear-Schrödinger one, ∼ sech p −2*µz* . For *δg*/*g* < 0, the *N*(*µ*) dependence for the soliton family carries over into the following form,

$$N(\mu)\Big|\_{\delta\mathcal{G}<0} = n\_0 \sqrt{\frac{2}{\mu\_0}} \left[ \sqrt{-\frac{\mu}{\mu\_0}} - \arctan\left(\sqrt{-\frac{\mu}{\mu\_0}}\right) \right],\tag{18}$$

which is an analytical continuation of Equation (15). This dependence also satisfies the VK criterion.

#### 3.1.2. The Plane-Wave Solution

The PW solution of Equation (12 )can be presented in a form *ψ*(*z*, *t*) = √ *n* exp (*iK*PW*z* − *iµt*) with wavenumber *K*PW and constant density *n*, which determine the corresponding chemical potential:

$$
\mu\_{\rm PW} = \frac{\delta \mathcal{g}}{\mathcal{g}} n - \frac{\sqrt{2}}{\pi} \sqrt{n} + \frac{1}{2} K\_{\rm PW}^2. \tag{19}
$$

The Galilean invariance of Equation (12) implies that any quiescent solution *ψ*0 (*z*, *t*) generates a family of moving ones, with arbitrary velocity *c*. Therefore, *K*PW may be canceled by means of transformation *ψc* (*z*, *t*) = exp *icz* − *ic*2 *t*/2 *ψ*0 (*z* − *ct*, *t*) with *c* = −*K*PW.

For given *µ*, Equation (19) produces two different branches of the density as a function of *µ* (here, *K*PW = 0 is set):

$$\sqrt{n^{\pm}(\mu)} = \frac{1}{\sqrt{2}\pi} \frac{\mathcal{g}}{\delta \mathcal{g}} \pm \sqrt{\frac{1}{2\pi^{2}} \left(\frac{\mathcal{g}}{\delta \mathcal{g}}\right)^{2} + \frac{\mathcal{g}}{\delta \mathcal{g}}\mu}. \tag{20}$$

For *δg*/*g* = 0.05, these branches are shown in Figure 1. As follows from Equation (20), they exist (for *δg*/*g* > 0) above a minimum value of *µ*: *µc* = −(2*π* 2 *δg*/*g*) −1 = (9/8)*µ*0, the respective density being

$$n\_{\varepsilon} = n^{\pm}(\mu\_{\varepsilon}) = \frac{1}{2\pi^2} \left(\frac{\text{g}}{\delta g}\right)^2 = \frac{9}{16} n\_0. \tag{21}$$

Values *µ* = *µc* and *n* = *nc* correspond to the spinodal point [5], and *n* +(*µ*0) = *n*0 (see Equation (13)). Note that the above-mentioned existence region of the soliton solution in terms of the chemical

potential, *µ*0 < *µ* < 0, lies completely inside that of the PW state, which is *µc* ≤ *µ*. Thus, the soliton always coexists with the PW (this fact is also obvious in Figure 1).

#### 3.1.3. Modulational Instability of the Plane Waves

Here, we aim to analyze the MI of PW solutions in the framework of the single-component GP Equation (12) and demonstrate how the development of the MI can help to generate QDs. We perform the analysis for the PWs with zero wavenumber *K*PW = 0, which is sufficient due to the aforementioned Galilean invariance of the underlying equation.

A small perturbation is added to the stationary PW state as

$$\psi(z,t) = \left[\sqrt{n} + \delta\psi(z,t)\right] \exp\left(-i\mu t\right). \tag{22}$$

The substitution of this expression in Equation (12) and linearization with respect to perturbation *δψ* leads to the corresponding Bogoliubov–de Gennes equation,

$$i\frac{\partial}{\partial t}\delta\psi = -\frac{1}{2}\frac{\partial^2}{\partial z^2}\delta\psi + \frac{\delta g}{g}n(\delta\psi + \delta\psi^\*) - \frac{\sqrt{n}}{\sqrt{2}\pi}(\delta\psi + \delta\psi^\*). \tag{23}$$

By looking for perturbation eigenmodes with wavenumber *k* and frequency Ω,

$$
\delta\psi = \zeta\cos(kz - \Omega t) + i\eta\sin(kz - \Omega t),
\tag{24}
$$

and real infinitesimal amplitudes *ζ* and *η*, Equation (23) yields a dispersion relation for the eigenfrequencies:

$$
\Omega^2 = \frac{k^4}{4} + \left(\frac{\delta g}{g}n - \frac{\sqrt{n}}{\sqrt{2}\pi}\right)k^2. \tag{25}
$$

The MI takes place when Ω acquires an imaginary part. As follows from Equation (25), this occurs when the density satisfies condition *n* < [2*π* 2 (*δg*/*g*) 2 ] −1 = *nc* see (Equation (21)), which corresponds to branch *n* − of the PW state. The instability region in terms of *k* is given by

$$k^2 < 4\left(\frac{\sqrt{n}}{\sqrt{2}\pi} - \frac{\delta\mathfrak{g}}{\mathfrak{g}}n\right) \equiv k\_0^2. \tag{26}$$

The MI gain *σ* ≡ |ImΩ| is plotted in Figure 2 versus |*k*| and *δg*/*g*, for given density *n* = 40 in panel (a), and versus |*k*| and *n*, for given *δg*/*g* = 0.05 in (b). It is easy to find from Equation (25) that the largest gain is attained at wavenumber

$$k\_{\text{max}} = \frac{k\_0}{\sqrt{2}},\tag{27}$$

with *k*0 defined as per Equation (26). Note that Figure 2a includes the case of the self-attractive cubic nonlinearity, with *δg*/*g* < 0, which naturally displays much stronger MI, as in this case it is driven by both the quadratic and cubic nonlinear terms. In fact, the extension of the MI chart to *δg*/*g* < 0 makes it possible to compare the MI in the present system and its well-known counterpart in the setting with the fully attractive nonlinearity.

Comparing parameter values at which the QD solutions are predicted to appear, and the MI condition for the PW with the corresponding density, the MI is expected to provide a mechanism for the creation of the QDs. This is confirmed by direct simulations of the GP Equation (12), as shown in Figure 3. The PW with *n* = 10 is taken as the input, so that it is subject to the MI for *δg*/*g* = 0.05, as seen in Figure 2b. As shown in Figure 3, small initial perturbations trigger the emergence of multiple-QD patterns (chains) at *t* ≥ 100. For these parameters, we get *k*max = 0.6508 and *σ* (*k*max) = 0.2118, which determines the wavelength of the fastest growing modulation, *λ* = 2*π*/*k*max ≈ 9.66, and the growth-time scale, *τ* = 2*π*/*σ* (*k*max) ≈ 30. The number of the generated droplets in Figure 3 is

consistent with estimate *L*/*λ* ' 10, where *L* = 100 is the size of the simulation domain. We have checked that the number of generated droplets is approximately given by *L*/*λ* for other values of parameters as well. This dynamical scenario is similar to those observed in other models in the course of the formation of soliton chains by MI of PWs [39,40]. The long-time evolution in Figure 3a shows that the number of the droplets becomes smaller due to merger of colliding droplets into a single one, which agrees with dynamical properties of 1D QDs reported in [20].

**Figure 2.** Color-coded values of the MI gain, *σ* = Im(Ω), are displayed for fixed *n* = 40 in (**a**), and for fixed *δg*/*g* = 0.05 in (**b**). Note that panel (**a**) covers both signs of the cubic nonlinearity, *δg* > 0 and *δg* < 0. Solid and dashed white curves represent the MI boundary (Equation (26)) and the peak value of the MI gain (Equation (27)), respectively.

To implement this mechanism of the generation of a chain of solitons in the experiment, i.e., make the density smaller than the critical value *nc*, one may either apply interaction quench (by means of the Feshbach resonance), suddenly decreasing *δg*/*g*, as was done in recent experimental works for different purposes [21–23,48]. Another option, which is specific to the 1D setting, is sudden decrease of density *n* by relaxing the transverse trapping.

#### *3.2. The Two-Component Gross–Pitaevskii Model*

In this section, we revert to the full two-component GP system Equation (7), aiming to explore the formation of QD states in it. The two-component setting may include parameter imbalance between the two components, as indicated theoretically [3] and observed experimentally [21–23,27]. Here, we present the analysis of asymmetric QDs in two cases: (i) the two-component GP system with different populations, *N*1/*N*2 6= 1, and equal intra-component coupling strength, *g*1 = *g*2 (i.e., *P* = 1, see Equation (8)), and (ii) the system with different intra-component coupling strengths, *g*1 6= *g*2 (i.e., *P* 6= 1). These options suggest a possibility to check the stability of the solutions of the symmetric system, reduced to the single-component form, against symmetry-breaking perturbations. That objective is relevant because, in the real experiment, scattering lengths of the self-interaction in the two components are never exactly equal [21–24]. We address, first, an asymmetric single-droplet solution, and, subsequently, MI of the PW states in the two-component system.

Because, as said above, solutions with mutually proportional components (written as *φ*1 = *Kφ*2) are possible solely in the strictly symmetric setting, asymmetric QDs cannot be found in an exact analytical form. As shown in Appendix B (see Equations (A6)–(A12)), asymptotic analytical solutions can be obtained for strongly asymmetric states, with one equation replaced by its linearized version. In this section, we chiefly rely on numerical solution of Equation (7).

**Figure 3.** A typical example of the MI development, starting from an unstable PW state, with density *n* = 10 and *δg*/*g* = 0.05, which is subject to the MI, pursuant to Figure 2. In (**a**), the spatiotemporal pattern of the evolution of the condensate density is shown. In the right-hand panels, cross sections of the density profiles are displayed at *t* = 100 (**b**), *t* = 120 (**c**), and *t* = 140 (**d**). The simulations were performed in domain −50 < *z* < +50 with 2500 grid points and periodic boundary conditions.

3.2.1. Asymmetric QDs with Unequal Populations (*N*1 6= *N*2) for *g*1 = *g*2 (*P* = 1)

In the system with *P* = 1 (see Equation (8)), we calculated the droplet states as stationary solutions of Equation (7) by means of the imaginary-time-evolution method with the Neumann's boundary conditions, under the constraint that the norm is fixed in the first component, R +∞ −∞ *dz*|*ψ*1(*z*)| 2 = *N*1, while chemical potential *µ*2 is fixed in the other one, allowing its norm *N*2 to vary.

Figure 4 displays essential features of weakly asymmetric droplets for *δg*/*g* = 0.05 and fixed *N*1 = 100. The symmetric (completely overlapping) solution with *N*1 = *N*2 is found at *µ*1 = *µ*2 = −0.88878. When *µ*2 deviates from this value, profiles of the two components become slightly different, as shown in Figure 4a. The profiles of the droplet solution hardly change for different values of *µ*2, but Figure 4b demonstrates that, at *µ*2 → −0, *ψ*2 develops small-amplitude extended tails, which are absent in *ψ*1. Due to the contribution of the tails, the approach of *µ*2 < 0 towards zero leads to the increase of norm *N*2, as seen in Figure 4c. Note that the growth of *N*2(*µ*2) at *µ*2 → −0 is opposite to the decay of the QD's norm in the single-component model at *µ* → −0, cf. Equation (15). At *µ*2 ≥ 0, the *ψ*2

component undergoes delocalization, with its tails developing a nonzero background at |*z*| → ∞, as seen in the density profile displayed in Figure 4b for *µ*2 = 0, and norm *N*2(*µ*2) diverging at *µ*2 → −0 in Figure 4c.

**Figure 4.** (**a**) Stationary weakly asymmetric (with respect to the two components) solutions of Equation (7), obtained for *µ*2 = −0.4 with fixed *N*1 = 100. Dashed and solid curves display density profiles of the first (*n*1 ) and second (*n*2) components, respectively. (**b**) The semi-log plot of the density profiles of *n*2 for *µ*2 = −0.4, −0.04, and 0 at *z* > 0. (**c**) Dependences of *N*2 (black dots: the left vertical axis) and asymmetry parameter *δ*21, defined as per Equation (28) (the red dashed line pertaining to the right vertical axis), on *µ*2 for fixed *N*1 = 100. The parameters are *P* = 1 (*g*1 = *g*2) and *δg*/*g* = 0.05. The symmetric point with *N*1 = *N*2 = 100 and *δ*21 = 0 corresponds to *µ*1 = *µ*2 = −0.88878.

In Figure 4c, we also plot the parameter of the asymmetry between the two components, defined as

$$\delta\_{21} = \frac{n\_2(z=0) - n\_1(z=0)}{n\_2(z=0) + n\_1(z=0)}. \tag{28}$$

It increases almost linearly with *µ*2, although its absolute value does not exceed 0.02. Thus, the droplet tends to keep a nearly symmetric profile, with respect to the two components, in the symmetric system, even if the population imbalance is admitted. In fact, this circumstance makes the analysis self-consistent, as the use of the GP system with the LHY correction implies that the MF intra- and inter-component interactions nearly cancel each other, which is possible only if shapes of the two components are nearly identical.

#### 3.2.2. Asymmetric QDs in the System with *P* 6= 1 (*g*1 6= *g*2)

Next, we consider the QDs for *P* 6= 1, setting *P* > 1 without loss of generality. Then, the MF energy is minimized for *n*2 > *n*1; the situation with *n*1 > *n*2 can be considered too, replacing *P* by *P* −1 .

Following the procedure similar to that employed in Section 3.2.1, we produce QD solutions for *δg*/*g* = 0.05, *N*1 = 100, and several different values of *P*, varying *µ*2. In Figure 5a, we plot density profiles for three different values of *P*. Naturally, the difference of the two components increases with the increase of *P*. In Figure 5b we display *N*2 and parameter *δ*21 (see Equation (28)) of the asymmetric QDs for *P* = 1.25 and 1.67. All these states have been checked to be stable in time-dependent simulations.

The density difference at the center of the droplet can be determined by the condition of the existence of the liquid phase in the free space. This condition is obtained by minimizing the grand-potential density E1D − *µ*1*ρ*1 − *µ*2*ρ*2 [5,14], which leads to the zero-pressure condition,

$$p(\rho\_1, \rho\_2) = -\mathcal{E}\_{\rm ID} + \sum\_{j=1,2} \left(\frac{\partial \mathcal{E}\_{\rm ID}}{\partial \rho\_j}\right) \rho\_j$$

$$\equiv -\mathcal{E}\_{\rm ID} + \mu\_1 \rho\_1 + \mu\_2 \rho\_2 = 0. \tag{29}$$

*Symmetry* **2020**, *12*, 174

From this, we obtain relation

$$\frac{\left(\sqrt{\mathcal{G}\_1}\rho\_1 - \sqrt{\mathcal{G}\_2}\rho\_2\right)^2}{2} + \frac{\mathcal{g}\delta\mathcal{g}\left(\sqrt{\mathcal{G}\_2}\rho\_1 + \sqrt{\mathcal{G}\_1}\rho\_2\right)^2}{(\mathcal{g}\_1 + \mathcal{g}\_2)^2} - \frac{\sqrt{m}}{3\pi\hbar}(\mathcal{g}\_1\rho\_1 + \mathcal{g}\_2\rho\_2)^{3/2} = 0,\tag{30}$$

which can be rewritten in the scaled form as

$$\frac{P+GP^{-1}}{2}n\_1^2 + \frac{P^{-1}+GP}{2}n\_2^2 + (G-1)n\_1n\_2 = \frac{1}{3\pi} \left(Pn\_1 + \frac{n\_2}{P}\right)^{3/2}.\tag{31}$$

For given *n*1, we solved Equation (31) to find the respective value of *n*2, which is shown in Figure 6 for *δg*/*g* = 0.05 and several values of *P*. There are two branches of the solutions, that enclose the negative-pressure region, in which QDs may exist. The maximum value of *nj* at the tip of the negative-pressure region corresponds to the density in the droplet's FT segment. The ascending negative-pressure region for each *P* nearly follows relation *n*2 = *Pn*1, which is derived by the minimization condition for the dominant first term in Equation (30) It is seen that a larger difference in the profiles of the two components occurs for larger *P*, as expected. Also, for given *n*1, the negative-pressure region becomes wider with respect to *n*2 for larger *P* (note that the figure displays a log–log plot).

**Figure 5.** (**a**) Stationary solutions of Equation (7), obtained for *δg*/*g* = 0.05 and *N*1 = 100. From the left panel to the right one, the parameter in Equation (8) is *P* = 1.25, 1.67, and 2.5, and the chemical potential for the second component is *µ*2 = −0.018, −0.011, and −0.006, respectively, just below the threshold above which the tails of *ψ*2 extend to infinity. Dashed and solid curves represent the density of the first (*n*1 ) and second (*n*2) components. (**b**) Dependences of *N*2 (black dots: the left vertical axis) and asymmetry parameter *δ*21, defined as per Equation (28) (the red dashed line pertaining to the right vertical axis), on *µ*2 for fixed *N*1 = 100 and *P* = 1.25 or *P* = 1.67.

As the QDs have a finite norm, it is relevant to characterize the asymmetry in terms of the norm, rather than density. Here, we aim to find a largest value of the norm difference,

$$
\Delta\_{21} = (N\_2 - N\_1) / N\_{\rm T} \tag{32}
$$

where *N*T = *N*1 + *N*2 is the total norm, which admits the existence of the QDs. For given *N*1, we obtain the upper bound for *N*2 above which the solution becomes delocalized, and calculate the corresponding critical value of ∆21. The results are shown in Figure 7. For the system with *P* = 1 , the curve demonstrates an empirical dependence ∆21 ∝ *N* −*α* T with exponent *α* ≈ 0.58. Accordingly, the asymmetry tends to vanish asymptotically for very "heavy" droplets, at *N*T → ∞. As the system becomes slightly asymmetric, with *P* = 1.25, exponent *α* is significantly reduced for small *N*T, and converges to a certain finite value at *N*T → ∞. Thus, it is again confirmed that values *P* > 1 maintain conspicuous asymmetry between the QD's components. Finally, strongly asymmetric non-FT (Gaussian-shaped [20]) solutions can be obtained in an approximate analytical form for any value of *P*, as shown in Appendix C.

**Figure 6.** The negative-pressure region in the (*n*1 , *n*2) plane for *δg*/*g* = 0.05 and values of asymmetry parameter in Equation ( 8) *P* = 1 (the solid curve), 1.25 (dashed), 2.5 (dashed-dotted), and 10 (dotted). Boundaries are determined by the zero-pressure condition, as given by Equation (31). The negative pressure, at which localized states may exist, occurs inside the boundaries. Thin lines represent relation *n*2 = *Pn*1 .

**Figure 7.** The inverse of the largest relative norm difference ∆21, up to which the asymmetric droplets exist (see Equation (32)), shown as a function of the total number, *N*T, at different values of asymmetry parameter of Equation (8). Here we set *δg*/*g* = 0.05.

#### 3.2.3. The MI of the Asymmetric PW States

The MI of two-component asymmetric PWs is a relevant subject too. Such solutions are written as *ψj*(*z*, *t*) = √*nj e* −*iµj t* , (*j* = 1, 2). The substitution of this in Equation (7) yields

$$\mu\_1 \quad = (P + GP^{-1})n\_1 + (-1 + G)n\_2 - \frac{P}{\pi} \sqrt{P n\_1 + \frac{n\_2}{P}},$$

$$\mu\_2 \quad = (P^{-1} + GP)n\_2 + (-1 + G)n\_1 - \frac{1}{\pi P} \sqrt{P n\_1 + \frac{n\_2}{P}}.\tag{33}$$

Accordingly, in the symmetric system with *P* = 1, densities of the asymmetric PW state are expressed in terms of the chemical potentials as

$$m\_{\overline{l}} = \frac{1}{4} \left[ \frac{1}{\pi^2 G^2} + \frac{\mu\_1 + \mu\_2}{G} + (-1)^{j+1} (\mu\_1 - \mu\_2) \right] \pm \frac{\sqrt{1 + 2\pi^2 G (\mu\_1 + \mu\_2)}}{4\pi^2 G^2}.$$

We introduce the perturbation around the PW states as

$$
\psi\_{\dot{\jmath}}(z,t) = \left[\sqrt{n\_{\dot{\jmath}}} + \delta\psi\_{\dot{\jmath}}(z,t)\right]e^{-i\mu\_{\dot{\jmath}}t},\tag{34}
$$

$$\delta\psi\_{\rangle} = \mathbb{Z}\_{\rangle}\cos(kz - \Omega t) + i\eta\_{\rangle}\sin(kz - \Omega t),\tag{35}$$

with infinitesimal amplitudes *ζj* and *ηj* , cf. Equation (24). The substitution of this in Equation (7) and the linearization with respect to *ζj* and *ηj* yields the dispersion equation for the perturbation:

$$
\Omega\_{\pm}^{2} = \frac{k^{2}}{4} \left[ k^{2} + 2(P\_{1} + P\_{2} - Q\_{1} - Q\_{2}) \right] \pm \frac{k^{2}}{2} \sqrt{(P\_{1} - P\_{2} - Q\_{1} + Q\_{2})^{2} + 4(R - S)^{2}} \tag{36}
$$

where

$$\begin{aligned} P\_1 &= (P + GP^{-1})n\_1, \quad P\_2 = (P^{-1} + GP)n\_2, \\ Q\_1 &= \frac{P^2 n\_1}{2\pi\sqrt{P n\_1 + P^{-1} n\_2}}, \quad Q\_2 = \frac{P^{-2} n\_2}{2\pi\sqrt{P n\_1 + P^{-1} n\_2}}.\end{aligned} \tag{37}$$

$$\mathcal{R} = (-1 + \mathcal{G})\sqrt{n\_1 n\_2}, \quad \mathcal{S} = \frac{\sqrt{n\_1 n\_2}}{2\pi\sqrt{P n\_1 + P^{-1} n\_2}}.$$

For *P* = 1 and *n*1 = *n*2, these results reproduce Equation (25) for the Ω− branch. A parameter region in which at least one squared eigenfrequency Ω2 ± is negative gives rise to the MI of the two-component state.

#### 3.2.4. The MI for *P* = 1

In Figure 8, we plot the gain spectrum *σ* = Im(Ω) for the asymmetric PWs in the symmetric system with *P* = 1 and *δg*/*g* = 0.05, in the plane of wavenumber *k* and density ratio *n*12 = *n*2/*n*2. For the consistency with the single-component situation displayed in Figure 3, we here fix the total density as (*n*1 + *n*2)/2 = 10. For given *n*12, the MI occurs at |*k*| < *k*0, and the gain attains its maximum at *k* = *k*max = *k*0/ √ 2. The largest gain is obtained at equal densities, *n*12 = 1. Both the *k*-band of the instability and magnitude of the gain slowly decrease as the deviation of *n*12 from unity increases. This means that the MI occurs in the PW states with a large density difference, thus giving rise to the formation of solitons with large asymmetry even for equal intra-component MF interaction strengths, *P* = 1 (see Equation (8)).

**Figure 8.** Color-coded values of the MI gain, *σ* = Im(Ω), for asymmetric PWs, as calculated from Equation (36) in the plane of wave number |*k*| and density ratio *n*12 = *n*1/*n*2, are displayed for (**a**) *P* = 1 and (**b**) *P* = 1.25 with fixed *δg*/*g* = 0.05 and (*n*1 + *n*2)/2 = 10. The solid and dashed white curves represent the MI boundary *k* = *k*0 and the peak value of the MI gain at *k* = *k*max = *k*0/ √ 2, respectively. In (**c**), we plot *σ*(*k*max) (circles) and *n* max 12 (triangles) versus *P*.

In Figure 9 we display typical examples of the numerically simulated development of the MI in the symmetric two-component system with *P* = 1 and population imbalance. Figure 9a shows the evolution of central-point values of the density of the first component, *n*1(*z* = 0), for different values of the density ratio, *n*12 = *n*1/*n*2. Time required for the actual onset of the instability increases with the increase in *n*12, as is clearly shown by the density-plot evolution in Figure 9d,e for *n*12 = 1 and Figure 9f,g for *n*12 = 9. This observation can be understood in terms of the MI gain *σ*, as shown in Figure 8c, where *σ* at *k* = *k*max becomes smaller with increasing *n*12.

Spatial profiles at fixed time, which are plotted in Figure 9b,c for these two cases, show fragmentation of the profiles into sets of localized structures. The decrease in the number of fragments with the increase of *n*12 is explained by the decrease of *k*max, see Figure 8a. For *n*12 = 1, the results are the same as in the single-component case, as coinciding profiles in the two components of the symmetric system are stable against spontaneous symmetry breaking. On the other hand, when *n*12 6= 1 an in-phase two-component localized structure appears, keeping the initial density imbalance. Since one can select an arbitrary ratio of densities of the two components for the initial PW state, a highly asymmetric structure, like the one displayed in Figure 9c, may emerge even for *P* = 1, as a result of the MI-induced nonlinear evolution.

#### 3.2.5. The MI for *P* 6= 1

Figure 8b represents the MI gain for *P* = 1.25 and a fixed total density, (*n*1 + *n*2)/2 = 10, in the case of slightly different strengths of the intra-component repulsion. The peak value of the MI gain is attained at *n*12 = *n* max 12 = 0.577, below the equal-densities point *n*12 = 1. This is consistent with the fact that, at *P* > 1, unequal values *n*1 < *n*2 are suitable to the formation of an asymmetric soliton structure, as seen in Figure 5a. In Figure 8c, we plot the peak MI gain, *σ*(*k*max), along with the respective value of the density ratio, *n* max 12 , as a function of *P*. Value *n* max 12 monotonously decreases as a function of *P*, while the peak gain attains a minimum at *P* = 1.

In Figure 10, we present the development of the MI in the two-component system for *P* = 1.25 and a fixed total density, (*n*1 + *n*2)/2 = 10. Figure 10a displays the evolution of the central-point density of the first component, *n*1(*z* = 0), for different values of the density ratio, *n*12 = *n*1/*n*2. It shows that time required for the development of the MI increases with the increase in the asymmetry of the density. This is also made evident by the density plots of the temporal evolution of the first component, shown in Figure 10e–g. This result is consistent with Equation (36), which shows a decrease of the MI gain with the increase of the asymmetry even for *P* 6= 1. Spatial profiles at fixed time, displayed in Figure 10b–d, show fragmentation of the profiles. Figure 10c clearly indicates that,

even for *n*12 = 1, the MI generates asymmetric droplet-like structures similar to Figure 5a, where the complete overlapping of the two densities does not occur.

**Figure 9.** Numerically simulated development of the MI of asymmetric PW states in the two-component system, with *P* = 1 and *δg*/*g* = 0.05 . The initial PW states are taken with fixed total density, (*n*1 + *n*2)/2 = 10. (**a**) The evolution of the central density of the first component, *n*1(*z* = 0), for different density ratios in the two components, *n*12 = *n*1/*n*2. (**b**,**c**) Snapshots of density profiles for the cases of (**b**) *n*12 ≡ *n*1/*n*2 = 1 at *t* = 200 and (**c**) *n*12 = 9 at *t* = 400. Panels (**d**,**e**) and (**f**,**g**) are top views of the spatiotemporal evolution of the densities, *n*1 (*z*, *t*) and *n*2(*z*, *t*), for *n*12 = 1 and *n*12 = 9, respectively. Simulations were performed in the domain −50 ≤ *z* ≤ +50 with 2048 grid points, subject to periodic boundary conditions. In this figure and in Figure 10, the scaled time unit corresponds to ∼ 1 µs in physical units.

**Figure 10.** Numerically simulated development of the modulational instability in the two-component system with *δg*/*g* = 0.05 and *P* = 1.25. The initial PW states are taken with a fixed total density, (*n*1 + *n*2)/2 = 10. (**a**) The evolution of the central density of the first component, *n*1(*z* = 0), for different density ratios in the two components, *n*12 = *n*1/*n*2. (**b**–**d**) Snapshots of density profiles for the cases of (**b**) *n*12 ≡ *n*1/*n*2 ∼ 0.1 at *t* = 300, (**c**) *n*12 = 1 at *t* = 200 and (**d**) *n*12 = 9 at *t* = 600. Panels (**e**–**g**) represent the top view of the spatiotemporal evolution of the densities, *n*1 (*z*, *t*), corresponding to (**b**–**d**), respectively (the evolution of *n*2 (*z*, *t*) shows similar patterns). Simulations were performed in the domain −50 ≤ *z* ≤ +50 with 2048 grid points, subject to periodic boundary conditions.

#### **4. Conclusions**

The main purpose of this work is to associate the MI (modulation instability) of plane waves (PWs) to the mechanism of the creation of QDs (quantum droplets) in the system described by the coupled GP (Gross–Pitaevskii) equations including the LHY (Lee–Huang–Yang) terms in the 1D setting. This system is the model of weakly interacting binary Bose gases with approximately balanced interactions between the intra-component self-repulsion and the inter-component attraction. We have investigated, analytically and numerically, the MI of the lower branch of PW states in both symmetric (effectively single-component) and asymmetric (two-component) GP systems, and ensuing formation of a chain of droplet-like states. In particular, numerical solution for QDs which are

asymmetric with respect to the two components are obtained, both in the system with equal repulsion strengths but unequal populations in the two components, and in the one with different self-repulsion strengths. The results corroborate that the previously known symmetric states are robust against symmetry-breaking disturbances.

These predictions can be tested experimentally by preparing uniform binary Bose gases with equal or different densities of two components, and suddenly reducing the strength of the effective MF (mean-field) interaction by means of the Feshbach-resonance quench, in order to enhance the relative strength of the LHY terms. In particular, for typical values of physical parameters, an estimate of the characteristic time of the modulation instability growth for typical values of the physical parameters is ∼1 µs. This time is much smaller than a typical lifetime of the droplet, which is &100 µs [21–23,27], thus making the observation of the MI feasible. The present analysis being restricted to the 1D setting, effects of the tight transverse confinement and crossover to the 3D configuration [12,44,45] deserves further consideration.

**Author Contributions:** Conceptualization and investigation, T.M., A.M., K.K., B.A.M., and A.K.; analytics, A.K. and B.A.M.; numerics, T.M., K.K., and A.M.; writing—original draft preparation, T.M.; writing—review and editing, A.M., K.K., B.A.M., and A.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We appreciate valuable comments received from M. Modugno. T.M. acknowledges support from IBS (Project Code IBS-R024-D1). A.M. acknowledges support from the Ministry of Education, Science, and Technological Development of Republic of Serbia (project III45010) and the COST Action CA 16221. The work of K.K. is partly supported by the Japan Society for the Promotion of Science (JSPS) Grant-in-Aid for Scientific Research (KAKENHI Grant No. 18K03472). B.A.M. appreciates support from the Israel Science Foundation through grant No. 1287/17. A.K. thanks the Indian National Science Academy for the grant of INSA Scientist Position at Physics Department, Savitribai Phule Pune University).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Other Exact Solutions for the Single-Component GP Equation**

Here we briefly list other types of exact solutions of the single-component Equation (12), in addition to the FT and PW solutions in Equations (13) and (22) which were considered in detail above (solutions to Equation (12) in the form of dark and anti-dark solitons were reported in [46]).The stability of a majority of these solutions is not addressed here, as it should be a subject for a separate work.

#### *Appendix A.1. δg*/*g* > 0

In the case of comparable quadratic self-attraction and cubic repulsion in Equation (12) with *δg* > 0, exact spatially-periodic solutions with odd parity can be expressed in terms of the Jacobi's elliptic sine, whose modulus *q* is an intrinsic parameter of the family:

$$\psi(z,t) = \exp\left(-i\mu\_{\text{sn}}t\right) \left[A\,\text{sn}(\beta z, q) + B\right],\tag{A1}$$

where

$$B = \frac{\sqrt{2}}{3\pi} \frac{\mathcal{g}}{\delta g} > 0, \quad A = \sqrt{\frac{2}{1+q^2}} B > 0, \quad \mu\_{\text{sn}} = -2\frac{\delta g}{\mathcal{g}} \mathcal{B}^2 < 0, \quad \mathcal{\beta}^2 = \frac{2}{(1+q^2)} \frac{\delta g}{\mathcal{g}} \mathcal{B}^2.$$

In the limit of *q* → 1, Equation (A1) goes over into the kink (the same as found in [46]),

$$\psi(z,t) = \exp\left(-i\mu\_{\text{kink}}t\right) \left[A\tanh(\beta z) + B\right],\tag{A2}$$

with parameters

$$A = B = \frac{\sqrt{2}}{3\pi} \frac{\text{g}}{\delta g} > 0, \quad \mu\_{\text{kink}} = -2 \frac{\delta \text{g}}{\mathcal{S}} B^2, \quad \mathcal{S}^2 = \frac{\delta \text{g}}{\mathcal{S}} B^2.$$

#### *Appendix A.2. δg*/*g* < 0

In the case when the inter-species MF attraction is stronger than the intra-species repulsion, resulting in *δg* < 0, spatially-periodic solutions are expressed in terms of even Jacobi's elliptic functions, dn(*x*, *q*) and cn(*x*, *q*). First, it is

$$\psi(z,t) = \exp\left(-i\mu\_{\mathrm{dn}}t\right) \left[A\,\mathrm{dn}(\beta z, q) + B\right],\tag{A3}$$

with the elliptic modulus taking all values 0 < *q* < 1, other parameters being

$$B = \frac{\sqrt{2}}{3\pi} \frac{\mathcal{g}}{\delta g} < 0, \quad A = -\sqrt{\frac{2}{2 - q^2}} B > 0, \quad \mu\_{\text{dn}} = -2B^2 \frac{\delta g}{\mathcal{g}} > 0, \quad \mathcal{J}^2 = -\frac{2}{(2 - q^2)} \frac{\delta g}{\mathcal{g}} B^2.$$

The second solution is expressed in terms of the elliptic cosine, with *q* 2 > 1/2:

$$\psi(z,t) = \exp\left(-i\mu\_{\rm cn}t\right) \left[A \operatorname{cn}(\beta z, q) + B\right],\tag{A4}$$

$$B = \frac{\sqrt{2}}{3\pi} \frac{\mathcal{g}}{\delta g} < 0, \quad A = -\sqrt{\frac{2}{2q^2 - 1}} B > 0, \quad \mu\_{cn} = -2\frac{\delta g}{\mathcal{g}} B^2 > 0, \quad \mathcal{J}^2 = -\frac{2}{(2q^2 - 1)} \frac{\delta g}{\mathcal{g}} B^2.$$

In the limit of *q* → 1, both solutions in Equations (A3) and (A4) carry over into a state of the "bubble" type [49], which changes the sign at two points (the same solution was reported as an "W-shaped soliton" in [46]):

$$\psi(z,t) = \exp\left(-i\mu\_{\text{bulk}}t\right) \left[A\text{sech}(\beta z) + B\right],\tag{A5}$$

with parameters

$$B = \frac{\sqrt{2}}{3\pi} \frac{\text{g}}{\delta \text{g}} < 0, \quad A = -\sqrt{2}B > 0, \quad \mu \text{ (bubble } = \beta^2 = -2\frac{\delta \text{g}}{\mathcal{S}}B^2 > 0.1)$$

#### **Appendix B. Analytical Solutions for Strongly Asymmetric Fundamental and Dipole States**

Here we consider analytical solutions of Equation (7) with strong asymmetry, *N*1 *N*2, which can be found under small-amplitude conditions, *n*1(*z* = 0) *n*2(*z* = 0) *n*0. Then, cubic terms may be neglected in Equation (7), and approximation p *P*|*ψ*1| 2 + *P*−1 |*ψ*2| 2 ≈ *P* −1/2 |*ψ*2| is used to simplify Equation (7) to the following equations for stationary states in Equation (9):

$$
\mu\_1 \phi\_1 = -\frac{1}{2} \frac{d^2 \phi\_1}{dz^2} - \frac{\sqrt{P}}{\pi} \phi\_2 \phi\_1 \tag{A6}
$$

$$
\mu\_2 \phi\_2 = -\frac{1}{2} \frac{d^2 \phi\_2}{dz^2} - \frac{1}{\pi P^{3/2}} \phi\_2^2. \tag{A7}
$$

Although this case is somewhat formal, in terms of the underlying concept of the quantum droplets, which is essentially based on the competition of residual MF and LHY terms, it is interesting to consider it too.

The soliton solution of Equation (A7) is obvious,

$$\phi\_2(z) = \frac{3\pi}{2} \left( -\mu\_2 \right) \frac{P^{3/2}}{\cosh^2 \left( \sqrt{-\mu\_2/2} z \right)} \tag{A8}$$

where the solution in Equation (13) takes essentially the same form in the limit of |*µ*| *µ*0. Then, the substitution of Equation (A8) in Equation (A6) makes it tantamount to the linear Schrödinger

equation with the Pöschl-Teller potential [50]. The ground-state (GS) solution of Equation (A6) for *φ*1, with arbitrary amplitude *φ* (0) 1 ,

$$\left(\phi\_1(z)\right)\_{\rm GS} = \frac{\phi\_1^{(0)}}{\left[\cosh\left(\sqrt{-\mu\_2/2}z\right)\right]^{\gamma}}\tag{A9}$$

exists with

$$\gamma = \frac{1}{2} \left( \sqrt{24P^2 + 1} - 1 \right),\tag{A10}$$

and eigenvalue

$$(\mu\_1)\_{\rm GS} = \left(\sqrt{24P^2 + 1} - 1\right)^2 \frac{\mu\_2}{16}.\tag{A11}$$

In this case, the QD solutions are quasi-Gaussian objects [20]. Note that, in the symmetric system with *P* = 1, Equations (A10) and (A11) yield *γ* = 2 and (*µ*1)GS = *µ*2, i.e., the eigenmode and eigenvalue coincide with their counterparts in the soliton solution in Equation (A8), while they are different in the asymmetric system, the GS level lying below or above the chemical potential of soliton in Equation (A8) at *g*1 > *g*2 and *g*1 < *g*2, respectively.

In Figure A1 we compare a typical asymptotic solution given by Equations (A8) and (A9) with a numerically obtained GS solution for the same values of the parameters. It is seen that the analytical and numerical results match well.

**Figure A1.** Comparison of the asymptotic analytical solutions, given by Equations (A8) and (A9), with their numerically obtained counterparts. The density of the first (*n*1 ) and second (*n*2) components are displayed in top and bottom panels, respectively. Solid blue lines represent the numerical results, while dashed red lines represent the analytical solution. Here, parameters are *δg*/*g* = 0.05, *N*1 = 0.0001067, *N*2 = 0.0148044 and (*µ*2)GS = *µ*2 = −0.005.

Further, it is also possible to produce the first excited state of Equation (A6) in the form of the dipole (antisymmetric) mode with an arbitrary amplitude:

$$\Phi\left(\phi\_1(z)\right)\_{\text{dip}} = \frac{\phi\_1^{(0)} \sinh\left(\sqrt{-\mu\_2/2}z\right)}{\left[\cosh\left(\sqrt{-\mu\_2/2}z\right)\right]^\gamma},\tag{A12}$$

where *γ* is the same as in Equation (A10), the respective eigenvalue being

$$\left(\left(\mu\_1\right)\_{\text{dip}} = \left(\sqrt{24P^2 + 1} - 3\right)^2 \frac{\mu\_2}{16}\right) \tag{A13}$$

which is obviously higher than its GS counterpart in Equation (A11) (at *P* = 1, Equation (A13) yields (*µ*1)dip = *µ*2/4, and (*µ*1)dip falls below *µ*2 for *P* > √ 2). Unlike the GS, the dipole mode exists not at all values of *P*, but only for *P* > √ 1/3. Exactly at *P* = √ 1/3, one has (*µ*1)dip = 0, and the dipole mode in Equation (A12), with *γ* = 1, is a delocalized one, ∼ tanh p −*µ*2/2*z* .

Linear Schrödinger Equation (A6) with the Pöschl-Teller potential may give rise to higher bound states of integer order *ν* as well, with eigenvalues

$$(\mu\_1)\_\nu = \left(\sqrt{24P^2 + 1} - (1 + 2\nu)\right)^2 \frac{\mu\_2}{16} \tag{A14}$$

where *ν* = 0 and 1 correspond to Equations (A11) and (A13), respectively, the *ν*-th spate existing at *P* 2 > *ν* (*ν* + 1) /6. The number of such solutions is always finite.

Unlike solutions considered in Appendices A and C, the stability of the solutions given by Equations (A8)–(A14) is obvious.

#### **Appendix C. Other Exact Solutions in the Case of** *N*1 *N*2

Here we provide periodic solutions to the semi-linear system of Equations (A6) and (A7) in terms of Jacobi elliptic functions. In the limit of *q* → 1, they go over into solutions given in the main text, in the form of Equations (A8), (A9), and (A12).

#### *Appendix C.1. Solution of Equation (A7)*

An exact periodic solution of Equation (A7) with the quadratic nonlinearity is

$$\phi\_2 = A \left[ \text{dn}^2(\beta z, q) + p \right], \tag{A15}$$

with

$$\beta^2 = -\frac{\mu\_2}{2\sqrt{1-q+q^2}}, \quad A = -\frac{3\pi\mu\_2 P^{3/2}}{2\sqrt{1-q+q^2}}, \quad p = \frac{-(2-q) + \sqrt{1-q+q^2}}{3}.\tag{A16}$$

In the limit of *q* → 1, the solution in Equation (A15) goes over into the solution in Equation (A8). Note that *p* is vanishing in this limit, according to Equation (A16).

#### *Appendix C.2. Solutions of Equation (A6)*

We now show that, with *φ*2 given by Equation (A15), linear Equation (A6) *φ*1 has several particular solutions depending on the value of *P*.

**Solutions For** *P* 2 = 1/3

Appendix C.2.1. Solution I

It is easy to check that

$$
\phi\_1 = \phi\_1^{(0)} \text{dn}(\beta z, q) \tag{A17}
$$

is an exact solution to Equation (A6), provided that

$$P^2 = \frac{1}{3}, \quad \mu\_1 = \left(\frac{\mu\_2}{12}\right) \frac{2 - q + 2\sqrt{1 - q + q^2}}{\sqrt{1 - q + q^2}}.$$

Appendix C.2.2. Solution II

$$\boldsymbol{\phi}\_{1} = \boldsymbol{\phi}\_{1}^{(0)} \mathbf{c} \mathbf{n}(\boldsymbol{\beta} \boldsymbol{z}, \boldsymbol{q}) \tag{A18}$$

is an exact solution to Equation (A6), provided that

$$P^2 = \frac{1}{3}, \quad \mu\_1 = \left(\frac{\mu\_2}{12}\right) \frac{2q - 1 + 2\sqrt{1 - q + q^2}}{\sqrt{1 - q + q^2}}.$$

In the limit of *q* → 1, solutions I and II go over into the solution Equation (A9) with *γ* = 1 and *µ*1 = *µ*2/4.

Appendix C.2.3. Solution III

$$\phi\_1 = \phi\_1^{(0)} \text{sn}(\beta z, q) \tag{A19}$$

is an exact solution to Equation (A6), provided that

*P* 2 = 1 3 , *µ*1 = *µ*2 12 2 p 1 − *q* + *q* 2 − (1 + *q*) p 1 − *q* + *q* 2 .

In the limit of *q* → 1, solution III goes over into the solution Equation (A12) with *γ* = 1 and *µ*1 = 0. **Solutions for** *P* 2 = 1**.**

Appendix C.2.4. Solution IV

It is easy to check that

$$\phi\_1 = \phi\_1^{(0)} [\text{dn}^2(\beta z, q) + p] \tag{A20}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \mu\_2 \dots$$

Appendix C.2.5. Solution V

$$\phi\_1 = \phi\_1^{(0)} \text{cn}(\beta z, q) \text{dn}(\beta z, q) \tag{A21}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \left(\frac{\mu\_2}{2}\right) \frac{q + \sqrt{1 - q + q^2}}{\sqrt{1 - q + q^2}}.$$

In the limit *q* = 1, solutions IV and V go over into solution Equation (A9) with *γ* = 2 and *µ*1 = *µ*2. Appendix C.2.6. Solution VI

$$\phi\_1 = \phi\_1^{(0)} \text{sn}(\beta z, q) \text{dn}(\beta z, q) \tag{A22}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \left(\frac{\mu\_2}{4}\right) \frac{\Im(1-q) + \sqrt{1-q+q^2}}{\sqrt{1-q+q^2}}.$$

Appendix C.2.7. Solution VII

$$\boldsymbol{\phi}\_{1} = \boldsymbol{\phi}\_{1}^{(0)} \operatorname{sn}(\beta z, \boldsymbol{\eta}) \operatorname{cn}(\beta z, \boldsymbol{\eta}) \tag{A23}$$

is an exact solution to Equation (A6), provided that

$$P^2 = 1, \quad \mu\_1 = \left(\frac{\mu\_2}{4}\right) \frac{2\sqrt{1-q+q^2} - (2-q)}{\sqrt{1-q+q^2}} \frac{1}{4}$$

In the limit of *q* → 1, solutions VI and VII go over into Equation (A12), with *γ* = 2 and *µ*1 = *µ*2/4.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).


**Meet the editor**

Faik Atroshi, MScBM, Licentiate, PhD clinical biomarkers, Docent, senior researcher in pharmacology and toxicology, Adjunct professor in clinical genetics and clinical nutrition, visiting professors at different international universities; Dip Genetic engineering; Dip Molecular genetics; Dip Analytical clinical pharmacology; Director, education & research satellite centre in

trace elements for UNESCO; President of the Global Society for Nutrition, Environment and Health (GSNEH.com); Vice president of the Finnish society for biological medicine, chief editor of International Journal of Trends in Medicine. Research interests and publications are in the field of biomedicine and health; Free radicals and inflammation, oxidative stress and diseases, antioxidants and DNA repair.

Contents

**Preface IX**

Jerzy Majkowski

**Human Diseases 43** Vincent S. Gallicchio

**Thalassemia 119**

Ruhi Cakir

A.V. Skalny

Chapter 1 **Epilepsy Treatment and Nutritional Intervention 1**

Chapter 4 **Antioxidants as Complementary Medication in**

Chapter 5 **General Aspects of Ozone Therapy 159**

**Constituents and Sleep 191**

Somdet Srichairatanakool and Suthat Fucharoen

Chapter 6 **Dietary Aspects in Cancer Prevention — A Mini-Review 179**

Chapter 7 **Nutrition, Sleep and Sleep Disorders – Relations of Some Food**

Chapter 8 **Bioelements and Bioelementology in Pharmacology and Nutrition: Fundamental and Practical Aspects 225**

Markku Partinen, Tuomas Westermarck and Faik Atroshi

Mohamed Abdulla and Sangeeta Shukla

Chapter 2 **Use of Trace Elements and Halotherapy in the Treatment of**

Chapter 3 **Oxidative Stress and Dietary Interventions in Autism: Exploring the Role of Zinc, Antioxidant Enzymes and Other Micronutrients in the Neurobiology of Autism 97** Maija L. Castrén, Tuomas Westermarck and Faik Atroshi


**1. Introduction**

Early diagnosis of tuberculosis and drug resistance improves survival and by identifying in‐ fectious cases promotes contact tracing, implementation of institutional cross-infection pro‐ cedures, and other public-health actions. There have been many advances in methodology for tuberculosis diagnosis [1-3].

For every stages of diagnosis, there are new approaches. New tests are available by level of laboratory and phase of application.


Imprisonment is not only removal of freedom, but also a stigma. The penalty is a state-sanctioned deprivation of individual glory. When a member of the elite breaks the law, it can result in a problematic situation for others within the elite. An "elite" human's crime can threaten the elite's chastity in a manner that offences committed by "sub-humans" do not.

While a white-collar criminal's actions can threaten the legitimacy of the existing hierarchy of high and low positions, between the honest and virtuous ones on the one side and the dishonest and dishonorable on the other, a social security fraudster's actions will confirm and reinforce the same hierarchy. Could this be a reason why representatives of the social elite are more likely to strike hard and visibly more forcefully on welfare fraud, while they rarely hold up white-collar criminals as a social problem, despite the fact that white-collar crime probably costs us much more than social security fraud? Is this the reason why petty thieves, drug addicts, and NAV fraudsters dominate everyday life in Norwegian courtrooms, since they are considered major problems in society, while the corruption defendant from the executive suite is an unfortunate single case and moreover did not mean it anyway?

This book has presented two illustrative examples of how a social security fraudster and a white-collar criminal respectively are treated in the courts. We know of more major cases that NAV has highlighted, characterized by systematic fraud in networks. Such cases do naturally have a different severity than the case with the fortune teller. But, at the same time, we know of financial crime cases that were far more serious than the CEO case that we have presented.

So, what should the police prioritize in their investigations—social security fraud or white-collar crime? It is probably pretty obvious by now what we think about this matter. We have shown that the damage from white-collar crime is at least as major as the damage from social security fraud in an economic sense. We have shown the characteristics of the disadvantaged social security scammer versus the gifted white-collar criminal. We believe combating economic crime starts at the top, not at the bottom. Why on earth should the bottom people follow laws and regulations when on a daily basis they see in the media that the elite gets away with all kinds of misconduct.

With 1200 to 1500 individuals reported for social security fraud every year, substantial police and prosecution resources are taken up with minor criminal cases. The question is whether this is an appropriate use of resources, or whether the police should instead dismiss the majority of cases and only legally pursue the largest and most organized fraudsters. Although white-collar cases can be far more costly than social security matters for law enforcement, it is important to have a form of criminal justice. There is supposed to be equality before the law.

Most white-collar crime cases never end up in police hands. Our experts believe that 59 percent of all white-collar criminals are never detected. The detected 41 percent are not all followed up: 19 percent are not investigated, 8 percent are not prosecuted, 5 percent are not convicted, while only 9 percent are convicted. When 8 percent are not prosecuted, it is often because the police failed to provide adequate evidence. Thus, according to our experts, 19 + 8 = 27 percent—or one in four felons—go free because of police failure to prioritize financial crime committed by the elite in society.

### References

© The Author(s) 2018 139

Medievitenskap (Master Thesis in Media Science). Universitetet i Oslo (University of Oslo).


**2. Vitamin D deficiency**

Vitamin D deficiency is found in 30–50% of the general population, and prevalence estimates suggest that more than 1 billion people worldwide are vitamin D insufficient or deficient [4, 11].

Vitamin D deficiency is indicated by serum levels of 25(OH)D < 20 ng/mL [4]. Serum levels >30 ng/mL are likely optimal for bone health, but some studies have shown benefits with lesser values. Parathyroid hormone (PTH) suppression appears to plateau at levels between 30 and 40 ng/mL [12]. There has been no agreement on optimum 25(OH)D levels required for purported health benefits beyond skeletal health. One study suggested that 25(OH)D levels below 11–14 ng/mL signify increased CVD risk [13]. Levels in the range of 21–29 ng/mL are considered by some as insufficient, a definition that would label the majority of the U.S. population vitamin D insufficient [4].

The Endocrine Society does not recommend screening for vitamin D deficiency in individuals who are not at risk for it [12]. Many of the risk factors for vitamin D deficiency have been identified. Some of these include inadequate cutaneous synthesis stemming from insufficient sun exposure or dark skin pigmentation and inadequate dietary intake. Other noted risk factors include aging, obesity, renal disease, liver disease, disorders that affect fat absorption (e.g. celiac disease, inflammatory bowel diseases, types of bariatric surgery), increased catabolism due to medications (e.g. steroids, anticonvulsants etc.), and other hereditable (e.g. rickets) and acquired disorders (e.g. hyperthyroidism) [4].


*3.2. LCOH with Historical Price Duration Curve*

Comparing the stochastically optimized LCOH to one optimized using historical PDC data is useful for understanding the implications of the distribution tails on cost. The stochastic optimization case outputs the expected LCOH under a wide range of possible synthetic PDC. LCOH optimization using the historical dataset gives an example of the LCOH found on a PDC that is slightly skewed to higher prices.

As with the stochastic case, Figure 9 shows that hydrogen is not dispatched at low hydrogen prices and that there is a range in which a small amount of energy is dispatched for hydrogen production despite the inability to recover the IS capital cost. In this case, when hydrogen prices fall below approximately 80 JPY/m3, the system does not dispatch hydrogen. At 80–98 JPY/m3, a small amount of hydrogen is dispatched. At 98.1 JPY/m3, the ΔNPV equals zero, thus representing the LCOH.

Figure 10 shows the IS unit's utilization rate plotted against the hydrogen price. The utilization rate is zero hours when the hydrogen price is low. At an LCOH of 98.1 JPY/m3, the IS unit produces hydrogen for 637 h. Figure 10 shows the utilization rate of the IS under different hydrogen price conditions.

**Figure 9.** ΔNPV for various hydrogen prices using the historical PDC as input. The red dot represents the breakeven LCOH.

**Figure 10.** Utilization rate of the IS unit plotted against the hydrogen price. As the hydrogen price rises, hydrogen deployment becomes increasingly more economically advantageous than electricity sale, so the number of hydrogen production hours increases. The red dot represents the breakeven LCOH.

Table 5 summarizes the dispatch parameters found at an LCOH of 98.1 JPY/m3 using historical price inputs.

**Table 5.** Dispatch values for the system at a levelized cost of hydrogen of 98.1 JPY/m3.

#### **4. Discussion**

The reported LCOH values should not be relied on as a basis for making investment decisions. Rather, they help us understand the implications of different inputs so that when economic competitiveness is evaluated, the correct breadth of input data can be applied.

The effects of price distribution can be viewed by comparing the stochastic optimization case to the historical PDC case. The hydrogen dispatch is driven by two factors: the hours having the lowest electricity prices and the price of hydrogen. By raising the price of hydrogen, selling hydrogen becomes more profitable during more hours. Lower electricity prices and more incidences of low electricity prices also make hydrogen more economically advantageous than electricity.

The lowest-priced hours of the electricity price distribution are what dictate system profitability, since the capacity is fixed and the hydrogen price varied. The stochastic optimization case uses synthetic price histories in an attempt to produce the expected LCOH. On average, the synthetic histories showed lower electricity prices at the tail than did the historical price distribution. This led to a lower LCOH than in the historical case.

This lowest-priced hour distribution phenomenon is illustrated in Table 6. The lowest 500 h of electricity prices from the year are averaged and compared with the LCOHs for several different synthetic histories. The average electricity price over the year is also provided. The LCOH shares a stronger correlation with the bottom-hour average than with the total yearly price average. By way of comparison, the LCOH at the stochastic mean and historical points are \$0.64/m3 and \$0.78/m3, respectively at an exchange rate of 106 JPY/USD.

**Table 6.** Impacts of the distribution tail on LCOH.

1 The historical scenario is the only scenario not produced by sampling the ARMA model.

Synthetic data produced using the ARMA method outputted cheaper bottom-500 h price averages, as well as overall average prices that were lower than the historical averages. This meant that the distribution of PDCs was slightly more favorable to hydrogen dispatch than the historical PDC. As such, the LCOH was lower in the stochastic case than in the historical case.

This analysis demonstrates that careful consideration should be taken when applying PDCs to this type of economic dispatch problem. The breakeven price of hydrogen highly depends on the PDC input. Stochastic optimization helps reduce uncertainty, but care should still be taken to produce PDCs that are meaningful with regard to the chosen timeframe of analysis. For example, using a 2020 PDC to predict the 2030 LCOH would be inappropriate. A projection of 2030 prices would be acceptable, but the best practice would be to use a host of projected possibilities to produce an expected LCOH.

The results from this study also show that lower overall electricity prices and more incidences of low prices would provide greater economic incentives for hydrogen production. This means that NPPs in locations with depressed electricity prices due to factors such as zero- or negative-bid renewable energies, mild climates, or low electricity demand could provide hydrogen at a lower price yet still break even or potentially turn a profit.

Several other pathways exist for reducing the LCOH. Reducing capital expenditure would depress the LCOH. The effects and sizes of potential storage options could be explored in more detail. Additional cashflows generated by the NPP's ability to participate in other areas of the electricity market would lower the ΔNPV and thus the LCOH, as well. Before investment decisions are made, each of these sensitivities should be investigated to better understand their feedback.


Advanced Analytics and Artificial Intelligence Applications

#### Table 5.

Summary with the best results of similarity.

Figures 9 and 10 show a comparison of results obtained from the algorithm with 100, 200, and 500 iterations of the instance of 1027 services. It can be seen that the best values of f(x) are produced with 500 iterations; however the result of the number of groups involved is more stable and better with 200 iterations.

This chapter describes a hybrid algorithm for Web services clustering; the approach is based on the ABC optimization algorithm combined with the K-means and Consensus. The clustering process starts from the extraction and processing of WSDL documents, then the calculation of semantic similarities between pairs of Web services operations. The semantic similarity measurements are based on WordNet, combining corpus-based and taxonomic approaches. As a result of the calculations, eight matrices are generated which are the input data set to the hybrid

Biologically inspired algorithms offer advantages over conventional clustering methods, such as the ability to find a near optimal solution by updating the candi-

The clustering algorithm designed is based on the behavior of the bees but was improved by incorporating K-means and Consensus so that the algorithm adjusts itself in each iteration. A series of experiments were conducted with the hybrid ABC algorithm, using optimal values for each adjustable parameter. The experiments were carried out with three collections with 50, 647, and 1027 Web services, and the algorithm was executed with variants of 100, 200, and 500 iterations. The hybrid

As future work, more combinations of semantic similarities, as well as incorporating more information of the description of the services and the data types incor-

date solutions iteratively and have self-organizing behavior.

Comparison of the number of groups in the best case for 1027 services.

Bio-Inspired Hybrid Algorithm for Web Services Clustering

DOI: http://dx.doi.org/10.5772/intechopen.85200

ABC algorithm has shown good results for Web services clustering.

8. Conclusions

Figure 10.

bio-inspired algorithm.

23

porated in the XML service definitions.

#### Table 6.

Summary with the best results of clusters.

Table 5 shows the results obtained with configurations: the best, worst, average, variance, and standard deviation of f(x) over the 20 executions of the algorithm (Table 6).

Figure 9. Comparison of similarity between the best solutions found with 1027 services. Bio-Inspired Hybrid Algorithm for Web Services Clustering DOI: http://dx.doi.org/10.5772/intechopen.85200

Figure 10. Comparison of the number of groups in the best case for 1027 services.

Figures 9 and 10 show a comparison of results obtained from the algorithm with 100, 200, and 500 iterations of the instance of 1027 services. It can be seen that the best values of f(x) are produced with 500 iterations; however the result of the number of groups involved is more stable and better with 200 iterations.

#### 8. Conclusions

Table 5 shows the results obtained with configurations: the best, worst, average, variance, and standard deviation of f(x) over the 20 executions of the algorithm

Comparison of similarity between the best solutions found with 1027 services.

Collection-iterations Best Worst Mean Variance Deviation –100 0.58936208 0.29274988 0.45901481 0.00593354 0.0770295 –200 0.54853117 0.44904535 0.50052478 0.00088301 0.02971549 –500 0.57163535 0.46554965 0.51820976 0.00090081 0.03001343 –100 0.54947584 0.2153186 0.44401451 0.01027133 0.10134758 –200 0.54778745 0.20257995 0.40923046 0.01554565 0.12468218 –500 0.54881541 0.22262872 0.41911288 0.01280278 0.11314939 –100 0.55531806 0.27378002 0.44488931 0.00847871 0.09207991 –200 0.56846422 0.28919387 0.47529505 0.00772078 0.08786795 –500 0.66112025 0.20302145 0.4674717 0.00977301 0.09885856

Advanced Analytics and Artificial Intelligence Applications

Collection-iterations Best Worst Mean Variance Deviation –100 2 10 4.31578947 6.22105263 2.49420381 –200 2 4 3 0.89210526 0.94451324 –500 2 6 3.94736842 2.36578947 1.53811231 –100 3 49 18.0526316 221.207895 14.8730594 –200 2 67 18.9473684 305.831579 17.488041 –500 2 73 23.2105263 426.368421 20.6486905 –100 2 110 31.2631579 1066.82895 32.6623475 –200 2 69 22.1052632 355.884211 18.8648936 –500 2 123 31 1179.21053 34.3396349

(Table 6).

Figure 9.

22

Table 6.

Table 5.

Summary with the best results of similarity.

Summary with the best results of clusters.

This chapter describes a hybrid algorithm for Web services clustering; the approach is based on the ABC optimization algorithm combined with the K-means and Consensus. The clustering process starts from the extraction and processing of WSDL documents, then the calculation of semantic similarities between pairs of Web services operations. The semantic similarity measurements are based on WordNet, combining corpus-based and taxonomic approaches. As a result of the calculations, eight matrices are generated which are the input data set to the hybrid bio-inspired algorithm.

Biologically inspired algorithms offer advantages over conventional clustering methods, such as the ability to find a near optimal solution by updating the candidate solutions iteratively and have self-organizing behavior.

The clustering algorithm designed is based on the behavior of the bees but was improved by incorporating K-means and Consensus so that the algorithm adjusts itself in each iteration. A series of experiments were conducted with the hybrid ABC algorithm, using optimal values for each adjustable parameter. The experiments were carried out with three collections with 50, 647, and 1027 Web services, and the algorithm was executed with variants of 100, 200, and 500 iterations. The hybrid ABC algorithm has shown good results for Web services clustering.

As future work, more combinations of semantic similarities, as well as incorporating more information of the description of the services and the data types incorporated in the XML service definitions.

Also the incorporation of other bio-inspired algorithms (swarm intelligence) for the classification and clustering of Web services, such as Particle Swarm Optimization (PSO), Bat Algorithm (BA), and Bird Swarm Algorithm (BSA).

References

[1] Liang Q, Li P, Hung PC, Wu X. Clustering web services for automatic categorization. In: SCC'09: IEEE International Conference on Services Computing. IEEE; 2009. pp. 380-387

DOI: http://dx.doi.org/10.5772/intechopen.85200

Bio-Inspired Hybrid Algorithm for Web Services Clustering

[9] Fellbaum C. WordNet: An Electronic Lexical Database. Cambridge, MA: MIT

disambiguation using machine readable dictionaries: How to tell a pine cone from an ice cream cone. In: Proceedings

[11] Banerjee S, Pedersen T. An adapted

[12] Wu Z, Palmer M. Verbs semantics and lexical selection. In: Proceedings of the 32nd Annual Meeting on Association

[13] Resnik P. Using information content to evaluate semantic similarity. In: Proceedings of the 14th International

[10] Lesk M. Automatic sense

Lesk algorithm for word sense disambiguation using WordNet. In: Computational Linguistics and Intelligent Text Processing. Berlin Heidelberg: Springer; 2002. pp. 136-145

for Computational Linguistics. Association for Computational Linguistics; 1994. pp. 133-138

Joint Conference on Artificial Intelligence. Montreal; 1995.

[14] Jiang JJ, Conrath DW. Semantic similarity based on corpus statistics and lexical taxonomy. In: Proceedings of International Conference on Research in Computational Linguistics, Taiwan. 1997

[15] Lin D. An information-theoretic definition of similarity. In: Proceedings of the 15th International Conf. on Machine Learning. San Francisco, CA: Morgan Kaufmann; 1998. pp. 296-304

[16] Hirst G, St-Onge D. Lexical Chains as Representations of Context for the

Combining local context and WordNet

Detection and Correction of Malapropisms; 1998. pp. 305-332

[17] Leacock C, Chodorow M.

pp. 448-453

of SIGDOC '86. 1986

Press; 1998

[2] Platzer C, Rosenberg F, Dustdar S.

multidimensional angles as proximity measures. ACM Transactions on Internet Technology. 2009;9(3):11

Computation: Practice and Experience.

[5] Wu J, Chen L, Zheng Z, Lyu MR, Wu Z. Clustering web services to facilitate service discovery. Knowledge and Information Systems. 2014;38(1):

[6] Prakash J, Singh PK. Evolutionary and swarm intelligence methods for partitional hard clustering. In: 2014 International Conference on

Information Technology (ICIT). IEEE;

[7] Sahoo G. A two-step artificial bee colony algorithm for clustering. Neural Computing and Applications. 2017;

[8] Kotekar S, Kamath SS. Enhancing web service discovery using metaheuristic CSO and PCA based clustering. In: Progress in Intelligent Computing Techniques: Theory, Practice, and Applications. Singapore: Springer; 2018.

[4] Du YY, Zhang YJ, Zhang XL. A semantic approach of service clustering and web service discovery. Information Technology Journal. 2013;12(5):967-974

Web service clustering using

[3] Pop CB, Chifu VR, Salomie I, Dinsoreanu M, David T, Acretoaie V, et al. Biologically-inspired clustering of semantic web services. Birds or ants intelligence? Concurrency and

2012;24(6):619-633

207-229

2014. pp. 264-269

28(3):537-551

pp. 393-403

25