J Syväranta
email: jari.syvaranta@jyu.fi

J Cucherousset

D Kopp

A Crivelli

R Céréghino

F Santoul

Dietary breadth and trophic position of introduced European catfish Silurus glanis in the River Tarn (Garonne River basin), southwest France

INTRODUCTION

Fishes are among the most frequently introduced freshwater organisms (Copp et al. 2005), and fish invasions are recognised as a serious threat to freshwater ecosystems (Dudgeon et al. 2006). Introduced predatory fish species can have important implications for native species and food webs (Vander Zanden et al. 1999, Eby et al. 2006). Large-bodied predators (i.e. larger than typical native predators) can destabilise natural communities by occupying the highest trophic position and introducing a new predatory pressure on lower trophic levels. Indeed, body size typically plays an important role in determining trophic positions (Cohen et al. 1993, Roy et al. 2002, Woodward et al. 2005, but see also Layman et al. 2005) due to morphological constraints in feeding related to limited gape size of predators (Cohen et al. 1993, Forsman 1996).

ABSTRACT: Although being a widely introduced and successfully established species, the European catfish Silurus glanis L. (the world's third largest freshwater fish) remains poorly studied in its introduced areas. Here we studied the trophic ecology of non-native European catfish in a large river system in south-western France using stomach content and stable isotope analyses (SIA). We used fin samples for SIA of catfish and hence tested the validity of using fin tissue as a proxy for muscle in SIA. The mean δ 15 N and δ 13 C values analysed from fin tissues did not differ from those analysed from muscle tissue and reflected strong and consistent relationships (r 2 = 0.95 for carbon and r 2 = 0.98 for nitrogen). The δ 15 N values varied almost 5 ‰ among the analysed catfish individuals, while δ 13 C values varied > 5 ‰. Total length of these catfish ranged from 200 to 2240 mm and was correlated with δ 15 N and especially with δ 13 C values. Although catfish length and δ 15 N values were positively correlated, the mean trophic positions of catfish increased only slightly from smaller individuals to larger ones (4.3 to 4.7). However, larger catfish were considerably 13 C-enriched in their δ 13 C values compared to smaller individuals and had up to 4 ‰ higher δ 13 C values than their expected aquatic prey. This might indicate frequent consumption of mammals and/or non-aquatic birds by the larger sized individuals, which were found in the catfish stomachs.

KEY WORDS: Stable isotopes • Introduced species • Top predator • River • Food web • Fin clips

European catfish Silurus glanis L (hereafter catfish) is the world's third largest (and the largest European) species (Stone 2007). Native to eastern Europe and western Asia, catfish has been widely introduced outside its native range due to its popularity as a sport fish and in aquaculture (Copp et al. 2009), making catfish an intriguing model to study the trophic ecology of large-bodied, non-native top predators. The species has now established self-sustained populations in many countries of western Europe. A recent review of the environmental biology of catfish by Copp et al. (2009) identified a serious lack of knowledge concerning this species, particularly in its introduced range and in riverine ecosystems. Such information is often a prerequisite for assessing potential risks to native species and ecosystems in areas where non-native species are established. The ecology of catfish outside its native range and interactions with native biota are poorly known, and this is especially true for the larger individuals (i.e. >1 m;Wysujack & Mehner 2005, Carol 2007). In its introduced areas, adult catfish can be at least twice as large as the native piscivorous fishes, but its precise position in the food web and potential impacts on other species are not thoroughly understood. In the few published studies, catfish have been reported to be a poor biomanipulator species (Wysujack & Mehner 2005), having only vague top-down impacts on the underlying fish community. But in areas with only a few or no native predatory fish species (such as in Spain), introduced catfish have had strong adverse impacts on the native fish community (Copp et al. 2009). Considering the low number of published studies on catfish trophic ecology, the practical and technical difficulties in studying these large-bodied species are no wonder, especially in large riverine ecosystems.

Stable isotope analysis (SIA) has been identified as a potentially valuable tool in studies on introductions of non-native species (e.g. Vander Zanden et al. 1999, 2004, Gorokhova et al. 2005, Cucherousset et al. 2007, Olsson et al. 2009), potentially including those examining the trophic position of catfish (Copp et al. 2009). The aim of the present study was to assess the validity of SIA to study the trophic ecology of catfish and then to assess the trophic position of introduced catfish within the food web of a large river in southwestern France. Our specific objectives were to: (1) compare stable isotope values from muscle and fin tissues to assess the suitability of fin material as a non-lethal proxy for muscle, (2) assess the trophic position and potential ontogenetic changes in trophic position of the catfish using both SIA and stomach contents, and (3) evaluate the potential impact of European catfish on the recipient food web. This should help in evaluating their potential ecological impacts on the ecosystem and on other species.

MATERIALS AND METHODS

Study area.

The study was conducted in 2007 and 2008 in the Tarn River (15 400 km 2 drainage area), a tributary of the Garonne River, in the locality of Villemur-sur-Tarn (43°52' N, 1°30' E, southwestern France). The studied stretch, located between 2 impassable weirs, was ~8 km long. The mean discharge was ~145 m 3 s -1 , ranging annually from 33 to 267 m 3 s -1 . The mean width of the studied stretch was ~81 m, with a mean depth of 3.6 m (maximum depth of 9 m). The catfish was introduced to the Tarn River in 1983. The species has been present in the study stretch since 1987, and a self-sustaining population has now been established. The Tarn River fish assemblage was composed of 24 other species, including 9 non-native species. Three cyprinid fishes dominated the fish community: bleak Alburnus alburnus L., roach Rutilus rutilus (L.) and bream Abramis brama (L.). Other predatory fish included 3 native species: pike Esox lucius L., perch Perca fluviatilis L. and European eel Anguilla anguilla L.; and 2 non-native species: largemouth bass Micropterus salmoides Lacepède and pikeperch Sander lucioperca (L.).

Sample collection. Stomach contents and SIA samples of catfish were collected from local recreational anglers. Fishing was done using a reel and rod, and most fish were caught at dusk. Samples for SIA were collected from mid-April to the end of June both years, which corresponds to the most intensive feeding and growth periods of catfish (Orlova & Popova 1976, Copp et al. 2009). Stomach contents were collected over a similar period, with additional samples being collected in August and September 2008 to increase the sample size. SIA samples of other predatory fish species and of the most abundant potential prey fish for catfish were also collected. In addition, the mussels Dreissena sp. and Corbicula sp., spinycheek crayfish Orconectes limosus Rafinesque, red swamp crawfish Procambarus clarkii Girard, edible frog Rana esculenta L., pigeon Columbia livia Gmelin and coypu Myocastor coypus Molina were collected from local angling and hunting societies to represent other potential prey which have been previously reported to be consumed by catfish (Orlova & Popova 1976, Carol et al. 2009, Copp et al. 2009). Mussels were also used as δ 15 N baseline indicators for calculating trophic positions of species further up the food chain following Vander Zanden et al. (1997). Total length (TL, mm) was recorded from all fish, and a sample of white muscle tissue was dissected for SIA from prey organisms (posterior to the dorsal fin in fish). Small samples of fin tissue were collected from larger predatory fish as recreational anglers typically catch and release these predatory species. Fin clipping was selected because stable isotope values of fins are known to correlate closely with those of muscle tissue in small salmonids and other small fish, allowing a non-lethal sampling for SIA (Jardine et al. 2005, Kelly et al. 2006, Sanderson et al. 2009). However, the validity of fin tissue as a proxy for muscle tissue for much larger fish, such as catfish, has not been previously tested. Larger-sized fins may contain greater proportions of bone material, leading to increases in inorganic carbon, thereby potentially affecting the δ 13 C values of fin samples and fin-muscle comparisons. Therefore, we collected additional fin and muscle samples of catfish from 3 other locations (Loire River in western France, Camargue [Rhône delta area] in southern France and a local fish hatchery) to compare δ 13 C and δ 15 N values between muscle and fin tissue in 28 individuals (TL = 60 to 1340 mm).

Stable isotope and data analyses. Stable isotopes, especially those of carbon and nitrogen, are used to examine consumers' trophic ecology from small tissue samples, providing a time-integrated measure of trophic position and energy sources. Stable carbon and nitrogen isotope ratios in consumer tissues reflect those in their diet (Peterson & Fry 1987), depending on the relative proportions of prey species assimilated over a long time period. Nitrogen stable isotope ratios ( 15 N/ 14 N, expressed relative to a standard as δ 15 N, in parts per thousand) in consumers are typically enriched in the heavier ( 15 N) isotope by from 2 to 4 ‰ per trophic level (Minagawa & Wada 1984, Peterson & Fry 1987), making δ 15 N values useful in defining trophic positions of consumers. The carbon isotope ratios ( 13 C/ 12 C, δ 13 C, in parts per thousand) fractionate to a lesser extent (0 to 1 ‰) and are typically used to define diet compositions or sources of energy (primary production). All samples for SIA were oven dried (60°C for 48 h) and ground into a fine homogeneous powder using a mixer mill (Retsch MM 200). Approximately 0.2 mg of sample was accurately weighed into small tin cups, and stable isotope ratios of carbon and nitrogen were analysed in a Carlo Erba NC2500 elemental analyser coupled to a Thermo Finnigan Mat Delta XP isotope ratio mass spectrometer. Isotope ratios are reported in delta notation as per international standards: Peedee belemnite carbonate for δ 13 C and atmospheric nitrogen for δ 15 N. Data were corrected using working standards (bass muscle, bovine liver, nicotinamide; SD < 0.2 ‰ for both δ 13 C and δ 15 N) that were previously calibrated against International Atomic Energy Agency (IAEA) standards. All stable isotope analyses were performed at the Stable Isotopes in Nature Laboratory, University of New Brunswick, Canada.

Linear regression and paired t-tests were used to examine differences between muscle and fin δ 13 C and δ 15 N values. In addition, linear regression analysis was also applied to examine potential ontogenetic shifts in catfish by comparing δ 13 C and δ 15 N values and TL of individuals. Trophic positions (TP) for all prey and predatory fish were calculated following Vander Zanden et al. (1997): TP = (δ 15 N fish -δ 15 N mussel / ∆) + 2, where ∆ is the trophic fractionation factor for δ 15 N values and 2 is the trophic position for mussels. Considerable variation exists in the literature around the mean trophic fractionation for δ 15 N (e.g. Martinez del Rio et al. 2009), the most often used values being between 2.8 and 3.4 ‰. We used a ∆ value of + 3 ‰, and using any value within the proposed range would not dramatically affect the conclusions drawn, when comparing organisms within the same food web. Catfish were also grouped into 3 main categories according to their TL; small-bodied (< 600 mm) and medium-bodied (600 to 1300 mm) individuals that correspond to typical body sizes of the other piscivorous fish species, whereas large-bodied catfish (>1300 mm) correspond to sizes not reached by any other fish species in the study area. Differences in mean isotope values and trophic positions between catfish size categories and pike and pikeperch (principal top predators within the system) were examined using analysis of variance (ANOVA) followed by Tukey's multiple comparisons. All statistical analyses were performed using SPSS for Windows 13.0 software package (SPSS, Inc.).

RESULTS

Altogether, 49 individuals of catfish Silurus glanis were caught and sampled for SIA from the Tarn River. Eighteen of these were categorised as small bodied ranging in TL from 200 to 580 mm, whereas 15 individuals were categorised as medium bodied (600 to 1300 mm) and 16 individuals as large bodied (1400 to 2240 mm; Table 1 also caught and analysed for stable isotopes, but was not categorised, or used in any of the statistical analyses. Although nearly 100 individual catfish were sampled for stomach contents by the local anglers, only 31 had prey in their stomach. Among these individuals, cyprinid fish species were the most frequently observed prey in all catfish size classes, whereas molluscs and crayfish became less frequent, and terrestrial prey became more frequent in large-bodied individuals (Table 2). Both δ 13 C and δ 15 N values between catfish muscle and fin tissues were highly correlated and consistent along the range of analysed isotope values, and the values did not differ between the 2 tissue types (Fig. 1). The mean ± SE value of δ 13 C fin tissues was only slightly elevated compared to muscle tissue (+ 0.27 ± 0.2 ‰, paired t-test: t 27 = -1.45, p = 0.158), whereas the mean δ 15 N value was 0.09 ± 0.1 ‰ lower in fin tissue compared to muscle (paired t-test: t 27 = -0.76, p = 0.457). Since no difference was found between catfish fin and muscle tissue δ 13 C and δ 15 N values and the values were highly correlated, we did not make any correction on the catfish fin tissue δ-values for comparison with muscle tissue of other food web components.

A clear relationship was found between catfish TL and both δ 13 C and δ 15 N values (Fig. 2), with larger individuals being both significantly 13 C and 15 N enriched (F 1,48 = 55.8, p < 0.001, r 2 = 0.54 for δ 13 C and F 1,48 = 19.5, p < 0.001, r 2 = 0.29 for δ 15 N). The δ 13 C values increased rather consistently along the range of TL, whereas the increase in δ 15 N values with catfish TL was more modest and accompanied with considerable variation. The same was evident in the trophic positions of catfish, which changed only slightly, albeit significantly, from small-bodied individuals to medium-and large-bodied individuals (Table 1). Medium-and large-bodied individuals had significantly elevated trophic positions compared to smallbodied individuals (F 2,48 = 10.2, p < 0.001), but did not differ between each other. The trophic position of pike did not differ from those of any catfish size category, whereas pikeperch had a significantly (p < 0.001) elevated trophic position compared to pike and all catfish sizes.

The individuals (-27.2 to -23.2 ‰) and by 2.5 ‰ for largebodied individuals (-24.6 to -22.1 ‰; Fig. 3). Similarly, the δ 15 N values for all individuals ranged by 4.7 ‰ (14.0 to 18.7 ‰); those of small-bodied individuals by 3.2 ‰ (14.0 to 17.2 ‰), by 1.4 ‰ for medium-bodied (16.2 to 17.6 ‰) and by 3.3 ‰ (15.4 to 18.7 ‰) for large-bodied individuals. The 3 TL groups of catfish differed significantly in their mean δ 13 C and δ 15 N values (F 2,48 = 15.46, p < 0.001 for δ 13 C and F 2,48 = 10.17, p < 0.001 for δ 15 N). Large-bodied catfish had a mean ± SE δ 13 C value of -23.3 ± 0.2 ‰, which was significantly enriched relative to the mean values of both small-bodied (-25.1 ± 0.2 ‰) and medium-bodied individuals (-24.5 ± 0.3 ‰).

In contrast, the mean ± SE δ 15 N value of small-bodied catfish (16.1 ± 0.2 ‰) was significantly lower than those of medium-bodied (16.8 ± 0.1 ‰) and large-bodied (17.3 ± 0.2 ‰) individuals. Most catfish individuals, except the juvenile one, had δ 15 N values consistently 2 to 4 ‰ higher than their putative prey sources. In contrast, the juvenile catfish individual had a δ 15 N value ~6 ‰ lower than most larger catfish (Fig. 3). The δ 13 C values of most smalland medium-bodied catfish were in accordance with their principal aquatic prey, as indicated by the stomach contents. However, larger catfish individuals were considerably 13 C-elevated (up to 4 ‰ higher δ 13 C values) than their expected fish prey.

DISCUSSION

Stomach contents of catfish Silurus glanis from the Tarn River indicated a dominance of cyprinid fishes and benthic prey (molluscs and crayfish) in their diet, particularly in the smaller size groups. These are also the commonly reported prey items for catfish in the literature (e.g. Carol et al. 2009, Copp et al. 2009). However, sampling for catfish stomach contents revealed the often encountered difficulty with primarily piscivorous fish. High frequency of empty stomachs is particularly troublesome when studying fish species that are difficult to capture (e.g. low density, inaccessible habitats) and handle (e.g. endangered, large body size), and are also preferably released alive after capture. This fact highlights the need for alternative methods to study the trophic ecology of large-bodied predators, especially methods that allow nonlethal sampling. Samples from the Tarn River were caught by local anglers, who, after measuring TL and collecting a small fin clip, released all individuals into the river. Fin and muscle tissue δ 13 C and δ 15 N comparisons in salmonids and sculpins have produced consistent relationships (Jardine et al. 2005, Kelly et al. 2006, Sanderson et al. 2009), thereby enabling the use of fin tissue as a proxy for muscle tissue. However, not much has been published on larger fish species or on anything even remotely similar to catfish. Therefore, it was necessary to test if the fin samples from catfish could be used in SIA with confidence. Our results clearly indicate that the values analysed from fin tissues provide a suitable alternative for muscle tissue from catfish and can be used to study the trophic ecology of catfish. The size of catfish did not affect the fin-muscle comparisons.

While the mean δ 13 C and δ 15 N values between fin and muscle tissue were not significantly different, muscle had a tendency for lower δ 13 C and higher δ 15 N values, as also reported by Sanderson et al. (2009).

Comparisons of catfish TL and stable isotope values revealed significant size-related relationships (Fig. 2). δ 13 C values in particular increased consistently with catfish TL, indicating a gradual shift in diet composition with increasing size (and age). Although the δ 15 N values and TL of catfish were positively correlated, the δ 15 N values and, hence, TP only slightly increased during ontogeny. The mean difference in δ 15 N values between the smallest and largest size groups was only 1.2 ‰. Although statistically significant, this translates only to < 0.5 trophic levels (Table 1). Such subtle changes in δ 15 N values can be further masked behind the potential inherent variability in stable isotope values (Barnes et al. 2008). It is also possible that δ 15 N values increase along increasing animal age irrespective of diet (Minagawa & Wada 1984, Overman & Parrish 2001), and the largest catfish in our samples are likely to be very old (> 20 yr). However, the δ 13 C values also indicated a change in diet composition that could partly explain the slight change in δ 15 N values, albeit no meaningful ontogenetic changes in trophic position were evident. The 1 juvenile catfish (60 mm) in our samples differed substantially from all other catfish samples in both δ 15 N and δ 13 C, suggesting that any notable ontogenetic shift occurs within the gap in our size range (60 to 200 mm). The δ 13 C values measured from Tarn River catfish spanned over 5 ‰, indicating both considerable individual-level differences in diet composition (Syväranta & Jones 2008) and a wide use of available resources (Fig. 3). This is consistent with the conclusions of Copp et al. (2009) who stated that the dietary spectrum of catfish is considerably greater than that for many other piscivorous fish, such as pike or pikeperch.

Catfish did occupy a trophic position typical of piscivorous fish and equal to that of pike (Table 1), al-though significantly lower than the trophic position of the non-native pikeperch as also reported by Kopp et al. (2009). However, the variation observed in catfish δ 13 C values was considerably higher than that for pike or pikeperch, suggesting that catfish is a much more flexible predator in terms of diet than the 2 other top predators in this ecosystem. This would also indicate that catfish does not have a fixed role as an efficient top predator in the food web and is perhaps unlikely to have any serious adverse impacts on the underlying fish community.

Trophic position is often expected to increase from smaller-bodied species to larger-bodied ones (e.g. Cohen et al. 1993, Roy et al. 2002, Woodward et al. 2005). Similarly, one might expect the trophic position to increase within a single predatory species during ontogeny, especially if body size increases as considerably as in catfish. Such increases in trophic position could result from changes in (feeding) behaviour and increased gape size, allowing the consumption of larger prey from potentially higher trophic levels. Therefore, it is perhaps surprising that catfish in the Tarn River reflected only negligible changes in trophic position with increasing body size. However, catfish are known to prey upon fish considerably smaller than what could be expected based on their gape size (Wysujack & Mehner 2005), and they are also considered an opportunistic species (Copp et al. 2009). This can actually mean that their diet composition is not expected to go through very drastic changes according to increasing body size for the larger individuals, at least in terms of feeding from higher trophic levels, which is consistent with the present results.

The δ 13 C values of large-bodied catfish were heavily enriched compared to those of most of the smaller individuals and to those of any analysed aquatic prey item, reflecting a significant dietary contribution from a relatively 13 C-enriched source. Such elevated δ 13 C values were only analysed from pigeons and coypu, but the δ 15 N values of this prey source were considerably lower than those in any other prey source and ~10 ‰ lower than the δ 15 N values of catfish. Therefore, any catfish heavily relying on this prey source ought to show corresponding decrease in δ 15 N values, which was not evident. Catfish are known to feed on many anadromous species, which would reflect marine δ 13 C values when migrating from marine to freshwater habitats, and feeding on prey with such distinct δ 13 C values can quickly elevate the δ 13 C values of catfish (MacAvoy et al. 2000, Syväranta et al. 2009). But migrations of anadromous fish are effectively prevented by damming in the Tarn River, so any dietary contribution from anadromous fish with elevated δ 13 C values is unlikely. Some terrestrial vertebrates, mainly birds and small mammals, are known to contribute to the diets of catfish (Copp et al. 2009). These items were also observed in the stomach contents of large-bodied catfish in our study (Table 2), and catfish in the Tarn River are frequently seen catching pigeons from rocks in the river. However, comprehensive sampling of this prey source is not always straightforward, and, hence, our stable isotope samples of potential terrestrial prey for catfish consisted only of 3 pigeons and 3 coypu. Therefore, it is possible that the inconsistency between the δ 13 C and especially δ 15 N values between the larger catfish and their prey may result from insufficient sampling of terrestrial prey sources and potentially biased estimates for the δ 15 N values.

In conclusion, the combined use of stomach contents and SIA provided new insights into the trophic ecology of catfish outside of its native range. Our first original objective was to assess the suitability of fin clips as a proxy for muscle tissue, and our results clearly indicate that fin clips provide a suitable and non-lethal proxy for muscle tissue in SIA. Also, we aimed to assess the trophic positions and ontogenetic changes in trophic ecology of catfish in the Tarn River using both SIA from fin clips and gut contents. While the changes in trophic position with increasing size were rather small in the size range we studied, we demonstrated that the ontogenic changes in body size may still be associated with a shift towards a prey less exploited by other predators. This mechanism might facilitate the establishment of this non-native species outside its native range, as it might reduce both intraspecific and interspecific competition by providing a partly vacant niche for larger catfish. In conclusion to our third objective, our results also indicate that catfish in the Tarn River ecosystem is unlikely a voracious top predator, having strong adverse impact on the native fish community, but rather a scavenging predator able to exploit a wide variety of prey sources. However, future studies should determine in greater detail whether this might alter the structure and function of the food web by comparing similar ecosystems were the species is absent and present.

Acknowledgements. Local fishing and hunting societies are thanked for providing us with many samples. Three anonymous reviewers are thanked for their valuable comments on an earlier version of this manuscript. J.S. and D.K. were supported by the CNRS and Region Midi-Pyrénées post-doctoral fellowships.

Fig. 1 .

 1 Fig. 1. Silurus glanis. Comparisons of δ 13 C (upper panel) and δ 15 N (lower panel) values between European catfish fin and muscle tissues. Solid thin lines indicate the fitted regressions (δ 13 C: y = 0.89x -3.24, r 2 = 0.95, n = 28; δ 15 N: y = 0.84x + 2.19, r 2 = 0.98, n = 28) and dashed lines indicate 95% confidence limits. The thick solid lines indicate a 1:1 fit

Fig. 2 .Fig. 3 .

 23 Fig. 2. Silurus glanis. Relationships between European catfish total length (TL) and δ 13 C (upper panel) and δ 15 N (lower panel) values. Solid lines indicate fitted regressions (δ 13 C: y = 0.001x -25.8, r 2 = 0.54, n = 49; δ 15 N: y = 0.001x + 15.9, r 2 = 0.29, n = 49)

Table 1 .

 1 ). One juvenile (60 mm) catfish was Mean ± SE total length (TL), minimum and maximum length, and mean ± SE trophic position (TP) for catfish Silurus glanis and each analysed fish species from the Tarn River

	Species	n	TL (mm) TL (min./	TP
				max., mm)
	European catfish 49 1000 ± 85	200/2240 4.5 ± 0.0
	Small bodied	18	401 ± 27	200/580 4.3 ± 0.1
	Medium bodied 15	946 ± 68	600/1300 4.5 ± 0.0
	Large bodied	16 1724 ± 65 1400/2240 4.7 ± 0.1
	Crayfish	7	-		3.0 ± 0.1
	Bleak	9	101 ± 3.4	85/108 3.4 ± 0.1
	Common bream	15	159 ± 9.6 110/220 3.7 ± 0.0
	Roach	9	143 ± 14	110/210 3.6 ± 0.0
	European eel	5	418 ± 68	230/640 4.2 ± 0.1
	Largemouth bass	3	213 ± 8.8 200/230 4.2 ± 0.1
	Pike	8	654 ± 42	500/770 4.5 ± 0.1
	Pikeperch	11	643 ± 17	550/750 5.1 ± 0.0

  Macko SA, McIninch SP, Garman GC (2000) Marine nutrient contributions to freshwater apex predators. Oecologia 122:568-573 Martinez del Rio C, Wolf N, Carleton SA, GannesLZ (2009) Isotopic ecology ten years after a call for more laboratory experiments. Biol Rev Camb Philos Soc 84:91-111 Minagawa M, Wada E (1984) Stepwise enrichment of 15 N along food chains: further evidence and the relation between δ 15 N and animal age. Geochim Cosmochim Acta 48:1135-1140 Olsson K, Stenroth P, Nyström P, Graneli W (2009) Invasions and niche width: Does niche width of an introduced crayfish differ from a native crayfish? Freshw Biol 54:1731-1740 Orlova EL, Popova OA (1976) The feeding of predatory fish, the sheatfish, Silurus glanis, and the pike, Esox lucius, in the Volga Delta following regulation of the discharge of the river. J Ichthyol 16:75-87 Overman NC, Parrish DL (2001) Stable isotope composition of walleye: 15 N accumulation with age and area-specific differences in δ 13 C. Can J Fish Aquat Sci 58:1253-1260 Peterson BJ, Fry B (1987) Stable isotopes in ecosystem studies. Annu Rev Ecol Syst 18:293-320 Roy K, Jablonski D, Valentine JW (2002) Body size and invasion success in marine bivalves. Ecol Lett 5:163-167 Sanderson BL, Coe HJ, Pelekis V, Steel EA, Reichert WL (2009) Nonlethal sampling of fish caudal fins yields valuable stable isotope data for threatened and endangered fishes. Trans Am Fish Soc 138:1166-1177 Stone R (2007) The last of the Leviathans. Science 316: 1684-1688 Syväranta J, Jones RI (2008) Changes in feeding niche widths of perch and roach following biomanipulation, revealed by stable isotope analysis. Freshw Biol 53:425-434 Syväranta J, Cucherousset J, Kopp D, Martino A, Céréghino R, Santoul F (2009) Contribution of anadromous fish to the diet of European catfish in a large river system. Naturwissenschaften 96:631-635 Vander Zanden MJ, Cabana G, Rasmussen JB (1997) Comparing trophic position of freshwater fish calculated using stable nitrogen isotope ratios (δ 15 N) and literature dietary data. Can J Fish Aquat Sci 54:1142-1158 Vander Zanden MJ, Casselman JM, Rasmussen JB (1999) Stable isotope evidence for the food web consequences of species invasions in lakes. Nature 401:464-467 Vander Zanden MJ, Olden JD, Thorne JH, Mandrak NE (2004) Predicting occurrences and impacts of smallmouth bass introductions in north temperate lakes. Ecol Appl 14: 132-148 Woodward G, Ebenman B, Emmerson M, Montoya JM, Olesen JM, Valido A, Warren PH (2005) Body-size in ecological networks. Trends Ecol Evol 20:402-409 Wysujack K, Mehner T (2005) Can feeding of European catfish prevent cyprinids from reaching a size refuge? Ecol Freshwat Fish 14:87-95

	➤	➤ ➤	Carol J, Benejam L, Benito J, García-Berthou E (2009) Growth and diet of European catfish (Silurus glanis) in early and
	➤ ➤	➤ ➤	late invasion stages. Fundam Appl Limnol 174:317-328 Cohen JE, Pimm SL, Yodzis P, Saldana J (1993) Body sizes of animal predators and animal prey in food webs. J Anim Ecol 62:67-78
	➤ ➤ ➤	➤ ➤ ➤ ➤	Copp GH, Bianco PG, Bogutskaya NG, Ero "s T and others (2005) To be, or not to be, a non-native freshwater fish? J Appl Ichthyol 21:242-262 Copp GH, Britton R, Cucherousset J, García-Berthou E, Kirk R, Peeler E, Stake • nas S (2009) Voracious invader or benign feline? A review of the environmental biology of European catfish Silurus glanis in its native and introduced ranges.
	➤ ➤	➤	Fish Fish 10:252-282 Cucherousset J, Aymes JC, Santoul F, Céréghino R (2007) Stable isotope evidence of trophic interactions between
			introduced brook trout (Salvelinus fontinalis) and native
			brown trout (Salmo trutta) in a mountain stream of south-
			west France. J Fish Biol 71(Suppl D):210-223
		➤ ➤	Dudgeon D, Arthington AH, Gessner MO, Kawabata ZI and others (2006) Freshwater biodiversity: importance, threats,
			status and conservation challenges. Biol Rev Camb Philos
			Soc 81:163-182
		➤ ➤	Eby LA, Roach WJ, Crowder LB, Stanford JA (2006) Effects of stocking-up freshwater food webs. Trends Ecol Evol 21:
			576-584
		➤	Forsman A (1996) Body size and net energy gain in gape-limited predators: a model. J Herpetol 30:307-319
		➤ ➤	Gorokhova E, Hansson S, Höglander H, Andersen CM (2005) Stable isotopes show food web changes after invasion by
			the predatory cladoceran Cercopagis pengoi in a Baltic
			Sea bay. Oecologia 143:251-259
		➤	Jardine TD, Gray MA, McWilliam SM, Cunjak RA (2005) Stable isotope variability in tissues of stream fishes. Trans
			Am Fish Soc 134:1103-1110
		➤	Kelly MH, Hagar WG, Jardine TD, Cunjak RA (2006) Non-lethal sampling of sunfish and slimy sculpin for stable
			isotope analysis: how scale and fin tissue compare with
			muscle tissue. N Am J Fish Manage 26:921-925
		➤ ➤	Kopp D, Cucherousset J, Syväranta J, Martino A, Céréghino R, Santoul F (2009) Trophic ecology of the pikeperch
			(Sander lucioperca) in its introduced areas: a stable
			isotope approach in southwestern France. C R Biol 332:
			741-746
		➤	Layman CA, Winemiller KO, Arrington DA, Jepsen DB (2005) Body size and trophic position in a diverse tropical food
			web. Ecology 86:2530-2535
	LITERATURE CITED Barnes C, Jennings S, Polunin NVC, Lancaster JE (2008) The importance of quantifying inherent variability when inter-MacAvoy SE, ➤ ➤ ➤ ➤ ➤ ➤ ➤
	preting stable isotope field data. Oecologia 155:227-235
	Carol J (2007) Ecology of an invasive fish (Silurus glanis) in Catalan reservoirs. PhD dissertation, University of Girona,	➤
	Girona


A V Kavokin

M Vladimirova

B Jouault

T C H Liew

J R Leonard

L V Butov

Ballistic spin transport in exciton gases

Keywords: numbers: 71.35.-y, 03.75.Kk, 03.75.Mn, 73.63.Hs, 78.55.Cr

Traditional spintronics relies on spin transport by charge carriers, such as electrons in semiconductor crystals. The challenges for the realization of long-range electron spin transport include rapid spin relaxation due to electron scattering. Scattering and, in turn, spin relaxation can be effectively suppressed in excitonic devices where the spin currents are carried by electrically neutral bosonic quasi-particles: excitons or exciton-polaritons. They can form coherent quantum liquids that carry spins over macroscopic distances. The price to pay is a finite life-time of the bosonic spin carriers. We present the theory of exciton ballistic spin transport which may be applied to a range of systems supporting bosonic spin transport, in particular, to indirect excitons in coupled quantum wells. We describe the effect of spin-orbit interaction for the electron and the hole on the exciton spin, account for the Zeeman effect induced by external magnetic fields and long range and short range exchange splittings of the exciton resonances. We also consider exciton transport in the non-linear regime and discuss the definitions of the exciton spin current, polarization current and spin conductivity.

I. INTRODUCTION

Excitons are electrically neutral and have finite lifetimes. These are two obstacles which make the development of excitonic spintronics, or spin-optronics challenging. How possibly one can explore the current, which is carried by neutral particles, and whose amplitude changes with distance and time? -is a fair question to ask. While electrons and holes have been considered as perfectly valid spin carriers, and exotic effects like the spin Hall effect 1 have been intensively studied for them [2][3][4][5][6] , the spin currents carried by excitons [7][8][9] and exciton-polaritons 10,11 over tens or even hundreds of micrometers remained relatively less explored. There existed a huge imbalance of theoretical works on fermionic and bosonic spin transport. This is changing now. A number of phenomena have been observed and studied in the field of bosonic spin currents recently [7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26] . To summarize tens of publications in one sentence: bosonic systems bring new quantum coherent effects to the physics of spin transport. For instance, stimulation [23][24][25] and amplification 26 of spin currents are possible in exciton and exciton-polariton Bose gases. Bosonic spintronics or spin-optronics operates with electrically neutral spin carriers which makes control of spin currents carried by excitons a non-trivial task. Fortunately, the exciton density replaces charge in many aspects: the density currents may be efficiently controlled by stationary or dynamic potential gradients as demonstrated in recent works 27,28 . Combined with evident advantages of bosonic amplification and low dephasing, this makes spin-optronics a valuable alternative to fermionic spintronics. Besides bosonic effects, exciton spin transport has another important specific feature: it is dissipative by its nature, as the spin carriers have a finite (and short for excitons in regular materials) life-time. In continuous wave optical experiments stationary spin textures can appear: excitons are injected in the structure, they propagate ballistically or diffusively, and eventually disappear by radiative recombination. Their polarization properties and spin are inherited by the emitted photons, that is why the polarization patterns observed in near field photoluminescence experiments directly characterize exciton spin currents in the plane of the structure.

The goal of this work is to define what the exciton spin, magnetization and polarization currents are, and to explain how they can be described within the most frequently used spin density matrix (DM) approach and mean-field 29,30 approximation. We consider a specific system, namely a planar zinc-blend semiconductor structure containing quantum wells, where excitons can be formed. This choice is motivated by recent experimental results in GaAs/AlGaAs based coupled quantum wells. We limit the scope of this paper to heavy-hole excitons, however, our approach can be easily extended to lighthole excitons or excitons in quantum wells of a different symmetry. We do not speak here about the large variety of recent experimental results and application of the formalism presented here to the description of one particular experiment, as this would make this paper too long and too specific. For a direct comparison of theoretical simulations with the experimental data we address the reader to Ref. 8. The approaches formulated here are suitable for the description of a variety of excitonic spin effects in quantum wells.

The paper is organized as follows. In Section II we introduce the spin DM formalism accounting for the different mechanisms of spin re-orientation and the relation to electron and hole spin currents. In Section III we present numerical results obtained within the spin DM formalism and analyze them. In Section IV we study the non-linear spin dynamics of propagating excitons using the Gross-Pitaevskii (GP) equations. The next three sections of the paper are devoted to exciton spin currents and polarization currents. Conclusions are perspectives are given in Section VIII.

II. THE SPIN MATRIX FORMALISM FOR PROPAGATING EXCITONS

In zinc-blend semiconductor quantum wells (e.g. in the most popular GaAs/AlGaAs system), the lowest energy exciton states are formed by electrons with spin projections on the structure axis of +1/2 and -1/2 and heavy holes whose quasi-spin (sum of spin and orbital momentum) projection to the structure axis is +3/2 or -3/2. Consequently, the exciton spin defined as the sum of the electron spin and heavy-hole quasi-spin may have one of four projections on the structure axis: +1,-1,+2,-2 [START_REF] See | Superlattices and other heterostructures[END_REF] . These states are usually nearly degenerate, while there may be some splitting between them due to the short and long-range exchange interactions. Only the states with quasi-spin projections ±1 are coupled to the light, these are so-called bright states. The states with quasispin projections ±2 are called dark states.

It is important to note that the present formalism addresses the spin part of the exciton wavefunction, which is a product of electron and hole spin functions. For example, the probability to find the exciton in the spin state +1 is given by a product of probabilities to find an electron in the spin state -1/2 and the heavy hole in the spin state +3/2. The four component exciton wave-function is:

Ψ = (Ψ +1 , Ψ -1 , Ψ +2 , Ψ -2 ) = ( Ψ e,-1 2 Ψ h,+ 3 2 , Ψ e,+ 1 2 Ψ h,-3 2 , Ψ e,+ 1 2 Ψ h,+ 3 2 , Ψ e,-1 2 Ψ h,- 3 2 ) (1)
where Ψ e,+ 1 2 and Ψ e,-1 2 are the components of the electron spinor wavefunction; Ψ h,+ 3 2 and Ψ h,-3 2 are the components of the heavy-hole spinor wavefunction.

To describe the dynamics of the system we will first define the Hamiltonian (Section II A) describing the different physical mechanisms of spin evolution. We then introduce the spin DM (Section II B) for the description of exciton spin states. We relate its components to the observable Stokes' vectors of light emitted by bright excitons and use the Liouville equation to describe its evolution in time. From the DM one can describe the polarization state of excitons and thus exciton spin currents. It is also instructive to consider the relationship with electron and hole spin currents (Section II C).

A. Hamiltonian

Here we derive the exciton Hamiltonian in the basis of +1,-1,+2,-2 states, accounting for the spin-orbit interaction (Dresselhaus and Rashba effects) [START_REF] Rashba | [END_REF]37 , long-and short-range exchange interactions [START_REF] Ivchenko | Optical spectroscopy of semiconductor nanostructures[END_REF] and Zeeman effect, but neglecting exciton-exciton interactions, which will be discussed in the Section IV and neglecting magnetic field effect on center-of-mass motion and internal structure of exciton [34][35][36] . We consider excitons propagating ballistically in the plane of a quantum well. We shall characterize them by a fixed wave-vector, k ex . We represent the full exciton Hamiltonian Ĥtot ex as a sum of three parts describing the spin-orbit and Zeeman effects on electrons, Ĥe ex , and holes, Ĥh ex , and the exchange induced splittings of exciton states, Ĥex ex :

Ĥtot ex = Ĥe ex + Ĥh ex + Ĥex ex . (2)

Dresselhaus and Zeeman terms

We recall that the Rashba-Dresselhaus effect is a momentum-dependent splitting of spin bands in twodimensional semiconductor systems. It originates from a combined effect of the atomic spin orbit coupling and asymmetry of the potential in the direction perpendicular to the two-dimensional plane. This asymmetry comes either from the applied bias (which is described by the Rashba term in the Hamiltonian) or from the intrinsic asymmetry of the crystal lattice (described by the Dresselhaus term in the Hamiltonian). We shall separately consider both the Dresselhaus term (in this sub-section) and the Rashba term (in the next sub-section).

In order to build the 4×4 matrix Hamiltonian for excitons, we start with simpler 2 × 2 Hamiltonians describing the spin-orbit and Zeeman effects for electrons and holes.

The electron Hamiltonian in the basis of (+1/2,-1/2) spin states is:

Ĥe = β e (k e,x σx -k e,y σy ) - 1 2 g e µ B Bσ. (3)
Here g e is the electron g-factor, µ B is the Bohr magneton, B is a magnetic field, σ is the Pauli matrix vector, and β e is the Dresselhaus constant describing spin-orbit interactions of electrons. The Pauli matrix operators are:

σz = [ 1 0 0 -1 ] , σy = [ 0 -i i 0 ] , σx = [ 0 1 1 0 ] . (4
where φ is the angle between k ex and the chosen x-axis. The exciton Hamiltonian needs to be written in the basis of (+1,-1,+2,-2) exciton spin states, which correspond to (-1/2, +1/2, +1/2, -1/2) electron spin states. The electron spin-flip couples +1 and +2 states as well as -1 and -2 states. For each of these two couples of states we apply the Hamiltonian [Eq. ( 5)], which results in the following electronic contribution to the 4 × 4 exciton Hamiltonian:

Ĥe ex =     1 2 g e µ B B 0 β e k e e -iφ 0 0 -1 2 g e µ B B 0 β e k e e iφ β e k e e iφ 0 -1 2 g e µ B B 0 0 β e k e e -iφ 0
1 2 g e µ B B     (6)

Rashba terms

We note that another possible spin-orbit contribution to the Hamiltonian may come from the Rashba effect, which takes place in biased quantum wells. The Rashba term to be added in Eq. ( 3) is α e (σ x k e,y -σy k e,x ), where α e is a constant proportional to the Rashba field. The contribution of the Rashba term to the electron Hamiltonian, in the basis of (+1/2,-1/2) electron spin states, can be re-written:

Ĥ′ e = [ 0 iα e k e e -iφ -iα e k e e iφ 0 ] . (7)
Using the same procedure as for Dresselhaus terms, this gives an additional contribution to the exciton Hamiltonian, in the basis of (-1/2, +1/2, +1/2, -1/2) electron spin states:

Ĥe′ ex =     0 0
-iα e k e e iφ 0 0 0 0 iα e k e e -iφ iα e k e e -iφ 0 0 0 0 -iα e k e e iφ 0 0

    (8)
Unless stated explicitly, we will for simplicity omit the Rashba terms in the rest of this paper and consider only the Dresselhaus terms. Note that there are no linear in wave-vector Rashba terms for heavy holes in zinc-blend quantum wells grown along the (001)-axis.

Heavy-hole contribution (Faraday geometry)

The heavy hole contribution to the Hamiltonian can be calculated with the reasoning similar to the electron case. The heavy-hole Hamiltonian written in the basis of (+3/2,-3/2) states is:

Ĥh = β h (k h,x σx + k h,y σy ) - 1 2 g h µ B B σz . ( 9
)
Here g h is the heavy-hole g-factor and β h is the Dresselhaus constant for heavy holes [START_REF] Rashba | [END_REF]37 . Note that the Dresselhaus Hamiltonian is different for heavy holes formed by p-orbital states and for conduction band electrons formed by s-orbital electronic states in a zinc-blend crystal lattice. The resulting from Dresselhaus coupling effective magnetic fields acting upon electron and heavy hole spins are oriented differently as well. Re-writing Eq. ( 9), we obtain:

Ĥh = [ -1 2 g h µ B B β h (k h,x -ik h,y ) β h (k h,x + ik h,y ) 1 2 g h µ B B ] = [ -1 2 g h µ B B β h k h e -iφ β h k h e iφ 1 2 g h µ B B ] , (10)
The hole spin-flip couples +1 and -2 states as well as -1 and +2 states. For each of these two couples of states we apply the Hamiltonian (10), which results in the following hole contribution to the 4 × 4 exciton Hamiltonian:

Ĥh ex =     -1 2 g h µ B B 0 0 β h k h e -iφ 0 1 2 g h µ B B β h k h e iφ 0 0 β h k h e -iφ -1 2 g h µ B B 0 β h k h e iφ 0 0 1 2 g h µ B B     (11)
4. In-plane magnetic field (Voight geometry)

If the magnetic field is applied in the plane, it splits electron and hole states polarized in the plane of the quantum wells. Suppose that the field is applied in the x-direction. In the electron and hole basis the Zeeman Hamiltonian is in this case:

Ĥe,h = - 1 2 g e,h µ B B σx = [ 0 -1 2 g e,h µ B B -1 2 g e,h µ B B 0
] .

(12) Note that the hole g-factor in plane of the quantum well is different from the g-factor in Faraday configuration, in general. This maps into the (+1, -1, +2, -2) exciton basis as a Zeeman Hamiltonian of the form:

ĤZ = - µ B B 2      0 0 g e g h 0 0 g h g e g e g h 0 0 g h g e 0 0      (13)

Exchange terms

Besides the contributions from electron and hole spin orbit interactions and Zeeman splitting, there may be a purely excitonic contribution to the Hamiltonian, which is composed from the Hamiltonian for bright excitons written in the basis (+1,-1):

Ĥb = E b Î -δ b σx = [ E b -δ b -δ b E b ] (14)
and the Hamiltonian for dark excitons written in the basis (+2,-2):

Ĥd = E d Î -δ d σx = [ E d -δ d -δ d E d ] , (15)
where Î is the identity matrix. The terms with δ b and δ d describe the splittings of bright and dark states polarized along x and y axes in the plane of the structure due to the long-range exchange interaction. The structural anisotropy is virtually inevitable even in the best quality epitaxially grown quantum wells. It arises from the reduced symmetry of heterointerfaces, from local strains and from islands of quantum well with fluctuations elongated in certain crystallographic directions. E b -E d is the splitting between bright (+1 and -1) and dark (+2 and -2) exciton states due to the short-range exchange interaction. In microcavities, this splitting is additionally enhanced due to the vacuum field Rabi splitting of exciton-polariton modes formed by bright excitons and a confined optical mode of the cavity [START_REF] See | Microcavities[END_REF] .

The origin of Eqs. ( 14) and ( 15) can be easily seen from the exciton Hamiltonian written in the basis of linear x and y polarizations. For example, for the bright excitons:

ĤXY = [ E b -δ b 0 0 E b + δ b ] (16) Ĥb = Ĉ-1 ĤXY Ĉ, (17)
where:

Ĉ = 1 √ 2 [ 1 1 i -i ] , Ĉ-1 = 1 √ 2 [ 1 -i 1 i ] ( 18
)
are the transformation matrices from the linear to circular polarization basis and vice versa [START_REF] See | Principles of optics[END_REF] . The same reasoning can be applied to the dark excitons as well.

The sum of Hamiltonians Ĥb and Ĥd , written in the 4 × 4 exciton spin basis is:

Ĥex ex =    E b -δ b 0 0 -δ b E b 0 0 0 0 E d -δ d 0 0 -δ d E d    (19)
Now, the full exciton Hamiltonian can be written as:

Ĥtot ex = Ĥe ex + Ĥh ex +H ex ex =     E b + 1 2 (g e -g h ) µ B B -δ b β e k e e -iφ β h k h e -iφ -δ b E b -1 2 (g e -g h ) µ B B β h k h e iφ β e k e e iφ β e k e e iφ β h k h e -iφ E d -1 2 (g e + g h ) µ B B -δ d β h k h e iφ β e k e e -iφ -δ d E d + 1 2 (g e + g h ) µ B B     . (20)
For the translational motion of an exciton as a whole particle the exciton momentum is given by P ex = (m e + m h ) v ex , where m e and m h are in-plane effective masses of an electron and of a heavy hole, respectively; v ex is the exciton velocity. Having in mind that the exciton translational momentum is a sum of electron and hole translational momenta given by P e,h = m e,h v e,h , where v e and v h are the electron and hole velocity, respectively, one can easily see that v h = v e = v ex . Having in mind that P ex = k ex and P e,h = k e,h , we have k ex = k h + k e , with k e = me me+m h k ex and k h = m h me+m h k ex . Thus, Ĥtot ex depends on the exciton center of mass wave-vector k ex and on the angle φ between this angle and one of the structure axes (e.g. (100)-axis).

It should be noted that in this consideration the wavevectors k ex , k h , k e are related to the translational motion of the exciton as a whole particle, with hole and electron as its constituents. The wave-vector of relative motion of the electron and hole "inside" the exciton is zero on average but may be important for each given moment of time. Recently, the effect of relative electron-hole motion on the spin-orbit effects of excitons has been analyzed by Vishnevsky et al [START_REF] Vishnevsky | [END_REF] . Their analysis confirms the presence of linear in k ex spin orbit terms in the exciton Hamiltonian introduced above.

B. Spin Density Matrix

Having constructed the Hamiltonian for excitons propagating with a wavevector k ex , we now consider the description of their spin state. We shall use the spin density matrix, ρ = |Ψ⟩ ⟨Ψ|, where Ψ = (Ψ +1 , Ψ -1 , Ψ +2 , Ψ -2 ) are the components of the exciton wavefunction projected onto the four spin states, (

|Ψ +1 ⟩ , |Ψ -1 ⟩ , |Ψ +2 ⟩ , |Ψ -2 ⟩).

Relation to Stokes' vectors and polarization degrees of light

The exciton spin DM is given by:

ρ = |Ψ⟩ ⟨Ψ| =    Ψ * +1 Ψ +1 Ψ * -1 Ψ +1 Ψ * +2 Ψ +1 Ψ * -2 Ψ +1 Ψ * +1 Ψ -1 Ψ * -1 Ψ -1 Ψ * +2 Ψ -1 Ψ * -2 Ψ -1 Ψ * +1 Ψ +2 Ψ * -1 Ψ +2 Ψ * +2 Ψ +2 Ψ * -2 Ψ +2 Ψ * +1 Ψ -2 Ψ * -1 Ψ -2 Ψ * +2 Ψ -2 Ψ * -2 Ψ -2 ,    (21)
The elements of the upper left quarter of this DM are linked to the intensity of light emitted by bright exciton states,

I = Ψ * +1 Ψ +1 + Ψ * -1 Ψ -1
, and to the components of the Stokes' vector, S x , S y and S z of the emitted light:

ρ 11 = I 2 + S z , (22)
ρ 12 = S x -iS y , (23)
ρ 21 = S x + iS y , (24)
ρ 22 = I 2 -S z . (25)
These expressions can be summarized more succinctly using the Pauli matrices as:

[ ρ 11 ρ 12 ρ 21 ρ 22 ] = I 2 Î + S.σ, (26)
where S = (S x , S y , S z ) is the Stokes' vector and we recall that Î is the identity matrix. Note that the trace of the spin density matrix is a number of particles in the system, which is not conserved because of the finite lifetime, in contrast with the full quantum optical density matrix which has the trace equal to unity. Often when studying the polarization structure of fields with non-uniform intensity, it is useful to compare the polarization degrees of emitted light, which can be given by normalizing the Stokes' vectors to the light intensity. The circular polarization degree is:

ρ c = 2S z I = ρ 11 -ρ 22 ρ 11 + ρ 22 . (27)
The horizontal-vertical linear polarization degree is:

ρ l = 2S x I = ρ 12 + ρ 21 ρ 11 + ρ 22 . (28)
The linear polarization degree measured in the diagonal axes (also referred to as a diagonal polarization degree) is given by:

ρ d = 2S y I = i ρ 12 -ρ 21 ρ 11 + ρ 22 . ( 29
)

Liouville equation

The dynamics of the DM is given by the quantum Liouville equation:

i dρ dt = [ Ĥtot ex , ρ] , (30)
where the Hamiltonian is composed from the electron, hole and exciton contributions given by Eqs. ( 6), ( 11) and ( 19) (considering the Faraday magnetic field configuration).

So far, we have neglected all relaxation or scattering processes in the system. The commonly used way to account for these processes is through the introduction of a phenomenological Lindblad superoperator to the Liouville equation:

i dρ dt = [ Ĥtot ex , ρ] -L (ρ) , (31)
where the Lindblad superoperator is introduced as:

L (ρ) = i    ρ 11 /τ b ρ 12 /τ b ρ 13 /τ c ρ 14 /τ c ρ 21 /τ b ρ 22 /τ b ρ 23 /τ c ρ 24 /τ c ρ 31 /τ c ρ 32 /τ c ρ 33 /τ d ρ 34 /τ d ρ 41 /τ c ρ 42 /τ c ρ 43 /τ d ρ 44 /τ d    . ( 32
)
τ b is the bright exciton decoherence time, τ d is the dark exciton decoherence time and τ c is the characteristic decoherence time of processes between dark and bright excitons. Note that dissipation may be crucial in the description of exciton spin currents in realistic systems. In particular, within this formalism, in the presence of dissipation, the current conservation and flux conservation conditions become valid if completed by exciton generation and decay in the continuity equation. In the rest of this manuscript we shall neglect dissipation to simplify the model system and to clarify the physical mechanisms which govern the characteristics of exciton spin currents. Namely, we shall assume L (ρ) = 0. We stress that, in the experiment, the magnitude of predicted spin currents may be reduced by dissipation. The formalism described so far, in sections II A and this section (II B), has been successfully applied in the description of spin transport in gases of cold excitons in coupled GaAs/AlGaAs quantum wells 8 . In this work, cold excitons are generated within localized spots and then fly away ballistically in radial directions. The elements of the DM, ρ ij , are dependent on the distance from the excitation spot r = v ex t at time t and the polar angle, φ. The propagation speed v ex = k ex / (m e + m h ). In this work, when solving the Liouville equation introduced above, we shall refer to the experimental configuration of Ref. 8. In particular, this implies a specific choice of the initial conditions for Eq. ( 30): we shall assume that at zero time excitons are not moving. We shall assume that they populate the eigenstates of the exciton Hamiltonian Ĥtot ex taken with k ex = 0 following a thermal distribution with a temperature T . We shall assume that once created in the equilibrium state, the excitons start moving apart in the radial direction. Thus, implicitly, we account for a non-linear effect: dipole-dipole repulsion of excitons which makes them acquire a certain inplane velocity v ex . This non-linearity is crucial to move the system out of equilibirum. The rest of exciton propagation and spin dynamics is modelled using the linear equation (30), where the Hamiltonian Ĥtot ex contains now the off-diagonal terms proportional to k ex .

C. Electron and hole spin currents

We shall normalize exciton, electron and heavy-hole functions to unity, namely:

Ψ * +1 Ψ +1 + Ψ * -1 Ψ -1 + Ψ * +2 Ψ +2 + Ψ * -2 Ψ -2 = 1 (33) Ψ * e,+ 1 2 Ψ e,+ 1 2 + Ψ * e,-1 2 Ψ e,-1 2 = 1 (34) Ψ * h,+ 3 2 Ψ h,+ 3 2 + Ψ * h,-3 2 Ψ h,-3 2 = 1 (35)
Now, the exciton spin DM ( 21) can be represented in terms of electron and hole wavefunctions as:

ρ =      Ψ * e,-1 2 Ψ e,-1 2 Ψ * h,+ 3 2 Ψ h,+ 3 2 Ψ * e,+ 1 2 Ψ e,-1 2 Ψ * h,-3 2 Ψ h,+ 3 2 Ψ * e,+ 1 2 Ψ e,-1 2 Ψ * h,+ 3 2 Ψ h,+ 3 2 Ψ * e,-1 2 Ψ e,-1 2 Ψ * h,-3 2 Ψ h,+ 3 2 Ψ * e,-1 2 Ψ e,+ 1 2 Ψ * h,+ 3 2 Ψ h,-3 2 Ψ * e,+ 1 2 Ψ e,+ 1 2 Ψ * h,-3 2 Ψ h,-3 2 Ψ * e,+ 1 2 Ψ e,+ 1 2 Ψ * h,+ 3 2 Ψ h,-3 2 Ψ * e,-1 2 Ψ e,+ 1 2 Ψ * h,-3 2 Ψ h,-3 2 Ψ * e,-1 2 Ψ e,+ 1 2 Ψ * h,+ 3 2 Ψ h,+ 3 2 Ψ * e,+ 1 2 Ψ e,+ 1 2 Ψ * h,-3 2 Ψ h,+ 3 2 Ψ * e,+ 1 2 Ψ e,+ 1 2 Ψ * h,+ 3 2 Ψ h,+ 3 2 Ψ * e,-1 2 Ψ e,+ 1 2 Ψ * h,-3 2 Ψ h,+ 3 2 Ψ * e,-1 2 Ψ e,-1 2 Ψ * h,+ 3 2 Ψ h,-3 2 Ψ * e,+ 1 2 Ψ e,-1 2 Ψ * h,-3 2 Ψ h,-3 2 Ψ * e,+ 1 2 Ψ e,-1 2 Ψ * h,+ 3 2 Ψ h,-3 2 Ψ * e,-1 2 Ψ e,-1 2 Ψ * h,-3 2 Ψ h,-3 2 ,      (36)
This representation allows us to obtain useful links between the elements of exciton, electron and hole density matrices, in particular:

ρe = |Ψ e ⟩ ⟨Ψ e | = [ Ψ * e,+ 1 2 Ψ e,+ 1 2 Ψ * e,-1 2 Ψ e,+ 1 2 Ψ * e,+ 1 2 Ψ e,-1 2 Ψ * e,-1 2 Ψ e,- 1 2 ]
= [ ρ 22 + ρ 33 ρ 24 + ρ 31 ρ 13 + ρ 42 ρ 11 + ρ 44 ] (37) ρh
= |Ψ h ⟩ ⟨Ψ h | = [ Ψ * h,+ 3 2 Ψ h,+ 3 2 Ψ * h,-3 2 Ψ h,+ 3 2 Ψ * h,+ 3 2 Ψ h,-3 2 Ψ * h,-3 2 Ψ h,-3 2 ] = [ ρ 11 + ρ 33 ρ 14 + ρ 32 ρ 23 + ρ 41 ρ 22 + ρ 44 ] (38)
We know that the components of electron and hole density matrices are linked with the projections of electron and hole spins as: ρe =

[ 1 2 + S e,z S e,x -iS e,y S e,x + iS e,y

1 2 -S e,z ] (39) ρh
= [ 1 2 + S h,z S h,x -iS h,y S h,x + iS h,y 1 2 -S h,z ] , (40)
where, for the heavy hole we have assigned spin +1/2 to the state +3/2 and spin -1/2 to the state -3/2 accounting for the orbital momentum of these states of +1 and -1, respectively. The z-component of the spin polarization carried by electrons can now be expressed as:

S e,z = (ρ 22 + ρ 33 -ρ 11 -ρ 44 ) /2. ( 41
)
Similarly, the z-component of the spin polarization carried by holes can now be expressed as:

S h,z = (ρ 11 + ρ 33 -ρ 22 -ρ 44 ) /2. ( 42
)
The in-plane component of electron and hole spins can be extracted from the off-diagonal elements of the DM. Namely, the x-component of electron spin is given by:

S e,x = (ρ 13 + ρ 31 + ρ 24 + ρ 42 ) /2, (43)
while the x-component of the hole spin is given by:

S h,x = (ρ 14 + ρ 23 + ρ 32 + ρ 41 ) /2. (44)
The y-component of electron spin is given by: S e,y = i (-

ρ 13 + ρ 31 + ρ 24 -ρ 42 ) /2, (45)
while the y-component of the hole spin is given by:

S h,y = i (ρ 14 -ρ 23 + ρ 32 -ρ 41 ) /2, (46)

III. NUMERICAL RESULTS IN THE DENSITY MATRIX FORMALISM

Figure 1 shows the numerical results obtained within the DM formalism for a model system with the same parameters as those of coupled double quantum wells studied in Ref. 8. The parameters are summarized in Table I. I. (e) Same but at kex = 5µm -1 . The source area was taken circular with a radius of 4µm.

The dispersion of bright and dark exciton modes obtained by diagonalisation of the Hamiltonian Ĥtot ex (20) is shown in Fig. 1(a) by green solid lines. The momentum has been chosen along the x direction, but the anisotropy of the band structure remains small. The initial splittings of dark and bright states makes these dispersion curves qualitatively different from those presented by Vishnevsky et al [START_REF] Vishnevsky | [END_REF] . Note also, that Ref. ( 40) accounts for the exciton kinetic energy which we neglect in Ĥtot ex , leaving only the spin-dependent contributions to the energy. The numerically calculated linear and circular polarization degrees are shown in Figs. 1(c) and(d), respectively. To stay close to the experimental conditions of Ref. 8 , we have chosen as initial condition the circular source area with a radius of 4 µm, where cold excitons are generated within localized spots. In Fig. 1(c), four lobes are unambiguously observable in the pattern of the linear and diagonal (not shown) polarizations. This pattern is a consequence of the Dresselhaus spin-orbit coupling for electrons, and it is characteristic of the chosen initial state at the source: four split eigen-states with zero inplane wave-vector are occupied. Relative occupation of these states corresponds to the Boltzmann distribution at temperature T = 0.1 K. This choice of initial conditions leads to a variety of polarization patterns observed in experiment 8 . In the absence of the damping this pattern is periodic and infinite in the radial direction.

In order to reveal the mechanism of formation of the polarization patterns it is instructive to consider a simplified version of the Hamiltonian 20. In the absence of magnetic field, zero Dresselhaus effect for holes, zero splitting of bright excitons and zero splitting E b -E d be-tween bright and dark excitons this Hamiltonian can be rewritten as:

Ĥ = -δ d      0 0 ξe -iφ 0 0 0 0 ξe iφ ξe iφ 0 0 1 0 ξe -iφ 1 0      (47)
where ξ = -β e k e /δ d . In other words, here we only take into account linear splitting of dark exciton states which inevitably results from structural anisotropy even in the best quality samples, and the Dresselhaus field acting on the spin of electron, bound to the hole. We will see that at sufficiently low temperature these two ingredients provide the in-plane asymmetry that ultimately results in the formation of linear polarisation patterns. Indeed, the eigenvalues of this Hamiltonian can be obtained analytically :

E = ± 1 2 δ d ± 1 2 δ d √ 1 + 4ξ 2 .
This corresponds to two dispersion branches at low energy ∝ ±ξ 2 and two branches at high energy ∝ ±(1 + ξ 2 ) for ξ ≪ 1. These branches are shown by red squares in Fig. 1(a).The eigenvectors, starting from the lowest energy and taking δ d < 0, can be approximated for small ξ by:

     ξe -iφ -ξe iφ -1 1      ,      e -iφ e iφ -ξ -ξ      ,      e -iφ -e iφ ξ -ξ      ,      ξe -iφ ξe iφ 1 1      . (48)
As mentioned already, at low temperatures k B T ≪ |δ d |, the lowest energy state with zero momentum is given by [0, 0, -1, 1] and is a linearly polarized dark exciton. By linearly polarised dark exciton we mean the dark state which has a dipole moment oriented in a certain way, and which has a zero spin projection to the grows axis of the structure. After an acceleration due to dipole repulsion, this initial state is no longer an eigenstate and oscillates as a function of time between the two eigenstates of the Hamiltonian with which the initial state is not orthogonal, namely the first and the third eigenstates listed in Eqs. 48. Among these two states, only the third gives a significant contribution to the observed polarization as it is essentially "bright" (has large projections to +1 and -1 exciton states). The linear and diagonal polarizations originating from this state are readily given by -cos(2φ) and -sin(2φ) respectively. This reproduces the essential features of the numerical results for the linear polarization pattern, as one can see comparing the images calculated with the reduced Hamiltonian (Fig. 1(b)) and the full Hamiltonian (Fig. 1(c)). In the particular case considered here (initial state formed essentially by linearly polarized dark excitons) the Dresselhaus spin-orbit term for electrons leads to formation of the linear polarization vortex: the polarization plane is always perpendicular to the wave vector direction. The linear polarization vortex has been observed experimentally by High et al 8,9 . Rapid oscillations in radial direction are due periodical change in the occupation of mainly dark and mainly bright states, indicated by double-ended arrows in Fig. (1(b)). We should emphasize that with a proper reordering of the basis vectors, the simplified Hamiltonian 47 is analytically equivalent to the Hamiltonian of bilayer graphene. In its simplest expression, the Hamiltonian of bilayer graphene can be written as 41 :

       0 v F ke -iφ 0 0 v F ke iφ 0 t ⊥ 0 0 t ⊥ 0 v F ke -iφ 0 0 v F ke iφ 0        . ( 49
)
where v F is the velocity, k the amplitude of the momentum, φ the angle between the momentum and the x axis, and t ⊥ is the main coupling term between the two graphene layers, which is given by hopping between two carbon atoms that are superimposed. The four coefficients of the associated wavevectors correspond to the probability amplitudes on the two independent sublattices of the two graphene layers. By permutation of the basis vectors:1 → 1, 2 → 3, 3 → 4, 4 → 2, the Hamiltonian is rewritten as:

       0 0 v F ke -iφ 0 0 0 0 v F ke iφ v F ke iφ 0 0 t ⊥ 0 v F ke -iφ t ⊥ 0        . ( 50
)
which is equivalent to Eq. 47 with t ⊥ ≡ -δ d and v F k/t ⊥ ≡ ξ. The 4 × 4 Hamiltonian of bilayer graphene is often restricted to the subspace of the two low energy bands. From Eq. 48, the corresponding eigenstates are:

[e iφ , ±e -iφ ]. ( 51
)
We can then define a pseudospin vector, which represents the relative phase between the two components of the wavevectors. From Eq. 51, the pseudospin rotates two times when the particle wavevector undergoes one full rotation.

Adapting the terminology used for graphene, the polarization pattern observed in Fig. 1(b,c) would be nothing else than the fingerprint of the "pseudospin" rotation in the exciton system. In the context of the exciton system under consideration, the phase φ corresponds to the angle between the exciton polarization vector and the chosen x-axis of the structure.

The build up of circular polarization requires introduction in the model of the splitting between linearly polarized bright exciton states δ b . This splitting acts as an effective magnetic field applied to the Stokes vector of light emitted by bright excitons, S = (S x , S y , S z ). In particular, if the x-polarized exciton state has a lower energy than the y-polarized exciton state, it creates an effective magnetic field in the x-direction which rotates the Stokes vector in the yz plane. This converts the diagonal polarization to the circular polarization and leads to the appearance of right-and left-circularly polarized sectors in the polarization map of exciton emission. Separation of spins due to linear-to-circular polarization conversion is known in exciton-polariton systems as the Optical spin Hall effect. This effect was theoretically predicted in Ref. 12 and experimentally observed in polaritonic 10 and excitonic systems 8 . For microcavity polaritons it may be described in terms of beats between TE and TM polarized polariton modes, while in the exciton system studied here the effect is more complex due to the mixture of four nearly degenerate (dark and bright) exciton states. A recent theoretical paper [START_REF] Vishnevsky | [END_REF] predicts the skyrmion formation in this case. Here we concentrate on the circular polarization patterns appearing due to the beats between linearly polarized exciton states mixed by the exchange interaction. Note that the circular polarization pattern is strongly sensitive to the chosen exciton wave-vector, which governs the energies of the four involved eigen-states. Figure 1(d) is calculated assuming k ex = 15.3µm -1 . This corresponds to the crossing point of the dispersion branches associated to the first and second exciton eigen-states (see Fig 1(a)). For comparison, Fig. 1(e) shows the circular polarization pattern calculated with k ex = 5µm -1 . One can see that the four lobes pattern of circular polarization is washed out, because the lowest energy state of the system remains essentially dark. Rapid oscillations that show up, have the period determined by the splitting between the two lowest states.

Finally, let us underline that the dispersion curves shown in Fig. 1(a) do not take into account the kinetic energy of excitons. The kinetic energy would shift all curves up by K = 2 k 2 ex 2mex . Note that this does not affect the splittings between exciton eigenstates and would not affect the spin dynamics of excitons. In the spin density formalism developed above we assign to all excitons the same kinetic energy, K. In a realistic system, the kinetic energy may be spread, in which case averaging of the obtained polarization patterns over k ex may be needed. This averaging would smooth the fast oscillations seen on the images Fig. 1 (c,d,e,f). In the next section devoted to the non-linear spin dynamics we will fully take into account the kinetic energy of propagating excitons.

The distribution of in-plane projections of electron and hole spins for the same choice of parameters as above is shown in Figs. 2(a-h). The left panels show electron spins and the right panels show the hole spins. The direction of in-plane spin component is shown by arrows, while the length of each arrow is proportional to the computed value of the transverse spin component. The upper panels, Figs. 2(a,b), show the spin distributions in the absence of a magnetic field. In this case the electron and hole spins are oriented along the effective Dresselhaus fields which are oriented differently for electrons and heavy holes, as we have discussed in the previous section. I, at kex = 15.3µm -1 .The source area was taken circular with a radius of 4µm.

The decrease of the in-plane spin component upon propagation corresponds to the build-up of the z-component of electron and hole spins (Figure not shown), due to rotation of the exciton spin around the effective magnetic field. The magnetic field strongly changes the spin distribution in real space. The spin textures become strongly anisotropic in the case of in-plane (x-or y-oriented) magnetic field. Note, that the in-plane isotropy in the system is broken by the splitting between x-and y-polarized exciton states, which is why switching of the magnetic field between x and y axes strongly affects the distribution of electron and hole spins. It should be noted also that electron and hole in-plane spin textures can hardly be observed directly in optical experiments. However, they can be deduced from fitting the exciton polarization maps, e.g., using the formalism described above. I, kex = 15.3µm -1 . The source area was taken circular with a radius of 4µm.

Figure 3 shows how the magnetic field affects spatial patterns of linear (a-d) and circular (e-h) polarization. Switching the magnetic field orientation between the x-, y-and z-axes one can dramatically affect the polarization patterns. Having in mind that the exciton polarization patterns can be directly observed in near-field photoluminescence experiments, fitting of these patterns to the experimental data would allow extracting the Dresselhaus constants and exciton exchange splittings, which, in turn, allow to restore electron and hole spin textures 8 .

Figure 4(a-d) illustrates a peculiar regime where the Dresselhaus fields for electrons and holes are taken to be zero and there is no magnetic field applied, but electrons are subjected to the Rashba field (the Rashba field for heavy holes is zero). Figure 4(a) shows the electron spin distribution in space, where the spins are clearly aligned along the Rashba field force lines. Interestingly, the hole spins become polarized as well, while no effective field acts on them (Figure 4(b)). This is an illustration of the exciton effect: bound in excitons by Coulomb interaction and subject to the exchange induced exciton effects, the holes acquire in-plane spin polarization. The non-zero spin polarization of heavy holes is possible due to the exciton exchange effects. Figure 4(c,d) shows the resulting linear and circular exciton polarization patterns. One can see that the Rashba effect induces polarization patterns strongly different from those induced by the Dresselhaus effect, which is why from the shape of polarization patterns one can conclude on the nature of spin-orbit coupling in the system.

IV. NON-LINEAR SPIN DYNAMICS OF PROPAGATING EXCITONS AND EXCITON-POLARITONS

In the previous section we operated with a spin density matrix which is very convenient for the description of partially coherent and partially polarized exciton gases. The quantum Liouville equation ( 31) is a very efficient tool for the description of effects linear in the exciton density. On the other hand, one cannot straightforwardly incorporate non-linear interaction terms in this equation. The treatment of non-linear effects in a partially coherent system is a non-trivial task. Much simpler is the treatment of non-linear effects in a perfectly coherent system, such as a condensate at zero temperature. In this case, the ensemble of excitons can be described by a single 4-component wave function Ψ = (Ψ +1 , Ψ -1 , Ψ +2 , Ψ -2 )

T . The linear dynamics of this wave-function for ballistically propagating excitons having a wavevector k ex is described by the Schrödinger equation:

i d dt |Ψ⟩ = Ĥ |Ψ⟩ , ( 52
)
where the Hamiltonian is the same as in Eq. ( 20). This equation represents a set of four coupled linear differential equations for four exciton spin components. Nonlinear effects lead to the condensate evolution in real and reciprocal space.

From now on we shall consider the exciton spin dynamics in real space (2D), so that the wave function Ψ will become coordinate-dependent and will not be restricted to one single value of k ex . The non-linear interaction terms for multi-component exciton gases are introduced and discussed in detail in Ref. 42. Here we expand Eq. ( 52) by introducing the kinetic energy (to describe the real space dynamics) and the interaction terms. On the other hand, we neglect the magnetic field, for simplicity. This results in a system of four non-linear Schrödinger or GP equations [42][43][44] :

i dΨ +1 dt = - 2 ∇2 2m ex Ψ +1 + β e m e m ex ( kx -i ky ) Ψ +2 + β h m h m ex ( kx -i ky ) Ψ -2 + α 1 |Ψ +1 | 2 Ψ +1 + α 2 |Ψ -1 | 2 Ψ +1 + α 3 |Ψ +2 | 2 Ψ +1 + α 4 |Ψ -2 | 2 Ψ +1 + W Ψ * -1 Ψ +2 Ψ -2 , ( 53
) i dΨ -1 dt = - 2 ∇2 2m ex Ψ -1 + β e m e m ex ( kx + i ky ) Ψ -2 + β h m h m ex ( kx + i ky ) Ψ +2 + α 1 |Ψ -1 | 2 Ψ -1 + α 2 |Ψ +1 | 2 Ψ -1 + α 3 |Ψ -2 | 2 Ψ -1 + α 4 |Ψ +2 | 2 Ψ -1 + W Ψ * +1 Ψ +2 Ψ -2 , ( 54
)
i dΨ +2 dt = - 2 ∇2 2m ex Ψ +2 + β e m e m ex ( kx + i ky ) Ψ +1 + β h m h m ex ( kx -i ky ) Ψ -1 + α 1 |Ψ +2 | 2 Ψ +2 + α 2 |Ψ -2 | 2 Ψ +2 + α 3 |Ψ +1 | 2 Ψ +2 + α 4 |Ψ -1 | 2 Ψ +2 + W Ψ * -2 Ψ +1 Ψ -1 , ( 55
) i dΨ -2 dt = - 2 ∇2 2m ex Ψ -2 + β e m e m ex ( kx -i ky ) Ψ -1 + β h m h m ex ( kx + i ky ) Ψ +1 + α 1 |Ψ -2 | 2 Ψ -2 + α 2 |Ψ +2 | 2 Ψ -2 + α 3 |Ψ -1 | 2 Ψ -2 + α 4 |Ψ +1 | 2 Ψ -2 + W Ψ * +2 Ψ +1 Ψ -1 . (56)
Here kx,y = -i ∇x,y , m ex = m e + m hh . To make this system more compact we have omitted the terms describing exchange induced exciton splittings given by the Hamiltonian (20). We do not discuss here the nature and value of the interaction constants α 1,2,3,4 and W . In the system of indirect excitons in coupled GaAs/AlGaAs quantum wells, as a zeroth approximation, one can take

α 1 = α 2 = α 3 = α 4 .
Note also that in microcavities, where the lower exciton-polariton mode is strongly decoupled from dark excitons, the dark exciton states may be almost empty at low temperatures. If this is the case, the spin dynamics of the exciton-polariton condensate is given by the first two of the four GP equations (53 and 54) with α 3,4 = W = 0. The remaining constants α 1,2 have been widely discussed in literature 45 .

The GP equations are widely used for the description of coherent propagation of exciton-polaritons in microcavities 46 . They allow for the studying of interesting topology effects such as: quantum vortices [47][48][49][50][51][52] ; halfquantum vortices [53][54][55] ; bright 56 and dark 51,57-60 solitons.

The polarization of light emitted by an exciton or exciton-polariton condensate can be obtained as:

ρ c = 2S z I = |Ψ +1 | 2 -|Ψ -1 | 2 |Ψ +1 | 2 + |Ψ -1 | 2 , (57)
ρ l = 2S x I = 2Re { Ψ * +1 Ψ -1 } |Ψ +1 | 2 + |Ψ -1 | 2 , (58)
ρ d = 2S x I = - 2Im { Ψ * +1 Ψ -1 } |Ψ +1 | 2 + |Ψ -1 | 2 . (59)
These expressions easily follow from the definition of the spin density matrix. A significant limitation of the GP equations as a theoretical tool is that they assume a coherent state of the system. If one is interested in the spin structure of the zero temperature ground state of excitons in a Bose-Einstein condensate 43 , then this assumption is fulfilled by definition. However, in real systems there is an incomplete coherence that, strictly speaking, requires a description of statistical mixtures, perhaps involving density matrices. Furthermore, Eqs. (53-56) have been written assuming an infinite lifetime for the particles (be they excitons or exciton-polaritons), which is never the case of real nonequilibrium systems. Often pumping and radiative decay terms are introduced into Eqs. (53-56) phenomenologically 30,44 . While in the case of a resonant coherent pump, one can imagine that the exciton/exciton-polariton distribution inherits coherence directly, it is less obvious how an incoherent pump can be modelled. A phenomenologi-cal model introduced by Wouters and Carusotto 61 of incoherent pumping has allowed the modelling of the first order coherent fraction observed in many experimental configurations based on condensation 26,[62][63][64][65] . While the GP model cannot model the phase transition during formation of a condensate and/or superfluid, it can offer a suitable description of spin currents once spatial coherence has formed.

Here we will consider a localized source of the fourcomponent indirect exciton system, as corresponds to the localized bright spot sources generating exciton condensates 9 . We focus our attention on the possible spin polarization textures of coherent excitons propagating away from the source. We do not attempt to describe the partially coherent state within the source, noting that in exciton-polariton systems spin currents have been generated from both coherent 66 and incoherent 11 tightly focused spots utilizing the optical spin Hall effect 12 in a similar way. Since the Gross-Pitaevskii equations are only valid for coherent excitons, we restrict the exciton wavefunction to lie outside of the source area. The effect of the source is then characterized by the chosen boundary condition along the edge of the source area. Given that dark excitons have lower energy than bright excitons [due to the exchange splitting of Eq. ( 19)], it is reasonable to expect the source to provide linearly polarized dark excitons. By fixing the values of the exciton wavefunction along the edges of the source area, which is assumed circular, to such a distribution the boundary condition acts as an effective source for the exciton wavefunction outside the source area.

The indirect excitons are known to have very long lifetime, typically in the range of 10ns to 10µs 67 . This allows them to cover distances of a few hundreds of µm with negligible loss 8 . Consequently, when we focus on the behaviour of excitons in a small (10 × 10µm 2 ) area around the source, the main loss of excitons is caused by their escape from the area of interest rather than their decay (recombination). To model the spin currents we thus employ an absorbing boundary condition to allow the solution of Eqs. (53-56) in a finite area. This allows a balance between source and loss to achieve a steady state of the non-equilibrium system (for a continuous pump), where both source and loss appear as boundary conditions. 8 The steady state solution of the system is independent of the initial condition.

Exciton intensity and polarization distributions calculated within the GP approach are shown in Figs. 5 and7. Note that all the plotted quantities are spatially averaged over 1.5µm to account for the typical resolution of experimental setups 8 . In principle excitons can display features on the scale of the de Broglie wavelength, λ = 2π/|k ex |. Such features are on the sub-micron scale and are far beyond experimental resolution. We note also that while in the DM approach we could consider excitons having a fixed radial velocity, in the GP approach we necessarily cover the whole range of wavevectors and propagation velocities (the dispersion obtained from the GP approach is shown in Fig. 6). Also, exciton-exciton interactions modify the exciton dispersion. For all these reasons, the DM and GP approaches cannot give identical results, in principle. The DM approach is suitable for the description of both coherent and partly coherent exciton gases, while the GP approach better catches the dispersive propagation features and accounts for nonlinear effects. Both approaches are complementary, and it is instructive to compare the results, obtained within these two models.

The exciton density decreases as excitons travel away from the source (Fig. 5a). This is not due to exciton recombination, which is expected to be very slow, but more simply due to the spreading out of excitons in all directions. The intensity spread need not be perfectly circularly symmetric due to the presence of the spin-orbit (Dresselhaus terms), which can introduce a directional dependence of the exciton velocity. Fig. 5b shows the exciton brightness degree, defined as:

ρ b = |Ψ +1 | 2 + |Ψ -1 | 2 -|Ψ +2 | 2 -|Ψ -2 | 2 |Ψ +1 | 2 + |Ψ -1 | 2 + |Ψ +2 | 2 + |Ψ -2 | 2 . ( 60
)
This quantity represents the degree to which the bright exciton density exceeds the dark exciton density. There is a conversion of dark to bright excitons as they spread out from the source, which can be expected from the presence of Dresselhaus coupling terms. 60)] in the absence of a magnetic field. The parameters were the same as those used for the DM calculations, given in Table I, with: B = 0; W = 0.2α; mex was taken as 0.21 of the free electron mass. The source area was taken circular with a radius of 1µm. The images are presented with spatial averaging over 1.5µm. The scattering parameter α and intensity at the source center were chosen such that the interaction energy, α

( |Ψ+1| 2 + |Ψ-1| 2 + |Ψ+2| 2 + |Ψ-2| 2 )
= 1µeV (being comparable to the other energy scales in the system, we are in a nonlinear regime). The absorbing boundary condition used in calculations appears outside of the plotted range, at a radius of 15µm from the source center.

The polarization distribution is shown in Fig. 7 and can be significantly influenced by non-linearity in the system. The left-hand plots show the results for negligible nonlinearity (α = 0; W = 0), which is equivalent to a weak pump intensity. Here the polarization distributions are qualitatively similar to those calculated in the density matrix formalism. In analogy to the (intrinsic) spin Hall effect 2,3 and the optical spin Hall effect 10,12 the presence of the spin-orbit coupling terms introduces a directional dependence of the polarization. The patterns of the polarization degrees divide into quadrants. Some quantitative differences with the DM calculations appear due to the presence of different wavevectors.

The right-hand plots show the case of a moderate nonlinearity, with interaction strength comparable to the other energy scales of the system. The most drastic effect is on the circular polarization degree, which becomes higher and each quadrant of circular polarization divides further giving an eight-lobed pattern to the polarization degree. The interaction terms that we have introduced are all spin conserving and it can be noted that we have considered the spin isotropic case. Even the W nonlinear interaction term, which allows the interconversion of bright and dark exciton pairs does not appear to directly change the spin polarization, conserving both circular and linear polarizations upon scattering. Still, the nonlinear interaction terms can have a drastic effect on the polarization structure. This is because they are able to shift (renormalize) the dispersion branches in the system. Given that the potential energy of excitons is fixed by their interaction energy at the source and that this energy is converted into kinetic energy at distances away from the source, any shifts in the dispersion branches can change the wavevector of propagating excitons. Even if the nonlinear induced shifts of the dispersion branches were not polarization dependent, a change in the wavevector of an exciton can allow it to experience a different effect from the k-dependent spinorbit coupling terms. In this way, richer structures can appear in the nonlinear regime. Note, that the build-up of circular polarisation clearly seen in Fig. 7 (e,f) would not yield 100% circularly polarised excitons, as it might be expected in the ideal case of the optical spin Hall effect for exciton-polaritons 12 . In our system, the precession of electron and hole spins has different frequencies, and the interplay between dark-to-bright and linear-tocircular polarisation conversion prevents formation of a purely circularly polarised state. The current of each of these densities is given by a product of the exciton speed and the corresponding density:

j a = k ex m ex ρ jj (61)
with j = 1, 2, 3, 4 for a = +1, -1, +2, -2, respectively. Experimentally, one can measure the magnetization current associated with the exciton density current. The magnetization carried by propagating excitons can be found as:

M z = - µ B 2 [(g h -g e ) (ρ 11 -ρ 22 ) + (g h + g e ) (ρ 33 -ρ 44 )] (62
) This expression is obtained having in mind that an electron with a spin projection on the z-axis of ±1/2 contributes to the magnetization projection on the z-axis ∓ µB 2 g e , and a heavy hole with the spin projection of ±3/2 contributes to the magnetization ∓ µB 2 g h . Hence, the magnetization (spin) current produced by the excitons having a wave-vector k ex will be given by:

j M (k ex ) = - µ B k ex 2m ex [(g h -g e ) (ρ 11 (k ex ) -ρ 22 (k ex )) + (g h + g e ) (ρ 33 (k ex ) -ρ 44 (k ex ))] ( 63
)
The total magnetization current in the exciton gas can be obtained by integration over all wave-vectors:

j tot M = - A (2π) 2 ∫ j M (k ex )dk ex ( 64
)
Here A is the area of the sample. This current may be detected, for example, by spatially resolved Kerr rotation spectroscopy. Figure 8 shows the spin current density j M / (2πr) calculated for the system of indirect excitons which we considered above in the absence of external magnetic field (a) and in the presence of a magnetic field of 5T oriented normally to the plane of the structure (b) and along xand y-axes (c,d). In all cases the current intensity decreases as one moves away from the excitation spot, as the exciton density decreases inversely with the radius r. One can see that the total spin (magnetization) of propagating excitons experiences oscillations and has a strong angular dependence. The spin currents are suppressed by an in-plane magnetic field, which is not surprising: j M describes propagation of the normal-to plane spin component, which is strongly reduced by an in-plane magnetic field. The images show the total spin carried by both bright and dark excitons. They do not directly correspond to the polarized photoluminescence map for two reasons: first, dark excitons do not contribute to the photoluminescence; second, the polarization degree of photoluminescence does not experience the 1/r decay characteristic of the total spin density. On the other hand, the images presented in Figure 8 do correspond to the signal of spatially resolved Kerr or Faraday rotation, which is sensitive to the normal-to-plane magnetization.

VI. SPIN CURRENTS IN EXCITON CONDENSATES

The approach formulated above can be extended to the description of spin currents in coherent exciton (or exciton-polariton) condensates accounting for particleparticle interactions. In this case we need to replace the momentum k ex by a momentum operator p = -i ∇ and the diagonal components of the DM ρ 11 , ρ 22 , ρ 33 , ρ 44 by the exciton densities

|ψ +1 | 2 , |ψ -1 | 2 , |ψ +2 | 2 , |ψ -2 | 2 ,
respectively, in the expressions 61 and 64. In this case the density currents become:

j α = -i m ex Ψ * α ∇Ψ α , (65)
and the total magnetization current can be expressed as:

j tot M = iµ B 2m ex [ (g h -g e ) ( Ψ * +1 ∇Ψ +1 -Ψ * -1 ∇Ψ -1
)

+ (g h + g e ) ( Ψ * +2 ∇Ψ +2 -Ψ * -2 ∇Ψ -2 )] (66)
= - µ B 2 [(g h -g e ) (j +1 -j -1 ) + (g h + g e ) (j +2 -j -2 )] (67)
The distribution of the spin density current, j +1 , is shown in Fig. 9 (a). The current density propagates outward from the source in all directions, decreasing in intensity. The apparent rotation of the current density is a nonlinear effect coming from the interactions in the system. The other spin density currents, j -1 , j +2 and j +2 , display a similar behaviour.

The magnetization current is shown in Fig. 9 (b). The current is stronger closer to the source, where the intensities are stronger. The magnetization current is predicted to rotate around the source. One can also introduce the spin conductivity tensor linking the components of the density current ( 65) with the gradient of potential acting upon each of the exciton spin components:

j α,l = σ l,m α,β ∇U β,m , (68)
where l = x, y and m = x, y indicate the in-plane projections of the current and potential gradient, respectively. One can see that σ l,m α,β is a 64-component tensor in the general 2D case. The origin of the potential gradient ∇U β,m needs to be discussed separately. ∇U β,m can originate from the gradient of the quantum well width, gradient of the barrier height, or it can be induced by excitons themselves due to e.g. dipole-dipole repulsion. Indirect excitons have built-in dipole moments, the laterally modulated external electric field in the z-direction can create an in-plane potential landscape and, in turn,∇U for them. This was used in studies of transport of indirect excitons in various electrostatic potential landscapes including potential energy gradients 28,68,69 , circuit devices [70][71][72] , traps 73 lattices 74,75 , moving latticesconveyers 76 , and narrow channels 72,77,78 .

VII. POLARIZATION CURRENTS

Spatially resolved measurements of the polarization degrees ρ c , ρ l and ρ d of light emitted by excitons give access to the exciton polarization currents. In terms of the DM formalism, they can be defined as products of the exciton speed and the corresponding polarization degree: The polarization currents in an exciton condensate can be found from the GP equations (53-56) as:

j c = - i m ex ( Ψ * +1 ∇Ψ +1 -Ψ * -1 ∇Ψ -1
)

|Ψ +1 | 2 + |Ψ -1 | 2 (73)
j l = - i m ex ( Ψ * +1 ∇Ψ -1 + Ψ * -1 ∇Ψ +1 ) |Ψ +1 | 2 + |Ψ -1 | 2 (74)
j d = -m ex ( Ψ * +1 ∇Ψ -1 -Ψ * -1 ∇Ψ +1 ) |Ψ +1 | 2 + |Ψ -1 | 2 (75)
The distributions of the polarization currents calculated within GP approach are shown in Fig. 10. A striking nonuniform structure appears due to the presence of the spin-orbit coupling terms. Close to the source spot, there is a strong circularly polarized current that rotates around the source. This can be attributed to the rotating circular polarization degree already observed in Fig. 7. Away from the source, the circular polarization current decays, which can be expected due to the decay of the spin density current observed in Fig. 9. Along the vertical axis (x = 0), a strong circular polarization current remains due to the particularly fast change of the circular polarization degree in this region. The linearly polarized current can be stronger further away from the source than at closer distances. This is attributed to an increasing linear polarization degree further from the source. In addition, one can recall that while spin density currents are generally weaker further from the source, there is some compensation due to the conversion between dark and bright excitons (as shown in Fig. 5b, the bright exciton fraction increases further from the source).

VIII. CONCLUSIONS

Bosonic spin transport is a young and promising area of solid-state physics. The theories of mesoscopic transport of charge carriers and quantum transport are among the most interesting chapters of modern physics. Substitution of fermions by bosons and of a scalar electric charge by a spin vector cannot be formally done in these theories. Basically, all mesoscopic and quantum transport effects need to be reconsidered if we speak about electrically neutral bosonic spin carriers like excitons or excitonpolaritons. This is why the area of "spin-optronics" essentially remains terra incognita. Experimentally, direct measurements of transport of indirect excitons and exciton-polaritons in time-resolved imaging experiments have become possible in recent years. In this work, we have demonstrated that exciton polarization currents are inseparably connected with electron and hole spin currents. The intensity and direction of exciton polarization currents and electron and hole spin currents is governed by an interplay of spin-orbit effects, Zeeman effects and exciton exchange effects. In the non-linear regime, the pattern of spin currents may also be affected by spindependent exciton-exciton interactions.

We have developed two complementary approaches to the description of exciton spin currents and textures. The DM formalism allows for description of the spin transport effects in both classical exciton gases and condensates of non-interacting excitons, while the GP equations describe propagation of exciton condensates. We predict non-trivial topologies of interacting exciton spin in condensates, and suggest tools of their control, such as external magnetic and electric fields, and source intensity. We have demonstrated, that ballistic propagation of excitons may result in a build up of polarization patterns, which may be observed in near-field photoluminescence spectra.

This work has been supported by the EU FP7 ITN INDEX, EU FP7 PodiTrodi, the E.P.S.R.C., the EU FP7 Marie Curie EPOQUES, ANR-2011-NANO-004-06, and DOE. J. R. L. acknowledges Chateaubriand Fellowship from the Embassy of France in the United States. The authors are deeply grateful to T. Ostatnicky, Y.G. Rubo, M.M. Glazov, I.A. Shelykh, and A. Bramati for many useful discussions on the peculiarities of bosonic spin transport.

  ) Rewriting Eq. (3) and retaining only z-component of the magnetic field, which corresponds to the Faraday geometry, one can obtain: µ B B β e (k e,x + ik e,y ) β e (k e,x -ik e,y ) µ B B β e k e e iφ β e k e e -iφ 1 2 g e µ B B ] ,

FIG. 1 :

 1 FIG. 1: (color online) (a) Dispersion of the excitonic states calculated using the set of parameters from Table I (green lines) and reduced set of parameters with E b -E d = 0, δ b = 0, β h = 0 (red symbols). Black arrows indicate the values of the wavevectors used in (b)-(d) kex = 15.3µm -1 and (e) kex = 5µm -1 . Double-ended arrows indicate the energies of oscillations between the eigenstates which appear in spatial polarization patterns. (b) Linear polarization degree along xaxis calculated with the simplified set of parameters. (c) Same as (b) but for the full set of parameters from the Table I. (d) Circular polarization degree with parameters from TableI. (e) Same but at kex = 5µm -1 . The source area was taken circular with a radius of 4µm.

FIG. 2 :

 2 FIG.2:(color online) Spatial distribution of electron and hole spin components in the plane of the quantum well structure, calculated at zero magnetic field and at B = 5 T at three different orientations. The parameters are given in TableI, at kex = 15.3µm -1 .The source area was taken circular with a radius of 4µm.

FIG. 3 :

 3 FIG. 3: (color online) Spatial distribution of exciton linear polarization ((a)-(d), length of the bars maps the intensity in arbitrary units) and color maps of the exciton circular polarization degree ((e)-(f)), same color code as in Fig. 1). Magnetic field B = 0 in (a), (e), and B = 5T along z in (b),(f), x in (c),(g), y in (d),(h). Parameters are given in TableI, kex = 15.3µm -1 . The source area was taken circular with a radius of 4µm.

FIG. 4 :

 4 FIG. 4: (color online) Spatial distribution of electron (a) and hole (b) in-plane spin component calculated at B=0 in the absence of the Dresselhaus field βe = 0; β h = 0, but including Rashba field αe = 2.7µeVµm;. Other parameters are given in Table I, kex = 15.3µm -1 . The source area was taken circular with a radius of 4µm. Corresponding patterns of linear (c) and circular (d) exciton polarization degree are shown with the same color code as in Fig. 1.

FIG. 5 :

 5 FIG. 5: (color online) Spatial distribution of the total exciton density (a) and exciton brightness degree (b) [see Eq. (60)] in the absence of a magnetic field. The parameters were the same as those used for the DM calculations, given in TableI, with: B = 0; W = 0.2α; mex was taken as 0.21 of the free electron mass. The source area was taken circular with a radius of 1µm. The images are presented with spatial averaging over 1.5µm. The scattering parameter α and intensity at the source center were chosen such that the interaction energy, α( |Ψ+1| 2 + |Ψ-1| 2 + |Ψ+2| 2 + |Ψ-2| 2 )= 1µeV (being comparable to the other energy scales in the system, we are in a nonlinear regime). The absorbing boundary condition used in calculations appears outside of the plotted range, at a radius of 15µm from the source center.

FIG. 6 :

 6 FIG. 6: (color online) Dispersion relation obtained from the Gross-Pitaevskii equations, corresponding to Figs. 5 and 7b,d and f. The dispersion is obtained by Fourier transform of the wavefunctions in space and time, from which the grayscale map of intensity is obtained. The curves show the bare dispersion, obtained from diagonalization of Ĥtot ex + 2 k 2 2m .

FIG. 7 :

 7 FIG.7:(color online) Spatial distribution of the polarization state of excitons in the absence of a magnetic field: horizontallinear polarization degree (a, b), diagonal polarization degree (c, d) and circular polarization degree (e, f) for bright excitons, which corresponds to the near field emission pattern of light. The left-hand plots show results in the absence of nonlinear interactions (α = 0; W = 0), while the right-hand plots show the case of a moderate nonlinearity. The parameters were the same as in Fig.5

FIG. 8 :

 8 FIG. 8: (color online) Spatial distribution of the exciton spin current calculated in arbitrary units using the parameters summarized in Table I at kex = 15.3µm -1 at B = 0 (a) and B = 5T (b)-(d) along the z, x and y-axis, respectively. The source area was taken circular with a radius of 4µm.

FIG. 9 :

 9 FIG.9:(color online) Spatial structure of the spin density current j+1 (a) and the total magnetization current j tot M (b). The arrows show the directional dependence of the vector fields in space, while the colour code illustrates the intensity. The parameters were the same as in Fig.5. Arbitrary units are used for both the spin density and total magnetization current.

2 ∫

 2 j c (k ex ) = k ex m ex ρ c = k ex m ex ρ 11 -ρ 22 ρ 11 + ρ 22 (69) j l (k ex ) = k ex m ex ρ l = k ex m ex ρ 12 + ρ 21 ρ 11 + ρ 22 (70) j d (k ex ) = k ex m ex ρ d = k ex m ex ρ 12 -ρ 21 ρ 11 + ρ 22(71)The total polarization currents can be obtained integrating the expressions (69-71) over reciprocal space: j c,l,d tot = -dk ex j c,l,d (k ex )(72)

FIG. 10 :

 10 FIG.10:(color online) Spatial structure of the polarization currents jc (a) and j l (b). The arrows show the directional dependence of the vector fields in space, while the colour code illustrates the intensity. The parameters were the same as in Fig.5. Arbitrary units are used for the polarization currents.

TABLE I :

 I Parameters for numerical calculations.

	Electron mass	me	0.07 m0
	Heavy hole mass	m h	0.16 m0
	Electron Dresselhaus coupling	βe	2.7 µeVµm
	Heavy-hole Dresselhaus coupling	β h	0.92 µeVµm
	Bright exciton XY splitting	δ b	0.5 µeV
	Dark exciton XY splitting	δ d	-13 µeV
	Bright-dark exciton splitting E

b -E d 5 µeV