Nils Petermann

Niklas Stein

Gabi Schierning

Ralf Theissmann

Benedikt Stoib

Christian Hecht

Christof Schulz

Hartmut Wiggers
email: hartmut.wiggers@uni-due.de

Plasma synthesis of nanostructures for improved thermoelectric properties

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The continuous and increasing discussion regarding possibilities for energy harvesting, power efficiency and sustainability has led to an intensive (re)search for further development and optimization of energy conversion processes. This covers all areas from improvement in efficiency over improved insulation techniques to waste heat recovery. The latter one becomes more and more interesting as new approaches in thermoelectrics pave the way for its commercialization.

Thermoelectric materials are solids which allow for a direct conversion of temperature differences into electrical energy. Thermoelectric generators have proven to show long-term stability and reliability even under extreme conditions but up to now they are mainly used in niche applications like space probes [START_REF] Bennet | presented at the 4th International Energy Conversion Engineering Conference and Exhibit (IECEC)[END_REF] or at remote locations [START_REF] Congress | [END_REF] . Despite of still relatively low over-all performance with respect to costs and efficiency, thermoelectrics are highly interesting for heat recovery, especially at hightemperatures such as waste heat from processing plants or from automobile exhaust, containing about 30% of the overall energy delivered from the fuel [START_REF] Heywood | Internal Combustion Engine Fundamentals[END_REF] . Therefore, thermoelectric generators mounted at the exhaust pipe of combustion engines have the potential to significantly improve the fuel economy through conversion of waste heat to electric energy which is increasingly required to drive the electrical components of modern cars and will be even more important as an energy source in hybrid cars.

Thermoelectrics

Thermoelectric generators are operated based on the Seebeck effect. It describes the emergence of a thermo voltage, when the ends of an electrically conductive material are exposed to different temperatures. Equation 1S is the Seebeck-coefficient or thermo power, U the thermo voltage, and ∆T the temperature difference. The origin of the thermo voltage is a thermally driven diffusion of charge carriers from the hot to the cold end of the chunk. Thus, an electrical field is generated with the polarity of the majority charge carriers at the cold side. Hence, semiconductors can be used to transform temperature gradients into usable electrical energy. The efficiency of a thermoelectric generator is limited to a fraction of the Carnot efficiency, determined by the dimensionless figure of merit ZT defined as Equation 2where σ is the electrical conductivity, while κ el and κ lattice are the electronic and the lattice part of the thermal conductivity, respectively. The numerator S 2 σ is called the power factor.

For decades, the main approach towards advanced thermoelectric materials was focused on the development of degenerate semiconductor materials like Bi 2 Te 3 , Sb 2 Te3, PbTe, and alloys such as . Si x Ge 1-x . This led to reasonable values for ZT with peak values typically in the range of 0.8 to 1.1.

From the mid 1990s, efforts towards an improved performance of thermoelectric materials have concentrated to reduce thermal conductivity. Note that the electron thermal conductivity is proportional to the electrical conductivity, so ZT can not be increased by decreasing the electron thermal conductivity

The lattice thermal conductivity -responsible for heat transport by phonons/lattice vibrations -is directly related to the mean free path of phonons. From detailed microscopic investigations it was found that high efficiency thermoelectrics feature structures in the nanometer range. As a result, scientists started to further enhance the efficiency of thermoelectric materials by introduction of artificial nanostructures to synthesize materials with a significantly lowered thermal conductivity.

According to ISI Web of Science, the number of both, publications and citations with keywords nano* and thermoelectric* have experienced an impressive boom. The number of publications rose from 8 in 1996 to 290 in 2009 and the number of citations from 22 in 1996 to 4093 in 2009.

In this contribution, we focus on the formation of artificial nanostructures and grain boundaries as a highly promising and powerful tool to introduce enhanced phonon scattering mechanisms into the material. The goal is to lower its thermal conductivity while simultaneously preventing a decrease in electronic conductivity. Hence, the structure size has to be adapted between the electron wavelength and the phonon mean free path. While the electron mean free path is in the order of some nanometers, typical phonon mean free paths are in the range of up to several 100 nm. Especially low-frequency acoustic phonons, which carry most heat due to their high group velocity, show mean free paths up to some microns [START_REF] Singh | ASME/Pacific Rim Technical Conference and Exhibition on Packaging and Integration of Electronic and Photonic Systems, MEMS and NEMS (InterPACK 09)[END_REF] . Thus, tailoring nanostructures in an appropriate size can effectively scatter phonons while having less impact on the charge carrier transport properties.

Nanometer-scaled materials may be formed by compacting nanoparticles under preservation of the nanostructure which is the approach of this work. Other possibilities include intrinsic structural modulations at the nanometer scale [START_REF] Hsu | [END_REF]6 , the incorporation of precipitates and inclusions [START_REF] Hsu | [END_REF]7 , or the use of other obstacles for phonons like structural domains 8 . Further strategies to tune the thermal conductivity of thermoelectric materials are the incorporation of rattling atoms with large vibrational amplitudes at partially filled structure sites 9,10 and the transition to complex crystal structures where the number of optical phonon modes is increased and -as a consequence -the heat transport by the acoustic phonon branches reduced 11 . Alloying produces disorder at the atomic scale so that the thermal conductivity is at the 'alloy limit' 12 and may be combined with nano structuring in order to affect the complete phonon spectrum [13][14][15] . Remarkably high values of the figure of merit can be achieved in isolated nanostructures 16,17 and multilayer structures 18 which combine low values of the thermal conductivity and optimized electronic structure. For those multilayer structures, expensive fabrication techniques are employed like lithography or molecular beam epitaxy, and possible routes towards a large scale device fabrication are a long way off. Nevertheless, plasma assisted processes like plasma-enhanced chemical vapor deposition (PECVD), atomic layer deposition (PEALD), and sputtering are highly promising techniques that might also be suitable for commercial multilayer synthesis.

With respect to waste heat harvesting from combustion engines, the materials of interest have to fulfill several requirements: i) the thermal stability of possible materials has to cover a temperature range from 300 to 900 °C; ii) the materials should consist of abundant and -preferably -non-toxic materials; iii) the devices should consist of p-and n-type doped materials with similar thermal expansion coefficient for high thermal stability. According to 19 the most promising materials in the required temperature range are CoSb 3 , (GeTe) 0.85 (AgSbTe2) 0.15 (TAGS), CeFe 4 Sb 12 , Yb 14 MnSb 11 , and Si x Ge 1-x . To our opinion, the most promising high temperature material with respect to the requirements mentioned above is Si x Ge 1-x and it has been shown, that nanostructuring of Si and Si x Ge 1-x is a successful way towards enhanced thermolectric properties 20,21 .

The thermoelectric performance of bulk silicon is poor, as it exhibits diamond structure and therefore, its bulk thermal conductivity is quite high (about 140 W m -1 K -1 at ambient temperature 22 ). The high value of κ limits the efficiency of bulk silicon as a thermoelectric material, but this material will become an interesting candidate as a basis for widely-spread thermoelectric applications if its thermal conductivity could be significantly reduced as described above.

Decades ago, Dismukes et al. reported a thermal conductivity of about 5 W m -1 K -1 for n-type Si 0.8 Ge 0.2 alloys at room temperature with a peak ZT = 1 for the n-type alloy at 1200°C 23 . Alloys of silicon-germanium have proven to show a much lower thermal conductivity than pure silicon, which is related to alloy scattering of phonons, resulting in a much better thermoelectric performance 24 .

Microstructured SiGe alloys are in use in radioisotope thermoelectric generators (RTG) to power several space probes [START_REF] Bennet | presented at the 4th International Energy Conversion Engineering Conference and Exhibit (IECEC)[END_REF] . Alloying silicon with germanium is thus one way to lower thermal conductivity and to achieve good ZT values, but germanium is a rare and expensive material. Nevertheless, using germanium in low concentration could be an optimum with respect to costs and efficiency. Zhu et al. 14 reported about nanostructured Si 0.95 Ge 0.05 alloys with a significantly enhanced power factor, especially at room temperature, and a low thermal conductivity resulting in a similar figure of merit as found for the aforementioned Si 0.8 Ge 0.2 bulk alloys. Additionally, many other materials with promising thermoelectric properties consisting of highly available and sustainable compounds also contain silicon as a main component such as iron-, manganese-and magnesium silicides. Composite materials consisting of SiGe with embedded nanocrystalline silicides are discussed as high performance materials. In summary, nanostructuring of silicon-based materials has the potential to be used for high temperature thermoelectrics enabling for a reasonable and cost-efficient utilization.

Synthesis of nanosized materials

There are two main strategies to synthesize the required nanostructured materials, the top-down and the bottom-up approach. The top-down method starts from bulk materials by grinding 20,21 or highly specified methods like spark erosion or laser ablation 15,25 . Generally, these methods end up in nanostructured powder with a very broad particle size distribution and in most cases, these methods exhibit a specific lower limit concerning the particle size, which is in the range of a few ten nanometers. While the scattering rate of phonons is directly related to the number of grain boundaries, it is necessary to achieve a very small particle size to maximize the reduction of phonon thermal conductivity and methods other than top-down ones are required. As the formation of high performance thermoelectrics requires degenerate semiconductors with a very high content of dopants, these dopants must also be added as bulk material during formation of the nanosized powder, which might be difficult or even impossible for top-down methods.

The bottom-up approaches start from liquid or gaseous precursor materials 26,27 . These precursors react or decompose resulting in nucleation, coagulation and growth of nanoparticles. Depending on the reactor concept, thermodynamic as well as kinetic control of particle formation can be realized. As the time constants in liquid phase synthesis are usually much higher than in gas phase synthesis, liquid phase usually leads to thermodynamically controlled products, while gas phase synthesis enables for kinetic control. Moreover, liquid phase synthesis requires the addition of stabilizing, surface-active compounds while gas phase synthesis usually leads to high-purity products. Additionally, it is favored for the formation of non-oxidic nano materials because of its flexibility and scalability 28 . It is an established method for the formation of silicon nanoparticles 29,30 and has also been used recently for the formation of germanium nanoparticles 31 . High purity can be obtained as gases are available with a very low degree of impurities and particle synthesis can be carried out under high purity conditions. Control over stoichiometry, doping and residence time is attained as long as vaporizable precursor materials are available while mass flow controllers allow for a precise metering of the process gases. In gas phase synthesis, thermal decomposition and/or chemical reaction of the precursor forms a supersaturated vapor, which is thermodynamically less favored than a condensation and nucleation of single vapor atoms to clusters of several molecules. Depending on the characteristic times for collision and fusion, collision-limited and sintering-limited particles are formed, respectively 32 . Regarding the specific properties required for thermoelectric materials -small particles in the 10 nm rangesynthesis conditions have to be adapted in such a way that the characteristic time constants for nucleation, coalescence and sintering are much lower than the characteristic collision time. According to Tsantilis et al. 33 , a third parameter, the cooling rate downstream the reaction chamber, also has an important impact to particle morphology. Agglomeration and formation of aggregates ("hard agglomerates", see figure 2a) starts, when the sintering time exceeds the collision time 34 , which always will happen during the cooling process but can be suppressed at very high temperature gradients downstream the reaction zone leading to single, spherical particles with very small point contacts ("soft agglomerates", see figure 2b). As a result, the degree of agglomeration can be tuned by the adjustment of temperature and temperature gradient on the one hand and the adjustment of pressure and precursor concentration/vapor density on the other hand. Low precursor concentration a b

and low pressure in combination with high temperature and high temperature gradient result in softagglomerated particles, while high pressure and high concentration along with low temperature and shallow temperature gradient leads to hard agglomerates.

The morphology of raw materials used for the synthesis of bulk thermoelectrics is of high importance because of the fact that formation of dense bulk samples is required for high performance. The densification is done by sintering processes which always implicate restructuring and coarsening leading to bigger crystallites. This is in conflict with the requirement of small nanocrystals for low thermal conductivity. Especially materials with covalent bonds like silicon must be heated to very high temperatures just below their melting temperature for successful densification, which is always related to strong coarsening. Two main possibilities related to particle morphology help to overcome this disadvantage: i) reasonable content of very small nanoparticles sizing in the range of a few nanometer which show a significantly reduced melting temperature 35,36 and ii) addition of small amounts of homogeneously distributed sintering aids 37 . While the first one can be tuned by the particle formation process itself, the latter one is also directly associated to the particle size as it is strongly surfacerelated.

Plasma synthesis of nanoparticles

A versatile tool for synthesis of soft-agglomerated nanoparticles of all kinds is the plasma process, as it can provide high reactivity, high temperature, high temperature gradients and also high purity. It enables for production of all oxide and non oxide materials, for which respective gaseous or liquid precursors are available. Most important, it allows for a synthesis of nano materials in a continuous process with kinetic control of the particle formation, which makes thermodynamically instable materials accessible.

Plasma synthesis for thermoelectric materials is a developing technique [38][39][40][41] due to the fact that it enables for the formation of highly specific particles in the 10 nm range. A lot of investigation towards the formation of silicon and germanium nanoparticles within a plasma is carried out from the Kortshagen group 31,42,43 . They utilize a 200 W radio frequency (13.56 MHz) generator to ignite a capacitively coupled plasma within a quartz glass tube operated at a pressure of a few mbar resulting in non-thermal plasma. Depending on precursor concentration, mass flow and pressure, crystalline particles from a few to a few ten nm are available. Gorla et al. 44 investigated an argon/silane inductively coupled plasma also operated at 13.56 MHz at a maximum power of 1200 W. This type of non-thermal plasma can be used to either form silicon or germanium thin films by plasma-assisted chemical vapor deposition or to synthesize spherical silicon and germanium nanoparticles, depending on operation conditions. Nozaki et al. 45 present a microplasma reactor based on a very-highfrequency (VHF, 144 MHz) process similar to the setup used by the Kortshagen group. The VHF generator is connected to two electrodes to ignite a non-thermal plasma within a glass tube. The apparatus is operated at atmospheric pressure with a power of 15 to 20 W. Silicon nanoparticles sizing from 3 to 15 nm in diameter can be produced from a mixture of argon, hydrogen and SiCl 4 .

Direct current (DC) arc discharges are widely used to ignite thermal plasma such as the one used by Rao et al. to synthesize silicon nanoparticles from silicon tetrachloride 46 . Due to thermal decomposition, the precursors induce a supersaturated vapor with subsequent particle formation. The gas/particle mixture is expanded through a nozzle into a water-cooled reactor vessel, quenching the particle formation process. Calorimetric energy balance measurements determined temperatures in the range of a few thousand K and due to the immediate cooling process, a quench rate greater than 5 x 10 6 K s -1 was found. Sophisticated sampling techniques have been developed for near real time inline-measurements of particle size distribution and it is found that nanoparticles sizing between a few and hundred nm in diameter can be achieved depending on pressure and precursor concentration 47 . DC arc discharge enables for the production of a few grams per hour of spherical, nanosized silicon particles, while the steep temperature gradient prevents them from aggregation. It is a well-known and stable technology that is widely used in industry, especially in gas phase processes. Nevertheless, the DC arc based plasma reactor has a strong disadvantage regarding the purity when used for the synthesis of nanoparticles. Due to erosion from the electrodes, traces of the electrode material are always present in the desired product.

With respect to high and controlled purity of materials used for thermoelectrics, we favor a contamination-free reactor design. It should combine the possibilities for high purity synthesis with a high production rate. To this goal, we have implemented a technology based on a microwave-induced plasma. This technology enables to focus the microwave power into a very small volume within a quartz glass tube. The plasma can be stabilized in the center of the tube and reaches gas temperatures considerably higher than 1000 K. This is done by focusing the microwave within a cavity, resulting in a small spot inside a quartz glass tube by constructive interference of the incident waves. The energy transferred from the oscillating electrical field to charged particles is inverse proportional to their mass and to the square of the frequency 48 . Hence, free electrons are accelerated by the electrical field of the microwave, whereas the mass of the heavy ions greatly reduces their acceleration. Thus, electrons with their small mass have considerably more kinetic energy than the heavy particles. This makes microwave plasmas non-thermal plasmas, but especially at higher pressure above a few mbar an elevated gas temperature is found due to collisions of accelerated electrons with neutral gas atoms, molecules and their respective ions, transferring a part of their energy.

Simulation of nanoparticle formation and growth within plasmas is a challenging but increasingly used and highly promising tool for a better understanding of the growth process. The plasma energy coupled into a reactive flow creates very complicated non-equilibrium conditions resulting in various difficulties for modeling attempts. The gas phase is partially ionized, i.e. free electrons and ions are present in the system. The heavy and light species can have different temperatures and their detailed reaction kinetics and diffusional behavior is not known. We have developed a very simple model ignoring all the difficulties described above and considering the absorption of microwave energy as a heat source in the flow resulting in an associated temperature increase 49 . This assumption led to reasonable results for the initial particle formation. A very similar attempt was chosen from Shigeta et al. also neglecting different temperatures and charges 50 . They developed a more suitable model for a numerical analysis of the entire growth of binary nanoparticles in thermal plasma synthesis including homogeneous nucleation and heterogeneous condensation processes. Therefore, it is possible to model nanopowders of any composition and particle size distribution also taking into account the freezing point depression of the nanoparticles. For verification of their model, they made a comparison of the calculation and TEM analysis of nanoparticles from different experiments on silicon-based materials 51, 52 and found a surprisingly good agreement between model and experiment. Despite these encouraging results, up to now there is no suitable model taking into account the plasma-specific properties mentioned above, not to talk of much more complicated topics such as plasma chemistry.

Therefore, the detailed modelin thermal plasmas, is still a challen

Microwave plasma reactor

The Note that the plasma riginates from the hot F) measurements have ve plasmas 53,54 . It is orescence 55 . For the LIF) thermometry is a formation of silicon nanoparticles (see Figure 2). As can be seen from NO-LIF measurements performed to the afterglow of a microwave plasma by means of an excimer laser and an intensified charge coupled device (ICCD) camera, temperatures in the range of 1000 to 2000 K are found in agreement with results from the literature.

Figure 3: NO-LIF temperature measurements behind the microwave antenna exit at a microwave power of 500 W and a pressure of 60 mbar (measurements were performed without particle formation)

As in most flow reactor setups, the carrier gases argon and hydrogen as well as the precursor silane are injected through a central nozzle into a quartz glass tube usually leading to some inhomogenities in the gas temperature as can be seen from Figure 3.

In-situ doping of nanoparticles

For electronic as well as thermoelectric applications, doping of nanoparticles is essential to tailor electrical properties. Although it can be assumed, that solubility of dopants in nanoparticles remain the same as in a bulk material 57,58 , doping of nanoparticles is not trivial. Dopant incorporation, location and efficiency depend on the specific properties of the elements chosen. N-type doping of siliconbased nanoparticles in gas phase synthesis can be realized by adding phosphine (PH 3 ) to the gas mixture, while p-type doping can be achieved by adding diborane (B 2 H 6 ). The maximum solubility of phosphorous in Si is 2.4 at.% at 1180°C 59 , while the solubility of boron in silicon is found to vary between 0.21 and 3 at.% 60 .

Complete decomposition of the respective dopant precursor during plasma synthesis is observed and quantitative incorporation within the nanocrystalline powder is found 58,61 . Nevertheless, only a fraction of the added dopant was homogeneously embedded within the crystalline particles while a reasonable amount segregates at the particle surface. Up to now there is no final understanding concerning the incorporation of dopants at high levels much above those usually used for semiconductor application. Self-purification -assisted by short diffusion lengths within the nanoparticles -larger formation energies for substitutional incorporation, strain and mismatching size are some of the points under discussion 57,58,[61][62][63][64] . However, all investigations dealing with the localization of dopants in nanosized particles from gas phase synthesis find that the mean concentration matches the intended value, while fluctuations are found at the scale of the nanoparticles. Therefore, an almost homogeneous distribution of dopants can be assumed from plasma synthesis of small, nanoscale materials. This finding might be due to the fact, that plasma processes are very fast and usually performed under kinetic control.

Plasma-made nanoparticles

All plasma gas phase synthesis methods described above are capable to produce mostly spherical and soft-agglomerated nanoparticles. As mentioned above, this is the preferred morphology with respect to further processing by sintering. A typical image of an as-prepared material can be seen from Figure 1b. Despite the utilization of low melting point according to small particle size, we also want to take advantage of the often unwanted formation of a native oxide on the silicon surface. It is a well-known fact, that almost all materials when handled in air show the occurrence of chemically bond hydroxyl groups and oxides on their surface. Due to their high surface-to-volume ratio, this is an important process in surface termination of nanomaterials and can result in a double-digit percentage of concentration with respect to the main composition of the nanomaterial. Silicon when kept in air usually forms a silicon oxide layer with a thickness of 1-2 nm. Figure 4 confirms this finding as it clearly indicates that the particle surface is covered with a thin silicon oxide while the core consists of crystalline silicon. It is also found that areas with an initial point contact do not oxidize as much as the free surface and the native oxide forms a common envelope.

Thermoelectrics from plasma-made Si-NPs

The formation of bulk samples for thermoelectric applications can be done by densification of the doped silicon nanoparticles obtained from the gas phase plasma process. After precompaction to a density of about 50% and subsequent field assisted sintering 65 , dense pellets can be obtained. As the maximum voltage drop is expected to occur at the grain boundaries between the interconnecting particles, heating and subsequent sintering starts at the particle/particle contacts. Typical heating rates are in the range of 100 K min -1 and densification can be achieved within minutes. The density of the resulting pellets was found to vary between 95% and 98% of the bulk density of crystalline bulk silicon. The field assisted sintering process takes advantage of two materials properties that are directly related to the particle formation process. i) As is known from the literature the sintering temperature of nanosized silicon is dramatically reduced compared to bulk silicon 35,36 . It was observed that small silicon nanoparticles from the gas-phase sizing 20 nm and below partially melt and recrystallize immediately starting from about 800 °C while forming slightly bigger crystallites compared to their initial size. This process is especially observed for the very small crystallites that are "eaten" by the bigger crystallites. ii) As described above, plasma synthesis of the silicon nanoparticles results in spherical particles with weak point contacts due to van der Waals forces leading to a softagglomerated powder. The handling of these materials under ambient conditions results in a thin native oxide shell surrounding the crystalline silicon core. During sintering, this native oxide strongly supports the densification of the pressed pellets at surprisingly low temperatures due to a viscous flow sintering of the oxide [START_REF] Brinker | Sol-gel science: the physics and chemistry of sol-gel processing[END_REF] . The silicon oxide precipitates while a three-dimensional interconnecting silicon network is formed as can be seen from Figure 5. X-ray diffraction of as-prepared and sintered materials (not shown) reveals that the mean crystallite diameter is only little affected during sintering and remains almost constant. As a result, almost dense pellets are received at a temperature 350 °C below the melting temperature of silicon which is surprising all the more considering the fact, that silicon shows very poor sintering due to the covalent nature of its bonds. Additionally, we are able to synthesize thermoelectrics by keeping the initial crystallite size that can be adjusted by the plasma process parameters.

Thermoelectric properties

The samples used for the characterization of the thermoelectric properties were made from silicon nanoparticles with a nominal dopant concentration of 5*10 20 cm -3 . The Seebeck coefficient was measured by a reference method from MMR Technologies (Micro Miniature Refrigerators) using a constantan wire as reference material while the specific electrical conductivity used to calculate the charge carrier density was determined in van der Pauw geometry 67 . Thermal conductivity was investigated using a laser flash method with LFA 457 MicroFlash from NETZSCH Gerätebau GmbH, while ZT was measured applying a modified Harman method 68 .

The investigations concerning the characterization of charge carrier density reveal that after sintering about one third of the phosphorous atoms is electrically activated within the silicon, while the main amount is found to be accumulated in the silicon oxide forming a phosphorous glass. For all samples investigated, typical values found for the Seebeck coefficient are around 100 to 150 µV K -1 , and the electrical conductivity was measured between 200 and 900 S cm -1 , depending on the density of the samples. This indicates that the Seebeck coefficient and electrical conductivity are similar to literature values found for single crystalline silicon and are only weakly affected by the nanostructure of the samples as expected. In contrast, the thermal conductivity of our nanostructured samples is found to be typically between 10 W m -1 K -1 and 20 W m -1 K -1 , and therewith about one order of magnitude lower than that of single crystalline bulk silicon and about five times lower compared to polycrystalline silicon, see red dots in Figure 6a. Being slightly above 10 W m -1 K -1 , the thermal conductivity of our samples is on the same level as undoped polycrystalline layer silicon shown in Figure 6a where a comparison of thermal conductivities of different types of silicon after a review by Cahill et al. is presented 69 . To derive values for the figure of merit ZT, a modified Harman method was applied to our field assisted sintered samples at temperatures between 300 and 700 K. As can be seen from Figure 6b, our samples can compete with samples presented by Bux et al. prepared from ball-milled silicon and phosphorous powder 21 . Besides using a different technique for sample preparation, they also used a more than tenfold higher dopant concentration at around 5*10 21 cm -3 within their samples, which may be the main origin of the difference in ZT between their and our samples.

How to proceed?

As has been shown, gas phase plasma processing of gaseous precursors towards soft-agglomerated nanopowders is a reasonable and successful way to nanostructured thermoelectrics and enables for a one-step formation of raw materials required for the production of thermoelectric devices.

Nevertheless, a combination of the technologies described here with established plasma powder processes such as restructuring approaches as described by Shigeta et al. are highly promising 52 .

Regarding the formation of nanocomposites e.g. consisting of SiGe with embedded nanocrystalline silicides as proposed by the group of Shakouri 13 , one might think about the combination of pure gas phase processes (for the formation of highly doped Si/SiGe) with plasma spray technologies based on the reconstruction of (micro) crystalline silicides leading to the respective composite. To our opinion, plasma spray technologies are highly favored for such an approach due to the fact that they have been adapted to a lot of different atmospheres including inert conditions. As nanoparticle formation as well as reconstruction are well-established plasma spray technologies operated at production scale [70][71][72] , combination with gas phase microwave plasma reactors in the 100 kW range is highly promising.

Nevertheless, a lot of technological challenges have to be overcome such as clogging, contamination and particle deposition due to thermophoresis.

Figure 1 :

 1 Figure 1: Nanocrystalline silicon from gas phase synthesis produced under thermodynamic control

  experimental setup used for Cyrannus plasma source produc generator. Typical operation para 0.01 to 0.150 slm of silane (SiH output from 0.7 to 10 g h -1 of sof be tuned between 4 and 50 nm by

Figure 2 :

 2 Figure 2: a) Schematic sketch of

Figure 4 :

 4 Figure 4: High resolution energy filtered transmission electron microscopy of silicon nanoparticles

Figure 5 :

 5 Figure 5: Energy filtered TEM of nanocrystalline silicon, sintered at 1060 °C making use of the

Figure 6 :

 6 Figure 6: Comparison of the thermoelectric properties of our field assisted sintered nanocrystalline

Acknowledgements

Financial support by the German Research Foundation (DFG) within the Priority Program on nanoscaled thermoelectric materials, SPP 1386, and within the Collaborative Research Center on nanoparticles from the gas phase, SFB 445, is gratefully acknowledged. Financial support by the European Union and the Ministry of Economic Affairs and Energy of the State North Rhine-Westphalia in Germany in the frame of an Objective 2 Programme (European Regional Development Fund, ERDF) and a Young Investigator Grant is gratefully acknowledged.


Peter Mason

Amandine Aftalion

A classification of the ground states and topological defects in a rotating two-component Bose-Einstein condensate

Keywords:

We classify the ground states and topological defects of a rotating two-component condensate when varying several parameters: the intracomponent coupling strengths, the intercomponent coupling strength and the particle numbers. No restriction is placed on the masses or trapping frequencies of the individual components. We present numerical phase diagrams which show the boundaries between the regions of coexistence, spatial separation and symmetry breaking. Defects such as triangular coreless vortex lattices, square coreless vortex lattices and giant skyrmions are classified. Various aspects of the phase diagrams are analytically justified thanks to a non-linear σ-model that describes the condensate in terms of the total density and a pseudo-spin representation.

I. Introduction

Bose Einstein condensates (BECs) provide an excellent environment to study experimentally and theoretically a rich variety of macroscopic quantum phenomena. In a rotating single component condensate, topological defects to the order parameter often manifest themselves as vortices that correspond to a zero of the order parameter with a circulation of the phase. When they get numerous, these vortices arrange themselves on a triangular lattice. In fact, vortices were first observed in twocomponent BEC's [1]. Since then two-component BECs and the topological excitations within have been experimentally realised in a number of configurations: a single isotope that is in two different hyperfine spin states [1][2][3][4][5][6][7], two different isotopes of the same atom [8] or isotopes of two different atoms [9][10][11].

When a two-component condensate is under rotation, topological defects of both order parameters are created which lead to more exotic defects such as singly or multiply quantised skyrmions [7,[12][13][14][15][16]. Analogy with the Skyrme model from particle physics is often invoked to represent the defects [START_REF] Manton | Topological Solitons[END_REF][START_REF] Cho | [END_REF]. The singly quantised skyrmions contain a vortex in one component which has the effect of creating a corresponding density peak in the other component. These singly quantised skyrmions are often referred to as coreless vortices. Once numerous, these vortices and peaks arrange themselves in (coreless vortex) lattices, that can be either triangular or square. Other defects such as vortex sheets or giant skyrmions can also be observed [19][20][21]. The aim of this paper is to classify the type of defects according to the different parameters of the problem restricting ourselves to twodimensions. While this paper will only concern itself with magnetically trapped two-component BEC's, there is ac-tive research into spinor condensates (for a recent review see [23]).

In the mean-field regime, the two-component Bose-Einstein condensate at zero temperature is described in terms of two wave functions (order parameters), Ψ 1 and Ψ 2 , respectively representing component-1 and component-2. The two-component condensate is placed into rotation about the z-axis with Ω = Ωẑ where Ω is the rotation frequency assumed to apply equally to both components. The Gross-Pitaevskii (GP) energy functional of the rotating two-component, two-dimensional BEC is then given by

E[Ψ 1 , Ψ 2 ] = k=1,2 2 2m k |∇Ψ k | 2 + V k (r)|Ψ k | 2 -ΩΨ * k L z Ψ k + U k 2 |Ψ k | 4 + U 12 |Ψ 1 | 2 |Ψ 2 | 2 d 2 r (1)
where r 2 = x 2 + y 2 , and V k (r) = m k ω 2 k r 2 /2 are the harmonic trapping potentials with trapping frequencies ω k , centered at the origin. Here m 1 and m 2 are the masses of the bosons in component-1 and component-2 respectively, and is Planck's constant. The angular momentum is in the z-axis and is defined as L z = i[ẑ •(r ×p)] for linear momentum p. The energy functional (1) contains three interaction constants: U k representing the internal interactions in component k, and U 12 representing the interactions between the two components.

The time-dependent coupled Gross-Pitaevskii (GP) equations are obtained from a variational procedure, i ∂Ψ k /∂t = δE/δΨ * k , on the energy functional (1). Let ω = (ω 1 + ω 2 )/2 be the average of the trapping frequencies of the two components and introduce the reduced mass m 12 such that m -1 12 = m -1 1 + m -1 2 .

The GP energy and the coupled GP equations can be non-dimensionalised by choosing ω-1 , ω and r 0 = /(2m 12 ω) as units of time, energy and length respectively. On defining the non-dimensional intracomponent coupling parameters g k = 2U k m 12 / 2 and the intercomponent coupling parameter g 12 = 2U 12 m 12 / 2 (≡ g 21 ), the dimensionless coupled GP equations read

i ∂ψ k ∂t = - m 12 m k (∇ -iA k ) 2 ψ k + m k 4m 12 ω 2 k ω2 -Ω 2 r 2 ψ k + g k |ψ k | 2 ψ k + g 12 |ψ 3-k | 2 ψ k , (2)
for Ω = Ω/ω and where

A k = 1 2 m k m 12 Ω × r, (3)
for Ω = (0, 0, Ω) and r = (x, y, 0). The ground state of the energy (1) or the GP equations ( 2) is determined by preserving the normalisation condition which in this paper is either taken to be

|ψ k | 2 d 2 r = N k , (4)
where N k is the total particle number of the kthcomponent, or

|ψ 1 | 2 + |ψ 2 | 2 d 2 r = N 1 + N 2 . (5)
The following sections will only consider repulsive interactions so that g 1 , g 2 and g 12 are always non-negative. In order to separate the regimes of interest, a nondimensional parameter combining these g k and g 12 will appear:

Γ 12 = 1 - g 2 12 g 1 g 2 . (6)
Furthermore to simplify matters, we introduce ratio parameters η and ξ such that m 1 = ηm 2 and ω 1 = ξω 2 with {η, ξ} > 0.

The effective trapping potentials for each component coming from equation ( 2) are then, respectively,

V eff 1 (r) =(η + 1)

ξ 2 (ξ + 1) 2 - 1 4 Ω 2 r 2 , (8a)
V eff 2 (r) = (η + 1) η

1 (ξ + 1) 2 - 1 4 Ω 2 r 2 . ( 8b
)
It follows that the limiting rotational frequency for each component is Ω lim 1 = 2ξ/(ξ + 1) and Ω lim 2 = 2/(ξ + 1), so that it is necessary to consider a rotational frequency such that 0 < Ω < Ω lim = min{Ω lim 1 , Ω lim 2 }.

Experimentally, it is often the case that the U k , m k , ω k and N k are of the same order, so that much of the theoretical and numerical analysis concerning two-component condensates has assumed equality of these parameters. In the case where the intracomponent coupling strengths are equal, Kasamatsu et al. [19,20] produced a numerical phase diagram in terms of the rotation and the intercomponent coupling strengths: they found phase separation regions with either vortex sheets or droplet behavior and regions of coexistence of the components with coreless vortices, arranged in triangular or square lattices.

In this paper, we present a classification of the ground states and various types of topological defect when the intracomponent coupling strengths are distinct. Depending on the magnitude of the various parameters the components can either coexist, spatially separate or exhibit symmetry breaking. A richer phenomenology of topological defects is then found.

Much of the analysis carried out to investigate the ground states and topological defects will use a a nonlinear sigma model. It has been introduced previously in the literature [19,24] for η = ξ = 1 but we will generalise this to the cases η and ξ different from 1. This involves writing the total density as

ρ T = |ψ 1 | 2 + η|ψ 2 | 2 . (9)
A normalised complex-valued spinor χ = [χ 1 , χ 2 ] T is introduced from which the wave functions are decomposed as

ψ 1 = √ ρ T χ 1 , ψ 2 = ρ T /ηχ 2 where |χ 1 | 2 + |χ 2 | 2 = 1.
Then we define the spin density S = χσχ where σ are the Pauli matrices. This gives the components of S as

S x =χ * 1 χ 2 + χ * 2 χ 1 , (10a) S y = -i(χ * 1 χ 2 -χ * 2 χ 1 ), (10b)
S z =|χ 1 | 2 -|χ 2 | 2 , (10c)
with |S| 2 = 1 everywhere. We will write the GP energy (1) in terms of ρ T and S that will allow us to derive information on the ground state of the system. The paper is organised as follows. The different regions of the Ω -Γ 12 phase diagrams are outlined in Sect. II with a detailed analysis of the range of ground states and topological defects. Then the non-linear σ-model is developed in Sect. III to analyse the ground states in terms of the total density. Lastly Sect. IV takes the second normalisation condition (5) and presents an example Ω -Γ 12 phase diagram.

II. The Γ12 -Ω Phase Diagram

A. Numerical Parameter Sets

In this paper, three different configurations are considered, two of which relate directly to experimental config-urations. Firstly, we analyze a 87 Rb-87 Rb mixture with one isotope in spin state |F = 2, m f = 1 and the other in state |1, -1 . The complex relative motions between these two isotope components of rubidium were considered experimentally by Hall et al [2]. Here, the masses are equal (η = 1) and the transverse trapping potentials are equal, ω 1 = ω 2 (ξ = 1). The scattering lengths for each component are a 1 = 53.35 Å and a 2 = 56.65 Å. The intracomponent coupling strengths U k for a twodimensional model are defined as

U k = √ 8π 2 a k /[m k a zk ] (k = 1, 2)
, and the intercomponent coupling strength is defined as

U 12 = √ 2π 2 a 12 /[m 12 ãz ]
, where m 12 is the reduced mass, given by m -1 12 = m -1 1 + m -1 2 , a k , a 12 are the s-wave (radial) scattering lengths, a zk is the characteristic size of the condensate in the z direction, and ãz = (a z1 + a z2 )/2. In the experiment [2], the computations of approximate values of a z1 and a z2 can be made using a Thomas-Fermi approximation in the z direction as explained in [25], chapter 17. This relies on the assumption that N k a k /a zk are large. Then, it is reasonable to describe the experiments through a two-dimensional model. For our simulations, we choose g 1 = 0.0078 and g 2 = 0.0083, which are consistent with the experimental data and we set N 1 = N 2 = 10 5 . We denote this parameter set as 'Experimental Set 1' (ES1).

Secondly, we tackle a 41 K-87 Rb mixture with both isotopes in spin state |2, 2 , as was considered by Modugno et al [9]. Component-1 is identified with the 41 K isotope and component-2 with the 87 Rb isotope. In the experiment, the masses are different (η ∼ 0.48), but since

ηξ 2 = 1, then m 1 ω 2 1 = m 2 ω 2
2 so that both components experience the same trapping potential. Here the scattering lengths for each component are a 1 = 31.75 Å and a 2 = 52.39 Å and we choose intracomponent coupling strengths as g 1 = 0.0067 and g 2 = 0.0063. This set is denoted ES2.

Lastly, in order to make an analogy with previous theoretical studies, we consider the set (ES3) which contains a mixture chosen such that all the parameter groups are equal, i.e. here g 1 = g 2 = 0.0078 and N 1 = N 2 = 10 5 with equality of the m k and ω k (η = ξ = 1).

In each case, the features of the ground states will be explained as g 12 and Ω are varied. While the value of the g k are nominally fixed, the value of the intercomponent coupling strength, g 12 , can be altered by Feshbach resonance (see for instance [26][27][28][29]), which allows for an extensive experimental range in Γ 12 , given by (6) (keeping Γ 12 ≤ 1). Simulations are thus performed on the coupled GP equations (2) in imaginary time. In general, it is difficult to find the true minimizing energy state. But the use of various initial data converging to the same final state allows us to state that the true ground state will be of the same pattern as the one we exhibit.

We would like to note that while we have chosen these particular parameters, we have conducted extensive numerical simulations over a range of different parameters sets and find the sets presented here illustrate well the possible ground states.

B. Classification of the Regions of the Phase Diagram

The ground states can be classified according to 1. the symmetry properties 2. the properties of coexistence of the condensates or spatial separation.

When there is no rotation (Ω = 0), as illustrated in Fig. 1, the geometry of the ground state can either be • two disks (coexistence of the components and symmetry preserving state)

• a disk and an annulus (symmetry preserving with spatial separation of the components), which depends strongly on the fact that the g k , m k , ω k are not equal

• droplets (symmetry breaking and spatial separation).

Under rotation, the different ground states can be classified according to the parameters Γ 12 and Ω. When the condensates are two disks, or a disk and annulus, defects may break some symmetry of the system as Ω is increased. But the wave function retains some non trivial rotational symmetry. We will not refer to this as symmetry breaking since the bulk condensate keeps some symmetry. We will use the terminology symmetry breaking when the bulk completely breaks the symmetry of the system and is a droplet or has vortex sheets. We find that there are four distinct regions, determined by the geometry of the ground state.

Region 1: Both components are disk shaped, with no spatial separation. Above some critical velocity Ω, coreless vortices appear: a vortex in one component which has the effect of creating a corresponding density peak in the other component. They arrange themselves either on a triangular or on a square lattice. Figure 2 provides a typical example of the density profiles.

Region 2: Vortex sheets. Under the effect of strong rotation, the components spatially separate and completely break the symmetry of the system. Nevertheless, they are approximately disk shaped with similar radii. Many vortices are nucleated that arrange themselves into rows that can have various patterns: either they can be striped or bent and are often disconnected from each other. A vortex sheet in one component corresponds to a region of macroscopic density in the other component. These features can be seen in the density plots of Fig. 3(c). FIG. 1: Density (divided by 10 4 ) plots along y = 0 for component-1 (dashed lines), component-2 (dotted lines) and the total density (solid lines) in which the components (a) coexist both being disks (Γ12 = 0.1), (b) spatially separate to be a disk and an annulus (Γ12 = 0) and (c) have symmetry broken to be two droplets (Γ12 = -0.3). The parameters are g1 = 0.0078, g2 = 0.0083 and N1 = N2 = 10 5 with η = ξ = 1 (set ES1). The angular velocity of rotation is Ω = 0. Distance is measured in units of r0 and density in units of r -2 0 .

Region 3: Spatial separation preserving some symmetry. Here one component is a disk while the other component is an annulus and the disk fits within the annulus with a boundary layer region in which both components have microscopically small density as illustrated in Figure 3(b) and 4(b). Under rotation, the topological defects can either be coreless vortex lattices in the disk and (but not always) corresponding isolated density peaks in the annulus, and/or a giant skyrmion at the boundary interface between the two components, as will be described later.

Region 4: Rotating droplets. The components spatially separate and have broken symmetry such that the centres of each condensate are different and each component contains a single patch of density as illustrated in Figures 3(a) and 4(a). In the rotating droplets, vortices can appear provided the rotation is greater than some critical value. These features can be seen in the density plots of Fig. 5.

The Ω -Γ 12 phase diagrams corresponding to the experimental parameter sets introduced above (sets ES1, ES2 and ES3) are shown in Fig.'s 6, 7 and 8. The last one is of similar type as the one reported by Kasamatsu et al [20] (there g 1 = g 2 = 4000 and

N 1 = N 2 = 1/2).
New features can be observed in Fig 6 and7, such as isolated density peaks that eventually vanish as Γ 12 is made more negative and the multiply quantised skyrmions at the interface between the disk component and the annular component in region 3. When {η, ξ} = 1 (Fig. 7), no region 2 is found to exist. This absence is easily explained by two factors: the onset of region 3 for (large) positive values of Γ 12 and the lack of vortices nucleated in component-1 (and to some extent in component-2).

We will now analyze in more detail some of the features of regions 1, 3 and 4 of the phase diagrams. Vortex sheets (region 2) have been analysed (mainly in the case of equal g k , m k and ω k ) in [30]. 6), in (b) the components have spatially separated (but keep some symmetry) and have isolated density peaks (region 3 of Fig. 6) and in (c) there are vortex sheets present (region 2 of Fig. 6). Distance is measured in units of r0 and density in units of r -2 0 .

C. Analysis of the Features of the Phase Diagrams

Square lattices

It is a specific property of two components to stabilize square lattices. A requirement for the existence of square lattices is the nucleation of many vortices in both components, something not permitted in ES2 (Fig. 7), where only a small number of vortices are ever nucleated in component 1.

Square lattices generally occur at high rotational velocities and examples have been observed experimentally [7] and numerically [19,20,31]. Mueller and Ho [32] and later Keçeli and Öktel [33] have analysed the transition from triangular to square lattices as Γ 12 is varied in two-component condensates when g 1 ∼ g 2 and Ω is such that the condensate is in the lowest Landau level (LLL). Providing Ω is large, according to [32,33] and with η = ξ = 1 (set ES1) and Γ12 = -1.3 (which gives g12 = 0.0122). The angular velocity of rotation is Ω and it takes the values (a) 0 and (b) 0.9. In (a) the components are rotating droplets (region 4 of Fig. 6) while in (b) the components have spatially separated (but keep some symmetry) and there are no isolated density peaks (region 3 of Fig. 6). Distance is measured in units of r0 and density in units of r -2 0 .

the value of α = √ 1 -Γ 12 determines whether the vortex lattices are triangular (0 < α < 0.172), in transition to becoming square (0.172 < α < 0.373) or square (α > 0.373). Therefore, assuming that, in the example parameters of g 1 = 0.0078, g 2 = 0.0083 and η = ξ = 1 (set ES1), the g k are sufficiently close enough, the square lattice should first appear at α = 0.37, or equivalently Γ 12 = 0.86. Comparing this to the numerically obtained values, where the square lattice first appears for 0.8 < Γ 12 < 0.9 in Fig. 6, the agreement seems good. It will be interesting to see how the analysis of [32,33] has to be modified when g 1 becomes distant from g 2 .

Nevertheless, the square lattices are present for lower rotational velocities at which the LLL is not necessarily valid. In a forthcoming paper, using the nonlinear sigma model presented below, we hope to derive a vortex energy in terms of the positions of vortices. The ground state of this energy indeed leads to square lattices in some range of parameters (see also Section III.E). It turns out that this vortex energy is similar to that of Barnett et al [34,35].

Rotating droplets

In the parameter range of region 4, symmetry breaking with spatial separation occurs. When the condensate is not under rotation, the ground states of the two components will be 'half ball'-like structures, as can be seen in Fig. 4(a) (if g 1 = g 2 and η = ξ = 1, then the two components are exactly half-balls). The difference between the intercomponent strengths introduces a curvature to the inner boundary of each component with the component that corresponds to larger g k having positive curvature and the other component having negative curvature. This curvature depends on both Γ 12 and Ω. If Γ 12 is held constant, then the curvature increases as Ω increases. When Γ 12 goes to -∞, the droplets approach half-balls.

Vortex nucleation is also seen in region 4; see Fig. 5. In this figure there are four vortices in component-1 and three vortices in component-2. The number of vortices in each component will increase as Ω is increased but this increase will also increase the curvature of the inner boundaries of the components, thus preventing the vortices aligning themselves into vortex sheets. Examples of the rotating droplet ground states can be seen in Fig. 4(a), and further examples can be found in [22,[36][37][38][39].

Spatial separation preserving some symmetry: disk plus annulus

Region 3 is defined by the Γ 12 and Ω in which the condensates have spatially separated components but still possess some symmetry about the origin: the disk component-1 sits securely within the annular component-2 (see Fig. inner edge of the annulus.

Let us describe the onset of region 3 from region 4, captured in Fig. 3 (sub-image (a) to (b) or (c) to (b)). For a particular Γ 12 , as Ω is increased, then the curvature increases to such an extent that the components develop constant non-zero curvature, and are identified as a disk The boundary between region 1 and region 3 is also calculated analytically by Eq. (39) (dotted line). For these parameters there is no region 2. The unit of rotation is ω.

and an annulus. Conversely, if Ω is held constant and Γ 12 is pushed to more negative values, then the curvature decreases.

Under rotation, defects can be observed: coreless vortices and giant skyrmions. Coreless vortices sit inside the disk while giant skyrmions are observed in the boundary layer.

)). The solid lines indicate the boundary between different identified regions (determined by the geometry of the ground state) and the dashed lines the boundary between triangular and square lattices. For these parameters there is no region 3. The unit of rotation is ω.

Coreless Vortices in the disk plus annulus case

In the geometry of disk plus annulus, the vortices in the inside disk have corresponding isolated density peaks inside the annulus, hence in the microscopic density regions, as illustrated in Fig. 3(b). Thus any vortex lattice structure produced in component-1 is replicated by a peak lattice structure in component-2, i.e. there is the continuing presence of coreless vortex lattices in the spatially separated condensates [40].

As the number of vortices in the disk increases with increasing angular velocity, the number of density peaks inside the central hole of the annulus likewise increases. This has the effect that the central hole of the annulus can be masked at high angular velocities by the multiple occurrence of the density peaks. A recent analytical understanding of the interaction between vortices and peaks has been provided by [41]. It may be extended to the case of the disk and annulus, for which the average density in the central hole is very small for one component.

Pushing Γ 12 to lower values has the effect of reducing the size of the boundary layer between the disk and annular components, but also the isolated density peaks disappear; see Fig. 4(b). A recent work has analysed, from an energy perspective, the preference for the ground state to contain or not contain density peaks [42].

Figure 9 plots the maximum peak density of the density peaks in component-2 as a function of Γ 12 for the parameters of Fig. 6 and with Ω = 0.5 and Ω = 0.75. The disjoint region when Ω = 0.75 for -0.4 < Γ 12 < 0.1 on Fig. 9 is due to the absence of density peaks in component-2 as a result of the appearance of vortex sheets (region 4). We see from Fig. 9, that the maximum peak density occurs when Γ 12 ∼ 0.01 which is the value at which the components begin to spatially separate (the transition between region 1 and region 3) in set ES1. For Γ 12 and Ω that take values outside of region 1, the maximum density of the peaks decreases linearly as Γ 12 increases. The maximum density approaches microscopically small values for Γ 12 ∼ -1.1 when Ω = 0.5 and for Γ 12 ∼ -1 when Ω = 0.75. An example of the transition can be seen in

Giant skyrmion

A boundary layer between the overlap of component-1 and component-2 is present for all values of Γ 12 in region 3 but reduces in width as Γ 12 becomes more negative (indeed the boundary layer disappears only as Γ 12 → -∞). There are additional topological defects that occur in the boundary layer that are not discernible on a traditional density plot. These topological defects' presence can be observed in a phase profile, however a better visualisation is to use the pseudospin representation ( 9)- (10). One can plot the functions S x , S y and/or S z which reveal the presence of all the topological defects -the coreless vortices (singly quantised skyrmions) that have already been visualised on the plots, and the new defect, a multiply quantised skyrmion.

The distinct nature of the two types of topological defect can be clearly seen in Fig. 10(a,b) where an S x plot and an (S x , S y ) vectorial plot are shown respectively. A density plot of each component for these same parame- ters is seen in Fig. 3(b). The coreless vortices evident in Fig. 3(b) are again clearly evident in Fig. 10(a) and (b). A blow-up of the region close to one of (the three) coreless vortices in the (S x , S y ) vectorial plot, centred at (1, -1) and exhibiting circular disgyration is shown in Fig. 10(c). The texture of S can also exhibit crossand radial-disgyration [43]. Conversely, the multiply quantised skyrmion, not present in the density plots of Fig. 3(b), can clearly be seen in both Fig. 10(a) and (b). The multiply quantised skyrmion forms a ring along the boundary layer. A blow-up around (0, -3.5) for the (S x , S y ) vectorial plot again shows the multiply quantised nature of this defect. A multiply quantised skyrmion was evidenced by [21] who termed it a giant skyrmion, a terminology retained in this paper. Increasing the rotation results in an increase in both the number of coreless vortices and in the multiplicity of the giant skyrmion. Ref [21] gave a relationship, m = n + 1, between the multiplicity of the giant skyrmion n and the total circulation m in the central hole of component-2. However the numerical simulations here suggest that a more appropriate correlation is provided by m = n + l where l is the number of isolated coreless vortices. As an example, consider Fig. 10 where Ω = 0.5 and n = 8, l = 3 and m = 11, satisfying m = n+l. It should be noted that the relation-ship m = n + l still holds as the isolated coreless density peaks disappear from the density profile of component-2; this is because the phase imprints due to the vortices in component-1, which constitute l, are still present.

III. A Non-Linear σ-Model

The ground states of the rotating two-component condensate can be recovered from a 'non-linear σ-model' formulation of the energy functional in terms of the total density ρ T and spin vector S (Eq.'s ( 9)-( 10)). We will analyze in detail the possible regimes in this section and show how a Thomas Fermi approximation can hold in the case Γ 12 > 0 or can be generalized and provide relevant information in the case Γ 12 < 0.

A. Energy Functional Representation

We write the energy functional of the two wave functions (Eq. ( 1), non-dimensionalised) as

E[ψ 1 , ψ 2 ] = E KE + E PE + E I where E KE = 1 η + 1 ∇ - i 2 (η + 1)Ω × r ψ 1 2 + η ∇ - i 2η (η + 1)Ω × r ψ 2 2 d 2 r, (11a)
E PE = 2j 1 r 2 |ψ 1 | 2 + 2j 2 r 2 |ψ 2 | 2 d 2 r, (11b)
E I = 1 2 g 1 |ψ 1 | 4 + 1 2 g 2 |ψ 2 | 4 + g 12 |ψ 1 | 2 |ψ 2 | 2 d 2 r. (11c)
with

j 1 = 1 2 (1 + η)ξ 2 (1 + ξ) 2 - 1 8 (1 + η)Ω 2 , (12a)
j 2 = 1 2 (1 + η) η(1 + ξ) 2 - 1 8η (1 + η)Ω 2 . ( 12b
)
It is assumed that the g k , m k , ω k and N k are distinct so that a weighted total density can be defined as (9) with ψ 1 = √ ρ T χ 1 and ψ 2 = ρ T /ηχ 2 . The spin density vector S, which satisfies |S| 2 = 1 everywhere, has components given by Eq.'s (10). We then have

|ψ 1 | 2 = 1 2 ρ T (1 + S z ), |ψ 2 | 2 = 1 2η ρ T (1 -S z ). ( 13
)
We introduce the phases θ 1 and θ 2 defined by χ 1 = |χ 1 | exp(iθ 1 ) and χ 2 = |χ 2 | exp(iθ 2 ). The energy functional is expressed in terms of 4 variables: the total density ρ T , the component S z and the angles θ 1 , θ 2 . We see

that E KE = E ρT + E Sz + E θ1,θ2
where

E ρT = 1 (η + 1) (∇ √ ρ T ) 2 d 2 r, (14a)
E Sz = 1 4 ρ T (η + 1) (∇S z ) 2 (1 -S 2 z ) d 2 r, (14b)
E θ1,θ2 = 1 2 ρ T (1 + η) (1 + S z ) ∇θ 1 - 1 2 (1 + η)Ω × r 2 + (1 -S z ) ∇θ 2 - 1 2η (1 + η)Ω × r 2 d 2 r. (14c)
The other terms of the energy straightforwardly become

E PE = [(j 1 + j 2 /η) + (j 1 -j 2 /η)S z ] r 2 ρ T d 2 r, (15a)
E I = ρ 2 T 2 c0 + c1 S z + c2 S 2 z d 2 r, ( 15b
) with c0 = 1 4η 2 (η 2 g 1 + g 2 + 2ηg 12 ), (16a) c1 = 1 2η 2 (η 2 g 1 -g 2 ), (16b) c2 = 1 4η 2 (η 2 g 1 + g 2 -2ηg 12 ). (16c)
Thus the complete energy is

E = 1 (η + 1) (∇ √ ρ T ) 2 + ρ T 4(η + 1) (∇S z ) 2 (1 -S 2 z ) + ρ T 2(1 + η) × (1 + S z ) ∇θ 1 - 1 2 (1 + η)Ω × r 2 + (1 -S z ) ∇θ 2 - 1 2η (1 + η)Ω × r 2 + [(j 1 + j 2 /η) + (j 1 -j 2 /η)S z ] r 2 ρ T + ρ 2 T 2 c0 + c1 S z + c2 S 2 z d 2 r. (17)
Energy ( 17) is subject to the constraints (4) that can be rewritten as

ρ T (1 -S z )/2 d 2 r =N 2 η, (18a)
ρ T (1 + S z )/2 d 2 r =N 1 . (18b)
or equivalently

ρ T d 2 r =N 1 + N 2 η, (19a)
ρ T S z d 2 r =N 1 -N 2 η, (19b)
The minimization of the energy under the constraints (19) yields two coupled equations with two Lagrange multipliers, µ and λ. We will write them when the phase is constant (the gradient of the phases θ 1 and θ 2 can be ignored) and the effect of rotation is contained in an effective centrifugal dilation of the trapping potential included in j 1 , j 2 :

µ + λS z = - 1 (1 + η) ∆ √ ρ T √ ρ T + 1 4(1 + η) (∇S z ) 2 (1 -S 2 z ) + [(j 1 + j 2 /η) + (j 1 -j 2 /η)S z ] r 2 + ρ T (c 0 + c1 S z + c2 S 2 z ), (20)
and

λ = - 1 4(1 + η) ∆S z + ∇ρ T • ∇S z (1 -S 2 z ) + 1 2(1 + η) S z (∇S z ) 2 (1 -S 2 z ) + (j 1 -j 2 /η)r 2 + (c 1 + 2c 2 S z ) ρ T /2. ( 21
)
As pointed out in [14], in the general case where η and ξ are different from 1, we have seen from the expression E θ1,θ2 that the effective velocity in each component is different. Nevertheless, in the case when η = ξ = 1, it is possible to define an effective velocity shared by both components,

v eff = ∇Θ 2 + RS z 2(1 -S 2 z ) (22)
where Θ = θ 1 + θ 2 and R = S y ∇S x -S x ∇S y . We note the identity

(∇S z ) 2 (1 -S 2 z ) = (∇S) 2 - R 2 (1 -S 2 z ) , (23)
where (∇S

) 2 = (∇S x ) 2 + (∇S y ) 2 + (∇S z ) 2 .
Expansion of the square in E θ1,θ2 , substituting in the v eff and using the identity from above reduces the energy to the simple form found in [19,24].

E = 1 2 (∇ √ ρ T ) 2 + ρ T 8 (∇S) 2 + ρ T 2 (v eff -Ω × r) 2 + 1 2 r 2 (1 -Ω 2 )ρ T + ρ 2 T 2 c 0 + c 1 S z + c 2 S 2 z d 2 r, (24)
where c 0 , c 1 and c 2 are equal to c0 , c1 and c2 with η set equal to unity.

B. A Thomas-Fermi approximation

The typical Thomas-Fermi (TF) approximation to Eq.'s ( 20)-( 21) is to assume that derivatives in ρ T and S z are negligible in front of the other terms. If we apply the TF approximation we then get

µ + λS z = [(j 1 + j 2 /η) + (j 1 -j 2 /η)S z ] r 2 + ρ T (c 0 + c1 S z + c2 S 2 z ), (25)
and

λ = (j 1 -j 2 /η)r 2 + 1 2 (c 1 + 2c 2 S z ) ρ T . ( 26
)
The TF energy is then

E TF = [(j 1 + j 2 /η) + (j 1 -j 2 /η)S z ] r 2 ρ T + ρ 2 T 2 c0 + c1 S z + c2 S 2 z -µρ T -λρ T S z d 2 r. ( 27
)
It is important to point out that the reduction of this quadratic form in ρ T and ρ T S z yields

E TF = c2 2 ρ T S z + c1 2c 2 ρ T + 1 c2 (j 1 -j 2 /η)r 2 - λ c2 2 + 1 2 c0 - c2 1 4c 2 × ρ T - 1 (c 0 -c2 1 /4c 2 ) (µ -(j 1 + j 2 /η)r 2 ) + 1 2 c1 (c 0 c2 -c2 1 /4) (λ -(j 1 + j 2 /η)r 2 ) 2 d 2 r + constant terms. ( 28
)
The existence of a minimum for this quadratic form implies that c2 ≥ 0 and c0 -c2

1 /(4c 2 ) ≥ 0. Since c0 -c 2 1 /(4c 2 ) = g 1 g 2 Γ 12 /(4η 2 c2
), this implies in particular that Γ 12 > 0.

Therefore, if Γ 12 > 0, (which implies c2 > 0), a Thomas Fermi approximation can be performed as has previously been considered, generally taking an approximation on the individual component wave functions ψ 1 and ψ 2 [22,36,[46][47][48][49][50][51]. Nevertheless, as we will see below, if c2 < 0 and Γ 12 < 0, then (c 0 -c2 1 4c2 ) > 0 and a TF approximation can still be performed on ρ T , provided we keep gradient terms in S z .

Multiplying (26) by S z and subtracting (25) leads to

µ = (j 1 + j 2 /η)r 2 + ρ T c1 2 S z + c0 . ( 29
)
Simplification by ρ T S z with (26) yields, when

ψ 1 × ψ 2 = 0, ρ T = a 3 + a 4 r 2 g 1 g 2 Γ 12 , (30)
and

S z = a 1 + a 2 r 2 a 3 + a 4 r 2 , (31)
where

a 1 = 4η 2 (λc 0 -µc 1 /2), (32a)
a 2 = ηg 1 h 2 -g 2 h 1 , (32b)
a 3 = 4η 2 (µc 2 -λc 1 /2), (32c) a 4 = -(ηg 1 h 2 + g 2 h 1 ), ( 32d
)
and where we define

h k = 2 j k - g 12 g 3-k j 3-k . ( 33
)
When only one component is present (ψ 1 × ψ 2 = 0), this simplifies to

ρ T = 1 g1 µ + λ -2j 1 r 2 if ψ 2 = 0 η g2 η(µ -λ) -2j 2 r 2 if ψ 1 = 0 ( 34
)
since S z = +1 when ψ 2 = 0 and S z = -1 when ψ 1 = 0.

To begin, we note from Eq. ( 31) that r 2 = (a 1a 3 S z )/(a 4 S z -a 2 ), so that

ρ T = a 1 a 4 -a 2 a 3 g 1 g 2 Γ 12 (a 4 S z -a 2 ) . ( 35
)
Since the intracomponent coupling strengths are chosen to be distinct, the support of each component will not necessarily be equal. In order to lead the computations, we have to assume a geometry for the components: two disks, a disk and annulus or two half-balls (droplets) and a sign for Γ 12 .

C. Γ12 > 0

Both components are disks

We start when both components are disks. Without loss of generality, we can assume that the outer boundary of component-2 (at r = R 2 ) is larger than that of component-1 (at r = R 1 ). Therefore, S z = -1 and (34) holds in the annulus, while (35) holds in the coexisting region, which is the inside disk. The integrals from Eq.'s [START_REF] Cho | [END_REF] then give

N 1 2π = R1 0 ρ T (1 + S z ) 2 rdr = - (a 1 a 4 -a 2 a 3 ) 2 4g 1 g 2 Γ 12 -1 s0 (1 + S z ) (a 4 S z -a 2 ) 3 dS z , (36a)
where s 0 = S z | r=0 = a 1 /a 3 and

ηN 2 2π = R1 0 ρ T (1 -S z ) 2 rdr + R2 R1 ρ T rdr = - (a 1 a 4 -a 2 a 3 ) 2 4g 1 g 2 Γ 12 -1 s0 (1 -S z ) (a 4 S z -a 2 ) 3 dS z - g 2 4ηj 2 0 ρ d ρ T dρ T , (36b)
with

ρ d = a 3 + a 4 r 2 g 1 g 2 Γ 12 r=R1 = a 2 a 3 -a 1 a 4 g 1 g 2 Γ 12 (a 2 + a 4 ) , (37)
since

S z = -1 at r = R 1 .
Completion of these integrals and noting that {a

2 + a 4 , a 2 -a 4 } = {-2g 2 h 1 , 2ηg 1 h 2 } gives a 2 3 = 8N 1 g 1 g 2 2 Γ 12 h 1 π(s 0 + 1) 2 (38a)
which must necessarily be positive, i.e. h 1 Γ 12 > 0, and

N 1 (a 2 -s 0 a 4 ) 2 =2η 2 N 2 g 1 g 2 j 2 h 1 Γ 12 (1 + s 0 ) 2 + 2N 1 (1 + s 0 )ηj 2 g 1 Γ 12 × [ηg 1 h 2 (1 + s 0 ) + 2g 2 h 1 (s 0 -1)] (38b)
where in Eq. (38b) the expression for a 3 from Eq. (38a) has been substituted. This equation can always be solved in terms of s 0 when h

1 Γ 12 is positive since the discrim- inant is equal to 8N 1 η 2 g 2 2 g 1 j 2 h 1 Γ 12 (N 1 g 12 + N 2 g 2 ). We find that 1+s 0 = 2N 1 g 2 N 1 g 2 -ηN 1 g 12 + 2N 1 η 2 g 1 j 2 Γ 12 (N 1 g 12 + N 2 g 2 )/h 1 .
A similar calculation can be completed if the outer boundary of component-2 is larger than that of component-1. This then gives that h 2 Γ 12 > 0 (although expressions (38) change slightly).

If an annulus develops in

component 3 -k, it means that s 0 = 1 (s 0 = -1) for k = 1 (k = 2). Inputting this choice into Eq. (38b) gives ḡ12 = N k g k j 3-k 2(N 1 j 1 + N 2 j 2 ) + 1 2 N k g k j 3-k N 1 j 1 + N 2 j 2 2 + 4N 3-k g 1 g 2 j 3-k N 1 j 1 + N 2 j 2 1/2 (39)
as the critical g 12 at which an annulus forms in component-{3 -k}. In the case of Fig. 7, the curve (39) has been plotted in dashed lines and is close to the numerical curve. Notice that if ξ 2 = 1 (equal trapping frequencies for both components), then Eq. ( 39) becomes independent of Ω, as in the phase diagrams of Fig. 6, where it yields Γ 12 = 0.008.

Let us point out that before the transition to the disk plus annulus takes place, there is a subregion of region 1, where there are 2 disks, but in one component the wave function has a local minimum at the origin. For instance, in the case of Fig. 6, it corresponds to Γ 1 changing sign.

As a conclusion, in order for the ground state to be composed of two disks, assuming that component-k is the component with smaller support, we need that h k > 0, Γ 12 > 0 and g 12 < ḡ12 . These three conditions can be summarised as

g 12 < min j k j 3-k g 3-k , j 3-k j k g k , ḡ12 , √ g 1 g 2 if h 1 , h 2 > 0, (40a)
j 3-k j k g k < g 12 < min j k j 3-k g 3-k , ḡ12 , √ g 1 g 2 if h 3-k < 0.
(40b)

A disk and an annular component

In this case, we can assume a disk in component-1 and an annulus in component-2; there are three regions: an inner disk where only component-1 is present and (34) holds, an outer annulus where only component-2 is present and (34) holds and an inner annulus where both components coexist and (35) holds.

In order to use the TF approximation when one component is annular, computations similar to (38a)-(38b) lead to h 1 h 2 < 0 and g 12 > ḡ12 . This can be summarized (for an annulus in component-{3 -k}) as

max j 3-k j k g k , ḡ12 < g 12 < min j k j 3-k g 3-k , √ g 1 g 2 , (41)
with h 3-k < 0. Equation ( 41) also places the restriction that

g k <g 3-k j k j 3-k 2 =g 3-k Λ k , (42)
where Λ k = (j k /j 3-k ) 2 . Note that Λ 1 Λ 2 = 1 such that an annulus develops in component-2 (-1) if g 2 > (<) g 1 Λ 2 .

Orders of Intracomponent Strengths and Special Cases

The effect that changing the order of the intracomponent strengths and particle numbers has on ḡ12 , and thus on the phase diagrams, is now investigated. There are two cases to consider, depending on the relative orders of the particle numbers. Drawing aid from experimental values, it is always expected that min{N 1 , N 2 } ≫ max{g 1 , g 2 }. Throughout it will be assumed that g 2 > Λ 2 g 1 , so the annulus develops in component-2 and that j 1 and j 2 are of order unity.

In the first case when the particle number of component-1 is much greater than the particle number of component-2 (N 1 ≫ N 2 ), it follows that ḡ12 ∼ g 1 j 2 /j 1 . The boundary between region 1 and region 3 is then directly dependent on the value of the ratio

j 2 g 1 /[j 1 g 2 ].
Notice that if g 2 ≫ g 1 , Γ 12 evaluated at g 12 = ḡ12 tends to unity and as such an annulus would always be present in component-2, whatever the value of g 12 .

Conversely in the second case when the particle number of component-1 is much smaller than the particle number of component-2 (N 1 ≪ N 2 ), it follows that ḡ12 ∼ √ g 1 g 2 which implies that the annulus will only develop near Γ 12 = 0.

We can also look at some special cases -there are four that can be considered:

Case (i). Λ k g 3-k = g k . When Λ 2 (equivalently Λ 1
) is such that this equality is made, the two components are both disks and no annulus develops.

Case (ii). η

= ξ = 1. Then Λ k = 1 and R 2 ≷ R 1 ⇐⇒ N 2 g 2 Γ 1 ≷ N 1 g 1 Γ 1
when there are two disks and the annulus develops in the component which has the larger interaction strength. Case (iii). η = ξ = 1 and g 1 = g 2 ≡ g. When the intracomponent coupling strengths are equal, ḡ12 = g and there are always two disks with

R 2 ≷ R 1 ⇐⇒ N 2 ≷ N 1 .
Case (iv). η = ξ = 1, g 1 = g 2 ≡ g and N 1 = N 2 ≡ N . When the particle numbers and intracomponent coupling strengths are equal, ḡ12 = g and there are always two disks with R 2 = R 1 . A detailed analysis of this case, explored numerically in a phase diagram for all Γ 12 and analytically in the TF limit, was considered by [20].

Justification of the Thomas-Fermi approximation

The gradient terms in ( 20)-( 21) can be neglected if d hl , the characteristic length of variation of ρ T and S z is much smaller than d c , the characteristic size of the condensates. We have that 1/d 2 hl is of order of µ and λ, which are of order √ N k g k , N k g 3-k . Hence d hl is bounded above by the maximum of (N k g k ) -1/4 , (N k g 3-k ) -1/4 . From the expression of a 3 , the characteristic size of the condensate is of order the minimum of (g k N k Γ 12 ) 1/4 . Therefore, the Thomas Fermi approximation holds if

N k g k √ Γ 12 , N k g 3-k
√ Γ 12 are large. This requires the usual Thomas Fermi criterion that N k g k , N k g 3-k are large, but breaks down if Γ 12 is too small. For Γ 12 = 0, the equations lead to spatial separation: either the radii of the disks tend to 0 in the case of two disks or the outer radius of the inner disk tends to the inner radius of the annulus in the case of disk plus annulus so that ψ 1 × ψ 2 = 0 everywhere. Another analysis has to be carried out to understand the region of coexistence, which is of small size and has strong gradients. D. Γ12 < 0, beyond the TF approximation For negative Γ 12 , and if N k g k are large, the Thomas Fermi approximation can be extended provided some model takes into account the small region where the condensates coexist.

Indeed, if we go back to (28), and have both c2 < 0 and Γ 12 < 0, then the coefficient in front of the second square is positive and the optimal situation is to have the square equal to 0, which leads to the inverted parabola (30). On the other hand, the coefficient in front of the first square is negative, and the ground state involves derivatives in S z to compensate it. Under the assumption that the boundary layer where S z varies is small, we are going to derive a TF model with jump for ρ T . We will analyze it for the different geometries (disk plus annulus, droplets and vortex sheets) and show that it provides information consistent with the numerics.

Disk Plus Annulus

Assuming that the boundary layer is present only at some r = r s , then S z = +1 in the region in which component-2 is zero (r ≤ r - s ) and S z = -1 in the region in which component-1 is zero (r ≥ r + s ). The transition from S z = +1 to S z = -1 is not smooth, therefore creating the jump in density. Therefore, we are lead to minimise the integral

I = Br s 2j 1 r 2 ρ T + g 1 2 ρ 2 T d 2 r + B R\rs 2 j 2 η r 2 ρ T + g 2 2η 2 ρ 2 T d 2 r (43)
with respect to r s . Here B rs is a ball of radius r s and B R\rs is a torus with outer boundary at r = R and inner boundary at r = r s . Thus

2πr s 2j 1 r 2 ρ T + g 1 2 ρ 2 T r=r - s -2πr s 2 j 2 η r 2 ρ T + g 2 2η 2 ρ 2 T r=r + s = 0, (44)
which implies that

ρ -= 1 g 1 µ 1 -2j 1 r 2 r ∈ B rs (45a) ρ + = η g 2 µ 2 -2j 2 r 2 r ∈ B R\rs . (45b)
Then using the normalisation conditions ρ -d 2 r = N 1 and ρ + d 2 r = ηN 2 we get an outer radius

R = r 2 s + g 2 N 2 πj 2 1/2 , (46a)
and chemical potentials

µ 1 = N 1 g 1 πr 2 s + j 1 r 2 s , (46b)
µ 2 =2j 2 r 2 s + g 2 N 2 πj 2 . ( 46c
)
It remains to find r s . But, from (44), it follows that

r 2 s = 4N2j2 π -g 1 N1 πr 2 s - j1r 2 s g1 2 4 j 1 N1 πr 2 s - j1r 2 s g1 -2j 2 N2j2 πg2 . ( 47
)
The above leaves a quartic equation in r 2 s , the solution of which can be found numerically.

We can then complete the energy calculation to give

I = g 1 N 2 1 2πr 2 s +(N 1 j 1 +2N 2 j 2 )r 2 s - πj 2 1 6g 1 r 6 s + 4 3 N 3 2 g 2 j 2 π 1/2 .
(48) A check plot of the density profile with the parameters g 1 = 0.0078, g 2 = 0.0083 and N 1 = N 2 = 10 5 where η = ξ = 1 agrees well with this model where the boundary is calculated to be at r s = 3.62 (see Fig. 11(a)).

Droplets

This formalism can also be extended to study the droplet case. As before, assuming a thin boundary layer, one can allow for a jump in density. We thus need to minimise the Thomas-Fermi energy for two droplets described by B 1 and B 2 that exist in the regions 0 ≤ θ ≤ α and α ≤ θ ≤ 2π respectively. The energy is then

I = α R1 0 2j 1 r 2 ρ T + g 1 2 ρ 2 T rdr + (2π -α) R2 0 2 j 2 η r 2 ρ T + g 2 2η 2 ρ 2 T rdr. (49)
The expressions for the density in each domain (45) and the normalisation conditions allow completion of this integral to give

I = 4 √ 2 3 N 3 1 g 1 j 1 α 1/2 + N 3 2 g 2 j 2 (2π -α) 1/2 . ( 50
)
It remains to find the optimum α. This is achieved through the condition dI/dα = 0 which gives

(N 3 1 g 1 j 1 ) 1/2 α 3/2 = (N 3 2 g 2 j 2 ) 1/2 (2π -α) 3/2 ⇒ α = 2π N (ḡ j) 1/3 (1 + N (ḡ j) 1/3 ) (51)
where we have set N = N 1 /N 2 , ḡ = g 1 /g 2 and j = j 1 /j 2 . As expected, equality of the N k , k and setting η = ξ = 1 (⇒ j 1 = j 2 ) gives α = π and the condensate is then composed of two half-balls. Otherwise, a curvature is present. A check plot of the density profile with the parameters g 1 = 0.0078, g 2 = 0.0083 and N 1 = N 2 = 10 5 where η = ξ = 1 is given in Fig 11(b).

Finally we can note the energy for the droplets is

I = 4 3 √ π 1 + N (ḡ j) 1/3 -1/2 × (N 3 1 g 1 j 1 ) 1/2 + (N 1 N 2 2 (g 1 g 2 2 j 1 j 2 2 ) 1/3 ) 1/2 . (52)
This energy can be compared to the energy of the disk plus annulus (48) to determine which is the optimum geometry. Indeed, in the numerical cases studied before, the droplets are preferred states for small Ω.

Regions of Vortex Sheets

In the case of vortex sheets, we can assume that the global profile of the total density is TF-like, obeying Eq. ( 30) in the bulk of the condensate. By working with the total density, we do not require any information on the vortex sheets themselves (and consequently S z ).

We thus take the form of ρ T from Eq. ( 30) from which we note that the outer boundary at r = R satisfies R = -a 3 /a 4 and that completion of the normalisation condition (19a) gives

a 3 = - 2(N 1 + ηN 2 )g 1 g 2 Γ 12 a 4 π . (53)
This expression evidently requires a 4 Γ 12 < 0. We however expect the vortex sheets to be present only in the Γ 12 < 0 domain, thus we can be more specific on the condition, namely that a 4 > 0, or

g 12 > ηg 1 j 2 + g 2 j 1 (ηj 1 + j 2 ) . (54)
We point out that this critical number corresponds to Γ 12 = 0 in the cases studied above. With a 3 as above and a 4 given by the parameters of the system, the density profile is then fully accessible.

E. Analysis of defects

The advantage of the Thomas Fermi analysis of ρ T in the nonlinear sigma model is that it allows for analysis of defects as a perturbation calculation when N k g k , N k g 3-k and N k g 12 are large, in the spirit of what has been done for a single condensate in [44,[52][53][54] or for two condensates in [45]. We do not need to analyze specifically the peaks because they are taken into account in the S z formulation. This has been attacked in a more complicated fashion in the case of a single coreless vortex in [19]. We will develop our ideas in a later work but let us point out that the starting point is the energy ( 17) and if we call p i and q i the location of the vortex in each component, the main terms coming from the vortex contributions lead to:

• Each vortex core of each component providing a kinetic energy term proportional to ρ T (0) log d hl /(2(1 + η)) where d hl is the healing length, the characteristic size of the vortex core, here of order 1/ √ N k g k .

• the rotational term providing a term in -cΩ(ρ T (0)) 2 , where c is a numerical constant. The balance of these two terms allows us to compute the critical velocity for the nucleation of the first vortex.

• the kinetic energy then yielding a term in -log |p ip j | -log |q i -q j |.

• the rotation term yielding a term in Ω(

|p i | 2 + |q i | 2 ).
• the interaction term providing a term ρ 2 T S 2 z for the perturbation of S z close to a coreless vortex. The ansatz can be made in several ways leading to an interaction term in e -|pi-qi| 2 .

This leads to a point energy of the type

i a(|p i | 2 +|q i | 2 +be -|pi-qi| 2 )- i,j (log |p i -p j |+log |q i -q j |)
where a and b are related to the parameters of the problem. The ground state of such a point energy leads to a square lattice for a sufficient number of points and for some range of a and b.

IV. Phase Diagram Under Conserved Total Particle Number

In the experiment of Hall et al [2] a single component BEC of 87 Rb in the |1,-1 state was initially created, before a transfer of any desired fraction of the atoms from this |1,-1 state to the |2,1 state created the two-component BEC. Thus the ratio of particle numbers N 1 /N 2 is controllable experimentally, with the constraint that N 1 + N 2 is constant (in the case of the experiment of [2], N 1 + N 2 = 5 × 10 5 ). As such, experimentally, it is possible to keep the individual particle numbers constant (as in the normalisation condition (4)) or to keep the total particle number N 1 + N 2 constant, allowing N 1 and N 2 to vary (as in (5)). We produce an Ω -Γ 12 phase diagram for the parameters g 1 = 0.003, g 2 = 0.006 and η = ξ = 1 using the normalisation condition given by FIG. 12: (Color online) Ω-Γ12 phase diagram for parameters g1 = 0.003 and g2 = 0.006 with η = ξ = 1 using the normalisation condition given by ( 5). (a) Numerical simulations where triangles indicate that the lattice in both components is triangular, squares that the lattice in both components is square and diamonds where no vortices have been nucleated. Filled triangles, squares and diamonds are where the two components are disk-shaped and coexist, empty triangles, squares and diamonds are where only component-1 exists; those with a dot in the centre represent the appearance of coreless vortices in component-2 and those without a dot in the centre represent the complete disappearance of component-2. (b) A schematic representation of the numerical simulations. The solid lines indicate the boundary between different identified regions (determined by the geometry of the ground state) and the dashed lines the boundary between triangular and square lattices. The unit of rotation is ω.

(5) with N 1 + N 2 = 2.1 × 10 5 . Note that if g 1 = g 2 , m 1 = m 2 and ω 1 = ω 2 (η = ξ = 1), then the phase diagram would be identical to that of Fig. 8 (i.e. the normalisation condition in this case would not be important). The Ω -Γ 12 phase diagram is presented in Fig. 12. There are three distinct regions (determined by the geometry of the ground state):

Region 1. In the first region, both components are diskshaped. As before, the coreless vortices can either form a triangular or a square lattice depending on the values of Γ 12 and Ω. Figure 13 shows this case for Γ 12 = 0.7 and Ω = 0.65 (where a triangular lattice is present) and The parameters are g1 = 0.003 and g2 = 0.006 with η = ξ = 1 and Γ12 = 0.7 (which gives g12 = 0.0023) and normalisation taken over the total density (Eq. ( 5)). The angular velocity of rotation is Ω and it takes the values (a) 0.65 and (b) 0.9. There is a triangular lattice in (a) and a square lattice in (b). At these parameters the components are in region 1 of the phase diagram Fig. 12. Distance is measured in units of r0 and density in units of r -2 0 . Ω = 0.9 (where a square lattice is present). In region 1, both components have the same radii. Region 2. In the second region, only component-1 exists except for isolated density peaks that exist in component-2. These isolated density peaks occur at the same location as the vortices do in component-1, and are thus identical to the isolated coreless vortices described in detail in Sect. IV. Figure 14 shows this case Γ 12 = 0.3 with Ω = 0.5 and Ω = 0.9.

Region 3. In the third region, only component-1 exists (the density peaks in component-2 that were present in region 2 are no longer present). Furthermore, only triangular vortex lattices are observed in this effective one component condensate (in component-1).

Computations similar to Section III hold, except that now we have to set λ = 0. This leads to a 1 /a 3 = -c 1 /(2c 2 ). We see that this ratio (which is S z (0)), reaches 1 or -1 when c1 = ±2c 2 . In our numerical case, this leads to Γ12 = 0.5. We see clearly 3 regimes: Γ 12 > Γ12 , where the condensates are 2 disks, Γ 12 < 0, which is phase separation, in which case S z = 1 is the preferred state and 0 < Γ 12 < Γ12 , in which case the TF approximation leads to a computation of ρ T with coreless vortex lattices and variations in S z which improve the energy and lead to this intermediate state, still to be The parameters are g1 = 0.003 and g2 = 0.006 with η = ξ = 1 and Γ12 = 0.3 (which gives g12 = 0.0035) and normalisation taken over the total density (Eq. ( 5)). The angular velocity of rotation is Ω and it takes the values (a) 0.5 and (b) 0.9. At these parameters the components are in region 2 of the phase diagram Fig. 12. Distance is measured in units of r0 and density in units of r -2 0 .

studied in more detail.

V. Conclusion

We have presented phase diagrams of rotating two component condensates in terms of the angular velocity Ω and a nondimensionnalized parameter related to the coupling strengths Γ 12 = 1 -g 2 12 /(g 1 g 2 ). We have analyzed the various ground states and topological defects and have found four sets characterized by the symmetry preserving/symmetry breaking, coexistence or spatial separation of the components. When the geometry of the ground states is either two disks (coexistence of components, region 1) or a disk and an annulus (spatial separation keeping some symmetry, region 3), the topological defects are coreless vortex lattices (with possible stabilization of the square lattice) or giant skyrmions at the boundary interface between the two components. In the complete symmetry breaking case, we have found vortex sheets and droplets. The difference of masses or coupling strengths between the components can induce very different patterns.

We have introduced an energy (17) related to the total density and a pseudo spin vector. The minimization in a generalized Thomas Fermi approximation provides a lot of information on the ground states for general masses, trapping frequencies and coupling strengths. Some parts of the phase diagrams can be justified rigorously, both in the case Γ 12 > 0, which had been studied before, but also in the case Γ 12 < 0 with generalized models. This formulation of the energy should bring in the future more information on the defects.

FIG. 2 :

 2 FIG. 2: (Color online) A series of density plots for component-1 (left column) and component-2 (right column) in which the components coexist and are both disks. The parameters are g1 = 0.0078, g2 = 0.0083 and N1 = N2 = 10 5 and with η = ξ = 1 (set ES1) and Γ12 = 0.5 (which gives g12 = 0.0057). The angular velocity of rotation is Ω and it takes the values (a) 0.25, (b) 0.5, and (c) 0.75. Vortices in one component have a corresponding density peak in the other component (coreless vortices). In (c) the coreless vortex lattice is square. At these parameters the components are in region 1 of the phase diagram of Fig. 6. Distance is measured in units of r0 and density in units of r -2 0 .

FIG. 3 :

 3 FIG.3:(Color online) A series of density plots for component-1 (left column) and component-2 (right column) in which the components have spatially separated. The parameters are g1 = 0.0078, g2 = 0.0083 and N1 = N2 = 10 5 and with η = ξ = 1 (set ES1) and Γ12 = -0.3 (which gives g12 = 0.0092). The angular velocity of rotation is Ω and it takes the values (a) 0.1, (b) 0.5, and (c) 0.9. In (a) the components are rotating droplets (region 4 of Fig.6), in (b) the components have spatially separated (but keep some symmetry) and have isolated density peaks (region 3 of Fig.6) and in (c) there are vortex sheets present (region 2 of Fig.6). Distance is measured in units of r0 and density in units of r -2 0 .

FIG. 4 :

 4 FIG.4: (Color online) A series of density plots for component-1 (left column) and component-2 (right column). The parameters are g1 = 0.0078, g2 = 0.0083 and N1 = N2 = 10 5 and with η = ξ = 1 (set ES1) and Γ12 = -1.3 (which gives g12 = 0.0122). The angular velocity of rotation is Ω and it takes the values (a) 0 and (b) 0.9. In (a) the components are rotating droplets (region 4 of Fig.6) while in (b) the components have spatially separated (but keep some symmetry) and there are no isolated density peaks (region 3 of Fig.6). Distance is measured in units of r0 and density in units of r -2 0 .

FIG. 5 :

 5 FIG. 5: (Color online) A density plot for component-1 (left column) and component-2 (right column) showing examples of vortex nucleation in rotating droplets. The parameters are g1 = 0.0078, g2 = 0.0083 and N1 = N2 = 10 5 and with η = ξ = 1 (set ES1) and Γ12 = -9 and Ω = 0.5. At these parameters the components are in region 4 of the phase diagram of Fig. 6. Distance is measured in units of r0 and density in units of r -2 0 .

  's 3(b) and 4(b)). There is a boundary layer evidenced where the outer edge of the disk overlaps the

FIG. 6 :

 6 FIG.6: (Color online) Ω -Γ12 phase diagrams for parameters g1 = 0.0078, g2 = 0.0083, N1 = N2 = 10 5 with η = ξ = 1 (set ES1) with normalisation taken over the individual components (Eq. (4)). (a) Numerical simulations where triangles (squares) indicate that the vortex lattice in both components is triangular (square) and diamonds that no vortices have been nucleated. Filled triangles, squares and diamonds are where the two components are disk-shaped and coexist. The empty triangles and empty diamonds are where the two components have spatially separated: one component is a disk and the other an annulus with a giant skyrmion at the boundary layer; those triangles with a dot in the centre represent the appearance of coreless vortices inside the annulus. The crosses 'x' are where the two components have broken symmetry and are vortex sheets and the crosses '+' are rotating droplets. (b) A schematic representation of the numerical simulations. The solid lines indicate the boundary between different identified regions (determined by the geometry of the ground state) and the dashed lines the boundary between triangular and square lattices in region 1 and the boundary between peaks and no peaks in region 3. The boundary between region 1 and the others can be calculated analytically by Eq. (39) to give Γ12 = 0.008. The unit of rotation is ω.

  FIG. 7: (Color online) Ω -Γ12 phase diagrams for parameters g1 = 0.0067, g2 = 0.0063, N1 = N2 = 10 4 with η = 0.48 and ηξ 2 = 1 (set ES2) with normalisation taken over the individual components (Eq. (4)). (a) Numerical simulations where triangles indicate that the vortex lattice in component-2 is triangular, circles that component-2 contains rings of vortices and diamonds where no vortices have been nucleated. Filled triangles, circles and diamonds are where the two components are disk-shaped and coexist. The empty circles and empty diamonds are where the two components have spatially separated: one component is a disk and the other an annulus with a giant skyrmion at the boundary layer; those circles with a dot in the centre represent the appearance of coreless vortices inside the annulus. The crosses '+' are where the two components have broken symmetry and are rotating droplets. (b) A schematic representation of the numerical simulations. The solid lines indicate the boundary between different identified regions (determined by the geometry of the ground state).The boundary between region 1 and region 3 is also calculated analytically by Eq. (39) (dotted line). For these parameters there is no region 2. The unit of rotation is ω.

FIG. 8 :

 8 FIG.8:A schematic Ω-Γ12 phase diagram for g1 = g2 ≡ g = 0.0078, N1 = N2 = 10 5 where η = ξ = 1 (set ES3) and with the normalisation taken over the individual components (Eq. (4)). The solid lines indicate the boundary between different identified regions (determined by the geometry of the ground state) and the dashed lines the boundary between triangular and square lattices. For these parameters there is no region 3. The unit of rotation is ω.

2 FIG. 9 :

 29 FIG.9: Plots of maximum peak density in component-2 against Γ12 for Ω = 0.5 (dashed line) and Ω = 0.75 (solid line) when g1 = 0.0078, g2 = 0.0083, N1 = N2 = 10 5 and with η = ξ = 1 (set ES1) . The disjoint region when Ω = 0.75 for -0.4 < Γ12 < 0.1 is due to the absence of isolated density peaks in component-2 as a result of the appearance of vortex sheets. Density is measured in units of r -2 0 .

  Fig.'s 3(b) and 4(b).

FIG. 10

 10 FIG. 10: (Color online) Plots of (a) Sx and (b) a vectorial plot of (Sx, Sy) for Γ12 = -0.3 and Ω = 0.5 when g1 = 0.0078, g2 = 0.0083, N1 = N2 = 10 5 and with η = ξ = 1 (set ES1). The ring of skyrmions, at r = rs, traces the boundary layer and the coreless vortices exist for r < rs. Blow-ups (c) and (d) of (b) highlight the nature of the two types of topological defect; (c) a coreless vortex at (1, -1) and (d) a section of the multiply quantised skyrmion. The respective density plots for this case are shown in Fig. 3(b). Distance is measured in units of r0.

FIG. 11 :

 11 FIG.11: Total density profiles (divided by 10 4 ) obtained numerically (solid lines) and analytically (dashed lines) for two spatial separation cases with g1 = 0.0078, g2 = 0.0083, N1 = N2 = 10 5 , η = ξ = 1 (set ES1) and Γ12 = -10 and Ω = 0: (a) annulus plus disk, analytical estimate coming from (43) and (b) droplet, analytical estimate coming from(49). The inset in (a) shows the discontinuity of density at rs = 3.63. Distance is measured in units of r0 and density in units of r -2 0 .

  lattices in component-1 and density peaks in component-2.

FIG. 13

 13 FIG.13:(Color online) A series of density plots for component-1 (left column) and component-2 (right column). The parameters are g1 = 0.003 and g2 = 0.006 with η = ξ = 1 and Γ12 = 0.7 (which gives g12 = 0.0023) and normalisation taken over the total density (Eq. (5)). The angular velocity of rotation is Ω and it takes the values (a) 0.65 and (b) 0.9. There is a triangular lattice in (a) and a square lattice in (b). At these parameters the components are in region 1 of the phase diagram Fig.12. Distance is measured in units of r0 and density in units of r -2 0 .

FIG. 14

 14 FIG.14:(Color online) A series of density plots for component-1 (left column) and component-2 (right column). The parameters are g1 = 0.003 and g2 = 0.006 with η = ξ = 1 and Γ12 = 0.3 (which gives g12 = 0.0035) and normalisation taken over the total density (Eq. (5)). The angular velocity of rotation is Ω and it takes the values (a) 0.5 and (b) 0.9. At these parameters the components are in region 2 of the phase diagram Fig.12. Distance is measured in units of r0 and density in units of r -2 0 .

Acknowledgments

The authors wish to thank Thierry Jolicoeur for useful discussions that took place for the duration of this work. They are very grateful to the referee for his careful reading of the manuscript and his appropriate comments. We acknowledge support from the French ministry Grant ANR-BLAN-0238, VoLQuan.