Hélène Angellier

Luc Choisnard

Sonia Molina-Boisseau

Patrick Ozil

Alain Dufresne
email: alain.dufresne@efpg.inpg.fr

Optimization of the preparation of aqueous suspensions of waxy maize starch nanocrystals using a response surface methodology

Keywords: Design of experiments (DOE), waxy maize starch, acid hydrolysis, nanocrystals

Response surface methodology was used to investigate the effect of five selected factors on the selective H 2 SO 4 hydrolysis of waxy maize starch granules. These predictors were temperature, acid concentration, starch concentration, hydrolysis duration and stirring speed. The goal of this study was to optimize the preparation of aqueous suspensions of starch nanocrystals, i.e. to determine the operative conditions leading to the smallest size of insoluble hydrolyzed residue within the shortest time and with the highest yield. Therefore empirical models were elaborated for the hydrolysis yield and the size of the insoluble residues using a Central Composite Face Design (CCFD) involving 31 trials. They allowed us to show that it was possible to obtain starch nanocrystals after only 5 days of H 2 SO 4 hydrolysis with a yield of 15 wt% and having the same shape as those obtained from the classical procedure after 40 days of HCl treatment, with a yield of 0.5 wt%.

Introduction

There is currently a considerable interest in processing polymeric composite materials filled with nanosized rigid particles (essentially inorganic). This class of material attracting both scientific and industrial communities is called "nanocomposites". Because of the nanometric size effect, these composites have some unique outstanding properties with respect to their conventional microcomposite counterparts. Nowadays the application of nanoparticles and then the development of new nanocomposite materials are restricted by both their limited availability and their strong tendency to aggregate preventing their homogeneous dispersion within a continuous matrix, which is the key step required for high mechanical performances.

Starch is a natural polymer available in large amounts from several renewable plant sources and it is produced in abundance beyond available markets. Starch is the cheapest biopolymer and is totally biodegradable. These two main reasons lead to the growing interest in the non-food usage of starch-based products for applications in which synthetic polymers have traditionally been the materials of choice. It is well known that native starch granules contain more or less concentric "growth rings" that are readily visible by optical or electron microscopy. 1 Acid treatment is needed to reveal the concentric lamellar structure of starch granules. It has been shown that these lamellae, around 5000 Å thick, have subspacing of a few hundred Å. [START_REF] Frey-Wyssling | Ultrastructural Plant Cytology[END_REF] The purpose of this treatment using hydrochloric acid is to dissolve away regions of low lateral order so that the water-insoluble, highly crystalline residue may be converted into a stable suspension by subsequent vigorous mechanical shearing action. [START_REF] Battista | Microcrystal Polymer Science[END_REF] In previous works [START_REF] Dufresne | [END_REF]5 such starch nanocrystals obtained from potato starch granules were used as a reinforcing phase in a polymeric matrix and displayed substantially improved mechanical properties. The insoluble hydrolyzed residue obtained from waxy maize was found to be composed of crystalline nanoplatelets around 5-7 nm thick with a length of 20-40 nm and a width of 15-30 nm. 6 However, the main drawbacks for the more extensive use of such nanoplatelets in nanocomposite applications are the duration (40 days of treatment) and the yield (0.5 wt%) of the HCl hydrolysis step. [START_REF] Battista | Microcrystal Polymer Science[END_REF] Starch disruption by acid hydrolysis depends on many factors such as the botanic origin, [7][START_REF] Duprat | In Les polymères végétaux[END_REF][START_REF] Jane | [END_REF][10][11][12][13] namely the crystalline type, the relative proportion of amylose and amylopectin, and the granules morphology. It also depends on the conditions of acid hydrolysis, namely the acid type, 14,15 acid concentration, 7 temperature 7 and hydrolysis duration. 7 No deep interest has been brought to the influence of stirring and starch concentration. Suspensions were usually stirred manually every day and a large panel of starch concentration has been used, for instance 1.5 wt%, 10 1.67 wt%, 16,17 2 wt%, 13 2.5 wt%, 12 5 wt% [START_REF] Battista | Microcrystal Polymer Science[END_REF][START_REF] Dufresne | [END_REF]6,[START_REF] Jane | [END_REF] and 33 wt%. 18 If HCl hydrolysis has been largely studied, 7 equivalent works about sulfuric hydrolysis remain limited.

It is now accepted that the partial crystallinity of native starch granules is due to a clustered organization of amylopectin side chains 16 and that an increasing amylose content results in a decrease of the susceptibility of starch granule to acidic degradation. 19 Therefore amylopectinrich starch, i.e. waxy maize, was chosen for the present study. Furthermore, a previous work 20 has shown that H 2 SO 4 acid hydrolysis resulted in more stable suspensions than hydrochloric ones, which is an important characteristic for composite materials processing. Thereby, we chose to focus this work on optimizing the H 2 SO 4 hydrolysis of waxy maize starch granules taking into account the following parameters: hydrolysis duration, temperature, acid concentration, starch concentration and stirring speed.

Response surface methodology was carried out to investigate the effect of these selected factors. This classical method is largely used and well adapted to process optimization in the macromolecular science area. [21][22][23] Laser granulometry, yield calculation and transmission electron microscopy were used to characterize the insoluble hydrolysis residues.

Experimental

Acid hydrolysis. A given weight of native waxy maize starch granules (Waxylis, Roquette S.A.) was mixed with 250 mL of H 2 SO 4 solution at a known concentration in a 500 mL erlenmeyer. The suspensions were then put on a plat-form in a thermo-stated atmosphere and continuously stirred at a selected speed with an orbital shaking action. After various durations of hydrolysis, the suspensions were washed by successive centrifugations in distilled water until neutrality. They were stored at 4°C with several drops of chloroform.

The hydrolysis yield (wt%) was calculated as the ratio between the weight of freeze-dried hydrolyzed particles and the initial weight of native granules for an aliquot of 50 mL taken in the 250 mL of hydrolyzed suspensions. It was verified that these aliquots were representative of the entire volume of 250 mL.

Characterization. Transmission electron microscopy (TEM) observations were performed using a Philips CM200 microscope with a 80 kV accelerating voltage.

Laser granulometry measurements were carried out with a Malvern Mastersizer. The suspensions were characterized from the median particle size d 50 , which divides the population into two equal halves.

Design of experiments. The five following parameters were varied:

-hydrolysis temperature u 1 , -acid concentration u 2 , -initial starch concentration u 3 , -hydrolysis duration u 4 , -speed of stirring, u 5 .

Parameters u 3 and u 5 do not influence hydrolysis kinetic because of its catalytic nature.

These two parameters have been however taken into account in the present work because they

were varying in a large range in the literature.

The settings of the different parameters were determined for all the predictors from our preliminary experiments 20 and keeping in mind the two following intuitive rules. On the one hand, the optimal conditions have a lot of chance to be outside of the region of interest if selecting a too narrow variation range. On the other hand, the predictive power of the model risks to be poor if choosing a too large range. Table 1 shows the minimal (u i min), the midrange (u i mid) and the maximal (u i max) values used for each parameter, which respectively correspond to -1, 0 and +1 levels in terms of orthogonal variable x i defined as:

min) max ( min) .( 2
i i i i i u u u u x    Table 1.
Here the two responses under study were:

-The hydrolysis yield as previously defined, y yield (wt%), -the median size y size (m) of hydrolyzed particles treated by Ultra Turrax (13,000 rpm, concentration of starch of 1 wt%, treated volume of 60 mL) during 5 min measured by laser granulometry.

Since both response non linearities and interactions between factors were expected, the response model was postulated to be a quadratic one for each response  k . This model may be expressed in terms of orthogonal variables as:

           5 1 , 5 1 2 5 1 0 j i j i i i ii j i ij i i i k x b x x b x b b



This postulated model takes into account the linear effects (b i x i ) and the quadratic ones (b ii x i 2 ) as well as first-order interactions (b ij x i x j ) and its knowledge requires to estimate the

coefficients b i, b ii, b ij .
Experiments were conducted adopting a Central Composite Face Design (CCFD). This type of design was suitable for our objective, which was the optimization of a potential complex process. The CCFD (Table 2) involves a fractional factorial design 2 5-1 (trials 1 to 16), the face centers (trials 17 to 26 -all coordinates equal to zero except one equal to +1 or -1)) and five replications at the central point (trials 27 to 31 -all coordinates equal to zero). The run order of trials was randomized in order to prevent systematic errors.

The determination of the optimal conditions proposed by the MODDE software is based on the maximization of a desirability coefficient D defined as a weighted average of the individual response desirabilities d i :

24 2 1 1 1 1                 i i i i M i i M i i L T T y w M d M D
M being the number of responses and the experimenter having to provide the following data for each response y i : the weight w i (here, w i = 1 for all the responses), the desired response target T i and its worst acceptable value L i .

Results and discussion

Models and analysis. The responses measured for each trial are reported in Table 2.

Table 2.

Plotting the N-probability versus deleted studentized residuals has shown that there were no deviating experiments also called outliers, considering that the action limit was  4 standard deviations (not shown in this paper). Thereby, no experiment was performed again nor excluded from the analysis. Moreover the analysis of raw data through the Box-Cox transformation showed that no response transformation was useful to improve the models, considering a 0.95 confidence interval.

The responses were fitted owing to a multi-linear regression method (MLR) as quadratic models expressed in terms of orthogonal variables x i as said above and refined using a backward step by step technique based on a 0,95 confidence level. [START_REF] Goupy | In Plans d'expériences pour surfaces de réponse[END_REF][START_REF] Darlington | In Regression and linear models[END_REF] The final refined models ( yield and  size ) and their corresponding statistics (coefficient of regression R 2 , coefficient of determination R 2 adj and coefficient of prediction Q 2 ) are as follow: For both models, R 2 , R 2 adj and Q 2 statistics have quite high values close to unity with a difference R 2 and Q 2 lower than 0.2, which is an indication of suitable models. Here the adjusted R 2 statistics are greater than 0.85, so proving a good descriptive power of the models. This power can be illustrated besides by the plot of experimental responses versus the predicted ones (Figure 1) which should be ideally the diagonal line and is here characterized by an acceptable dispersion. Moreover the predictive power, as evidenced by Q 2 values greater than 0.5, is excellent for  yield and not so good for  size (Q 2 =0.751) but still acceptable. 21 Figure 1.

x x x x x x x x yield          R 2 = 0.940, R 2 adj = 0.922, Q 2 = 0.884
Moreover a classical analysis of variance (ANOVA) using F-tests allows to analyze the total response variation by identifying the parts corresponding to the sources of variation (regression model, pure experimental error) and to analyze the residuals in order to point out the possible lack of fit of the postulated model when replicates are available. Here we could conclude that there is no lack of fit for the model concerning the hydrolysis yield while some lack of fit exists for the median size. That is however not critical because of a Q 2 statistic greater than 0.5

(Q 2 =0.751).
The stirring speed (x 5 ) and overall the starch concentration (x 3 ) do not have any significant influence on the yield and the median size (Figure 2). Effectively, the probabilities of significance (p) of estimated coefficients b 3 and b 5 were superior to the critical limit of 0.05.

They were respectively p = 0.1598 and p = 0.0653 for  yield , and p = 0.7230 and p = 0.0356 for  size . However, these terms were kept in models according to the hierarchy principle. It is quite surprising that starch concentration did not have an effect on the yield of hydrolysis and the median size of insoluble residues. We could have thought that an increase of starch concentration would have slowed down the diffusion of acid in the suspension. Starch concentration did not have an effect on the yield of hydrolysis and the median size of insoluble residues. We could have thought that a decrease of starch concentration would have favored the action of the catalyst. This result proves that, in the chosen ranges, acid is largely in excess in the reactor.

Only one interaction between starch concentration and stirring has been detected. The response surfaces and the corresponding contour plots were drawn for both responses.

They allow illustrating the impact of two selected parameters by keeping the other three parameters at constant values (Figures 3 to 6).

Figures 3 to 6.

In order to study the reliability of the models and validate them, several experimental conditions were tested. We succeeded in producing controlled yields with a respective confidence level of 0.95. On the other hand, the experimental value of the median size was outside the predicted interval (difference between the measured sizes and the interval of prediction and the predicted values of 6.1% and 24%, respectively) so confirming a light lack of fit for the model.

Nevertheless the model will be used for optimizing the median size because of the suitable value of the Q 2 statistic greater than 0.5.

Optimization. The multivariable models obtained from the statically design of experiments were used for predicting the optimal conditions of H 2 SO 4 hydrolysis that should allow obtaining aqueous suspensions of waxy maize starch nanocrystals in the shortest time, with the smallest median size particle and the highest yield.

We assumed that the highest yield for the preparation of starch nanocrystal suspensions that could be reached was about 30-40% which is the relative crystallinity of native waxy maize starch granules. [START_REF] Zobel | [END_REF] Preliminary studies 20 have shown that the morphology of nanocrystals began to be observed after 7 days of acid hydrolysis with a corresponding yield of 26.9 wt% and that 40 days of HCl hydrolysis were needed to obtain a suspension of insoluble residues all having the shape of nanoplatelets, with a corresponding yield of 0.5 wt%. These results showed that selective acid hydrolysis takes time. Thereby, a reasonable goal was to obtain such suspensions in less than 7 days, with a yield of 20 wt%.

For optimization, the yield and the median size were considered as targets, with a constrained time value and all other factors varying freely in the considered experimental domain (Table 3). First, a target of 20 wt% for the yield ( yield ) and of 5 m for the median size ( size )

have been chosen, with a duration of hydrolysis of 7 days. These criteria corresponded to "optimization 1". Secondly, duration of hydrolysis of 5 days has been chosen, with the same target for the yield but a median size of 6 m, considering that we would not have been able to obtain the same median size than in 7 days. These criteria corresponded to "optimization 2".

Table 3.

The sets of conditions proposed by the MODDE software are given in Table 4. The given log D value is negative or near zero indicating that we are doing rather well. The predictions calculated with a confidence interval of 0.95 and the measured responses corresponding to these sets of conditions are given Table 5.

Table 5.

Measured yields were inside the predictive intervals but not the median sizes, which were smaller than the lower limits. Whereas it confirmed the lack of prediction of the model  size , this result is pleasantly surprising and interesting for our goal which was to prepare small residues, among others. Furthermore, contrary to all expectations, a smaller median size was obtained with the same yield after 5 days than after 7 days of acid treatment. Observations by transmission electron microscopy (TEM) were performed to verify the shape of the insoluble residues.

Characterization of optimized suspension of nanocrystals.

Observations by TEM clearly shown that the insoluble residue obtained after 5 days of optimized H 2 SO 4 hydrolysis (Table 5) have the shape of parallelepiped nanoplatelets (Figure 7a).

Nanoplatelets were generally observed in aggregates of 1 to 5 m (Figure 7a) or at best in barrets of several platelets (Figure 7b). Even if the parallelepipedic shape was the general shape that we observed, a lot of varying organizations (Figures 7b-d) were distinguished. Few stacks of nanoplatelets oriented edge-on were observed in a very little proportion (not shown here), what let assume that the platelets were well separated.

These observations lead to conclude that the optimization process was successful. We have shown that it was possible to obtain starch nanocrystals after 5 days of H 2 SO 4 hydrolysis with a yield of 15 wt% having the same shape than those obtained after 40 days of HCl treatment, with a yield of 0.5 wt%.

Conclusion

The statistical experimental design and the multi-linear regression analysis used in this study have proven to be very useful for establishing predictive models for both the yield of H 2 SO 4 hydrolysis of waxy maize starch granules and the median size of insoluble residues after acid treatment. We achieved to produce aqueous suspensions of starch nanocrystals after 5 days of

Figure 2 .

 2 Figure 2.

Figure 1 :

 1 Figure 1: Plot of measured responses versus predicted responses for (a)  yield and (b)  size .

Figure 2 :

 2 Figure 2: The main effects for (a) the yield ( yield ) and (b) the median size ( size ) for a confidence

Figure 3 :

 3 Figure 3: Response-surface plot (upper) and its contour plot of yield ( yield ): temperature versus

Figure 4 :

 4 Figure 4: Response-surface plot (upper) and its contour plot of median size ( size ): temperature

Figure 5 :

 5 Figure 5: Response-surface plot (upper) and its contour plot of yield ( yield ): acid concentration

Figure 6 :

 6 Figure 6: Response-surface plot (upper) and its contour plot of median size ( size ): acid

Figure 7 :

 7 Figure 7: TEM of negatively stained nanocrystals obtained after 3.16M H 2 SO 4 hydrolysis of

Table 4 .

 4 Achieving an accuracy of 4 decimals for experimental conditions settings was of course impossible. Run 1 was performed using 35.35 g of starch (14.14 wt%) mixed in 250 mL of 2.87 M H 2 SO 4 solution (161.8 mL of H 2 SO 4 95% and 838.2 mL of distilled water), at 40°C, 100 rpm and during 7 days. For run 2, 36.725 g of starch (14.69 wt%) were mixed in 250 mL of

	3.16 M H 2 SO 4 solution (178.1 mL of H 2 SO 4 95% and 821.9 mL of distilled water), at 40°C, 100
	rpm and during 5 days.

Table 1 :

 1 Setting levels of parameters u i

	3.16 M H 2 SO 4 hydrolysis at 40°C, 100 rpm and with a starch concentration of 14.69 wt% with a
	yield of 15.7 wt% and the same shape than those obtained after 40 days of HCl hydrolysis. It
	should allow considering preparation of starch nanocrystals for nanocomposite applications.

Table 2 :

 2 Results obtained for the complete set of 31 experimental points

	Trial	X 1	X 2	X 3	X 4	X 5	Y yield	Y size
							(%)	(m)
	1	-1	-1	-1	-1	+1	76.3	13.40
	2	+1	-1	-1	-1	-1	68.1	13.00
	3	-1	+1	-1	-1	-1	47.6	12.44
	4	+1	+1	-1	-1	+1	26.3	9.23
	5	-1	-1	+1	-1	-1	70.7	13.18
	6	+1	-1	+1	-1	+1	54.8	13.22
	7	-1	+1	+1	-1	+1	57.8	12.66
	8	+1	+1	+1	-1	-1	35.9	10.93
	9	-1	-1	-1	+1	-1	43.9	11.91
	10	+1	-1	-1	+1	+1	20.3	6.29
	11	-1	+1	-1	+1	+1	5.4	5.84
	12	+1	+1	-1	+1	-1	2.8	5.21
	13	-1	-1	+1	+1	+1	44.8	8.00
	14	+1	-1	+1	+1	-1	29.3	7.61
	15	-1	+1	+1	+1	-1	16.7	5.90
	16	+1	+1	+1	+1	+1	2.1	5.08
	17	-1	0	0	0	0	42.3	11.81
	18	+1	0	0	0	0	26.5	6.74
	19	0	-1	0	0	0	30.4	11.56
	20	0	+1	0	0	0	21.3	6.45
	21	0	0	-1	0	0	24.1	9.43
	22	0	0	+1	0	0	34.9	8.33
	23	0	0	0	-1	0	56.4	12.85
	24	0	0	0	+1	0	20.6	6.75
	25	0	0	0	0	-1	36.0	9.65
	26	0	0	0	0	+1	20.3	4.64
	27	0	0	0	0	0	37.6	7.56
	28	0	0	0	0	0	31.8	8.28
	29	0	0	0	0	0	28.3	8.43
	30	0	0	0	0	0	29.7	7.52
	31	0	0	0	0	0	28.9	7.79

Table 3 :

 3 Criteria for optimizations 1 and 2.

		Criteria	Target	Min	Max
	Optimization 1
	Yield ( yield )	Target	20	15	30
	Median size ( size )	Target	5	4.5	5.8
	u 1	Free		35	40
	u 2	Free		2.2	3.4
	u 3	Free		5	15
	u 4	Constant	7
	u 5	Free		0	100
	Optimization 2
	Yield ( yield )	Target	20	15	30
	Median size ( size )	Target	6	5.5	6.5
	u 1	Free		35	40
	u 2	Free		2.2	3.4
	u 3	Free		5	15
	u 4	Constant	5
	u 5	Free		0	100

Table 4 :

 4 Sets of conditions proposed for optimizations 1 and 2.

	Optimization	U 1	U 2	U 3	U 4	U 5	Log(D)
	1	39.9993	2.8728	14.1457	7	99.9942	0.0205
	2	39.9914	3.1636	14.6946	5	99.3155	-0.3706

Table 5 :

 5 Predicted and measured responses for sets of conditions 1 and 2.

		Predicted	Low limit	High limit	Measured
	yield
	Optimization 1	15.4881	10.0454	20.9315	16.63
	Optimization 2	16.6481	10.6395	22.6568	15.7
	size
	Optimization 1	5.9059	4.6442	7.1677	4.64
	Optimization 2	6.3172	4.9243	7.7101	4.40

Acknowledgements

The authors are grateful to Roquette (Lestrem, France) for the supply of waxy maize starch. Special thanks are given to Dr J.L. Putaux (CERMAV-CNRS) for his support in TEM observations and the discussions about acid hydrolysis.


S Laurent

P Parnaudeau

F Chevy

I Danaila

Nonlinear dynamics of coupled superfluids

Keywords:

Following recent experiments on ultracold dual superflows, we model in this work the dynamics of two harmonically trapped counterflowing superfluids. Using complementary analytical and numerical approaches, we study the shedding of elementary excitations triggered by the relative motion of the two species. We exhibit two different excitation mechanisms leading to distinct threshold velocities for the onset of dissipation: in addition to the parametric pair production present in homogeneous, galilean-invariant systems, we show that non-uniform motion and density inhomogeneities allow for a Landau-like decay mechanism where single excitations are produced.

I. INTRODUCTION

In recent years, progress in manipulation of ultracold boson/boson [START_REF] Myatt | Production of two overlapping Bose-Einstein condensates by sympathetic cooling[END_REF] or boson/fermion mixtures [START_REF] Ferrier-Barbut | A mixture of Bose and Fermi superfluids[END_REF][START_REF] Ikemachi | All-optical production of dual BoseEinstein condensates of paired fermions and bosons with 6 Li and 7 Li[END_REF][START_REF] Roy | Two-element mixture of Bose and Fermi superfluids[END_REF][START_REF] Yao | Observation of coupled vortex lattices in a mass-imbalance Bose and Fermi superfluid mixture[END_REF] has led to the observation of dual-superfluid gases. In these systems, the interplay between the two species leads to a rich physics, like the creation of spin domains [START_REF] Stenger | Spin domains in groundstate Bose-Einstein condensates[END_REF], dark-bright solitons [START_REF] Th | Dark-bright solitons in inhomogeneous Bose-Einstein condensates[END_REF][START_REF] Tylutki | Dark-bright solitons in a superfluid Bose-Fermi mixture[END_REF] or "liquid" droplets [START_REF] Petrov | Quantum mechanical stabilization of a collapsing Bose-Bose mixture[END_REF][START_REF] Cabrera | Quantum liquid droplets in a mixture of Bose-Einstein condensates[END_REF]. Dual-superfluid mixtures play a fundamental role going beyond the field of atomic physics and quantum gases. In neutron stars, kaon condensates are believed to coexist with neutronproton superfluids [START_REF] Paul | Kaon condensation in neutron star matter with hyperons[END_REF] while among the liquid Helium community the quest for a dual 3 He/ 4 He dual superfluid mixture has been considered as one of the holy grails of low-temperature physics ever since the observation of superfluidity in 3 He [START_REF] Rysti | Effective 3 He interactions in dilute 3 He-4 He mixtures[END_REF].

An intriguing possibility is the existence of superfluid counterflows where the two fluids move with different velocities. This multiple-fluid hydrodynamics was considered first theoretically for Helium mixtures [START_REF] Volovik | Theory of solutions of a superfluid Fermi liquid in a superfluid Bose liquid[END_REF][START_REF] Af Andreev | Three-velocity hydrodynamics of superfluid solutions[END_REF][START_REF] Nespolo | Andreev-Bashkin effect in superfluid cold gases mixtures[END_REF] and was observed experimentally in ultracold boson/boson [START_REF] Maddaloni | Collective oscillations of two colliding Bose-Einstein condensates[END_REF][START_REF] Smerzi | Macroscopic oscillations between two weakly coupled Bose-Einstein condensates[END_REF][START_REF] Hamner | Generation of dark-bright soliton trains in superfluid-superfluid counterflow[END_REF] and fermion/boson mixtures [START_REF] Yao | Observation of coupled vortex lattices in a mass-imbalance Bose and Fermi superfluid mixture[END_REF][START_REF] Delehaye | Critical velocity and dissipation of an ultracold Bose-Fermi counterflow[END_REF][START_REF] Wu | Coupled dipole oscillations of a mass-imbalanced Bose-Fermi superfluid mixture[END_REF]. These experiments raised the question of the critical relative velocity above which dissipation arises and suggested that superfluid counterflows were unstable against generation of pairs of excitations, a mechanism providing a Galileaninvariant generalization of Landau's celebrated scenario [START_REF] Castin | The landau critical velocity for a particle in a Fermi superfluid[END_REF][START_REF] Zheng | Quasiparticle lifetime in a mixture of Bose and Fermi superfluids[END_REF][START_REF] Abad | Counter-flow instability of a quantum mixture of two superfluids[END_REF].

In the experiments reported in [START_REF] Delehaye | Critical velocity and dissipation of an ultracold Bose-Fermi counterflow[END_REF], the counterflow was obtained by releasing a dual Bose/Fermi superfluid of 6 Li and 7 Li in a harmonic trap. Due to the mass difference between the two species, the motion of the two superfluid get out of phase after a few oscillation periods. This creates a relative motion between the two components and dissipation was only observed above a * Electronic address: Corresponding author: chevy@lkb.ens.fr critical velocity confirming the superfluid nature of the counterflow. However, the harmonic trapping confining the atoms breaks translational and Galilean invariance assumptions underpinning simple dissipation mechanisms. In the present work, we explore theoretically the stability of the counterflow generated by two harmonically trapped superfluids. We first consider the case of two weakly interacting Bose-Einstein condensates described using coupled Gross-Pitaevskii equations. These equations are solved numerically and the density profiles are characterized using a Principal Component Analysis (PCA) scheme [START_REF] Jolliffe | Principal Component Analysis[END_REF] that allows for a model-free identification of the modes triggered by the relative motion of the two clouds. In the second part of the paper we investigate the mode-coupling mechanism using an analytical hydrodynamic approach. This approximation is valid only for long wavelength excitations but it is complementary to the numerical simulation, since it is applicable to any kind of superfluid (bosons or fermions, weakly or strongly coupled, etc.). We show that the relative motion excites eigenmodes of the two clouds and we identify two different excitation mechanisms. First, a parametric process akin to [START_REF] Castin | The landau critical velocity for a particle in a Fermi superfluid[END_REF] leads to the formation of pairs of excitations in both superfluids. Second, the motion of the smaller cloud excites linearly the density profile of the larger one by the creation of single phonon-like modes, following a process similar to Landau's traditional scenario. We show in particular that this latter process is a direct consequence of the translation-symmetry breaking induced by the presence of the trap.

II. NUMERICAL SIMULATIONS

A. Coupled Gross-Pitaevskii equations

We first simulate numerically the counterflow of two harmonically trapped superfluids using a set of coupled Gross-Pitaevskii equations describing zero-temperature Bose-Einstein condensates (BECs) with mean-field interactions:

i ∂ψ 1 ∂t = - 2 2m 1 ∇ 2 + U (r) + U mf,1 (r) ψ 1 , (1)
i ∂ψ 2 ∂t = - 2 2m 2 ∇ 2 + U (r) + U mf,2 (r) ψ 2 , (2)
with N 1 and N 2 the number of particles of each species. The wavefunctions ψ 1 and ψ 2 are thus normalized to unity. To reproduce the setting used in recent experiments [START_REF] Ferrier-Barbut | A mixture of Bose and Fermi superfluids[END_REF], the two species are trapped using the same cigar-shaped harmonic external potential

U (r) = 1 2 m 2 ω 2 ⊥ (x 2 + y 2 ) + ω 2 z z 2 , ( 3
)
where ω ⊥ and ω z are the radial and axial trapping frequencies, respectively. Finally U mf,1,2 describe the meanfield interaction between the atoms and we have

U mf,i = j=1,2 N j g ij |ψ j | 2
, where the g ij are the s-wave coupling constants satisfying the symmetry condition g ij = g ji .

The parameters are chosen to reproduce two important features of the experiment described in [START_REF] Ferrier-Barbut | A mixture of Bose and Fermi superfluids[END_REF]. First, we consider an atomic mass ratio of m 1 /m 2 = 7/6, such that the bare trapping frequencies of the two clouds are related via ω z,1 = 6/7 ω z . Thus, by displacing the clouds by the same distance b along the z-axis, a relative motion between the clouds progressively sets in and a periodic counterflow between the two superfluids is created. In addition, the small atomic mass difference allows for a strong coherent energy exchange, as observed in [START_REF] Ferrier-Barbut | A mixture of Bose and Fermi superfluids[END_REF]. Second, we chose very different atoms numbers and interspecies coupling constants to create a large shape asymmetry between the two clouds. The i = 1 cloud is weakly interacting and small while the i = 2 cloud is strongly repulsive and is much broader (see density profiles in the top panel of Fig. 1).

We start the numerical simulation by generating the static density profiles of the two Bose-Einstein condensates. For this purpose, we numerically solve the coupled stationary Gross-Pitaevskii equations corresponding to (1)-(2) using an imaginary-time propagation method with Fourier spectral space accuracy [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF][START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF]. Then we shift the two stationary profiles and let the system evolve in time. The real dynamics described by (1)-( 2) is computed using a second-order Strang splitting method [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. Both stationary and real-time dynamics computations are performed using the state-of-the-art computational code GPS (Gross-Pitaevskii-Simulator) offering various modern numerical methods to solve GP equations on high-performance parallel computers [START_REF] Parnaudeau | GPS: an efficient and spectrally accurate code for computing Gross-Pitaevskii equation[END_REF]. A typical grid used for this study contained 1024 × 128 × 128 computational points along the z and x, y axes, respectively. To capture the oscillations of the system for long times, up to 400, 000 time steps were needed. Special care was devoted to the accuracy of the time and space numerical scheme, in order to conserve the mass and energy of the system during this long-time integration. More details on the numerical simulations are given in Appendix A). Supplemental Material is provided with animations depicting the oscillations of the clouds for several run cases [29].

The influence of the coupling parameters was explored by simulating more than 30 cases, for values in the range 0 < g 12 /g 22 < 0.3 and 2 ≤ b/a ho ≤ 8, with a ho = /mω z . A typical numerical result is displayed in Fig. 1. The response of the center of mass of the two clouds is illustrated by plotting their relative position z rel = z 2 (t) -z 1 (t). For the case g 12 /g 22 = 0.031, b = 2, low amplitude oscillations of the clouds are observed and z rel beats with the frequency difference ω z,1 -ω z,2 . 1 combined with a large atom number ratio N2/N1 = 30 induce a strong amplitude modulation of the small condensate oscillations (i = 1, in blue), while the oscillations of the large condensate are mostly not affected (i = 2, in red).

� � �� �� �� �� �� -�� -� � � �� � (�π/ω � ) � ��� (�) (� �� ) 5
• • • • • • • • • • • • • • • • • • • • ���� ���� ���� ���� ���� ���� ���� ���� ���� g �� /g �� ω ��� /ω � FIG.
2: Dipole-frequency shift (ωz,1/ωz) of the oscillations of the small cloud (i = 1) for different interspecies coupling g12/g22: numerical simulation (blue dots), theoretical prediction (black line) based on the perturbation theory Eq. ( 5) and coupled-oscillator prediction (orange line).

� � �� �� �� �� �� -� -��� � ��� � � (�π/ω � ) �(�)(����) � � �� �� �� �� �� -� -��� � ��� � � (�π/ω � ) �(�) (����)
FIG. 3: Examples of the two nonlinear modes observed in simulations and analysed using the Principal Component Analysis (PCA). Space structure (upper panels) and time evolution of the amplitude of the modes (lower panels) extracted using the PCA. Blue lines are for the small cloud (i = 1) and red lines for the large cloud (i = 2). Linear forced mode (left, g12/g22 = 0.291, b/a ho = 1) and parametric modes (right, g12/g22 = 0.003, b/a ho = 5). Animations depicting the oscillations of the clouds for these regimes are provided in the Supplemental Material [29].

B. Dipole-mode frequency shift

Numerical simulations similar to that displayed in Fig. 1, allowed us to follow the evolution of the oscillation frequency of the small cloud with the interspecies coupling constant g 12 /g 22 . Figure 2 confirms the shift of the dipole mode, as predicted by the sum-rule approach that was previously used in [START_REF] Ferrier-Barbut | A mixture of Bose and Fermi superfluids[END_REF] to measure the equation of state of an attractive Fermi gas in the BEC-BCS crossover or more recently in [START_REF] Roy | Two-element mixture of Bose and Fermi superfluids[END_REF] to measure the Lithium-Cesium scattering length and in [START_REF] Wu | Coupled dipole oscillations of a mass-imbalanced Bose-Fermi superfluid mixture[END_REF] in the case of LiK mixtures. In our case, the sum-rule associated with the shift of the two trapping potentials predicts that the oscillation frequencies of the two clouds are the eigenvalues of the susceptibility matrix M defined by

M =   m 1 (1 -N1 N2 χ 12 ) √ m 1 m 2 N1 N2 χ 12 √ m 1 m 2 N1 N2 χ 12 m 2 (1 -χ 12 )   (4)
where χ 12 = ∂z1 ∂b2 is the displacement of the center of cloud i = 1 after a shift by a distance b 2 of the potential trapping of cloud i = 2. In the weak-coupling limit, this model predicts a shift δω of the smaller cloud frequency given by

δω 1,z ω 1,z - g 12 2 ∂ρ 2 ∂µ 2 0 , (5)
where ρ and µ are the density and the chemical potential of the gas, respectively.

We observe in Fig. (2) that the simulation agrees with this asymptotic regime for a weak coupling. For intermediate coupling, we observe a slight departure of the simulation results with respect to Eq. ( 5) but this discrepancy can be cured using the exact diagonalization of M (orange line in Fig. 2).

C. Mode coupling

To get further insight on the dynamics of the system, we analyze the density profiles of the two clouds by performing a Principal Component Analysis (PCA) [START_REF] Jolliffe | Principal Component Analysis[END_REF]. This method allows us to identify the modes involved in the dynamics of system without any a priori assumption on their spacial structure [START_REF] Dubessy | Imaging the collective excitations of an ultracold gas using statistical correlations[END_REF]. When applied to cloud images (atomic density integrated along the y axis), the PCA extracts a set of 4000 modes together with their associated eigenvalues and their temporal evolutions. Only a limited number of modes play a non-negligible role in the fluctuations (at least 99% of the eigenvalues are 10 6 times smaller than the largest one). For each run, we typically limited our analysis to the 10-30 most populated modes given by the PCA. Among these modes, we could identify the following types of modes: parametric, linear forced, filtered dipole modes and harmonic modes. The first two types of nonlinear modes are illustrated in Fig. 3 and additional modes are depicted in in the Supplemental Material [29].

We monitor the mode coupling by plotting the evolution of the weight c ik (t) of the k-th mode of the superfluid i unveiled by the PCA (see Fig. 3). As a first check, we plot the Fourier spectrum of c ik to extract the frequencies of each mode. In Fig. [START_REF] Roy | Two-element mixture of Bose and Fermi superfluids[END_REF] we compare the PCA-extracted oscillations with the analytic prediction of frequencies in elongated traps [START_REF] Stringari | Dynamics of Bose-Einstein condensed gases in highly deformed traps[END_REF]. For µ ω ⊥ the cloud is hydrodynamic in all three directions and the frequency of the k-th mode is:

ω k = 1 2 k(k + 3)ω z . ( 6
)
This assumption is satisfied by the larger cloud, as illustrated in Fig. 4 by the agreement between the numerical results and the predictions of Eq. [START_REF] Stenger | Spin domains in groundstate Bose-Einstein condensates[END_REF]. By contrast, 3D hydrodynamics breaks down for the smaller cloud where the low atom number yields a lower chemical potential.

In this regime, since µ ω ⊥ , the transverse degrees of freedom are frozen and the collective dynamics is onedimensional, leading to the following dispersion relation

ω k = k(k + 3) 2 ω z , (7)
which is in very good agreement with the results of numerical simulations (see the upper panel of Fig. 4) . Likewise, we compare in the lower panel of Fig. 4 the spatial structure of the modes unveiled by the PCA and show that they agree with the Geigenbauer polynomials known to describe the low-lying modes of a harmonically trapped Bose-Einstein condensate [START_REF] Jin | Hydrodynamic response of a trapped superfluid to a periodic perturbation[END_REF].

III. HYDRODYNAMIC APPROACH A. General formalism

To explain the excitation mechanism of the modes observed in numerical simulations, we analyze their dynamics using the hydrodynamic approximation. This approach is rather standard and is summarized in the following section. The starting point of this second analysis is the classical Hamiltonian

H = d 3 r N i=1 2 2m i ρ i (r) (∇Φ i ) 2 + ρ i U i (r) + e(ρ 1 , ..., ρ N ) , (8)
that describes the dynamics of an ensemble of N superfluids. In this expression, Φ i and ρ i are the phase and the density of superfluid i, while e is the energy density of the system. Taking the phase and density as dynamically conjugate variables, Hamilton's equations of motion yield

m i ∂ t v i + ∇v 2 i /2 = -∇ ∂e ∂ρ i + U i (9)
∂ t ρ i + ∇ (ρ i v i ) = 0, (10)
� � �� �� � � � � � ���� ������ �� ����� ���� ��������� (ω � ) -� � � -� � � �/� � � � (�) -� � � -� � � �/� � � � (�) -� � � -� � � �/� � � �� (�) -� � � -� � � �/� � � �� (�)
FIG. 4: Analysis of the eigenmodes. Upper panel: Comparison between PCA-extracted mode frequencies and the theoretical eigenmodes of an elongated Bose-Einstein condensate (dashed line). Lower panels: space structure of modes k =5, 6, 10, and 14, and their fit using Geigenbauer polynomials.

where v i = ∇Φ i /m i is the local superfluid velocity and ∂ ρi e is the chemical potential of species i. In the limit of weak interspecies-coupling, we can expand the chemical potential versus the densities of the other components:

∂e ∂ρ i µ i (ρ i ) + j =i g ij ρ j ,
where µ i (ρ i ) is species i alone's zero-temperature equation of state. Since by definition we have g ij = ∂ 2 ρj ρi e, the interspecies coupling constants obey the symmetry relation g ij = g ji .

Consider first a single species (g ij = 0). In the stationary state, the phase varies as Φ i,0 = µ 0 i t/ and the density profile ρ i,0 (r) is time independent and satisfies a Thomas-Fermi equation

µ i (ρ i,0 (r)) + U i (r) = µ 0 i , (11)
where the index 0 indicates equilibrium quantities. Low-lying excitations are obtained by considering v i and δρ i = ρ i -ρ i,0 as small parameters. Expanding the hydrodynamic equation yields after some straightforward algebra

∂ 2 t δµ i - 1 m i ∂µ i ∂ρ i 0 ∇ (ρ i,0 ∇δµ i ) = 0. ( 12
)
with δµ i = (∂µ i /∂ρ i ) 0 δρ i . The eigenmodes of the superfluid are found using the ansatz δµ i (r, t) = u i,k (r)e -iω i,k t , where u i,k and ω i,k are solutions of the eigenproblem

ω 2 i,k u i,k = L i [u i,k ], (13)
with

L i [u] = - 1 m i ∂µ i ∂ρ i 0 ∇ [ρ i,0 ∇u] . ( 14
)
L i is a positive operator and is symmetric for the inner product •|• i defined by

u|v i = d 3 r ∂ρ i ∂µ i 0 u(r) * v(r). (15)
The u i,k can then be chosen as an orthonormal basis for this inner product.

In the literature, general properties of Eq. ( 13) were derived for simple cases:

• u i =cte is a solution to Eq. ( 13) for the eigenfrequency ω i = 0. This solution corresponds to a shift of the global chemical potential µ 0 i in Eq. ( 11) and implies that all other eigenstates u ik are orthogonal to constant functions for the inner product [START_REF] Nespolo | Andreev-Bashkin effect in superfluid cold gases mixtures[END_REF].

• For a homogenous system, the solutions to Eq. [START_REF] Volovik | Theory of solutions of a superfluid Fermi liquid in a superfluid Bose liquid[END_REF] are plane waves and describe first-sound propagation. The eigenfunctions are characterized by a wave-vector k i and obey the dispersion delation ω i = c i k i , where the sound velocity is defined by m i c 2 i = ρ i,0 (∂µ i /∂ρ i ) 0 . • The plane-wave structure also applies to cylindrical traps [START_REF] Stringari | Dynamics of Bose-Einstein condensed gases in highly deformed traps[END_REF]. In this case, the waves propagate along the symmetry axis of the cloud at a velocity c 1D = √ ρi,0 (∂µ i /∂ ρi ) 0 /m i ], where ρi is the density integrated over the transverse degrees of freedom.

• In a harmonic trap, we recover Kohn's theorem [START_REF] Kohn | Cyclotron resonance and de haas-van alphen oscillations of an interacting electron gas[END_REF] since the oscillation of the center of mass along the principal direction x α of the trap corresponds to the function u ik = x α and are associated with the eigenvalue ω i,α .

• The previous result can be generalized to any mode for polytropic equation of state µ i ∝ ρ γ i . In this case, the eigenmodes of L i are Gegenbauer polynomials of the spatial coordinates.

Furthermore, for a harmonic (or flat) potential, the dynamics of the cloud obeys Kohn's theorem [START_REF] Kohn | Cyclotron resonance and de haas-van alphen oscillations of an interacting electron gas[END_REF] and the previous results can be extended to the case where the superfluid is initially oscillating. More precisely, the hydrodynamic equations are invariant under the transformation r → r = r -R i (t), v i → v i = v i -V i , where R i and V i are the position and the velocity of the center of mass of the cloud and obey Newton's second law m i Vi = -∇U i (R i ). As a consequence, the eigenmodes u i,k (r, t) = u i,k (r -R i (t)) can also be used to describe the low-lying excitations of an oscillating superfluid [37].

We now consider an ensemble of coupled superfluids. As before, in absence of coupling, the density in the laboratory frame is given by ρ i (r, t) = ρ i,0 (r -R i (t)). The coupling affects the density profiles which can be expanded over the eigenmodes u i,k as

ρ i (r, t) = ρ i,0 (r -R i (t)) + k c i,k (t) ∂ρ i ∂µ i 0 (r -R i (t))u k,i (r -R i (t)). ( 16
)
To simplify the notation, we hereafter denote with a prime physical quantities evaluated in the moving frame. In other words, for any function F i (r), we define F i (r, t)

as F i (r, t) = F i (r -R i (t)).
Inserting this expansion in the hydrodynamic equations, we obtain the following set of coupled differential equations for the coefficients

c i,k ci,k + ω 2 i,k c i,k + ω 2 i,k j =i g ij A i,j,k (t) + k B ik,jk (t)c jk = 0, (17)
where

A ijk (t) = u i,k (r, t)|ρ j,0 (r, t) i (18)
B ik,jk (t) = u i,k (r, t)| ∂ρ j ∂µ j 0 u j,k (r, t) i [START_REF] Delehaye | Critical velocity and dissipation of an ultracold Bose-Fermi counterflow[END_REF] are time-dependent coefficients (note that by definition of the inner product, B ik,jk = B * jk ,ik ).

B. Homogeneous cloud

We first consider the special case of a homogeneous system. The static density profiles are uniform, and the eigenmodes u jk are therefore orthogonal to ρ j . Coefficients A ijk vanish and the dynamics is set by the coupling coefficients B ik,jk .

As stated earlier, the unperturbed superfluids move at constant velocities V i and the eigenmodes are plane waves characterized by wave-vectors k j=1..N . We then have

B ik,jk (t) = g i,j ∂ρ i ∂µ i ∂ρ j ∂µ j e iki(Vi-Vj )t δ ki,kj . (20)
Under these assumptions, the linearized hydrodynamic equations ( 17) can be written as

ci,k + ω 2 i,k c 1,k = -ω 2 i,k j =i g ij ∂ρ i ∂µ i ∂ρ j ∂µ j c j,k e ik•Vij t (21)
From ( 21), we infer that the inter-species interactions only couple same-momentum modes and that the coupling constant oscillates in time at a frequency k • (V i -V j ). This behaviour is reminiscent of a parametric oscillator characterized by a dynamical instability that we interpret as follows: In absence of coupling, the free solutions of Eq. ( 21) are c i,k ∝ exp(±iω i,k t). To the leading order of the perturbation, we can insert this solution in the right-hand side of Eq. ( 21) which is now equivalent to an ensemble of N harmonic oscillators driven at frequencies k • V ij ± ω j,k . We notice that when

|k • V ij | ω i,k + ω j,k , (22)
the drive is resonant and leads to the instability of the system. We prove the existence of this instability in the the appendix B for a mixture of N = 2 superfluids. For 3D phonons characterized by a velocity c j , excitations can propagate in any direction and the instability criterion ( 22) can be reformulated as V ij ≥ c i + c j , which is the extended Landau criterion [START_REF] Castin | The landau critical velocity for a particle in a Fermi superfluid[END_REF]. For an elongated harmonic trap, this argument is still valid but phonons can only propagate along the trap axis. In this case the instability criterion selects a velocity window centered on c i + c j (see for instance [START_REF] Abad | Counter-flow instability of a quantum mixture of two superfluids[END_REF] and Appendix B).

C. Trapped superfluids

We consider now the case of a harmonically trapped mixture of superfluids. The parametric mechanism discussed in the case of a homogeneous system is still present. Consider for simplicity a situation where the trapping frequencies of the two clouds are equal. In this case, the relative position and the coefficients B ik,jk oscillate with frequency ω z . After expanding these coefficients, the qualitative argument put forward for homogeneous systems predicts that the parametric amplification of pairs of spatially-matching modes occurs when their frequencies meet the resonance conditions

ω 1,k + ω 2,k = nω z , (23)
where n is an integer. In Apppendix C we present a toy model leading to this resonance condition. This parametric mechanism is illustrated in the right panel of Fig. 3 where we notice that the modes excited in the two clouds share approximately the same wavelength. Another specific feature of the trap is the inhomogeneity of the density profiles allowing for non-zero coupling coefficients A ijk . Qualitatively, this mechanism can be interpreted as an excitation of the superfluid i by the mean-field potential g ij ρ j,0 (r -R j (t)). Moreover, this coupling affects only the large cloud. Indeed, at the scale of the small cloud j, the density profile of a larger cloud i is approximately flat and the overlap between ρ i,0 and any mode u jk is therefore vanishingly small according to the general properties of the operator L.

To interpret the result of the simulation displayed in Fig. 3, we consider the case N = 2, with, for simplicity, B ik,jk = 0. We note that Eq. ( 17) is formally equivalent to that of a driven linear oscillator. Moreover, if we expand the driving term with the relative displacement

z ij = (R i (t) -R j (t)) • u z , we obtain [38] A ij,k = ∞ n=1 (-1) n z n ij u i |∂ n z ρ j i . (24)
In Fig. [START_REF] Yao | Observation of coupled vortex lattices in a mass-imbalance Bose and Fermi superfluid mixture[END_REF], we assess the accuracy of this model using the first two harmonics (n = 1 and n = 2) of the dynamics of the breathing mode (k = 2) obtained from numerical simulations. We consider the amplitudes and the frequencies of the breathing mode as free parameters and we find ω 2,2 = 1.580(3)ω z , in agreement with the expected value ω 2,2 /ω z = 5/2 1.581.

Finally, we note that in the case of the dipole mode k = 1, the excitation is resonant with the mode-frequency and triggers a slow divergence of the amplitude of the mode. In Appendix D we show that this resonant behaviour gives rise to the dipole-mode frequency shift, expressed in Eq. 5.

IV. CONCLUSION AND IMPLICATIONS FOR THE SUPERFLUID CRITICAL VELOCITY IN TRAPPED SYSTEMS

Both numerical and analytical approaches conclude to the existence of two mechanisms for the coupling of the center of mass motion to low-lying excitation modes. The parametric scenario is associated with the formation of pairs of excitations in both superfluids. By contrast, the linear mechanism can be interpreted as the resonant creation of excitation in the large cloud. This mechanism is the analogous to the motion of a potential in the larger superfluid and is consequently very similar to Landau's argument. Our zero-temperature approach predicts a coherent mode-coupling. At finite temperature, these modes are damped and will pump energy out of motion of the center of mass of the two clouds, as observed experimentally [START_REF] Wen | Dipole modes of a superfluid Bose-Fermi mixture in the bcsbec crossover[END_REF][START_REF] Kean Loon Lee | Phase separation and dynamics of two-component Bose-Einstein condensates[END_REF][START_REF] Edmonds | Nonequilibrium kinetic theory for trapped binary condensates[END_REF].

The resonance condition for the parametric instability can be reinterpreted in terms of relative velocity of the two clouds. By analogy with the homogeneous case, the instability is triggered when the relative velocity cross the instability window centered on c 1 + c 2 . This scenario is supported by the toy model presented in Appendix C and is validated by our numerical simulations (see Fig. 3, where we observe that the amplitude of the mode grows only inside narrow time windows). By contrast, the linear coupling predicts a critical velocity equal to the sound velocity of the large cloud. In an imbalanced system as in [START_REF] Delehaye | Critical velocity and dissipation of an ultracold Bose-Fermi counterflow[END_REF], c 1 c 2 and the two velocities are similar. Experimentally, the critical velocity reported in [START_REF] Delehaye | Critical velocity and dissipation of an ultracold Bose-Fermi counterflow[END_REF] is therefore compatible with both scenarios.

x

→ x x s , t → t t s , u 1 = ψ 1 x -3/2 s , u 2 = ψ 2 x -3/2 s . (A1) with t s = 1 ω , x s = a ho , a ho = m 2 ω z . ( A2
)
The non-dimensional form of Eqs. ( 1)-(2) becomes:

i ∂u 1 ∂t = - ∇ 2 2d 1 + U a (r) + β 11 |u 1 | 2 + β 12 |u 2 | 2 u 1 , (A3) i ∂u 2 ∂t = - ∇ 2 2d 2 + U a (r) + β 21 |u 1 | 2 + β 22 |u 2 | 2 u 2 , (A4)
with u 1 and u 2 normalized to unity:

R 3 |u 1 | 2 = 1, R 3 |u 2 | 2 = 1. (A5)
The non-dimensional trapping potential takes into account the initial shift b (in a ho units) of the clouds:

U a (r) = d 2 2 γ 2 ⊥ (x 2 + y 2 ) + γ 2 z (z -b(t)) 2 , ( A6
)
where γ ⊥ = (ω ⊥ /ω) , γ z = (ω z /ω) and b(t) = 0 for t > 0. Dimensionless parameters in Eqs. (A3)-(A4) are expressed by:

β 11 = 4π 1 d 1 N 1 a 11 a ho , β 12 = 2π d 1 + d 2 d 1 d 2 N 2 a 12 a ho , (
A7)

β 22 = 4π 1 d 2 N 2 a 22 a ho , β 21 = 2π d 1 + d 2 d 1 d 2 N 1 a 12 a ho . (A8)
Each run is identified by the value of the parameter g 12 /g 22 = β 12 /β 22 and the value of b, the initial shift of the clouds. The explored values of these parameters is depicted in Fig. 6. The long-time integration for these 30 cases requested over 50000 hours of CPU time and generated teraoctets of data.

FIG. 6: Summary of parameters for which simulations were performed (blue dots). The red dashed line shows the amplitude above which the critical velocity is expected to be reached in the limit of a vanishing coupling between the superfluids.

To avoid clouds coming close to the boundaries during oscillations, the dimensions of the computational domain were fixed as L x = L y = 8 a ho and L z = 64 a ho . This allowed us to capture amplitude oscillations up to 16 a ho along the z direction, which is more than enough to obtain a relative velocity larger than the sum of the central sound velocities of the condensates (typically reached for 5.5 a ho ). The grid resolution was δx = δy = δz = a ho /16, resulting in computational grids of 128×128×1024 points following x, y and z directions, respectively.

To capture both fast and slow dynamics present during the oscillations of the clouds, we used a refined time resolution δt = 5 • 10 -4 /ω z and a large number of steps N t = 4 • 10 5 . This allowed us to simulate at least 32 periods of oscillations for each run case. is therefore driven by the coupling coefficient where we assumed that the amplitude A of the motion is smaller than the size of the larger (i = 2) cloud. Using the fact that the size of the static cloud is larger than that of the moving one and that the trapping potentials are identical for the two species, we can further simplify this expression as (r = 0) z 1 t sin(ω z,1 t).

A 12,1 = d 3 r
A 12,1 -m 1 ω 2 z,

(D4)

The onset of this divergence in a perturbation expansion is usually the signature of a shift of the natural oscillation frequency of the system. Let us indeed assume that the oscillation is shifted by δΩ. The density profile of species i = 1 now evolves as ρ 1 (r, t) = ρ 1,0 (z -A cos((ω z,1 + δΩ)t)). Using the fact that in the stationary regime, the density profile of the unperturbed system is given by the LDA expression ρ 1,0 (r) = ρ 1 (µ (0) 1 -V (r)), expanding the density profile to first order in δΩ yields ρ 1 (r, t) = ρ 2,0 (z -A cos(ω z,1 t))

dρ 1 dµ 1 0 m 1 ω 2 z,1 zAδΩt sin(ω z,1 t)

(D5)
We recover the frequency shift expressed in Eq. ( 5) by comparing Eq. (D5) to the general expansion [START_REF] Maddaloni | Collective oscillations of two colliding Bose-Einstein condensates[END_REF], where c 1,1 is given by Eq. ( 16).

ahoFIG. 1 :

 1 FIG. 1: Example of low amplitude oscillations of the relative position z rel of centers of mass of the two superfluids. Case b/a ho = 2 and g12/g22 = 0.031. The small atomic mass difference (m1 -m2)/m11 combined with a large atom number ratio N2/N1 = 30 induce a strong amplitude modulation of the small condensate oscillations (i = 1, in blue), while the oscillations of the large condensate are mostly not affected (i = 2, in red).

FIG. 5 :

 5 FIG. 5: Spectral weight of the breathing mode (k = 2) of the large cloud. Numerical results (blue solid line) are fitted using a driven oscillator model (black dashed line).

  ∂ρ 1 ∂µ 1 0 ρ 2,0 (r + A cos(ω z,1 t)u z ) z z 1 -m 2 ω 2 z,2 A cos(ω z,1 t) d 3 r

	∂ρ 2 ∂µ 2 0	∂ρ 1 ∂µ 1	z 2 z 1 (D1) ,

  We note that the driving term A 12,1 oscillates at the resonance frequency ω z,1 of the dipole mode. The solution to the equationc11 + ω 2 z,1 c 11 = -g 12 ω 2 z,1 A 12,1(D3)is therefore characterized by a linearly divergent behaviour in long time, and we obtain in the steady state

			regime
			c 11 (t) ∼	g 12 2	m 1 ω 3 z,1 A	∂ρ 2 ∂µ 2 0
	1 A	∂ρ 2 ∂µ 2 0	(r = 0) z 1 cos(ω z,1 t)
			(D2)

Acknowledgments

The authors thank N. Proukakis, C. Salomon and ENS ultracold Fermi group for helpful discussions. F. Chevy acknowledges support from ANR (grant SpifBox), ERC (advanced grant CritiSup2) and Fondation del Duca. P. Parnaudeau and I. Danaila acknowledge financial support from the French ANR grant ANR-18-CE46-0013 QUTE-HPC. Part of this work used computational resources provided by IDRIS (Institut du développement et des ressources en informatique scientifique) and CRI-ANN (Centre Régional Informatique et d'Applications Numériques de Normandie).

Appendix A: Details of numerical simulations Numerical simulations are performed using dimensionless variables. We take the usual scaling:

Appendix B: Parametric instability criterion for a mixture of two homogeneous superfluids

We prove here the existence of parametric instability for a mixture of N = 2 superfluids. In this case, Eq. ( 21) reduces to a set of decoupled two-dimensional problems governed by the following equations

The amplification rate γ is thus solution to the eigenvalue equation

We recover the equations derived in [START_REF] Abad | Counter-flow instability of a quantum mixture of two superfluids[END_REF] showing that for a small coupling, the acoustic modes become dynamically unstable when the velocity satisfies the condition

with

Appendix C: A toy model for the parametric instability in a trap

We describe the parametric mode-coupling in a trap using the following assumptions. First, we approximate the space structure of the mode by local plane waves characterized by a wave-vector k = ω k /c i (this is true for high frequency mode at the center of the trap, following the WKB approximation [START_REF] Jin | Hydrodynamic response of a trapped superfluid to a periodic perturbation[END_REF]). Second, we assume that the trap frequencies of the two clouds are identical and that the initial displacement of the two clouds are different. As a consequence, the relative distance between the centers of the two clouds evolves as R cos(Ωt), where Ω is the common trapping frequency.

With these assumptions, the amplitudes of two matching modes will satisfy equations similar to Eqs. (B1)-(B2), where the relative displacement

We can repeat the qualitative argument developed for the homogeneous cloud by noting that exp(ikR cos(Ωt)) can be expanded as

where J n are Bessel functions of the first kind. As before, each harmonic of the sum will give rise to a parametric instability when the resonance condition nΩ = ω 1,k +ω 2,k is met and, in this case, the effective coupling is proportional to J n (kR). Consider now the case of a coupling with a high-frequency phonon mode, as in Fig. 3). In this case, Ω ω α,k and n 1. But, at a given kR, J n (kR) is vanishingly small for large n and, to maintain a significant coupling, we need to take kR n. Indeed, from Eq. (C3), we infer that

For large values of n, the exponential term oscillates rapidly, unless the phase is stationary. The integral will therefore be dominated by values of θ close to θ * , defined by the stationary phase condition kR sin θ * = n, implying that kR ≥ n. Since nΩ = ω k,1 + ω k,2 , this condition can be recast as

We therefore recover that, even for an oscillatory motion, dissipation occurs when the maximal relative velocity ΩR is above the sum of the sound velocities of the two superfluids.

Appendix D: Dipole mode frequency shift

We consider here the linear coupling to the mode k = 1 corresponding to the dipole mode along the direction z. As mentioned in the general properties of the linearized hydrodynamic equations, this mode is associated with the eigenvector u i,1 = z/ z i , with eigenfrequency ω i,1 = ω z,i . We assume that the larger superfluid (i = 2) is at rest while the smaller one (i = 1) oscillates with an amplitude A.

The time-evolution of the dipole-mode amplitude c 1,1


Jian Cheng
email: jiancheng@buaa.edu.cn

Tao Liu

Feng Shi

Ruiliang Bai

Jicong Zhang

Haogang Zhu

Dacheng Tao

Peter J Basser

On Quantifying Local Geometric Structures of Fiber Tracts

In diffusion MRI, fiber tracts, represented by densely distributed 3D curves, can be estimated from diffusion weighted images using tractography. The spatial geometric structure of white matter fiber tracts is known to be complex in human brain, but it carries intrinsic information of human brain. In this paper, inspired by studies of liquid crystals, we propose tract-based director field analysis (tDFA) with total six rotationally invariant scalar indices to quantify local geometric structures of fiber tracts. The contributions of tDFA include: 1) We propose orientational order (OO) and orientational dispersion (OD) indices to quantify the degree of alignment and dispersion of fiber tracts; 2) We define the local orthogonal frame for a set of unoriented curves, which is proved to be a generalization of the Frenet frame defined for a single oriented curve; 3) With the local orthogonal frame, we propose splay, bend, and twist indices to quantify three types of orientational distortion of local fiber tracts, and a total distortion index to describe distortions of all three types. The proposed tDFA for fiber tracts is a generalization of the voxel-based DFA (vDFA) which was recently proposed for a spherical function field (i.e., an ODF field). To our knowledge, this is the first work to quantify orientational distortion (splay, bend, twist, and total distortion) of fiber tracts. Experiments show that the proposed scalar indices are useful descriptors of local geometric structures to visualize and analyze fiber tracts.

Introduction

Diffusion MRI (dMRI) provides a unique window to non-invasively reveal anatomical connections (i.e., fiber tracts) and white matter tissue properties in human brain [START_REF] Basser | In vivo fiber tractography using DT-MRI data[END_REF]. In dMRI, a typical processing pipeline before statistical analysis is: 1) fit a diffusion model (e.g., diffusion tensor imaging (DTI) [START_REF] Peter | MR Diffusion Tensor Spectroscropy and Imaging[END_REF]) to measured diffusion signals in each voxel; 2) calculate various voxel-wise scalar indices from the model parameters; 3) calculate local fiber directions in each voxel based on the model parameters, and then perform tractography. Scalar indices from diffusion models in step 2 are to quantify different tissue properties at different levels. For example, at the voxel level, fractional anisotropy (FA), generalized FA, return-to-origin probability, mean-squared displacement, and orientation dispersion index are indices for tissue properties inside a single voxel in DTI, Q-ball imaging [START_REF] Tuch | Q-ball imaging[END_REF], spherical polar Fourier imaging [START_REF] Cheng | Model-free and Analytical EAP Reconstruction via Spherical Polar Fourier Diffusion MRI[END_REF], NODDI [START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF], etc. At the local neighborhood level, gradient norm [START_REF] Pajevic | A continuous tensor field approximation of discrete DT-MRI data for extracting microstructural and architectural features of tissue[END_REF], sheet probability index [START_REF] Chantal Mw Tax | Sheet Probability Index (SPI): Characterizing the geometrical organization of the white matter with diffusion MRI[END_REF], and orientational distortion [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF] are indices for tissue properties within a local neighborhood.

Fiber tracts estimated by tractography are curves densely distributed in the 3D space. The spatial geometric structure of fiber tracts is known to be very complex in human brain [START_REF] Basser | In vivo fiber tractography using DT-MRI data[END_REF]. In a local spatial region, typical geometric structures of fiber tracts include splay (aka diverge, converge, or fanning), bend, twist, crossing, kissing, etc. Grid and sheet structures were also reported in some areas in human brain [START_REF] Wedeen | The geometric structure of the brain fiber pathways[END_REF][START_REF] Chantal Mw Tax | Sheet Probability Index (SPI): Characterizing the geometrical organization of the white matter with diffusion MRI[END_REF]. See Fig. 1 for a demonstration of splay, bend, and twist of fiber tracts. Compared with various voxel-wise indices, only a few tract-based indices have been proposed. Savadjiev et al. [START_REF] Savadjiev | Fusion of white and gray matter geometry: a framework for investigating brain development[END_REF] proposed the total dispersion index for each point in a tract. Curvature and torsion based on the Frenet frame along the curve are two famous features of a single tract [START_REF] Peter | New histological and physiological stains derived from diffusiontensor MR images[END_REF][START_REF] Batchelor | Quantification of the shape of fiber tracts[END_REF]. However, the Frenet frame and its features are not designed for a set of curves, thus they cannot quantify local geometric structures of fiber tracts shown in Fig. 1.

Cheng et al. [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF] proposed a framework called director field analysis (DFA) to quantify orientational order, dispersion of spherical functions, and orientational distortion (splay, bend, twist, and total distortion) of a spherical function field at the local neighborhood level. Since the DFA in [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF] was developed for a spherical function field, we call it as voxel based DFA (vDFA). vDFA does not work for fiber tracts. In this paper, inspired by studies of liquid crystals [START_REF] Peter | Introduction to liquid crystals: chemistry and physics[END_REF], we propose tract-based director field analysis (tDFA) which generalizes vDFA in [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF] to fiber tracts. tDFA is not proposed to replace vDFA, but adopts methods and concepts in vDFA to fiber tracts. In tDFA, we define total 6 scalar indices at each point in fiber tracts, where orientational order and dispersion quantify the alignment and dispersion of fiber tracts, and splay, bend, twist and total orientational distortion quantify local orientational distortion of fiber tracts. tDFA is applied directly on fiber tracts after tractography. To our knowledge, this is the first work to quantify orientational distortion of fiber tracts (i.e., splay, bend and twist).

Method

Directors, Oriented Curves and Unoriented Curves

A director, borrowed from studies of liquid crystals [START_REF] Peter | Introduction to liquid crystals: chemistry and physics[END_REF], is defined as a vector v that is equivalent to its negative v [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF]. A director v can be represented as a dyadic tensor vv T without sign ambiguity. A local fiber direction in a voxel (i.e., a peak of the ODF in that voxel) is a director by definition. "Sign the opposite orientation, i.e., γ(0) = C(1) and γ(1) = C(0). Note that these two oriented curves actually represent the same unoriented curve in geometry, however, for a point C(t) (i.e., γ(1 t)) in the curve, the two tangent vectors in these two parameterization have opposite directions, i.e., T C(t) = T γ(1 t) . Thus, tangent vectors of fiber tracts are actually directors.

The sign ambiguity of fiber tracts is not a problem for some cases. For example, the curvature and torsion do not change under the above two parameterization. However, it is indeed a crucial problem if we would like to compare two tangent vectors. Thus, we will use tools of directors in [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF] for calculation to avoid this problem. The difference of two directors v 1 and v 2 in the director representation, denoted as

Diff d , is Diff d (v 1 ; v 2 ) = s 1 v 1 s 2 v 2 ,
where

s i = ¦1 such that s 1 s 2 v T 1 v 2 ! 0 [5].

Orientational Order and Dispersion of Fiber Tracts

Considering a set of fiber tracts fC i g which are curves densely distributed in R 3 , for each point x P C i , inspired by liquid crystals [START_REF] Peter | Introduction to liquid crystals: chemistry and physics[END_REF][START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF], we define orientational order index (OO) as

OO = ∑ y2Ω (x) w(y; x) 3(u 1 (y) T u 1 (x)) 2 1 2 ; ( 1
)
where y is a point in a curve C j in a spatial neighborhood Ω (x) of x, u 1 (x) and u 1 (y) are unit norm tangent vectors of curves at x and y, w(y; x) is a spatial weighting function (e.g., uniform or Gaussian weighting) and ∑ yPΩ (x) w(y; x) = 1. See Fig. 2. We have OO P [ 0:5; 1]. Then, we define the orientational dispersion index (OD) as OD = 1 OO. Thus, OD P [0;

1:5]. If u 1 (y)
is parallel to u 1 (x), Vy P Ω (x), i.e., the least dispersion case, then OO = 1 and OD = 0 at x. If u 1 (y) is orthogonal to u 1 (x), Vy P Ω (x), i.e., the most dispersion case, then OO = 0:5 and OD = 1:5 at x. Note that the sign ambiguity of tangent vectors u 1 (y) and u 1 (x) dose not change the values of OO and OD. In practice, we normally set Ω x as a ball centered at x with a radius r = 4 mm, considering typical isotropic voxel size of a DW image is 2 mm. Note that the total dispersion in [START_REF] Savadjiev | Fusion of white and gray matter geometry: a framework for investigating brain development[END_REF] is defined as the mean of differences between tangent vectors in a circle and u 1 (x), while our definition of OD is inspired by the order parameter in liquid crystals [START_REF] Peter | Introduction to liquid crystals: chemistry and physics[END_REF][START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF].

Local Orthogonal Frame for a Set of Unoriented Curves

Curvature and torsion are typical features of a single curve based on the Frenet frame [START_REF] Peter | New histological and physiological stains derived from diffusiontensor MR images[END_REF]. However, the Frenet frame has several limitations: 1) It requires a direction (i.e., orientation) of a curve; 2) It is defined only for a single curve, not for a set of curves; 3) For a straight line, the Frenet frame is not well defined because the curvature of straight lines is 0.

We propose a local orthogonal frame at each point x of fiber tracts. The local orthogonal frame has three directors. The first director, denoted as u 1 (x), is the tangent vector at x, and the other two directors are in the orthogonal plane of the first director. We project all tangent vectors of fiber tracts in a local neighborhood Ω x of x onto the orthogonal plane, and define the covariance matrix of the projected directors as

Q x = ∑ y2Ω x w(y; x)u 1;? (y)u T 1;? (y);
where u 1;?

(y) = u 1 (y) (u T 1 (y)u 1 (x))u 1 (x): (2)
Then we set the second director u 2 (x) as the eigenvector associated with the largest eigenvalue of Q x , i.e., the principal component in PCA. In this way, u 2 (x) indicates the direction of the largest change of u 1 (x) in the orthogonal plane. Then, the third director is the cross product of the first two directors. See Fig. 2. We prove in Proposition 1 that when Ω x tends infinitesimally small, the local orthogonal frame will converge to the Frenet frame.

Proposition 1.

Let fC i g be a set of differentiable curves. Let x P C 1 be a point in C 1 , and u 1 (x) is the tangent vector at x in C 1 . Assume x is not a crossing point. Given a local ball Ω x = fy j ky xk rg with a radius r, let the constructed local orthogonal frame be fu 1 (x); u 2 (x); u 3 (x)g. For x P C 1 , let the Frenet frame be fu 1 (x); v 2 (x); v 3 (x)g. The vectors in two frames all have unit norm. Then we have lim r30

(u 2 (x) T v 2 (x)) 2 = 1, lim r30 (u 3 (x) T v 3 (x)) 2 = 1.

Orientational Distortions (Splay, Bend, and Twist)

Based on distortion analysis of liquid crystals [START_REF] Peter | Introduction to liquid crystals: chemistry and physics[END_REF], there are 3 types of orientational distortions as showed in Fig. 1. 1) splay: bending occurs perpendicular to the director; 2) bend: bending is parallel to the director; 3) twist: neighboring directors are rotated with respect to one another, rather than aligned. Inspired by liquid crystals [START_REF] Peter | Introduction to liquid crystals: chemistry and physics[END_REF][START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF], after we obtain the local orthogonal frame for each point in tracts, we define at each point three scalar indices to describe the three types of local distortions, and a total distortion index:

Splay index: s = s (u T 2 ∂ u 1 ∂ u 2 ) 2 + (u T 3 ∂ u 1 ∂ u 3 ) 2 ;
(3)

Bend index: b = s (u T 2 ∂ u 1 ∂ u 1 ) 2 + (u T 3 ∂ u 1 ∂ u 1 ) 2 ; (4)
Twist index:

t = s (u T 2 ∂ u 1 ∂ u 3 ) 2 + (u T 3 ∂ u 1 ∂ u 2 ) 2 ; (5)
Total distortion index:

d = p s 2 + b 2 + t 2 ; ( 6
)
where ∂ u 1 ∂ u i , i = 1; 2; 3, is the directional derivative of u 1 (x) along u i (x), i.e.,

∂ u 1 ∂ u i % Diff d (u 1 (x + ku i ); u 1 (x ku i )) 2k ; (7)
where k is small, and Diff d means the difference of two directors in the director representation [START_REF] Cheng | Director Field Analysis (DFA): Exploring Local White Matter Geometric Structure in diffusion MRI[END_REF]. See Section 2.1. Since the director field is not continuous in practice, we use an interpolation method to estimate u 1 (x + ku i )

and u 1 (x ku i ) from local neighborhoods of x + ku i and x ku i . Let z = x ¦ku i , then we estimate u 1 (z) from its neighborhood tangent directors of tracts:

u 1 (z) = arg min u;λ ∑ y2I(z) kλuu T w s (y; z)w b (y; x)u 1 (y)u T 1 (y)k 2 ; (8)
where

I(z) is a neighborhood of z, w s (y; z) = 1
Zky zk 2 , w b (y; x) is the bundle probability weight, and Z is the normalization factor such that ∑ yPI(z) w s (y; z)w b (y; z) = 1. we use the inverse distance weight so that if z is a point in a tract, then w s = 1 and u 1 (z) is just the tangent vector of the tract. It can be proved that the interpolated u 1 (z) is the actually the principal eigenvector of ∑ yPI(z) w(y; z)u 1 (y)u T 1 (y). We set I(z) (i.e., I(x ¦ku i )) as a ball with a radius of 2k, and k is normally set as 1 mm, half of the typical voxel size of 2 mm. If one would like to calculate indices using all points y from all fiber bundels, then w b (y; x) = 1 as a constant. If one would like to calculate indices using points y only from the same fiber bundle as x, then we may set w b (y; x) as a pre-determined probability of fiber clustering, or we could simply set w b (y; x) = 1 only if the angle between u 1 (y) and u 1 (x) is smaller than a threshold θ 0 (e.g., 45 ), and w b (y; x) = 0, otherwise. If there is no specific mention, we set w b (y; x) = 1 only if the angle is smaller than 45 .

Experiments

Synthetic data experiments. We generated three synthetic fiber tracts (splay, bend, and twist data) which demonstrate these three types of distortions, and then calculated the proposed scalar indices for each point in the fiber tracts. These scalar indices were used to color the fiber tracts. See Fig. 3. We omit OO, because OO = 1 OD. The first three rows in Fig. 3 show that the proposed splay, bend, and twist indices completely separate these 3 datasets, where only one among these 3 values is non-zero in each dataset. We also combined fiber tracts in splay and bend data, then calculated these indices using interpolation with points from all fiber bundles and with points only from the same fiber bundle (angle threshold θ 0 = 45 ).

The last two rows in the figure show that in a crossing area, it is better to use fiber tracts from the same bundle to calculate bundle specific distortion indices. Note that increased OD happens in crossing areas, compared with single bundle areas, and calculation of OD does not require an interpolation.

Real data experiments. We performed tractography on a publicly available dataset with a single subject from DIPY (dipy.org).The data has a single shell with b = 2000 s/mm 2 and 150 directions on the shell. The corpus callosum (CC) was extracted based on diffusion ODF and deterministic tracking. Then, all six indices were calculated from the fiber tracts of CC. Fig. 4 shows the scalar indices as colors of tracts. In the right subfigure, splay, bend, and twist values are set as red, green and blue color channels, respectively. The bending areas of CC are mainly in green, which means that the bend index is higher than slay and twist indices in those bending areas. We can also see red (high splay values) in fanning areas of CC, and blue (high twist values) when tracts are twisting in the 3D space.

Conclusion

We propose a unified mathematical framework called tract-based director field analysis (tDFA) with six scalar indices to quantify local geometric structure of fiber tracts. OD and OO are useful to quantify the degree of alignment and dispersion of fiber tracts; The distortion indices (i.e., splay, bend, twist and total distortion) demonstrate good sensitivity in both synthetic datasets and a real dataset on CC. To our knowledge, this is the first work to quantify orientational distortion (splay, bend, and twist) of fiber

Fig. 1 :

 1 Fig. 1: Demonstration of three types of distortions of fiber tracts.

Fig. 2 :

 2 Fig.2: Demonstration of tract-based DFA at a point x, denoted as the red point. We use tubes, instead of traditional arrows, to denote tangent vectors because of the sign ambiguity of the unoriented curves and tangent vectors. Ω x is a local neighborhood of x, denoted as the grey sphere. The yellow points are within Ω x , and their tangent vectors are used to calculate OO and OD, and construct the second and third directors in the local orthogonal frame. The pink plane is the orthogonal plane of red director u 1 (x). The red, green, and blue tubes are the 3 directors fu 1 (x); u 2 (x); u 3 (x)g in the local orthogonal frame. The 6 purple points fx¦ku i g alone these 3 directors in the local orthogonal frame are used to calculate the distortion indices.

Fig. 3 :Fig. 4 :

 34 Fig. 3: Synthetic datasets of fiber tracts. Each row is a set of fiber tracts. The OD, splay, bend, twist and total distortion indices are used to color the fiber tracts. Low/high values are in blue/red color. The last two rows are crossing of bend and splay data with two interpolation strategies.

Acknowledgement: This work was supported by Australian Research Council Projects FL-170100117, DP-180103424, LP-150100671, Intramural Research Program of NICHD (ZIA-HD000266), the National Key R&D Program of China (2016YFF0201002), and the Beijing Municipal Science and Technology Commission (BJMSJY-160153).