Nicolas Houel

Marc Elipot

Frédéric André

Philippe Hellard

Influence of Angles of Attack, Frequency and Kick Amplitude on Swimmer's Horizontal Velocity During Underwater Phase of a Grab Start

Keywords: swimming starts, kinematics, undulations, performance

The underwater phase of starts represents an important part of the performance in sprint swimming's events. Kinematics variables that swimmers have to take into account to improve their underwater phase of starts are unknown. The aim of this study was to determine the kinematics variables that improve performance during the underwater phase of grab starts. A three-dimensional analysis of the underwater phase of ten swimmers of national level was conducted. Stepwise multiple linear regressions identified the main kinematics variables that influence the horizontal velocity of the swimmer each 0.5 m in the range of 5 to 7.5 m. The results show that the kinematics parameters change during the range of 5 to 7.5 m of the underwater phase of the starts. For this population of swimmers, the results enable proposals of four principles to improve the underwater phase: i) to be streamlined at the beginning of the underwater gliding phase, ii) to start the dolphin kicking after 6 m, iii) to generate propulsive forces using only feet and legs during underwater undulatory swimming, iv) to improve the frequency of underwater undulatory swimming.

In 50 and 100 m swimming races, performance has been strongly linked to starts performance [START_REF] Arellano | A device for quantitative measurement of starting time in swimming[END_REF][START_REF] Mason | What we can learn from competition analysis at 1999 pan Pacific Swimming Championships[END_REF]. Starts performance is defined as the time observed between the start signal and the moment when the swimmer's head reaches 10 m [START_REF] Arellano | A device for quantitative measurement of starting time in swimming[END_REF] or 15 m [START_REF] Issurin | Track start versus Grab Start: Evidence of the Sydney Olympic Games. Biomechanics and Medicine in Swimming[END_REF][START_REF] Mason | What we can learn from competition analysis at 1999 pan Pacific Swimming Championships[END_REF].

Swimmers currently use the grab start or the track start in national and international events. In the grab start, swimmer places both feet at the front of the block and hands grabbing the front edge of the block. In the track start, with the same hand placement, swimmer places one foot at the front edge of the block, the other foot is back on the starting block [START_REF] Holthe | Kinematic comparison of grab and track starts in swimming[END_REF][START_REF] Welcher | Front or rear weighted track start or grab start: which is the best for female swimmers[END_REF].

The start is defined as a sequence of three phases: the impulse phase on the starting block (including reaction time), the flight phase, and the underwater phase. The results of studies conducted since the 1970's defining the most effective start remain contradictory [START_REF] Issurin | Track start versus Grab Start: Evidence of the Sydney Olympic Games. Biomechanics and Medicine in Swimming[END_REF][START_REF] Vilas-Boas | Biomechanical analysis of ventral swimming starts: comparison of the grab start with two track start techniques[END_REF]. However, there is some agreement that horizontal and vertical velocities, as well as minimized hydrodynamic resistance when Nicolas Houel (Corresponding Author), Marc Elipot, Frédéric André, and Philippe Hellard are with the Département Recherche Fédération Française de Natation, Pôle Natation, INSEP, Paris, France.

entering in the water have direct impact on the underwater phase [START_REF] Holthe | Kinematic comparison of grab and track starts in swimming[END_REF][START_REF] Miller | A kinematic and kinetic comparison of the grab and the track starts in swimming[END_REF][START_REF] Elipot | Analysis of swimmers velocity during the underwater gliding motion following grab start[END_REF].

The underwater phase is described as gliding phase with underwater leg propulsion [START_REF] Elipot | Analysis of swimmers velocity during the underwater gliding motion following grab start[END_REF]. When entering the water, swimmer's velocity (around 3.61 m•s -1 ) is greater than the one that can be produced by underwater propulsion [START_REF] Elipot | Analysis of swimmers velocity during the underwater gliding motion following grab start[END_REF]. During the gliding phase, the streamline position [START_REF] Marinho | Hydrodynamic drag during gliding in swimming[END_REF] directly influences the hydrodynamic resistance. Swimmer propulsive movements should ideally be initiated when the underwater velocity reaches between 2.2 and 1.9 m•s -1 [START_REF] Lyttle | Bet forces during tethered simulation of underwater streamlined gliding and kicking techniques of freestyle turn[END_REF]. Based on this result, [START_REF] Elipot | Analysis of swimmers velocity during the underwater gliding motion following grab start[END_REF] showed that the optimal beginning of the propulsive movement appeared when the swimmer center of mass reached 5.61-6.01 m from the wall start. In this situation, underwater propulsion becomes more critical. The underwater leg propulsion phase is currently called "underwater undulatory swimming" [START_REF] Arellano | Usefulness of the strouhal number in evaluating human underwater undulatory swimming[END_REF][START_REF] Von Loebbecke | A comparaison of kinematics of the dolphin kick in humans and cetaceans[END_REF].

The underwater undulatory swimming of fish and human locomotion has been studied using an oscillating foil [START_REF] Anderson | Oscillating foils of high propulsive effi[END_REF]. These studies showed that the oscillation of the foils created a vortex wake associated with drag or thrust [START_REF] Anderson | Oscillating foils of high propulsive effi[END_REF]. Under specific conditions, an oscillating foil can modify the von Karman street of the wake and generate propulsive force or lift for maneuvering [START_REF] Anderson | Oscillating foils of high propulsive effi[END_REF]. These specific conditions are defined by the Reynolds number, the angle of attack (α) and the Strouhal number [START_REF] Anderson | Oscillating foils of high propulsive effi[END_REF]. The Reynolds number is a dimensionless number which divides forward speed (U) and characteristic length of the swimmer by kinematic viscosity. The Reynolds number represents the nature (laminar to turbulent) of the flow circulation around the body. The angle of attack is the angle between the flow velocity and the chord of the body. The angle of attack influences directly the drag and lift coefficient [START_REF] Rouboa | The effect of swimmer's hand/forearm acceleration on propulsive forces using computational fluid dynamics[END_REF]. The Strouhal number is a dimensionless number which divides strokes frequency ( f ) and amplitude (A) by forward speed [START_REF] Taylor | Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency[END_REF]. The Strouhal number represents the ratio of unsteady and steady motion. The Strouhal number is an important factor in determining propulsion forces and swimming efficiency [START_REF] Anderson | Oscillating foils of high propulsive effi[END_REF]. [START_REF] Taylor | Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency[END_REF] clearly defined the relation between the Strouhal number, the angle of attack and the aerodynamic or hydrodynamic coefficient which depends of the Reynolds number Re. For them: "The Strouhal number is bound to affect the hydrodynamic force coefficients and propulsive efficiency, because it defines the maximum angle of attack associated with modifications of vortices, which are the source of hydrodynamic force production" [START_REF] Taylor | Flying and swimming animals cruise at a Strouhal number tuned for high power efficiency[END_REF].

All these specific conditions can't be easily estimated during the start underwater phase because of the velocity decrease of the swimmer, as well as the modification of the swimmers' segment orientation and coordination. But some of these (U, α, A, f) can be described and studied during each part of the underwater phase. From what we know, each of the independent phases of gliding and underwater leg propulsion have been clearly investigated. But few researches clearly focus on the swimmers' strategy to optimize speed during the underwater phase of the start. The main goals of this research are: i) to estimate the velocity of the swimmer's center of mass and hip during the underwater phase of grab start and ii) to define the determinant factors in this underwater phase related to performance.

Method

The experiment took place at the INSEP swimming pool (French National Institute of Sport, INSEP, France). Ten volunteers' swimmers of national level (mean ± SD: age = 21.41 ± 4.54 years; height = 183.33 ± 4.88 cm; mass = 75.8 ± 5.09 kg) participated in the study. All were informed of the objectives and signed a consent form. All swimmers practiced grab start and underwater undulatory swimming on a regular basis and used them for competitive racing. For the purpose of the study, they were asked to perform a competitive grab start as efficiently as possible.

The underwater area between the wall (0 m) to 10 m was recorded using three camcorders (Panasonic NV-GS17 and Sony DCR-HC20E). Two camcorders (camera 1, 2) were fixed behind windows (Figure 1). These camcorders recorded a sagittal view of the swimmers' trajectories. The third camcorder was placed in a waterproof case (Sony HCG Sports pack waterproof case). This third video camera recorded a slanting view of the swimmer motion. The angles between the principal axis of camera 3 and the other cameras were between 60° and 70°. The underwater area was divided into two zones measuring 5 × 2 × 2 m (length, width, weight): The first zone was from the start wall to the 5 m distance. The second zone was from the 5 m to the 10 m distance. Each zone was recorded by two cameras. Cameras 1 and 3 recorded the first zone. Camera 2 and 3 recorded the second zone (Figure 1). To limit the effect of the image distortions (due to camera lens deformation) on reconstruction accuracy, only the points contained in the 2/3 center of the camera field have been reconstructed. To improve the reconstruction accuracy, the cameras were positioned to minimize optical refraction effects (Snell's law), the long distance between the cameras and the center of each zone [START_REF] Kwon | Object plane deformation due to refraction in two dimensional underwater motion analysis[END_REF] and the camcorders' optical axes were perpendicular to the air-water interface plane. Cameras were synchronized with a light flash. The sampling frequency was 25 Hz. The video was interlaced scan and the both odd and even fields were used. The gliding and underwater undulatory swimming phase was recorded from the instant when the swimmers were completely under water until the instant they began arm propulsion. To minimize errors during the digitizing process, the two sides of the swimmer's body were assumed to be symmetric. Only the right side was digitized. Nine anatomical landmarks were identified on the swimmers with Simimotion software: the toe, the lateral malleolus, the knee, the iliac spine, the acromion, a finger tip, the wrist, the elbow, and the center of the head. A modified double plane direct linear transformation method (inspired by [START_REF] Drenk | A 3D video technique for analysis of swimming in a flume[END_REF] was used to calculate the landmark coordinates in space. Space mean reconstruction accuracy, calculated as described by [START_REF] Kwon | Effects of light refraction on accuracy of camera calibration and reconstruction in underwater motion analysis[END_REF], was 6.2 mm. Data were filtered with a second-order Butterworth filter [START_REF] Winter | Biomechanics and motor control of human movement[END_REF]. Cut-off frequencies were between 5 and 7 Hz. The landmark's positions associated with Dempster's anthropometric data (1959) were used to determine the trajectory of the center of mass. Using the coordinates of the landmarks and the center of mass, the following variables were defined.

• The horizontal velocity of the center of mass (Vx g ) and hip (Vx h ) • The angle of attack of segments trunk (α tr ), thigh (α th ), leg (α le ), foot (α fo ) which were defined as the angle between each segment and the velocity vector of this proximal landmarks • The mean kick frequency ( f ) of the underwater undulatory swimming or the reverse period of the ankle motion • The mean kick amplitude (A) or the mean peak to peak amplitude of the ankle during the underwater undulatory swimming. These variables were calculated for each swimmer every 0.5 m between 5 m to 7.5 m. After 7.5 m, 30% (three swimmers) of the total number of swimmers began arm propulsion. For the range of distance, the normality of the variables was tested using Lilliefors test [START_REF] Lilliefors | On the Kolmogorov-Smirnov test for normality with mean and variance unknown[END_REF]. Mean frequency, kick amplitude and velocity of the swimmer's center of mass and hip were calculated at each normalized distance. The effect of the independent variables (α tr , α th , α le , α fo, f, A) on the dependent variables Vx g and Vx h at each normalized distance was analyzed using multiple stepwise linear regression. Correlation between independent variables (α tr , α th , α le , α fo, f, A) and dependant variables Vx g and Vx h at each normalized distance were determined. The limit of significance was set as p ≤ .05.

Results

Mean velocity of the swimmer's center of mass and hip decreased during the underwater undulatory swimming phase of a grab start (Figure 2). Multiple stepwise linear regression (Table 1 and2) and correlation (Table 3 and4) analysis showed that different independent variables (α tr , α th , α le , α fo, f, A) were ranked as best predictors of the underwater horizontal velocity Vx g and Vx h , at each distance between 5.5 and 7.5 m. These independent variables had different effect on the center of mass (Vx g ) and the hip velocities (Vx h ). At 5 m, mean velocities were Vx g = 2.46 ± 0.17 m•s -1 and Vx h = 2.59 ± 0.21 m•s -1 . No significant result was found at 5 m. At 5.5 m, mean velocities were Vx g = 2.18 ± 0.2 m•s -1 and Vx h = 2.14 ± 0.26 m•s -1 . Multiple stepwise regression analysis showed that the angle of attack of the trunk (α tr ) and amplitude (A) were the variables that present better Vx g (R 2 = .44; p = .036) and Vx h (R 2 = .57; p = .0012). At this distance, Vx g was significantly correlated with α tr (r = -0.64, p = .042), A (r = -0.66, p = .036) and α le (r = -0.63, p = .049). Vx h was significantly correlated with α tr (r = -0.75, p = .012), A (r = -0.69, p = .026), α le (r = -0.73, p = .016) and α th (r = -0.63, p = .046). At 6 m, mean velocities were Vx g = 1.99 ± 0.13 m•s -1 and Vx h = 2.05 ± 0.11 m•s -1 . Multiple stepwise regression analysis showed that

Kinetic Parameters r value P value

Angle of attack of trunk (α tr ) at 5.5 m -0.752 0.012 Amplitude (A) at 5.5 m -0.691 0.026 Angle of attack of leg (α le ) at 5.5 m -0.731 0.016 Angle of attack of thigh (α th ) at 5.5 m -0.638 0.046 Angle of attack of thigh (α th ) at 6.5 m -0.837 0.002 Angle of attack of trunk (α tr ) at 7 m -0.731 0.016 Frequency ( f ) at 7.5 m 0.826 0.011 Angle of attack of trunk (α tr ) at 7.5 m -0.738 0.036 the angle of attack of the foot (α fo ) was the variable that present better Vx g (R 2 = .43; p = .04). At this distance, Vx g was only significantly correlated with α fo (r = -0.65, p = .040). At 6.5 m, mean velocities were Vx g = 1.93 ± 0.14 m•s -1 and Vx h = 1.9 ± 0.19 m•s -1 . Multiple stepwise regression analysis showed that the angle of attack of the thigh (α th ), the angle of attack of the leg (α le ) and the frequency ( f ) were the best predictors of Vx g (R 2 = .86, p = .005). The angle of attack of the thigh (α th ) and the frequency ( f ) were the best predictors of Vx h (R 2 = .88, p = .0005). At this distance, Vx g was only significantly correlated to α th (r = -0.64, p = .042). Vx h was only significantly correlated to α th (r = -0.83, p = .002). At 7 m, mean velocities were Vx g = 1.74 ± 0.25 m•s -1 and Vx h = 1.73 ± 0.24 m•s -1 . Multiple stepwise regression analysis showed that the frequency ( f ) was the variable that presents better Vx g (R 2 = .45, p = .032). The angle of attack of the trunk (α tr ,) was the variable that presents better Vx h (R 2 = .53, p = .016). Vx g was significantly correlated to f (r = .67, p = .032). Vx h was significantly correlated to α tr (r = -0.73, p = .036). At 7.5 m, mean velocities were Vx g = 1.76 ± 0.17 m•s -1 and Vx h = 1.81 ± 0.17 m•s -1 . Mean kick frequency were f = 2.32 ± 0.22 Hz and mean amplitude were A = 70.8 ± 6.04 cm. Multiple stepwise regression analysis showed that the frequency ( f ) was the variable that presents better Vx h (R 2 = .68, p = .011). Vx h was significantly correlated to f (r = .82, p = .011) and α tr (r = -0.73, p = .036).

Discussion

This study showed that mean velocity of the swimmer decreases during the underwater phase of a grab start.

In this situation, different kinematics variables should be used by the swimmer who wants to limit the loss of his center of mass and hip horizontal velocities during underwater undulatory swimming. These kinematics variables are different at each distance of the underwater undulatory swimming.

The most important limitation of the results using multiple stepwise linear regression analysis is the low number (n) of subjects to predict independent variables (k) using stepwise regression [START_REF] Fonton | Comparaison empirique de methods de prediction en regression lineaire multiple[END_REF].

In our present study: n < 3×k and k < n. For R 2 > .4, the multiple stepwise regression analysis can only select the independent variables that present better the dependant variables. For R 2 > .8, the multiple stepwise regression analysis presents the independent variables as best predictors and explains more than 50% of the event variability [START_REF] Fonton | Comparaison empirique de methods de prediction en regression lineaire multiple[END_REF]. The regression results were based on the specific way that independent variables were entered on the regression equation.

This study showed that the horizontal velocity of the center of mass and the hip during underwater undulatory swimming were explained by different kinematics variables. This is in agreement with results observed in intracycle swimming where the hip motion can't be used as center of mass indicator to explain the swimmer's motion [START_REF] Barbosa | Comparative study of the reponse of kinemtical variables from the hip and the center of mass in butterfliers[END_REF][START_REF] Psycharakis | Validity of the use of a fixed point for intracycle velocity calculation in swimming[END_REF]. In the case of our present study, multiple stepwise regression analysis realized on the center of mass and the hip velocities helped to understand the swimmer's motion and how to optimize his underwater undulatory swimming.

The results of the current study indicated that the swimmer should stay in a streamlined position and limit the underwater undulatory swimming before he reaches 6 m. At this distance, mean velocities were Vx g = 1.99 ± 0.13 m•s -1 and Vx h = 2.05 ± 0.11 m•s -1 . These results are in agreement with the previous works [START_REF] Lyttle | Bet forces during tethered simulation of underwater streamlined gliding and kicking techniques of freestyle turn[END_REF][START_REF] Elipot | Analysis of swimmers velocity during the underwater gliding motion following grab start[END_REF]. In case the swimmer would have begun its underwater undulatory swimming too early, hydrodynamic resistances would have increased and limited the performance of the underwater phase of the start [START_REF] Elipot | Analysis of swimmers velocity during the underwater gliding motion following grab start[END_REF]. The decrease of the angle of attack of the trunk (α tr ) and the foot (α fo ) as selected variable to improve Vx g and Vx h , confirmed that the angles of attack directly influence drag and lift coefficients of the body [START_REF] Rouboa | The effect of swimmer's hand/forearm acceleration on propulsive forces using computational fluid dynamics[END_REF][START_REF] Naemi | Hydrodynamic glide efficiency in swimming[END_REF].

Decreasing the angle of attack of the thigh (α th ) at 6.5 m associated with increasing the angle of attack of the leg (α le ) and increasing the frequency ( f ) have improved the horizontal velocities Vx g and Vx h . So angles of attack of the human segments have different impacts on the swimming propulsion efficiency. This result is in agreement with studies realized on fish [START_REF] Von Loebbecke | A comparaison of kinematics of the dolphin kick in humans and cetaceans[END_REF]. As like for dolphins, the human underwater undulatory swimming tends to minimize the displacement of the drag produced by the forward parts of the body, and to maximize the displacement of the thrust produced by the legs and the feet. The displacement wave that travels the length of the body also has a small magnitude along the torso, and reaches a maximum at the toes.

At 7.5 m, the mean velocities of the swimmers (Vx g = 1.76 ± 0.17 m•s -1 ; Vx h = 1.81 ± 0.17 m•s -1 ), mean kick frequency (f = 2.32 ± 0.22 Hz) and mean amplitude (A = 70.8 ± 6.04 cm) were higher than the results observed in the study of [START_REF] Gavilan | Underwater undulatory swimming: study of frequency, amplitude and phase characteristics of the "body wave[END_REF][START_REF] Von Loebbecke | A comparaison of kinematics of the dolphin kick in humans and cetaceans[END_REF]. The result of the multiple stepwise regression analysis and the comparison with the results of the previous studies [START_REF] Gavilan | Underwater undulatory swimming: study of frequency, amplitude and phase characteristics of the "body wave[END_REF][START_REF] Von Loebbecke | A comparaison of kinematics of the dolphin kick in humans and cetaceans[END_REF] confirms that the swimmer can improve his velocity by increasing kick frequency and maintaining large kick amplitude.

Conclusion

This study tends to suggest that the swimmer should apply three principles to improve his underwater start phase: i) to be streamlined with linear adjustment of the trunk and the lower body segments at the beginning of the underwater gliding phase, ii) to start the underwater undulatory swimming after 6 m, iii) to generate propulsive kicking using only feet and legs during underwater undulatory swimming, iv) to improve the kick frequency of underwater undulatory swimming.

Figure 1 -

 1 Figure 1 -Experimental setup.

Figure 2 -

 2 Figure 2 -Variation of the mean velocity center of mass (Vx g ) and hip (Vx h ) with standard deviation.

Table 1 Regression equation with statistical coefficient (R 2 ) and P values of independent variables which present the horizontal velocity of the center of mass (Vx g )

 1

			Correlation
	Distance (m)	Equation	Coefficient, R 2	P value
	5	-	-	-
	5.5	Vx g = -0.03 A	R 2 = .438	P = 0.036
	6	Vx g = -0.06 α fo	R 2 = .428	P = 0.040
	6.5	Vx g = -0.01 α th + 0.004 α le + 0.33 f	R 2 = .862	P = 0.005
	7	Vx g = 0.81 f	R 2 = .455	P = 0.032
	7.5	-	-	-

Table 2 Regression equation with statistical coefficient and P values of independent variables which present the horizontal velocity of the hip (Vx h )

 2

	Correlation

Table 3 Significant correlation between kinetic parameters and the horizontal velocity of the center of mass (Vx g )

 3

	Kinetic Parameters	r value	P value
	Angle of attack of trunk (α tr ) at 5.5 m	-0.648	0.042
	Amplitude (A) at 5.5 m	-0.662	0.036
	Angle of attack of leg (α le ) at 5.5 m	-0.632	0.049
	Angle of attack of foot (α fo ) at 6 m	-0.654	0.040
	Angle of attack of thigh (α th ) at 6.5 m	-0.648	0.042
	Frequency ( f ) at 7 m	0.674	0.032

Table 4 Significant correlation between kinetic parameters and the horizontal velocity of the hip (Vx h )

 4

Acknowledgments

The authors wish to thank the "Ministère des Sports" and the "Fédération Française de Natation" for financing this study. The authors also wish to thank Frédéric Frontier, Frédéric Clerc, Isabelle Amaudry and the "INSEP" for logistic support, and Jean-Lyonel Rey, Stéphane Lecat, Eric Boissière and Yves Thomassin for their contributions.


Paul Best
email: paul.best@ircam.fr

Jean Bresson
email: jean.bresson@ircam.fr

Diemo Schwarz
email: diemo.schwarz@ircam.f

Musical Gesture Recognition Using Machine Learning and Audio Descriptors

Keywords: Machine learning, gesture, hidden Markov chains, audio descriptors

We report preliminary results of an ongoing project on automatic recognition and classification of musical "gestures" from audio extracts. We use a machine learning tool designed for motion tracking and recognition, applied to labeled vectors of audio descriptors in order to recognize hypothetical gestures formed by these descriptors. A hypothesis is that the classes detected in audio descriptors can be used to identify higherlevel/abstract musical structures which might not be described easily using standard/symbolic representations.

I. INTRODUCTION

Musical structures carry abstract perceptual features and characteristics that can be straightforward to identify by composers and/or listeners, but difficult or impossible to formally describe using the elements of standard score representations (e.g. identifying harmonic/melodic patterns etc.) Composers and authors often use the term of gesture to characterize these dynamic elements constituting musical forms, as an analogy with the idea of gesture in physical movements [START_REF] Farhang | Modelling a gesture: Tak-Sīm for string quartet and live electronics[END_REF], [START_REF]Musical Gestures: Sound, Movement, and Meaning[END_REF].

Physical movements can be analyzed using machine learning techniques: Ircam's XMM library [START_REF] Franc ¸oise | A Multimodal Probabilistic Model for Gesture-based Control of Sound Synthesis[END_REF], for instance, uses a combination of Gaussian mixture models (GMM) and hidden Markov models (HMM) in order to recognize performed gestures from parallel streams of descriptors extracted from motion capture (Cartesian coordinates, speed, etc.) HMM have shown good performance processing movements in such timeseries [START_REF] Lee | An HMM-based threshold model approach for gesture recognition[END_REF], [START_REF] Bevilacqua | Continuous realtime gesture following and recognition[END_REF], and require much smaller training sets than neural networks to build reliable models. We propose here that similar models could be used by composers working with audio descriptors to recognize and classify musical gestures.

This hypothesis is tested on a real, specific example, in the context a composer's musical research residency project at Ircam (Paris). We developed a binding of the XMM library in the OpenMusic composition and visual programming software [START_REF] Bresson | Next-generation Computer-aided Composition Environment: A New Implementation of OpenMusic[END_REF], associated with audio description and analysis tools, and used this framework to run preliminary experiments.

A general objective of this work is to provide music composers with tools and interfaces to work with machine learning in a computer-aided composition environment [START_REF] Bresson | From Motion to Musical Gesture: Experiments with Machine Learning in Computer-Aided Composition[END_REF]. HMM techniques fit well with this situation, where algorithms will be trained on data input by end-users, rather than on databases or other large datasets distributed for instance over the web. A crucial part of the work, however, resides in the strategy employed to efficiently test and select the relevant features and parameters to use in order to carry out a recognition task given some input data type and training set at hand.

II. PROTOCOL

The dataset preparation, model testing and validation are carried out and controlled in OpenMusic visual programs such as the one visible in Fig. 1, using the OM-XMM library.

A. Dataset

The data-set used in this experiment is a 9-minutes-long extract of the piece Zāmyād for cello by contemporary music composer Alireza Fahrang. 1 This piece is widely built around the concept of musical gestures. Indeed, some easily noticeable patterns occur -with variations -throughout the piece. The definition of a gesture in such musical context can be relatively subjective, and based both on perceptual and more abstract compositional considerations (gestures are for instance defined by the composers in terms like "an unstable glissando sound with reversed attack"). The composer identified 20 different gesture classes in the score, and annotated the audio extract accordingly, cutting it into 187 labeled segments (the approx. duration of a segment is between 1 and 3 seconds). The number of elements per class varies between 1 (a gesture that appears only once in the whole training set) and 22.

B. Audio descriptors analysis

Each audio sample is analyzed and turned into a vector of signal descriptors using the IAE/pipo libraries [START_REF] Schnell | IAE & IAEOU: The IMTR Audio Engine[END_REF], [START_REF] Schnell | PiPo, A Plugin Interface for Afferent Data Stream Processing Modules[END_REF]. The main audio descriptors considered in this study are the Mel-frequency cepstral coefficients (MFCCs), the Fundamental Frequency, Energy, Periodicity, First-order autocorrelation coefficient (AC1), Loudness, Centroid, Spread, Skewness, and Kurtosis. We use the first 12 MFCCs and therefore work with a total of 21 descriptors, sampled according to a given analysis window size and overlap factor. Not all of these descriptors are relevant: further in this paper, we discuss strategies to select and combine them efficiently.

C. Model validation

We use a 10-fold cross-validation to measure our model's performances: for 10 folds, we successively train the model with 90% of the data, and compute the accuracy of recognition for the training and remaining data sets. Eventually the average of each accuracy is computed. Our measure for the model's performance thus consists in two average accuracy values: the training-set accuracy and the test-set accuracy. The final purpose of this project being to recognize gesture on new audio pieces, we aim at improving the test set accuracy first and foremost.

Confusion matrices allow to identify among the different classes of segments, the ones that are problematic to recognize. Table I reproduces part of a matrix such as the one output by the OpenMusic visual program in Fig. 1. We can observe how some classes (such as C) are easily recognized, while some others are not (e.g. K, 45% times confused with L).

D. Hyperparameters optimization

Feature selection and the tuning of hyperparameters are key elements for the recognition model performance. In addition to the choice of audio descriptors, three main parameters can vary in the model:

• The number of hidden states of the Markov model.

A C G H K L ...

PREDICTED CLASS

• The two regularization coefficients (offset added to the covariance matrices of the Gaussian distributions at each re-estimation in the algorithm):

-Relative: offset relative to data-variance.

-Absolute: minimum offset value.

A simple format has been defined to describe the model features and hyperparameters in the validation process. For instance: <mfcc,descr> (0 2 15) 15 (0.1 0.05) means:

• Audio descriptors are extracted using the mfcc and descr analysis modules (which respectively generate 12 MFCC coefficients and the vector of 9 descriptors). • The subset of selected descriptors are the descriptors number 0, 2 and 15 (among 21 in total). 2• 15 hidden states are used in the Markov model.

• The relative and absolute regularization values in the model are 0.1 and 0.05. Because results can be highly dependent on the dataset used for validation, it is important to have an easy and automated way for users to optimize their model and fine-tune these hyperparameters. Different techniques have been applied to choose optimum values for these different parameters.

1) Brute force algorithm: This algorithm consists in an exhaustive test (10-fold cross-validation) of a whole list of hyperparameters sets. It gives a clear view of which parameters set works and which doesn't, but requires a pre-definition of all hyperparameters configurations to be tested beforehand (testing all possible combinations would take too much time).

2) Genetic algorithm: To avoid the test of an exhaustive set of hyperparameter configurations, a genetic approach has been developed: starting from a small number of hyperparameter configurations (in principle, already identified for fairly good performances), and applying small variations to them at each iteration.

The algorithm searches for optimal number of states and regularization values, and selects audio descriptors to use in the model. The range of random variations for each hyperparameter is set as follows (but adjustable for different cases and data-sets):

• Number of states: [ -4 , +4 ];

• Relative regularization: [ -0.05 , +0.05 ];

• Absolute regularization: [ -0.005 , +0.005 ]; At each iteration a new audio descriptor is also randomly chosen and added or removed in the descriptors list.

After applying these random variations to the hyperparameters configurations, the algorithm evaluates the newly formed configurations and keeps the 4 best-performing (in terms of test-set accuracy) among the old and new ones.

The main limitation of this algorithm is that the search is likely to get stuck in local optima: it is efficient in finetuning hyperparameter configurations already known for good performances, but is not reliable in a global search context.

E. Computation time

These search and evaluation procedures for hyperparameters take time to compute. Although efforts were made to optimize them, in average it takes approximately 5 minutes on a standard personal computer to train and cross-validate a model with a given configuration of hyperparameters on the 187 elements of our dataset.

III. RESULTS

In this section, we report elements of the model performance depending on selected hyperparameters, focusing on the selection of audio descriptors used for training and classifying gestures. Additional results and data can be found at http://repmus.ircam.fr/paco/cbmi18.

Our experiments have shown that the more descriptors are used in the model, the more hidden states are required for the HMM to perform well. We used first-order Hierarchical HMMs, with only one gaussian per state. In the following sections III-A, III-B, and III-C, the number of hidden state has been fixed manually. In section III-D, this number is searched using the genetic algorithm. In all reported results, relative and absolute regularization values were set to 0.1 and 0.05 respectively.

A. One descriptor

Table II reports the model performance when trained and run with a single audio descriptor. The model is tested with 10 hidden states. These initial results show that none of the descriptors is reliable enough by itself to accurately analyze and recognize gestures in the audio samples. We study the model performances using combined descriptors in the following sections.

B. Combination of two descriptors

Table III presents performances of models trained with selected combinations of two descriptors. Those models were tested with 15 hidden states.

Performances are visibly improved. We also note that the best combinations of two do not necessarily correspond to the combination of best-performing descriptors when tested alone (section III-A), which emphasize emerging, non-predictable aspects of the descriptors selection task. For example, Periodicity alone gave a 0.06 test-set accuracy, but performs a 0.44 test-set accuracy when combined with Spread.

D. Genetic algorithm results

Using the directed search of the genetic algorithm allowed to fine-tune relevant hyperparameter configurations found previously. Table V shows the four best hyperparameter sets obtained with the algorithm. For all of them, the relative and absolute regularization values converged to 0.1 and 0.05 respectively. The resulting selected sets of descriptors are far more complex than the previous combinations, and the performance in recognition is slightly improved. IV. CONCLUSION We tested machine learning techniques designed for motion processing on audio signal descriptors, in order to recognize abstract musical gestures in audio extracts: elements of a temporal form which carry a perceptual identity, but can not be easily described using standard score descriptions.

We also developed an integrated framework including the IAE and XMM libraries in the OpenMusic computer-aided music composition and visual programming environment, which allows users (composers) to prepare data-sets, train and test HMM models, and select relevant hyperparameters. 3The choice of hyperparameters still requires a significant involvement from the user, and is highly dependent on the datasets and types of gestures to be recognized. The work accomplished so far facilitates a workflow for such hyperparameter optimization, thanks to built-in tools such as crossvalidation and genetic algorithms.

We can not be fully satisfied yet with our model performance on the tested dataset (test-set accuracy of 0.65 at maximum) but no conclusion can be drawn from one single dataset: further work will focus on applying and improving these tools and workflow on other examples and datasets.

Fig. 1 .

 1 Fig. 1. Example of an OpenMusic patch and graphical interface (v. o7) including XMM model testing and validation tools applied to audio signal descriptors.

TABLE I A

 I CONFUSION MATRIX EXCERPT OUTPUT AFTER THE TEST PROCEDURE OF Zāmyād'S GESTURE RECOGNITION MODEL.

	GROUND TRUTH CLASS	A C G H K L	0.89 0.10	0.10 1 0.09 0.05	0.81 0.16	0.05 0.74 0.19	0.05 0.54 0.05	0.45 0.71

TABLE II ACCURACY

 II OF THE MODEL WITH SINGLE DESCRIPTORS.

	Descriptor	Test set accuracy	Training set accuracy
	MFCC 1	0.21	0.30
	MFCC2	0.24	0.29
	MFCC 3	0.24	0.31
	MFCC 4	0.17	0.27
	MFCC 5	0.11	0.24
	MFCC 6	0.12	0.23
	MFCC 7	0.19	0.30
	MFCC 8	0.12	0.12
	[...]	[...]	[...]
	Frequency	0.23	0.29
	Energy	0.04	0.05
	Periodicity	0.06	0.12
	AC1	0.04	0.06
	Loudness	0.24	0.30
	Centroid	0.21	0.35
	Spread	0.22	0.34
	Skewness	0.14	0.22
	Kurtosis	0.11	0.18

TABLE III ACCURACY

 III OF THE MODEL WITH COMBINATIONS OF 2 DESCRIPTORS.Table IV presents a selection of performance results from models trained and run with combinations of 3 descriptors. Those models were tested with 18 hidden states. As we can see, results are again sensibly better, and relevant combinations of descriptors become more salient (e.g. Spread, Centroid, and the first few MFCC coefficient).

	Descriptors	Test set accuracy	Training set accuracy
	MFCC 1, 4	0.44	0.66
	MFCC 2, Frequency	0.47	0.64
	MFCC 2, Loudness	0.47	0.63
	MFCC 2, Spread	0.46	0.65
	MFCC 3, Spread	0.45	0.65
	Frequency, Loudness	0.46	0.66
	Frequency, Spread	0.45	0.63
	Periodicity, Spread	0.44	0.56
	Centroid, Spread	0.52	0.69
	C. Combination of three descriptors

TABLE IV ACCURACY

 IV OF THE MODEL WITH COMBINATIONS OF 3 DESCRIPTORS.

	Descriptors	Test set accuracy	Training set accuracy
	MFCC 2, Centroid, Spread	0.57	0.79
	Periodicity, Centroid, Spread	0.57	0.77
	MFCC 3, 4, Spread	0.56	0.79
	MFCC 2, Frequency, Spread	0.56	0.75
	MFCC 2, 4, Centroid	0.56	0.76

TABLE V MODEL

 V ACCURACIES WITH HYPERPARAMETER SETS AND BESTCOMBINATIONS OF DESCRIPTOR AS FOUND BY THE GENETIC ALGORITHM.

	Descriptors	Hidden states	Test set accuracy	Training set accuracy
	MFCC 2, 3, 4, 5, 6, 7, 12,
	Frequency, Energy,	29	0.61	0.91
	Periodicity, AC1
	MFCC 2, 3, 4, 5, 12,
	Frequency, Energy,	21	0.60	0.87
	Periodicity, AC1, Loudness
	MFCC 2, 3, 4, 5, 12,
	Frequency, Energy,	22	0.60	0.86
	Periodicity, AC1
	MFCC 2, 3, 4 5, 12, Frequency, Energy, AC1	22	0.60	0.86

https://www.youtube.com/watch?v=loy6JaW7L2c

These numbers correspond to the 1st and

3rd MFCCs, and to the first-order autocorrelation coefficient descriptor (AC1).

The library is open-source and available for download at: https://github.com/openmusic-project/om-xmm.

ACKNOWLEDGMENT

This work is carried out within the PEPS I3A support framework of the French National Center for Scientific Research (CNRS). The authors thank composer Alireza Farhang for providing the annotated gesture dataset.


David S Dean

Pierre Le Doussal

Satya N Majumdar

Gregory Schehr

Grégory Schehr

Wigner function of noninteracting trapped fermions

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

I. INTRODUCTION

A. Background

The Heisenberg uncertainty principle, the basic cornerstone of quantum mechanics, tells us that the position and the momentum of a single quantum particle cannot be measured simultaneously. In position space, the squared wave function |ψ(x)| 2 is the probability density. Similarly, | ψ(p)| 2 [where ψ(p) is the Fourier transform of ψ(x)] gives the probability density in momentum space. Although the joint probability density function (PDF) cannot be defined in phase space (x,p), the closest object to such a joint PDF is the celebrated single-particle "Wigner function" [1] W 1 (x,p) = 1 2π h +∞ -∞ dy e ipy/ hψ * x + y 2 ψ x -y 2 .

(

) 1
By integrating W 1 over p one recovers the spatial PDF, |ψ(x)| 2 , and similarly, by integrating W 1 over x, one recovers the momentum PDF, | ψ(p)| 2 . However, in general W 1 (x,p) need not be positive, and hence does not have the interpretation of a joint PDF. Nevertheless, the Wigner function has been useful in numerous contexts [2,3], including quantum chaos and semiclassical physics [4,5], in quantum optics [START_REF] Walls | Quantum Optics[END_REF], in the modeling of optical devices [3], and in quantum information [START_REF] Douce | [END_REF]. It has also been used in the context of many-body (fermionic and bosonic) systems. For one-dimensional (1D) fermions at zero temperature and in the presence of a confining potential, the Wigner function was first studied by Balazs and Zipfel [8]. Interestingly, the same results were derived recently by using a completely different method by Bettelheim and Wiegmann, in the context of nonequilibrium dynamics of a perturbed Fermi gas [9]. These results were found to be useful in the context of quantum mirror curves [10]. The Wigner function has also been used in the implementation of numerical methods to study fermions [11] and also in the study of Bose-Einstein condensates [12].

As the Wigner function contains information both about the positions and the momenta of particles, its experimental measurement can be useful to probe subtle questions concerning the nature of the quantum state of the system. However, measuring the Wigner function experimentally is a considerable technical challenge. In most cold atom experiments, only either the position or the momentum distributions can be measured, but not both simultaneously due to the uncertainty principle. In certain special cases, e.g., for a Fermi gas trapped in a two-dimensional box potential, measuring the momentum distribution can indirectly provide the measure of the Wigner function [13], but this is not generic. However, the development of quantum tomography has paved the way for measuring the full Wigner function (and not just its position and momentum marginals) in a number of recent contexts, e.g., in quantum optics [14], in trapped atom setups [15][16][17], and for electrons in ballistic conductors [18].

Given that most experiments are performed at finite temperature, it is also important to provide a finite temperature generalization of the T = 0 Wigner function introduced in (1). In fact, one of the main goals of this paper is to characterize this temperature dependence of the Wigner function for noninteracting fermions trapped in a confining potential in d dimensions.

Recently, there has been considerable interest in trapped Fermi gases, both theoretically [19] and in cold atom experiments [20]. Even in the noninteracting limit this system displays rich and universal quantum and thermal fluctuations, as was demonstrated recently [21][22][23][24][25][26][27][28]. The case of the harmonic trap played a fundamental role because it is solvable and makes an important connection between trapped noninteracting fermions and the eigenvalues of a random matrix. Indeed, in one dimension (d = 1) and at zero temperature T = 0, the positions of the fermions are in one-to-one correspondence with the eigenvalues of the Gaussian unitary ensemble (GUE) of random matrix theory (RMT) [22,25,27,28]. Consequently, at T = 0, the quantum fluctuations of N fermions, characterized by the squared many-body ground-state wave function, | 0 (x 1 , . . . ,x N )| 2 , was shown to be identical to the joint PDF of the eigenvalues of a GUE random matrix. Similarly, the joint distribution of the momenta is given by | ˆ 0 (p 1 , . . . ,p N )| 2 , where ˆ 0 is the N -variable Fourier transform of 0 . In the case of the harmonic trap, because of the symmetry x ↔ p (in scaled units), it is identical to the joint PDF of the positions, i.e., ˆ 0 = 0 .

The squared many-body wave function | 0 (x 1 , . . . ,x N )| 2 (in real space) or its Fourier counterpart | ˆ 0 (p 1 , . . . ,p N )| 2 (in momentum space) encodes information about quantum fluctuations. For instance, by integrating

| 0 (x 1 , . . . ,x N )| 2 [respectively | ˆ 0 (p 1 , . . . ,p N )| 2 ]
over N -1 positions (respectively momenta), one obtains the average density of fermions in real space (respectively in momentum space). In the N → +∞ limit, from the mapping to the GUE, it is given (in scaled units and normalized to unity) by the Wigner semicircle law of RMT, ρ W (y) = π -1 2y 2 , with |y| √ 2 and zero elsewhere. Near the soft edge y = √ 2 (and similarly around y = -√ 2), the density gets smeared over a width w N ∼ N -1/6 , which defines the edge regime, and the density profile is described by a nontrivial scaling function F 1 , known in RMT [29,30]. These results extend to all n-point correlation functions, either in position or momentum space. In particular, the scaled PDF of the position of the rightmost fermion, x max = max 1 i N x i , is given [23,25] by the Tracy-Widom (TW) distribution of the GUE [31]. Interestingly, the (scaled) largest fermion momentum p max = max 1 i N p i , measurable in time of flight experiments [START_REF] See | Ultra-cold Fermi Gases[END_REF], is also distributed with the same TW distribution. This analysis has been recently extended to any spatial dimension d [24,26], to finite temperature, and beyond the harmonic oscillator for more general smooth potentials [25].

It is thus natural to ask which of these universal properties extend to the Wigner function for N noninteracting trapped fermions, to gain insight on the quantum fluctuations in the phase space. The N body Hamiltonian is ĤN = N i=1 Ĥ (x i , pi ), where the single-particle Hamiltonian for spinless fermions of mass m is given by

Ĥ = Ĥ (x, p) = p2 2m + V (x), (2)
with V (x) = 1 2 mω 2 x 2 in the case of the harmonic oscillator. The many-body Wigner function is defined at T = 0 as a generalization for any N and d of (1)

W N (x,p) = N (2π h) d +∞ -∞ dy dx 2 . . . dx N e ip•y h × * 0 x + y 2 ,x 2 , . . . ,x N 0 x - y 2 ,x 2 , . . . ,x N , (3)
which by construction satisfies

+∞ -∞ dp W N (x,p) = ρ N (x), +∞ -∞ dx W N (x,p) = ρN (p), (4)
+∞ -∞ dx dp W N (x,p) = N,
where ρ N (x) is the average density of fermions (here normalized to N ) and ρN (p) its counterpart in momentum space.

B. Main results

In this paper, we compute W N (x,p) exactly in the large N limit, both in the bulk and at the edge of a noninteracting Fermi gas trapped by a confining potential V (x). We perform the derivation in arbitrary dimension d, first at T = 0 and for the harmonic oscillator and then at finite temperature and for a large class of smooth potentials. Our results generalize, to d > 1 and T > 0, the results previously obtained for d = 1 and T = 0 [8,9].

Zero temperature T = 0. The result in the bulk is particularly simple,

W N (x,p) 1 (2π h) d (μ -E(x,p)), (5)
where

E(x,p) = p 2 2m + V (x) ( 6 )
is the classical energy in the phase space. Here (x) is the Heaviside unit step function and μ is the Fermi energy which is related to N via the normalization (4). Note that Eq. ( 5) is valid for large N (equivalently large μ) and for an arbitrary potential V (x). This result, which can be obtained by semiclassical methods such as the local density approximation (see e.g. [START_REF] Castin | Ultra-cold Fermi Gases[END_REF]), is obtained here through a controlled asymptotic analysis of an exact formula. Clearly, the form of W N (x,p), given in Eq. ( 5), vanishes beyond the surface parameterized by (x e ,p e ) where

p 2 e 2m + V (x e ) = μ. ( 7
)
Following Ref. [9], we will call this surface the "Fermi surf;" it is the semiclassical version of the Fermi surface in classical phase space (see Fig. 1). Note that by integrating (5) over momentum as in (4), one recovers the result for the average number den-

sity ρ N (x) = (2 h) -d [p e (x)] d /γ d with γ d = π d/2 (1 + d/2),
where p e (x) := √ 2m(μ -V (x)) is the Fermi momentum, which sets the typical inverse interparticle spacing ∝ h(|p e (x)|) -1 . A similar result can be obtained for the average momentum density ρN (p) by integrating (5) over x. Both the position as well as the momentum densities exhibit marginal edges x em (respectively p em ) beyond which they vanish. The number density ρ N (x) vanishes at x = x em , where x em satisfies V (x em ) = μ. Similarly, the average momentum density ρN (p) vanishes at p = p em , where one can show that

|p em | = max x |p e (x)| = √ max x [2m(μ -V (x))
]. In the

where (x,p) is a point in the phase space close to the Fermi surf and e N is an energy scale given by

e N = ( h) 2/3 (2m) 1/3 1 m (p e • ∇) 2 V (x e ) + |∇V (x e )| 2 1/3 . ( 9
)
We then show that the Wigner function W N (x,p), at T = 0 and in arbitrary d, can be expressed as a universal function of the dimensionless variable a as

W N (x,p) W(a) (2π h) d , ( 10
)
where the scaling function

W(a) = +∞ 2 2/3 a
Ai(u)du (11) is independent of the space dimension d. In Eq. ( 11), Ai(u) is the Airy function. The function W(a) has the asymptotic behaviors

W(a)∼ (8π ) -1/2 a -3/4 exp -4 3 a 3/2 , a → +∞, 1, a → -∞. (12)
In particular, the limit lim a→-∞ W(a) = 1 ensures a smooth matching with the bulk result (5). In the inset of Fig. 2 we show a plot of this function W(a). Note that in d = 1 our results coincide exactly with the one obtained using semiclassical methods [8,9]. Our results here provide a generalization of the d = 1 result to arbitrary d. Interestingly, exactly the same scaling function as in Eq. ( 11) also appeared recently in the context of electron quantum optics [START_REF] Ferraro | [END_REF]. Note that for the case of harmonic oscillator, where V (x) = (1/2)m ω 2 x 2 , the energy scale e N in (9) reduces to where r e = 2μ/(mω 2 ) and

e N = mω 2 r e w N , (13)
w N = 1 α √ 2 (μ/ hω) -1/6 , with α = mω/ h, (14)
represents the width of the edge region in the real space [24,25]. Hence, the energy scale e N for the harmonic oscillator reads

e N = ( hω) 2/3 μ 1/3 . ( 15
)
Furthermore, for the harmonic oscillator, the Fermi energy μ is related to N , for large N, via [24,25] μ ∼ hω(N (d

+ 1)) 1/d . ( 16
)
Consequently, the argument a of the scaling function W(a) in Eq. ( 8) reduces, in this case, to

a = 1 w N p 2 m 2 ω 2 + x 2 -r e . ( 17
)
Finite temperature T > 0. Next, we generalize our T = 0 results for the Wigner function to finite temperature T . As in the T = 0 case, there are two regimes, namely the bulk and the edge. The sharp bulk behavior at T = 0 in Eq. ( 5) is smeared out by thermal fluctuations at finite T and is replaced by

W μ(x,p) = 1 1 + e β p 2 2m +V (x)-μ , ( 18
)
where β = 1/T . The finite temperature chemical potential μ, in the canonical ensemble, can be determined as a function of β and N via the Fermi relation

N = k n k = k 1 1 + e β( k -μ) , ( 19
)
where the k denote the single-particle energy levels of the Hamiltonian Ĥ in (2). In the limit T → 0, μ → μ from Eq. ( 19) and Eq. ( 18) reduces to the T = 0 result in (5). This semiclassical finite-temperature bulk result in Eq. ( 18) was also derived by other methods [35].

Near the finite-temperature edge, where E(x,p) → μ, we show that the Wigner function has a universal scaling behavior for large N . This universal behavior emerges when the temperature T scales with N (or equivalently with μ) in a particular fashion, namely when temperature T ∼ e N , where e N is the energy scale defined in Eq. ( 9). Note that e N just depends on the Fermi energy μ, but not on the temperature. Hence, we set

β e N = b, ( 20
)
where the dimensionless parameter b = O( 1) is kept fixed in the limit of large N . For instance, for the harmonic oscillator in 1D, using Eq. ( 15) and μ ∼ hω N , one gets b = ( hω/T )N 1/3 . This is the same temperature scale that appears in the analysis of the spatial correlations near the edge in real space [25]. In this temperature regime, one can show that the finite-temperature chemical potential μ ∼ μ, indicating that the finitetemperature edge is the same as the zero-temperature edge. Hence, as in the T = 0 case (8), we consider the same dimensionless variable a. We show that, in this temperature regime characterized by the single dimensionless parameter b (20), the Wigner function takes a scaling form

W μ(x,p) ∼ W b (a) (2π h) d , a = 1 e N p 2 2m + V (x) -μ , ( 21
)
with e N given in Eq. ( 9). The scaling function W b (a) is given by

W b (a) = ∞ -∞ dy Ai(y) 1 + e a b e -b y2 -2/3 . ( 22
)
In Fig. 2

W b (a) ∼ e b 3 /12 e -a b , a → +∞, 1, a→ -∞. (23)
Note that, for any finite b, the right tail of W b (a), as a → ∞, decays exponentially with a. It is only exactly at T = 0, i.e., when b → ∞ limit, that the right tail decays faster than exponentially as in Eq. ( 12). Finally, we note that, as in the T = 0 case, this edge scaling function W b (a) is completely universal, i.e., independent of the dimension d as well as the confining potential V (x) as long as the potential is nonsingular. The rest of the paper is organized as follows. In Sec. II, we compute the Wigner function at zero temperature. Section II A contains the exact solution (for any finite N ) for the 1D harmonic oscillator, Sec. II B discusses the d-dimensional harmonic oscillator, while in Sec. II C we generalize these results for arbitrary smooth confining potentials. In Sec. III, these results are generalized to finite temperature T > 0. Finally, we conclude with a summary and discussion in Sec. IV. Some details are relegated to the Appendixes.

II. WIGNER FUNCTION AT ZERO TEMPERATURE

At T = 0, the quantum correlation functions of noninteracting fermions can be written as determinants constructed from a central object, the so-called kernel (see e.g. [25])

K μ (x,x ) = k (μ -k )ψ * k (x)ψ k (x ) (24)
in terms of the single-particle eigenfunctions ψ k (x) of ( 2) and their associated eigenenergies k , labeled by quantum numbers k. In (24) μ is chosen so that the sum contains exactly N levels.

It turns out that one can relate the Wigner function (3) to the kernel, using the following formula: 25) which follows from the property that 0 is a Slater determinant constructed from the ψ k 's (see the derivation in Appendix A). Comparing ( 25) and ( 3) we obtain

K μ (x,x ) = N +∞ -∞ dx 2 . . . dx N × * 0 (x,x 2 , . . . ,x N ) 0 (x ,x 2 , . . . ,x N ), (
W N (x,p) = 1 (2π h) d +∞ -∞ dy e ip•y h K μ x + y 2 ,x - y 2 .
(26) The scaling behavior of this kernel K μ has been well studied in the large μ and N limit [25]. One can then use these results in (26) to obtain information about the Wigner function as shown below.

A. Calculation in d = 1 for the harmonic oscillator

Let us first present an exact calculation for the d = 1 harmonic oscillator, using space, momentum, time, and energy dimensionless units, i.e., in units of

x 0 = 1/α, p 0 = hα, t 0 = 1/ω, e 0 = hω. ( 27
)
In these scaled units μ N for large N . The kernel reads

K μ (x,x ) = N-1 k=0 ψ k (x)ψ k (x ), (28)
where

ψ k (x) = ( 1 √ π 2 k k! ) 1 2 H k (x)e -1 2
x 2 and H k is the kth Hermite polynomial. Plugging it in Eq. ( 26), and specifying d = 1, we obtain in dimensionless units and at T = 0

W N (x,p) = 1 2π N-1 k=0 ∞ -∞ dy e ipy ψ k x + y 2 ψ k x - y 2 .
(29) Remarkably, the integral over y in Eq. ( 29) can be performed explicitly using an identity first derived by Groenewold [36]

+∞ -∞ dy e ipy ψ k x + y 2 ψ k x - y 2 = 2(-1) k L k (2(x 2 + p 2 ))e -x 2 -p 2 , ( 30
)
where L k (y) is the Laguerre polynomial of degree k, defined via its generating function

∞ k=0 z k L k (y) = 1 (1 -z) e -z y/(1-z) . ( 31
)
Substituting (30) in Eq. ( 29) gives the explicit result We have plotted W N (x,p) as a function of r = x 2 + p 2 for N = 200 in Fig. 3. From this figure one sees that there are three distinct regimes at large N: (i) the bulk regime where W N (x,p) oscillates around the bulk value 1/(2π ) = 0.159155 . . ., (ii) the edge regime around r = √ 2N where the Wigner function vanishes over a width of order N -1/6 (see below), and (iii) an "anomalous" regime near r = 0. This anomalous regime has been pointed out [4] for the single-particle case (N = 1); here we show that it persists for multiparticle systems. At r = 0 the Wigner function vanishes exactly for N even, and equals 1/π for N odd. At small r ∼ 1/ √ N there is a scaling regime describing the Wigner function near r = 0. We now study these three regimes in detail.

W N (x,p) = 1 π N-1 k=0 (-1) k L k (2(x 2 + p 2 )) e -x 2 -p 2 . ( 32
)
We start by multiplying Eq. ( 32) by z N and sum over N . We obtain

∞ N=1 z N W N (x,p) = 1 π ∞ N=1 z N N-1 k=0 (-1) k L k (2r 2 )e -r 2 , ( 33
)
with r 2 = x 2 + p 2 . To perform this double sum, we write z N = z N-k z k and perform the sums separately over k and m = Nk. This gives

∞ N=1 z N W N (x,p) = 1 π ∞ m=1 z m ∞ k=0 (-z) k L k (2r 2 )e -r 2 . ( 34
)
Using the generating function of Laguerre polynomials in Eq. ( 31), we get

∞ N=1 z N W N (x,p) = 1 π z 1 -z 2 e -1-z 1+z (x 2 +p 2 ) . ( 35
)
For convenience, we use z = e -s and get

W (x,p; s) := ∞ N=1 W N (x,p)e -Ns = e -s-(x 2 +p 2 ) tanh s 2 π(1 -e -2s ) . ( 36
)
Bulk behavior. To extract the bulk result we consider the limit s 1 with s(x 2 + p 2 ) fixed, leading to

W (x,p; s) e -s 2 (x 2 +p 2 ) 2πs , ( 37
)
which yields, after Laplace inversion,

W N (x,p) 1 2π (2N -(x 2 + p 2 )), (38)
in agreement with the general result (5). Edge behavior. We now show that the exact formula (35) can be used to derive the edge scaling function (10), (11) in d = 1 for the harmonic oscillator. Our starting point is the exact relation (35). Inverting the generating function using Cauchy's inversion formula,

W N (x,p) = 1 π c+i×∞ c-i×∞ ds 2iπ e s (N-1) (1 -e -2s ) exp -(x 2 + p 2 ) tanh s 2 , ( 39
)
where c is to the right of all singularities in the complex s plane. For large N , the most important contributions come from the vicinity of s → 0. We set, near the edge

x 2 + p 2 = √ 2N , x 2 + p 2 = √ 2N + 1 √ 2 N -1/6 a, ( 40
)
where a denotes the distance (on the scale of N -1/6 ) from the edge; see (8) in d = 1 and dimensionless units. Using tanh(s/2) = s/2s 3 /24 as s → 0, we find that the integral in Eq. ( 40) reduces for large N to

W N (x,p) ≈ 1 2π c+i×∞ c-i×∞ ds 2iπs exp -s N 1/3 a + s 3 N 12 .
(41)

Rescaling further by setting t = s N 1/3 , we finally get, near the edge,

W N (x,p) ≈ 1 2 π W(a), (42)
where

a = √ 2 N 1/6 ( x 2 + p 2 - √ 2N ), (43)
and the scaling function W(a) is given exactly by

W(a) = c+i×∞ c-i×∞ dt 2iπt exp -t a + t 3 12 . ( 44
)
Using the integral representation of the Airy function

Ai(x) = c+i×∞ c-i×∞ dτ 2iπ exp -τ x + τ 3 3 , ( 45
)
we find that this scaling function W(a) is indeed given by formula (11). One can easily check that it has the asymptotic behaviors given in Eq. (12). Anomalous behavior near r = x 2 + p 2 = 0. It is easy to see that the Wigner function vanishes at r = 0 for N even.

From the definition (29),

W N (x = 0,p = 0) = 1 2π ∞ -∞ dy N-1 k=0 ψ * k y 2 ψ k - y 2 = 1 π N-1 k=0 (-1) k = 1 π , N odd, 0, N even, (46)
where we used (i) the orthonormality of the single-particle wave functions ψ k (x) and (ii) the fact that for the harmonicoscillator potential ψ k (-x) = (-1) k ψ k (x). Note the property (46) extends to an arbitrary even potential

V (x) = V (-x) in d = 1.
As we show in Appendix D, near r = 0 in a regime where r ∼ 1/ √ N , the Wigner function has the following scaling behavior for large N:

W N (x,p) ∼ 1 2π -(-1) N F ( N (x 2 + p 2 ))
,

F (z) = 1 2π J 0 (2 √ 2z), (47)
where J ν (x) is the Bessel function with index ν. Interestingly, this parity dependence in Eq. ( 47) persists even for large N .

We have verified the scaling behavior in (47) by numerically evaluating W N (x,p) in Eq. ( 32).

B. d-dimensional harmonic oscillator

In dimension d > 1 it is more convenient to use a different method using the quantum propagator to calculate the kernel, and in turn the Wigner function via (26). In addition, as we show later this method is more versatile as one can treat more general potentials and demonstrate universal properties of the Wigner function. The method relies on the following representation of the kernel in d dimensions for arbitrary single-particle Hamiltonian [25]

K μ (x,x ) = C dt 2πit e μt/ hG(x,x ; t), ( 48
)
where C is the Bromwich contour in the complex plane and G(x,x ; t) = x |e -Ĥ t/h |x is the one-particle Euclidean quantum propagator associated to the Hamiltonian Ĥ in (2).

Let us apply this relation to the case of the harmonic oscillator in dimension d for which the exact propagator is known [START_REF] Feynman | Quantum Mechanics and Path Integrals[END_REF].

We work again here in the aforementioned dimensionless units in which the propagator reads

G(x,x ; t) = 1 (2π sinh t) d/2 × exp - (x-x ) 2 +(x 2 +(x ) 2 )(cosh t-1) 2 sinh t . ( 49
)
We first insert (49) into (48) and then use Eq. (26). Performing the Gaussian integration over y we obtain

W N (x,p) = 1 (2π ) d C dt 2πit e μt 1 cosh t 2 d e -(x 2 +p 2 ) tanh t 2 ,
(50) an exact formula. Note that the N dependence of W N (x,p) is only through μ. As in the d = 1 case we now analyze this formula both in the bulk as well as in the edge regime.

Bulk behavior. In the bulk one can show that the values of t which dominate the integral are O(1/μ). Hence, we need to expand the factor tanh(t/2) inside the exponential only to O(t). This leads to

W N (x,p) = 1 (2π ) d C dt 2πit e (μ-p 2 +x 2 2 )t . ( 51
)
The integral over t just gives a Heaviside θ function, hence establishing the result for the bulk in Eq. ( 5) in the case of the harmonic potential. Edge behavior. To analyze the edge behavior in the large μ limit, we need to expand the exponential up to order t 3 , as in [25] for the study of the edge in real space. We obtain [discarding terms of O(t 4 ) in the exponential] Let us now define t = t N τ , with t N = 2 2/3 μ -1/3 . Discarding the term ∝ t 2 which is subdominant, this leads to

W N (x,p) = C dt (2π ) d+1 it e μ-x
S = -a2 2/3 τ + τ 3 3 + O(μ -2/3 ). ( 54
)
One thus obtains, upon restoring the physical units, the scaling form (10) for the Wigner function at large μ as

W N (x,p) W(a) (2π ) d , W(a) = C dτ 2πiτ e -a2 2/3 τ + τ 3 3 , ( 55
)
which is precisely the integral representation given in Eq. ( 44) with the substitution 2 2/3 τ = t. Thus, we obtain the remarkable result that the edge scaling form W(a) of the Wigner function for the harmonic oscillator is completely independent of the space dimension d. The same result can also be obtained (see Appendix B), directly from the known scaling behavior of the kernel [25]. A natural question is whether this scaling form is also universal with respect to the details of the shape of the confining potential, as we will discuss below.

In d = 1, we have seen in the previous subsection that there is an additional anomalous regime when x 2 + p 2 = O(1/N ). In d > 1, similar anomalous regimes are likely to exist, though we have not investigated them in detail.

C. Wigner function for other smooth confining potentials:

Beyond the harmonic oscillator

The case of the harmonic oscillator, treated in the previous section, is special because x and p [appropriately rescaled in the dimensionless units (27)] play symmetric roles, and the Wigner function depends only on the single variable x 2 + p 2 .

For a general potential V (x) this is no longer the case and a different treatment is needed, as we now show.

Consider a more general smooth potential V (x) in d dimensions. In this case we show, using the above propagator method, that in the bulk the Wigner function is given by (5). Putting together formulas (48) and ( 26) one can express the Wigner function directly in terms of the Euclidean propagator

W N (x,p) = 1 (2π h) d C dt 2πit e μt h +∞ -∞ dy e ip•y h G x + y 2 ,x - y 2 ,t . ( 56
)
Let us recall that the Euclidean quantum propagator satisfies the Feynman-Kac equation G(x,x ,t)

∂ t G = -Ĥ G = h2 2m ∇ 2 x -V (x) G, ( 57
)
with G(x,x ,0) = δ d (xx ). As was shown in [25], and further detailed in Appendix C, the large N limit can be obtained from the small time t expansion of the quantum propagator, which up to O(t 3 ) reads

G(x,x ,t) = m 2π ht d/2 exp - m 2 ht (x -x ) 2 × exp - t h S 1 (x,x ) - t 2 2m S 2 (x,x ) + t 3 2m h S 3 (x,x ) , ( 58
)
where S 1 ,S 2 ,S 3 for an arbitrary potential V (x) are given explicitly in Eqs. ( 243)-(245) of [25].

To analyze the bulk behavior of the Wigner function we only need the leading O(t) term S 1 , which then reads

S 1 x - y 2 ,x + y 2 = 1 0 du V x + u - 1 2 y . ( 59
)
We substitute this expression of S 1 in (58) and keep only up to O(t) terms. Next we substitute this propagator in (56)

W N (x,p) = 1 (2π h) d C dt 2πit m 2π ht d/2 -∞ dy e ip•y h × exp - m 2 ht y 2 + t h μ- 1 0 du V x + u - 1 2 y . ( 60
)
One can show that the values of t which dominate the integral are O(1/μ). Hence, by rescaling y as shown in Appendix C, one can neglect the term (u -1/2)y in the argument of V , leading, after integration over y, to

W N (x,p) = 1 (2π h) d C dt 2πit e (μ-p 2 2m -V (x))t . ( 61
)
The integral over t just gives a Heaviside θ function, hence establishing the result for the bulk in Eq. ( 5). From Eq. ( 5) it is clear that the Wigner function vanishes beyond the boundary of a bounded support, which defines an edge in phase space, i.e., a surface parameterized by (x e ,p e ), which satisfies the equation

p 2 e 2m + V (x e ) = μ. ( 62
)
As for the case of the harmonic oscillator, for large but finite N the jump of the Wigner function described by ( 5) is smoothed over a scale w N which now explicitly depends on the potential. However, the appropriately centered and scaled Wigner function at any point of the edge surface is again universal and is given by W(a) in (11). More precisely the Wigner function takes the following edge scaling form for

p 2 2m + V (x) -μ ∼ μ 1/3 μ: W N (x,p) W(a) (2π h) d , ( 63
)
where the scaled variable a = O( 1) is now naturally expressed as the ratio of two energies

a = 1 e N p 2 2m + V (x) -μ , ( 64
)
where the energy scale e N is given:

e N = ( h) 2/3 (2m) 1/3 1 m (p e • ∇) 2 V (x e ) + |∇V (x e )| 2 1/3 . ( 65
)
For the derivation of these results, see Appendix C. For the harmonic oscillator V (x) = 1 2 mω 2 x 2 , one finds

e N = hω μ hω 1/3 = mω 2 r e w N = hω 2α 2 w 2 N . ( 66
)
We thus see that Eq. ( 64) is consistent with Eq. ( 8) and Eq. ( 65) is consistent with (13), as discussed earlier. For a general potential the above results are valid on any point (x e ,p e ) on the edge surface. In particular for p e = 0 we have

a = 1 e N (V (x) -μ) 1 e N ∇V (x e ) • (x -x e ), ( 67
) e N = ( h) 2/3 (2m) 1/3 |∇V (x e )| 2/3 . ( 68
)
Furthermore, to make contact with the 1D result derived in Ref. [9], let us consider a point in the (x,p) plane near the edge (x e ,p e ), where we set x = x e and p = p e + p with |p| 1 (see Fig. 1). In particular, if we focus on d = 1, the formulas (64) and (65) simplify a lot. In this case the formula (65) becomes

e N = ( h) 2/3 (2m) 1/3 p 2 e m V (x e ) + (V (x e )) 2 , ( 69
)
and Eq. ( 64) reduces, to leading order in p,

a = p e (x e ) e N p. (70)
Here, p e (x e ) is the point on the edge surface parameterized by the function

p 2 e (x) = 2m(μ -V (x)). ( 71
)
In Ref. [9], this p e (x) was called the "Fermi surf." By taking twice the derivative of (71), the expression for e N in Eq. ( 69)

simplifies to

e N = h2 2 1/3 |p e (x e )| 1/3 p e (x e ) m . ( 72
)
Consequently, from Eq. ( 70), one gets

a = κ p, κ = h2 2 p e (x e ) -1/3 . ( 73
)
Thus, the argument a of the scaling function W(a) in Eq. ( 63) reduces, in d = 1, to precisely the argument derived in Ref. [9].

III. WIGNER FUNCTION AT FINITE TEMPERATURE

We now extend our analysis to finite temperature. The Wigner function at temperature T = 1/β in the canonical ensemble can be defined from the many-body density matrix DN = e -β ĤN /Z N (β) as

W N,T (x,p) = N (2π h) d +∞ -∞ dy dx 2 . . . dx N ,e ip•y h × x + y 2 ,x 2 , . . . ,x N DN x - y 2 ,x 2 , . . . ,x N . ( 74
)
By decomposing on the basis of eigenvectors |E of ĤN one can equivalently write it as

W (x,p) = 1 Z N (β) E e -βE W N,E (x,p), ( 75
)
where Z N (β) = E e -βE is the canonical partition sum and

W N,E (x,p) = 1 (2π h) d +∞ -∞ dy e ip•y h K x + y 2 ,x - y 2 ; {n k } ,
(76) where the N-body eigenstate |E is labeled by a set of occupation numbers n k = 0,1 of the single-particle eigenstates, such that E = k n k k and N = k n k . The kernel K(x,x ; {n k }) has the expression

K(x,x ; {n k }) = k n k ψ * k (x)ψ k (x ), ( 77
)
where ψ k (x) is the single-particle eigenstate labeled by k. Note that in the T = 0 limit (77) reduces to ( 28) and (76) reduces to (26). One can also define the Wigner function in the grandcanonical ensemble with chemical potential μ as

W μ(x,p) = 1 Z gr (β, μ) N Z N (β)W N,T (x,p)e μβN , ( 78
)
where Z gr (β, μ) = N Z N (β)e μβN is the grand-canonical partition function. Substituting (75) in (78) and summing Eq. ( 76) over the eigenstates |E and N , i.e., over the independent variables n k 's, leads to

W μ(x,p) = 1 (2π h) d +∞ -∞ dy e ip•y h K μ x + y 2 ,x - y 2 , ( 79
)
where K μ(x,x ) is the kernel defined in the grand-canonical ensemble

K μ(x,x ) = k n k ψ * k (x)ψ k (x ), (80)
where

n k = 1 1 + e β( k -μ) (81)
is the mean occupation number of state k, over the Fermi distribution.

In the grand-canonical ensemble one can relate the finitetemperature Wigner function W μ(x,p) to the zero-temperature Wigner function W N (x,p). To see this, we consider first the zero-temperature kernel in (24). Taking a derivative with respect to the Fermi energy μ gives

∂ μ K μ (x,x ) = k δ(μ -k )ψ * k (x)ψ k (x ). (82)
Now we start with Eqs. ( 80) and ( 81) and rewrite them as

K μ(x,x ) = +∞ -∞ dμ 1 + e β(μ -μ) k δ(μ -k )ψ * k (x)ψ k (x ) = +∞ -∞ dμ 1 + e β(μ -μ) ∂ μ K μ (x,x ), (83)
where in the last line we used (82). This equation was previously derived as Eq. ( 240) in [25]. We can now apply this formula to the Wigner function using (79), leading to

W μ(x,p) = +∞ -∞ dμ 1 1 + e β(μ -μ) ∂ μ (W N (x,p)| μ=μ ), (84)
which is exact in the grand-canonical ensemble. We have denoted the integration variable by μ to avoid confusion with the variable μ which denotes the Fermi energy. Note that the zero-temperature Wigner function W N (x,p) depends implicitly on μ. This relation (84) will now allow us to derive the bulk and edge properties of the Wigner function at finite temperature.

The above results are obtained in the grand-canonical ensemble where μ is a given parameter. However, our main goal is to describe the Wigner function in the canonical ensemble where N is fixed. Indeed, in the large N limit, adapting the saddle-point method of [25] to the Wigner function, we expect that one can use the grand-canonical results for the canonical ensemble, provided we determine μ as a function of N by the relation

N = k n k = k 1 1 + e β( k -μ) . ( 85
)
Note that by definition of μ one also has

N = k n k = k (μ -k ), (86)
which is also the T = 0 limit of (85). This implies that μ is related to μ by equating the two relations (85) and (86). Bulk behavior. Substituting the result (5) for the bulk zerotemperature Wigner function in (84) we obtain

W μ(x,p) = 1 1 + e β( p 2 2m +V (x)-μ) . ( 87
)
This equation is valid for all x,p in the phase space where p 2 2m + V (x)μ = O(T ). In the canonical ensemble μ is related to N via (85), and using also (86), one sees that the bulk regime corresponds to scaling T ∼ μ and T ∼ μ. By integrating (87) over p (respectively over x) one obtains the finite-temperature bulk density ρ N (x) [respectively momentum density ρN (p)] [see, e.g., Eqs. ( 270)-(273) in [25]].

Edge behavior. For simplicity let us first focus on the harmonic oscillator. There it turns out that, in the edge regime, one needs to scale the temperature as T ∼ μ 1/3 in the limit of large μ. Consequently, one defines a reduced inverse temperature

b = βμ 1/3 , ( 88
)
with b = O(1) in this regime. In addition, in this regime μ μ [25]; hence we set μ = μ in the following. Our starting point is the integral in (84). For T ∼ μ 1/3 , this integral is dominated by the regime where μ is close to μ. In fact, by setting

β(μ -μ) = -bu, (89)
we see that the Fermi factor in (84) takes the dimensionless form 1/(1 + e -bu ). This suggests that the integral will be controlled by values of u which are of order unity; hence

μ -μ = O(μ -1/3 ). (90)
Therefore, at any point x,p in phase space close to the edge defined in (62), we can use the scaling form ( 10) and ( 11) of W N (x,p) inside the integral in (84). We note that the only dependence on μ of W(a) is through the scaling variable

a = a μ := √ 2 (μ ) 1/6 p 2 + x 2 -2μ . ( 91
)
Using the scaling form we have

∂ μ W N (x,p) W (a μ )∂ μ a μ . ( 92
)
We can now expand a μ in (91) around μ = μ as

a μ = a μ - 1 μ 1/3 (μ -μ). (93)
Using (89) and b = βμ 1/3 we obtain

a μ = a μ + u + O(μ -2/3 ). ( 94
)
Inside the integral in (84) the factor

∂ μ a μ = du dμ ∂ u a μ = du dμ (1 + O(μ -2/3 )). ( 95
)
We now rewrite the integral in (84) in terms of the u variable using (94) in the argument of W in the right-hand side (RHS) of (92). This gives us the finite-temperature Wigner function, which near the edge takes the scaling form

W N (x,p) W b (a) (2π h) d , a = 1 w N p 2 m 2 ω 2 + x 2 -r e , ( 96
)
where w N is given in (14). The finite-temperature scaling function, parametrized by b = βμ 1/3 , is (98)

W b (a) = +∞ -∞ 2 2/3 du 1 + e -bu
For a → -∞, this gives

W b (a) → ∞ -∞ dy Ai(y) = 1. (99)
In contrast, when a → +∞, we get to leading order

W b (a) ∼ C e -a b , ( 100
)
where the prefactor

C = ∞ -∞ dy Ai(y) e b y 2 -2/3 = e b 3 /12 . ( 101
)
As in the zero-temperature case discussed in Sec. II C, we can extend the above finite-temperature results to the case of arbitrary smooth potentials. As discussed in Sec. II C the zerotemperature scaling form W(a) of the Wigner function near the edge is identical to that of the harmonic oscillator. In contrast the scaling variable takes the form

a = a μ = 1 e N (μ ) p 2 2m + V (x) -μ , (102)
where the energy scale e N (μ ) is an implicit function of μ which can be obtained from formula (65). It thus depends nonuniversally on the shape of the potential V (x) and on the precise location of x e ,p e on the edge surface in the phase space. Following the same steps as for the harmonic-oscillator case, the reduced inverse temperature variable is now

b = β e N (μ). ( 103
)
Defining u as in (89) we can again expand a μ in (102) around μ = μ as

a μ a μ - 1 e N (μ) + a μ ∂ μ ln e N (μ) (μ -μ), (104)
which leads exactly to Eq. (94) to leading order in μ. The rest of the argument simply goes through, and we find that, for arbitrary smooth potentials, the finite-temperature Wigner function near the edge takes the form

W N (x,p) W b (a) (2π h) d , a = 1 e N p 2 2m + V (x) -μ , ( 105
)
with exactly the same scaling function W b (a) given in (96).

IV. CONCLUSION

In this paper we have studied the Wigner function W N (x,p) for N noninteracting fermions in a confining trap in d dimensions. At zero temperature and large N , we have shown that there are two main regimes for W N (x,p) in the phase space. A bulk regime where the Wigner function is flat over a finite support and vanishes outside. The edge of this support in the phase space is given by p 2 2m + V (x) = μ, where μ is the Fermi energy. Around this edge W N (x,p), appropriately centered and scaled, is described by a scaling function W. We have shown that this scaling function is universal, i.e., the same for a large class of confining smooth potentials, and, strikingly, is independent of the space dimension d. We then extended these results to finite temperature and found a one-parameter edge scaling function W b where b is the scaled temperature. The finite-temperature scaling function is also universal and independent of the space dimension d. In both the bulk and the edge regimes, the scaling functions are non-negative everywhere.

In addition to these universal features (bulk and edge) there appears to be anomalous small scale regimes around special points in the phase space where the Wigner function is rapidly varying. For the harmonic oscillator in d = 1 we have analyzed such an anomalous regime [4] near x = p = 0 in detail. Indeed, it is well known, already at the single-particle level, that the Wigner function is not guaranteed to be non-negative for arbitrary potentials V (x). However, for the scaling regimes at large N (bulk and edge), our results show that the Wigner function remains positive, for a large class of smooth potentials.

Finally, in this paper, we have focused on the standard Wigner function W N (x,p) of an N-body system, as defined in Eq. ( 4). This can be interpreted as the one-point probability density function in the phase space (in the semiclassical sense). Naturally, one can also investigate higher-order correlation functions in the phase space. At T = 0, this can be naturally done by introducing a generalized Wigner function [1]

W (N) N (x 1 ,p 1 , . . . ,x N ,p N ) = 1 (2π h) N d dy 1 • • • dy N e i h N i=1 p i •y i × * 0 x 1 + y 1 2 , . . . ,x N + y N 2 × 0 x 1 - y 1 2 , . . . ,x N - y N 2 , ( 106
)
where 0 is the ground-state many-body wave function.

From this generalized Wigner function, one can construct successively n-point correlation functions by integrating out Nn phase-space coordinates as follows:

C (N) n (x 1 ,p 1 , . . . ,x n ,p n ) = N ! (N -n)! dx n+1 dp n+1 × • • • dx N dp N W (N) N (x 1 ,p 1 , . . . ,x N ,p N ). (107)
For instance, for n = 1, Eq. (107) reduces precisely to the standard Wigner function W N (x,p) defined in Eq. ( 4), i.e., C (N) 1 (x,p) = W N (x,p). Investigations of these higher-order correlation functions with n > 1, both at T = 0 and T >0, would be interesting and will be studied in a future publication [START_REF] Dean | [END_REF].

where S(x,p,t) = ln exp -t

1 0 du V x + (itp + ỹ√ t) u - 1 2 + √ tB u -V (x) B,ỹ , (C9)
where the "average" over ỹ is over a unit Gaussian random variable, uncorrelated with B, i.e.,

. . . B,ỹ = +∞ -∞ d ỹ (2π ) d/2 exp - ỹ2 2 . . . B . (C10)
Since the averaging measure is even in ỹ and even in B u it is clear that S(x,p,t) starts at O(t 2 ). Hence the form (C8) is quite convenient to study the short-time expansion. The leading term, O(t) in the exponential, is obtained by setting S(x,p,t) to zero, which recovers the result Eq. ( 5) of the text for the Wigner function in the bulk. Since ỹ, B, and p do not depend on t, the t dependence of S(x,p,t) in (C9) is explicit, and its expansion in powers of t at small time is straightforward, although tedious. It is done by a gradient expansion of the argument of V around x x e . In doing so we also need to check that this gradient expansion is consistent with the expansion at large μ, i.e., that in the large μ limit all terms in the argument of V are small compared to x e . This is clearly the case for the terms proportional to ỹ and B u , for which the gradient expansion is an expansion in the parameter t 1/2 N /x e ∼ μ -1/3 x -2/3 e 1, using our above anticipated estimate for t N . For the term itp, the gradient expansion parameter is t N p e /x e ∼ μ -1/6 x -1/3 e 1. We now calculate S(x,p,t) up to O(t 3 ). For this it is easy to see that we need only the first two cumulants in (C9). The first cumulant is

S 1 (x,p,t) = -t 1 0 du V x + (itp + ỹ√ t) u - 1 2 + √ tB u B,ỹ -V (x) (C11) = -t 2 1 24 ỹj ỹk + 1 2 1 0 du B ju B ku B ∇ j ∇ k V (x) -t 3 - 1 8 1 0 du(1 -2u) 2 (p • ∇) 2 V (x) + 1 1920 ỹj ỹk ỹ ỹn ỹ + 1 16 ỹj ỹk ỹ 1 0 du(1 -2u) 2 B u B nu B + 1 24 1 0 du B ju B ku B u B nu B ∇ j ∇ k ∇ ∇ n V (x) + O(t 4 ). (C12)
Note that all terms proportional to an odd number of gradients vanish due to the symmetry u → 1u in the integral over u. Performing all the averages we obtain

S 1 (x,p,t) = -t 2 1 24 + 1 12 ∇ 2 V (x) -t 3 - 1 24 (p • ∇) 2 V (x) + 1 640 + 1 480 + 1 240 ∇ 2 ∇ 2 V (x) , ( C13
)
where we have purposely indicated the contribution of each term, and checked that the last term is identical to the result in Eqs. ( 246)-(248) in [25]. The second cumulant expanded to O(t 3 ) reads (C18) where in the cubic term it is consistent to replace x by x e . Performing the change of variable t = 2 2/3 τ/e N , with e N given by (C2), we obtain exactly the integral representation (55) of the function W(a), hence demonstrating the result (C3) with the scaling variable a defined in (C2). Furthermore, this confirms that t N ∼ 1/e N is the time scale which dominates the integral in the edge regime as anticipated above.

S 2 (x,p,t) = t 2 2 1 0 du 1 0 du (C14) × V x + (itp + ỹ√ t) u - 1 2 + √ tB u -V (x) V x + (itp + ỹ√ t) u - 1 2 + √ tB u -V (x) c B,ỹ (C15) = t 3 2 1 0 du 1 0 du u - 1 2 u - 1 2 ỹj ỹk + B ju B ku ∇ j V (x)∇ k V (x) + O(t 4 ) (C16) = t 3 24 |∇V (x)| 2 + O(t 4 ). ( C17
Estimating the term O(t 2 ) in (C13) to be of order t 2 N V (x e )/x 2 e ∼ μ -1/3 x -2/3 e 1, we see that it is indeed negligible, as was already the case for p e = 0 in [25]. Inside the t 3 term in (C13) we see that the second term is smaller than the first by a factor 1 p e x 2 e ∼ μ -1/2 x -2 e 1. The examination of terms O(t 4 ) and higher is very tedious and can be performed along the lines of Appendix A in [25]. We will not reproduce this analysis here.

In summary, the above shows that for a large class of smooth potentials, at a generic point of the edge surface in phase space, the universal edge form of the Wigner function holds. The analysis is a rather simple extension of the one in [25] (simple in the sense that the characteristic scales are not changed, apart from some prefactors). Note that a necessary condition is that e N does not vanish, which is true at a generic point, but could fail in some exceptional cases, for instance if p e and ∇V (x e ) vanish simultaneously. In this appendix, we show that for the 1D harmonic oscillator at T = 0 the Wigner function near (x = 0,p = 0) in the phase space has an anomalous behavior. Indeed, for r 2 = x 2 + p 2 = O(1/N ), we will show that W N (x,p) has the following behavior:

W N (x,p) ∼ 1 2π
-(-1) N F ( N(x 2 + p 2 )),

F (z) = 1 2π J 0 (2 √ 2z), (D1)
where J ν (x) is the Bessel function of index ν. Our starting point is the exact generating function in Eq. (35). Formally inverting this generating function using Cauchy's formula, we find

W N (x,p) = 1 2πi C 0 dz z N 1 π(1 -z 2 ) e -1-z 1+z r 2 , ( D2
)
where C 0 is the contour around the origin in the complex z plane, as shown in Fig. 4. For z such that Re(z) > -1, the integrand in Eq. (D2) has a simple pole at z = 1 and an

C 0 C 1 Γ C 0 -1 1 FIG. 4.
Contours in the complex z plane used to evaluate the integral in Eq. (D2) using the decomposition in Eq. (D3).

N th order pole at z = 0 and is analytic elsewhere. We can thus replace the contour integral on C 0 by three other contour integrals as follows (see Fig. 4):

C 0 = - C 1 - C -, ( D3
)
where the contours C 1 , C , and are shown in Fig. 4. We will eventually deform the contours C and such that is a straight vertical line passing infinitesimally close to the right of z = -1 and C will be eventually sent to infinity. Evaluating the simple pole around z = 1 gives

- C 1 = 1 2π . ( D4
)
The contribution from the contour integral C is exponentially small for large N when the contour C is sent to infinity. It remains to evaluate the contour integral over . It is clear that the leading contribution to this integral over comes from the vicinity of z = -1. Hence, it is natural to make the change of variable z = -1 + z and note that the dominant contribution comes from the regime close to z = 0, i.e., z = O(1/N ) for large N . Therefore, we set z = t/N and expand the integrand for large N . To leading order, we obtain

W N (x,p) ∼ 1 2π - ( -

FIG. 1 .

 1 FIG. 1. Schematic representation of the Fermi surf (x e ,p e ) (blue solid line) defined by Eq. (7) in the phase space (x,p). The gray shaded region represents the edge around the Fermi surf, while the white inner region represents the bulk.

FIG. 2 .

 2 FIG. 2. Plot of the universal scaling function W b (a), given in Eq. (22), for two different values of the scaled inverse temperature b = 2 and b = 5. In the inset, we show a plot of the zero temperature (i.e., b → ∞) scaling function W ∞ (a) ≡ W(a) given in Eq. (11). The oscillations become more pronounced as b increases.

  we show a plot of W b (a) for two different values of b = 2 and b = 5. In the T → 0 limit, i.e., b → ∞ limit, the function W b (a) reduces to W(a) given in Eq. (11), i.e., W ∞ (a) ≡ W(a). The asymptotic behaviors of W b (a) are given by

FIG. 3 .

 3 FIG.3. Plot of W N (x,p) as a function of r = x 2 + p 2 for the 1D harmonic oscillator, as given in Eq.[START_REF] See | Ultra-cold Fermi Gases[END_REF] for N = 200. Inset: zoom on the range r ∈ [0,0.8]. Note in particular that W N (0,0) = 0 here as N = 200 is even [see Eq. (46)].

1 +

 1 Ai(2 2/3 (u + a)). (97) It reduces to the T = 0 scaling form in the limit b → ∞, where the Fermi factor becomes a Heaviside θ function. The asymptotics of W b (a) can be computed easily. After making the change of variable 2 2/3 (u + a) = y, we obtain W b (a) e a b e -b y2 -2/3 .

  D5) immediately gives the result in Eq. (D1).

  Now remember that our goal is instead the large μ expansion of the Wigner function. If one can show that (i) all terms in S 1 except the term t3 2 (p • ∇) 2 V (x) are irrelevant in the edge regime and (ii) all terms O(t 4 ) or higher are also irrelevant, then we

				)
	The sum S(x,p,t) = S 1 (x,p,t) + S 2 (x,p,t) + O(t 4 ) together
	with (C13) and (C14) provides the exact short-time expansion up to O(t 3 ) of S(x,p,t), which enters the formula (C8) for the Wigner function [it gives in fact the short-time expansion up to O(t see that W N (x,p) 1 (2π ) d	C	dt 2πit	e μ-p 2

3

) of the Fourier transform of the Euclidean propagator

G(x -y 2 ,x + y 2 ,t)]. 2 -V (x) t+ t 3 24 (|∇V (x e )| 2 +(p e •∇) 2 V (x e )) ,

ACKNOWLEDGMENTS

We thank Paul Wiegmann for stimulating discussions about Ref. [9] and Marcos Mariño for pointing out Ref. [8]. This research was supported by ANR Grant No. ANR-17-CE30-0027-01 RaMaTraF.

APPENDIX A: DERIVATION OF EQ. (25)

We start from the expression of the ground state for N noninteracting fermions as a Slater determinant:

We now evaluate the integral on the RHS of Eq. (25) as

. . dx N * 0 (x,x 2 , . . . ,x N ) 0 (x ,x 2 , . . . ,x N )

where S N denotes the group of permutations over N elements and (σ ) is the signature of the permutation σ . In the middle line each integral over x i , i = 2, . . . ,N constrains the permutations in the double sum to be the same, i.e., σ (i) = τ (i), for each i between 2 and N , from orthonormality of the single-particle eigenfunctions. However, this also constrains σ (1) = τ (1) and we thus have σ = τ . The sum in the last line is obtained by setting σ (1) = i and summing over σ . Since we are dealing with the ground state, the k i 's correspond to the N lowest eigenstates, recovering (24).

APPENDIX B: EDGE KERNEL AND WIGNER FUNCTION

In [25] it was shown that, for a d-dimensional harmonic oscillator at T = 0, the kernel near a point on the edge in position space x em = r e takes the scaling form

where the scaling function

where a n = a • r e /r e and b n = b • r e /r e . In Eq. (B2) we have Ai

• r e w N r e .

(B3)

Fourier transforming with respect to y we obtain

from Eq. (11). Furthermore, we identify the argument which appears in the scaling function as 2 2/3 a, where a has been expanded around the point x = r e , p = 0 of the edge surface in the phase space. Indeed, to lowest order in (xr e ) and p 2 one has

using the relation w 3 N r e / h2 = 1/(2m 2 ω 2 ) valid for the harmonic oscillator. This provides an alternative derivation of W N (x,p), which is valid near the special point x e = r e , p e = 0 on the edge surface in the phase space. However, for the harmonic oscillator, since W N (x,p) is isotropic in the phase space (in dimensionless units) it clearly suffices to establish the result (10) for any point on the edge surface. This result is thus fully consistent with our derivation of the scaling behavior of the Wigner function for the harmonic oscillator.

Furthermore, using the fact that the scaling form (B1), (B2) is universal for a broad class of smooth potentials, one can similarly show that the zero-temperature Wigner scaling function is also universal around the special point x = r e , p e = 0, and matches with the formula for e N , given in Eq. ( 9) at this special point only. However, since for a general potential the isotropy of the Wigner function in phase space no longer holds, this method does not allow one to obtain the scaling form at a generic point on the edge surface. However, the method used in the text does not rely on isotropy and is valid anywhere on the edge surface.

APPENDIX C: SHORT-TIME EXPANSION AND UNIVERSALITY FOR A CLASS OF SMOOTH POTENTIALS

In this appendix, for simplicity, we use units such that m = h = 1 and restore the units in the text. Consider a generic point in phase space (x e ,p e ) on the edge surface defined by

We demonstrate (the result given in the main text) that for x,p near such a point, and for the scaling variable a = O(1) with a defined as

the Wigner function takes the form

This statement is valid in the limit of large μ, which can be studied using the short-time expansion of the Euclidean propagator, extending the calculation performed in Appendix A of [25] (see below). It is useful to anticipate the main idea of the proof. First, at a generic point (C1) one has |p e | ∼ μ 1/2 and V (x e ) ∼ μ. The two terms in the energy scale e N in (C2) are thus both of the same order, with e N ∼ |∇V (x e )| 2/3 . The typical time scale t = t N which will control the final integral over t [see below in Eq. (C4)] is t N ∼ 1/e N , which, in the particular case p e = 0, also agrees with the result given by Eq. ( 282) in [25]. For a potential V (x) ∼ |x| p at large |x| the estimate is e N ∼ (μ/x e ) 2/3 ∼ μ 2(p-1)/(3p) (with x e = |x e |) and t N ∼ (x e /μ) 2/3 ∼ μ -2(p-1)/(3p) , consistent for p = 2 with t N ∼ μ -1/3 obtained for the harmonic oscillator in the text. We will justify these statements below, but it is useful to keep them in mind for estimating the various terms. Before performing the short-time expansion, let us first derive some useful exact representations for the Wigner function. We use the relation (56) between the Wigner function and the Euclidean propagator

From Appendix A of [25] we can write, as an exact starting point, the following representation (using the symmetry of the Euclidean propagator):

where . . . B denotes an average over the d-dimensional Brownian bridge B u = {B iu } d i=1 on the interval [0,1], i.e., a Gaussian process with mean zero and correlation function

and consequently we have that B 0 = B 1 = 0. By first performing a cumulant expansion and then expanding in √ tB u , we generate the short-time expansion of Eq. (C5).

We now substitute (C5) into (C4). We note that, in the absence of a potential (or if we neglect the y dependence in the potential term), we have a Gaussian integral over y with a saddle point at y = itp. This suggests that it is natural to make the change of integration variable y = itp + ỹ√ t (C7) and rewrite

e (μ-p 2 2 -V (x))t+S(x,p,t) , (C8)