Clément Iphar

Aldo Napoli

Cyril Ray

Erwan Alincourt

David Brosset

Risk Analysis of falsified Automatic Identification System for the improvement of maritime traffic safety

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

INTRODUCTION

Crossroads of international issues, the maritime domain is facing growing human activities. The ocean has a central place in some particular domains such as goods transportation and energy transportation where 90% of the global traffic is done by sea. The traffic generated by maritime activities (including fishing, sailing and cruising) is important and still increases. This ever increasing traffic leads to navigation difficulties and risks in coastal and crowded areas where a large amount of vessels exhibit different movement objectives which can be conflicting.

The international traffic is globally dense, but some regions (such as Europe or South-Eastern Asia) and some locations (such as straits or canals) have a particularly dense traffic. For instance Malacca strait and the Suez canal have respectively a yearly traffic of circa 50,000 and 20,000 vessels. Moreover this traffic sometimes takes place in hazardous areas of the world. For the two former examples, as for 2016, the Malacca strait is subject to rampant piracy [START_REF] The | Malacca Strait rampant with pirates[END_REF] and the Suez canal is located on one side of the Sinai peninsula, which is subjected to an armed conflict [START_REF] Acled | Real-time Analysis of African Political Violence[END_REF].

For the enhancement of the security and safety of navigation, several electronical systems have been designed, such as radar or Long-Range Identification and Tracking system. One of those systems is the Automatic Identification System (AIS), which broadcasts messages from a ship or a coastal base station to all surrounding ships and coastal base station within the radio horizon range. In this paper we present a risk analysis of the AIS system supported by the EBIOS methodology and an anomaly typology. After the introduction, section II presents the Automatic Identification System, principles and weaknesses, then section III details the EBIOS risks analysis conducted on this system. Section IV presents the anomaly typology developed. Section V finally proposes the data processing needed for anomaly discovery and an exemplification of integrity assessment in the case of AIS messages.

THE AIS

The principles of the system

The AIS was put in place in 2000 by the International Maritime Organization and its characteristics and deployment schedule are defined in the Safety of Live at Sea convention [START_REF] Imo | International convention for the safety of life at sea[END_REF]. This convention, initiated two years after the sinking of the RMS Titanic in 1912, aims at defining the minimal requirements to which every vessel belonging to a signatory country should comply with. The convention deals with a large scope of subjects which range from the construction of vessels to the way radiocommunications shall be done.

One of the themes of the convention is the safety and security of maritime navigation, and it is in this scope that the system was created, as it provides a real-time spatiotemporal positioning of a vessel to every vessels and shore station located in its radio range of action.

Not all the vessels from the signatory countries are concerned with the AIS regulation, as the convention states that "All ships of 300 gross tonnage and upwards engaged on international voyages and cargo ships of 500 gross tonnage and upwards not engaged on international voyages and passenger ships irrespective of size shall be fitted with an automatic identification system" [START_REF] Imo | International convention for the safety of life at sea[END_REF].

The system uses transponders to send messages through Very High Frequency marine bandwidth on two worldwide dedicated frequencies: 161.975MHz and 162.025MHz. The transponder is linked with a Global Navigation Satellite System (GNSS) which computed its location as well as external sensors (electronic compass for instance) in order to fill in several data fields within the messages. The AIS can also transmit to satellite the messages (AIS-SAT), that makes possible to keep a track of a vessel even beyond the radio horizon. One kind of message (number 27) is shorter than all other messages and is especially dedicated to satellite reception.

The system, albeit being mainly initially designed for security purposes, the AIS has alternative uses such as the prevention of boarding (alarm triggering when a small closest point of approach is computed), investigation in case of accident, control of fishing fleets, cargo fleets, global traffic, traffic in specific hazardous areas, maritime safety (for a state), aid to navigation or search and rescue operations.

Indeed, some examples of the use of AIS for alternative purposes can be found in studies in several subjects such as accident investigation [START_REF] Wang | A spatial-temporal forensic analysis for inland-water ship collision using AIS data[END_REF], the detection of near miss collision between vessels [START_REF] Zhang | A method for detecting possible near miss ship collisions from AIS data[END_REF], the behaviour understanding in a waterway through traffic monitoring [START_REF] Xiao | Comparison study on AIS data of ship traffic behavior[END_REF] or the mapping of the fishing effort using AIS data [START_REF] Natale | Mapping Fishing Effort through AIS Data[END_REF].

Due to the variety of communications, there are 27 different kinds of messages. Position report messages, information messages and check messages are the three main groups of messages [START_REF] Tunaley | Utility of Various AIS Messages for Maritime Awareness[END_REF]. All 27 kinds of messages are standardized by the International Telecommunications Union (ITU) and have its particular outline with data fields, each one being allocated a certain number of bits. The information within each field can take several forms: boolean, text, number representing a physical quantity, number representing a choice in a given list. The content of the messages do vary largely from one message to another: in a position report message the data fields are the speed, the position and the cape, amongst others, while an aid to navigation message will display the type of beacon, its name or the location of a hazard.

A system with weaknesses

Three major cases of bad data quality can be distinguished: the errors (when false data in nondeliberately broadcasted), the falsifications (when false data is deliberately broadcasted) and the spoofing (when data is created or modified and broadcasted by an outsider) [START_REF] Ray | DeAIS project: Detection of AIS Spoofing and Resulting Risks In: The proceedings of OCEANS[END_REF]. Data contained in AIS messages can be erroneous, falsified or spoofed for several reasons: there is no strong verification of the transmission, the transmission is done using a non-secured channel, some pieces of information might not be well known by the crew or the crew may want to hide some data from other people's knowledge. Those operations modify and handicap the understanding of the maritime traffic.

The errors, by nature unintentional, can be caused by transponder deficiency, a wrong input of manual data, an input of manual data of poor quality, erroneous pieces of information that come from external sensors, and can have an impact on the name of the vessel, its physical characteristics, the position or the destination for instance. Those pieces of information can then be false, incomplete, impossible according to the norm or impossible according to the physics (for instance a latitude field value shall be inferior to 90°). According to (Harati-Mokhtari et al., 2007), circa 50% of the messages contain erroneous data.

A falsification is the fact to voluntarily degrade a message by the modification of a genuine value by a false value, or by stopping the broadcast of messages, made in order to mislead the outer world. Identity theft (The Maritime Executive, 2012), the disappearances [START_REF] Windward | AIS data on the high seas: an analysis of the magnitude and implications of growing data manipulation at sea[END_REF], the broadcast of false GNSS coordinates or the statement of a wrong activity (Katsilieris et al., 2013) are types of falsification. According to (Harati-Mokhtari et al., 2007), circa 1% of the vessels broadcast falsified data.

The spoofing of messages is done by an external actor by the creation ex nihilo of false messages and their broadcast on the AIS frequencies [START_REF] Balduzzi | A security evaluation of AIS automated identification system[END_REF]. Those spoofing activities are done in order to mislead both the outer world and the crews at sea, by the creation of ghost vessels, of false closest point of approach trigger, a false emergency message or even a false cape (in the case of a spoofed vessel).

AIS EBIOS RISKS ANALYSIS

The EBIOS

The EBIOS method (ANSSI, 2010) has been created by the ANSSI (French National Agency for the Security of Information Systems) and is used in both the public and private sectors, mainly in France but also out of France. It is an approach of risk evaluation which clarifies the entities of the system, their vulnerabilities, the inventoried threats, and contributes in the assessment of the right level of security needed by the specification the security means which must be put in place. The EBIOS analysis is compliant with ISO norms 27001, 27005 and 31000. The following figure shows the 5-module procedure of the EBIOS analysis.

The first module, "Context Study", is about the identification of the system. The knowledge about the system is gathered, the architecture of the Information System is studied, as well as the technical and regulation constraints, the equipment and software used, and the human organization around the systems and which interacts with it.

The second module is defined as the "study of dread events", and permits to estimate the potential events that could affect an essential good and to clarify the associated risks. The users of the Information System provide their needs in terms of security (confidentiality for instance) with respect to the impacts they consider as inacceptable (people's security or financial cost for instance).

The third module, "Study of the threat scenarios", takes an inventory of the threats according to the technical architecture of the Information System. A list of vulnerabilities and kinds of attacks is set with respect to the materials, the network architecture and the used pieces of software, would it be caused by accident or not, would it be caused by a human factor or not. Then the next module, "Study of risks", evaluates the security needs for the identified threats, in order to highlight the risks against which the Information System shall be protected. Objectives appear and form a bill of specification on security that will con-vey the level of importance given to the threats in proportion with the required level of security.

The last module, "Study of security measures" specifies the advised security measures, the expected level of demand (which compromise should be chosen: the acceptation, the reduction or the rejection of the risk), the planning of the implementation of the measures and evaluates the benefits in terms of risks of the then applied measures.

The analysis outcome

We conducted an EBIOS analysis of the AIS, in this analysis we compiled all known information about the system in order to obtain a complete analysis of it, with circa 350 threat scenarios identified.

The application of the EBIOS method on the AIS led to the construction of several tables that enable us to consider several risk levels and the importance to put in place security measures on certain areas which have been found out as particularly vulnerable according to threat scenarios and possible threat sources.

Context Study

In this part, essential goods are identified (Tab 1), they are the different kind of information within the AIS. Then, essential functions are identified (Tab 2), which correspond to elementary indispensable functions of the system. Support goods are defined (Tab 3), taking into consideration both the vessels and shore-based stations characteristics. The following table (Tab 5) displays the synthesis of dread events, with 20 events, 11 of them being of level 3 for intolerable seriousness, 5 of level 2 for major seriousness and 4 of level 1 for minor seriousness. The levels are defined with both seriousness and plausibility to occur scales, also gratuated with minor, major and intolerable levels. For instance, a major level of risk would be reached by either events of intolerable seriousness with both insignificant and minor plausibilities, or events of major seriousness with minor, major and intolerable plausibility, or events of minor seriousness with intolerable plausibility.

Study of the threat scenarios

The following table displays the most plausible threats scenarios. 8 scenarios are presented, all of intolerable seriousness (level 3). This analysis enables us to have a great overview of the risks linked to the use of the AIS, however this analysis must be strenghtened by the discovery of anomalies in actual data. Such a study needs a characterizaion and a classification of the different kinds of possible anomalies in the case of AIS messages, and it is the topic of the following part.

TYPOLOGY OF ANOMALIES

In the EBIOS analysis presented in the previous section, a large amount of the identified risks deal with integrity, and more particularly some deal with data integrity issues. Considering integrity as one of the data quality dimensions, the risk assessment through AIS data shall take into consideration the AIS data quality. It is the knowledge of anomalous data and the determination of outliers that will lead us to the data quality assessment. In the case of AIS messages, the anomaly assessment is quite complex as it involves several disjoint concepts. This is the subject of this section.

Anomalies are not all the same of a kind, their spectrum is wide and a classification in families and subfamilies is not trivial. In the scope of the study of AIS messages and according to the research presented in (Roy et Davenport, 2010) and [START_REF] Roy | Exploitation of Maritime Domain Ontologies for Anomaly Detection and Threat Analy-sis[END_REF], a classification in four main families has been chosen: the behaviour, the content, the lawfulness and the quality. It is presented in figure 3. By its size, the behavioural anomalies family is the largest. Kinematic anomalies are the main subfamily, with on the one hand the position-based (about either the destination or the area of location) and on the other hand the movement-based (about either the route or the speed, with or without engine on) anomalies. The other subfamily is route-based anomalies, including unexpected change of destination, illogical or non-understandable behaviour.

As for the content anomalies, two subfamilies are distinguishable: the anomalies in the content of the message itself that do not come under the vessel's behaviour (such as static data, data which usually do not change over time) and anomalies about the people on board (crew or passengers).

The lawfulness anomalies can be split in two subfamilies: criminal issues (terrorism or organised crime) and breach-level issues (such as forbidden behaviour, undeclared change of flag for instance).

About data quality anomalies are distinguished the unexpectedly changing data, the impossible data so as the piece of information is out of the possible scope for it, or impossible with respect to others pieces of information (when a comparison is possible), and missing data due to poor signal reception or voluntary lack of data providing.

EXEMPLIFICATION

Anomaly detection methods

Anomaly detection is used in the field of data analysis, with the three elementary steps that are: (1) the identification of the "normality" characteristics by computation and determination of data classical signatures (for instance trajectory clusters in the AIS messages case), (2) the determination of metrics for the computation of the distance of the studied behaviour from the standard behaviour and (3) the determination of threshold criteria for the distance to the standard behaviour, allowing the normality, the abnormality, the of normality and the magnitude of abnormality for a datum.

The determination of a normal behavioural pattern can be done using diverse methods such as statistical methods, neural networks or machine learning [START_REF] Martineau | Maritime Anomaly Detection: Domain Introduction and Review of Selected Literature[END_REF], that will depend on the nature of the data to study: the created patterns can be a sequence of events, a cluster or a statistical distribution. Statistical methods are fitted with the data in which extreme values are anomalous (the speed for instance), and do not work in the cases where the anomalies are evenly distributed. Neural networks are fitted for the discovery of hidden patterns with complex boundaries but their black box behaviour and their need of a learning phase. The principle of machine learning is to automatically learn complex structures and take data-based deci-sions. Amongst those methods are genetic algorithms [START_REF] Shapiro | Genetic Algorithms in Machine Learning[END_REF], Bayesian networks [START_REF] Friedman | Bayesian Metwork Classifiers[END_REF] or clustering algorithms [START_REF] Jain | Data Clustering: A Review[END_REF].

The distance determination will depend on the type of data and distances such as Euclidian distance in dimension 1 (in the case of speed), Euclidian distance in dimension 2 (in case of localization coordinates with small distance), orthodromic distance (in case of long distance localization coordinates), mean, Hausdorff or Fréchet distances (in the cases of computation of distance between trajectories), or edition distance (in the case of distance with textual data, for the field "destination" for instance). The choice of the right distance for each field is of major importance in the anomaly detection process.

The third step in this anomaly detection is the thresholding for outlier determination in different cases, among which those aforesaid. This will be the purpose of future work measures.

Message correctness assessment

We conducted a preliminary study between October 2014 and March 2015 with an antenna located in Brest, France. We received and gathered AIS messages and led the first analysis tasks. Over the six months, 20,544,654 messages were received by the antenna (it represents a mean value of more than one message per second), out of which 1,316,689 did not comply with the ITU standardized outline for the particular message. It represents a ratio of 6.41%. The message number 1 (standard position report for class A vessel) is the more transmitted message, with 62% of all messages, out of which 4.6% are not correct. High rates of not correct current messages (both high rate and current message is considered when prevalence is superior to 1%) are found for message number 9 (standard search and rescue aircraft position report) with 99.8% of not correct messages (2.6% prevalence) and message 21 (Aids-to-Navigation report) with 39.9% of not correct messages (2.0% prevalence).

Moreover, 13 of the 27 message types have a frequency less than 0.01%, and the relevance of a study on such a little amount of message is weak. 7 message types were received less than 10 times, out of which 2 were never received over the six months span of the study.

Message integrity assessments

The first and the simplest way to assess the integrity of an AIS message is to consider the value of their fields. All 27 messages have an outline defined by the ITU in (ITU, 2014), in which all bit is allocated to a specific field. All fields have a predetermined number of bits assigned, in a predetermined order, and the value within the field has a possible range of possibilities. For instance the field latitude shall refer to a value between -90° and +90°. As the latitude field in the message number 1 gives a value in ten thousandths of minute of arc, and given that negative values (for the southern hemisphere) are given with the 2's complement, the actual value of the field shall be either inferior to 54,000,000 or superior to 80,217,728. Any value between those two numbers (except the default value) is an anomaly.

The very same conduct can be used for every field of every message. We inventoried 183 different field value assessments over the 27 messages.

In addition, it is possible to compare the fields of a single message within themselves, the fields of several messages of the same kind (sent by the same vessel or by several vessels) between themselves and the fields of several messages of the different kinds (sent by the same vessel or by several vessels) between themselves. A total number of respectively 94, 114 and 278 different integrity assessments for the three kinds of comparison have been inventoried. Their study will be the purpose of future work.

CONCLUSION

This paper introduces the AIS principles and vulnerabilities that lead to errors and falsifications of the message. These errors and falsifications hamper the view of the maritime traffic, both for the mariners (for their immediate surrounding) and the local authorities (for the understanding of the activity off their coasts), causing a degradation of the maritime navigation safety. An EBIOS risks analysis has been conducted that highlighted the dreads events and the threat scenarios of the AIS when failing to transmit correct data. The concept of anomalies in the maritime environment of the broadcasting of AIS messages have been illustrated by a typology and anomaly detection techniques have been presented. An exemplification of AIS correctness and integrity as-sessment is proposed, showing the large scope of research that has to be done in order to assess the reliability of AIS, for the improvement of the maritime situation awareness and thus the improvement of the safety of navigation.

Figure 1 .

 1 Figure 1. The maritime traffic in Europe (AIS data)

Figure 2 .

 2 Figure 2. The EBIOS procedure

Figure 3 .

 3 Figure 3. Typology of anomalies in the AIS case

Figure 4 .

 4 Figure 4. Locations of all messages received during the study

Table 1 .

 1 The

essential goods Static AIS data: no evolution over time expected Dynamic AIS data: positioning and route-based data Other AIS data: text and binary messages, weather

Table 2. The essential functions Transmit AIS data Receive AIS data Visualize data by table or card Transponder configuration (possibility to change) Analyse, treat and store AIS data (for surveillance centre) Determine its position (whatever the positioning system)

Table 3 .

 3 List of support goods Study of dreads events A total of 17 threat sources have been identified, and presented in the following table(Tab 4)

	INTERNAL	EXTERNAL
	AIS chain	GNSS
	Surveillance tools	AIS data supplier (incl. AIS SAT)
	Vessel Organisation	Public or private IP network
		(incl. the Internet)
	Surveillance Centre
	Organisation
	Vessel premises
	Surveillance Centre
	premises

Table 4 .

 4 Identified threat sources

	Vessel's crew
	Rival vessel or ship-owner
	Pirate
	Operator of the collecting system
	Administrator of the collecting system
	Data supplier
	Contributor of the data supplier
	Manufacturer or installer of the material
	Hacker
	Hacktivist
	Criminal organisation
	State organisation
	Malicious code from unknown origin
	Natural phenomenon
	Natural or health disaster
	Animal activity
	Internal event

Table 5

 5

	. Synthesis of dread events
	3 Dynamic data are not available
	3 Dynamic data are false
	3 Other AIS data are false
	3 AIS data transmission is not working
	3 AIS data transmission is not reliable
	3 AIS data reception is not working
	3 AIS data reception is not reliable
	3 Position determination is impossible
	3 Position determination is incorrect
	3 Information visualisation is not available
	3 Information visualisation is incorrect
	2 Static data are accessible to an opponent
	2 Dynamic data are accessible to an opponent
	2 Other AIS data are accessible to an opponent
	2 Data analysis, treatment or storage is corrupted
	2 Transponder configuration is corrupted
	1 Static data are not available
	1 Static data are false
	1 Other AIS data are not available
	1 Data analysis, treatment or storage is not available

Table 6 .

 6 Most plausible threats scenarios

	Use of the computer for other tasks
	Use of the system for sending and receiving short
	messages (addressed or broadcasted)
	Branching of a non-legitimate input
	Identity data change on the transponder
	Wrong update of route-based and cargo data
	Unexpected case according to the norm
	Maritime traffic surveillance for villainous ends
	Injection of false data hijacking report tools
	3.2.4 Study of risks
	The following table (Tab 7) displays the 11 risks of
	intolerable level highlighted by the study.

Table 7 .

 7 Intolerable level risks identified Risks linked to the dynamic data availability Risks linked to the dynamic data integrity Risks linked to the other AIS data integrity Risks linked to the transmission function availability Risks linked to the transmission function integrity Risks linked to the reception function availability Risks linked to the reception function integrity Risks linked to the positioning function availability Risks linked to the positioning function integrity Risks linked to the visualisation function availability Risks linked to the visualisation function integrity

ACKNOWLEDGMENTS

Research presented in this paper is supported by The French National Research Agency (ANR) and cofunded by DGA (French Armaments Procurement Agency) under reference ANR-14-CE28-0028. The project is also labelled by French clusters Pôle Mer Bretagne Atlantique and Pôle Mer Méditerranée.


I Bouchoule

S S Szigeti

M J Davis

K V Kheruntsyan

Finite-temperature hydrodynamics for one-dimensional Bose gases: Breathing mode oscillations as a case study

Keywords: numbers: 03.75.Kk, 67.85.-d, 05.30.Jp

We develop a finite-temperature hydrodynamic approach for a harmonically trapped one-dimensional quasicondensate and apply it to describe the phenomenon of frequency doubling in the breathing-mode oscillations of the quasicondensate momentum distribution. The doubling here refers to the oscillation frequency relative to the oscillations of the real-space density distribution, invoked by a sudden confinement quench. By constructing a nonequilibrium phase diagram that characterises the regime of frequency doubling and its gradual disappearance, we find that this crossover is governed by the quench strength and the initial temperature, rather than by the equilibrium-state crossover from the quasicondensate to the ideal Bose gas regime. The hydrodynamic predictions are supported by the results of numerical simulations based on a finite-temperature c-field approach, and extend the utility of the hydrodynamic theory for low-dimensional quantum gases to the description of finite-temperature systems and their dynamics in momentum space.

Hydrodynamics is a powerful and broadly applicable approach for characterizing the collective nonequilibrium behavior of a wide range of classical and quantum fluids, including Fermi liquids, liquid helium, and ultra-cold atomic Bose and Fermi gases [1][2][3][4][5][6]. For ultra-cold gases, the hydrodynamic approach has been particularly successful in describing the breathing (monopole) and higher-order (multipole) collective oscillations of harmonically trapped three-dimensional (3D) Bose-Einstein condensates [2,6,7]. For condensates near zero temperature, the applicability of the approach stems from the fact that for long-wavelength (low-energy) excitations the hydrodynamic equations are essentially equivalent to those of superfluid hydrodynamics, which in turn can be derived from the Gross-Pitaevskii equation for the order parameter. For partially condensed samples at finite temperatures, the hydrodynamic equations should be generalized to the equations of two-fluid hydrodynamics, where the applicability of the approach to the normal (thermal) component of the gas is justified by fast thermalization times due to collisional relaxation [3,8].

In contrast to 3D systems, the applicability of the hydrodynamic approach to 1D Bose gases is not well established. Firstly, in the thermodynamic limit 1D Bose gases lack the long-range order required for superfluid hydrodynamics to be a priori applicable. Secondly, the very notion of local thermalisation, required for the validity of collisional hydrodynamics of normal fluids, is questionable due to the underlying integrability of the uniform 1D Bose gas model [9]. Despite these reservations, the hydrodynamic approach has already been applied to zero-temperature (T = 0) dynamics of 1D Bose gases in various scenarios [10][11][12][13][14][15] (for related experiments, see [16][17][18]). The comparison of hydrodynamic predictions with exact theoretical results is challenging. In Ref. [13], time-dependent density matrix renormalization group simulations of the collision of 1D Bose gases at T = 0 found reasonable agreement with the hydrodynamic approximation, although the latter failed to predict short wavelength dynamics such as shock waves. An alternative approximate approach, based on the conservation of Lieb-Liniger rapidities, has been applied to describe the free expansion dynamics of a T = 0 1D gas [15] and was able to reproduce the hydrodynamic results for both weak and strong interactions.

At finite temperatures finding exact predictions is extremely difficult, and thus developing a hydrodynamic approach is appealing, despite its lack of justification. Here, we develop a general finite-T hydrodynamic approach suitable for 1D Bose gases and specifically apply it to the breathing-mode oscillations of a harmonically trapped 1D quasicondensate. We find the predictions agree both with experimental observations [18] and numerical simulations of a finite-temperature c-field methodology [19,20]. More remarkably, our hydrodynamic approach not only adequately describes the dynamics of the density distribution of the gas (the standard observable of the hydrodynamic theory), but it can be also used to describe the dynamics of the momentum distribution. This is a key observable for quantum gas experiments, and has not previously been accessible from a hydrodynamic approach.

understanding the momentum-space properties of equilibrium quasicondensates [21][22][23][24][25][26], the said experimental data for dynamics were well-described by a simple zero-temperature classical hydrodynamic approach, wherein the frequency doubling was interpreted as a result of a self-reflection mechanism due to the mean-field interaction energy barrier. In contrast to this behavior, no frequency doubling was observed in the nearly ideal Bose gas regime, as expected for a noninteracting gas. The experimentally observed smooth crossover from the regime of frequency doubling to no doubling has so far not been explained theoretically.

Here, we explain this phenomenon within the hydrodynamic approach and construct, for the first time, the corresponding nonequilibrium phase diagram, showing that the frequency doubling crossover is governed by the quench strength and a nontrivial combination of the temperature and interaction strength. For small enough quenches, the crossover from frequency doubling to no doubling can lie entirely within the quasicondensate regime, and does not require an equilibriumstate crossover to the ideal Bose gas regime. Constructing and studying phase diagrams is an important goal in many areas of physics, and our findings here serve as an example where equilibrium and dynamical phase diagrams are not identical. We confirm our predictions by comparing the hydrodynamic results to those obtained numerically using finite-temperature c-field simulations (for a review, see [20]) based on the projected Gross-Pitaevskii equation (PGPE) [19].

1. Hydrodynamic equations and evolution of the density distribution.-The hydrodynamic approach relies on the local density approximation (LDA) and assumes that the 1D system can be divided into small locally uniform slices, each of which is in thermal equilibrium in the local moving frame. Moreover, one can assume that heat transfer between the slices is negligible for long-wavelength excitations [27], which implies that each slice undergoes isentropic (de)compression. The hydrodynamic description of this system is [1]

∂ t ρ + ∂ x (ρv) = 0, (1a)
∂ t v + v∂ x v = - 1 m ∂ x V (x, t) - 1 mρ ∂ x P, (1b)
∂ t s + v∂ x s = 0, (1c)
where ρ(x, t) is the local 1D density of the slice at position x, v(x, t) is the respective hydrodynamic velocity, s(x, t) is the entropy per particle, P (x, t) is the pressure, m is the mass of the constituent particles, and V (x, t) is the external trapping potential which for our case study is harmonic, V (x, t) = 1 2 mω(t) 2 x 2 , of frequency ω(t). We now apply the hydrodynamic approach to describe the post-quench dynamics induced by the following scenario. Initially the atomic cloud with density profile ρ 0 (x) is in thermal equilibrium at temperature T 0 in the trap of frequency ω 0 . Subsequently, at time t = 0, the trap frequency is suddenly changed to ω 1 . To characterize the ensuing dynamics in different regimes of the 1D Bose gas, we introduce the dimensionless interaction parameter γ 0 = mg/ 2 ρ 0 (0) and the dimensionless temperature t 0 = 2 2 k B T 0 /mg 2 [28,29], where g is the coupling strength of the pairwise δ-function interaction potential. The solutions of the HDEs (1a)-(1c) describing this harmonic-confinement quench depend only on the thermodynamic equation of state of the gas. In each of the following three cases, (i) -ideal gas regime (t 0 , γ 3/2 0 t 0 1), (ii) -strongly interacting or Tonks-Girardeau regime (γ 0 , 1/t 0 1), and (iii) -quasicondensate regime (γ 0 , γ 3/2 0 t 0 1), the solutions reduce to scaling equations of the form

ρ(x, t) = ρ 0 (x/λ(t))/λ(t), v(x, t) = x λ(t)/λ(t), (2)
T (t) = T 0 /λ(t) ν+1 , (3)
where the scaling parameter λ(t) [with λ ≡ dλ(t)/dt, λ(0) = 1, and λ(0) = 0] satisfies the ordinary differential equation,

λ = -ω 2 1 λ + ω 2 0 /λ 2ν+1 , (4)
with the value of ν in different regimes given below [START_REF]scaling solutions (2)-(3) remain applicable for arbitrary time[END_REF]. The hydrodynamic solution (3) for the temperature is one of the key results of this paper as it allows one to simply calculate the evolution of the temperature-dependent momentum distribution of the gas (see below).

(i) Ideal gas regime (t 0 , γ 3/2 0 t 0 1): In this case ν = 1, and the validity of the above scaling solutions can be demonstrated using a dimensional analysis of the equation of state (see Ref.

[31]), which we note is also applicable to an ideal Fermi gas. Equation (4) in this regime has an explicit analytic solution,

λ(t) = 1 + sin 2 (ω 1 t). (5)
This corresponds to harmonic oscillations of the mean squared width of the density profile, occurring at frequency ω B = 2ω 1 , with ≡ (ω 0 /ω 1 ) 2 -1 characterizing the quench strength. This result coincides with that for a noninteracting gas obtained from a single-particle picture. The fact that the hydrodynamic approach, which a priori assumes sufficient collisions to ensure local thermal equilibrium, agrees with the results for a noninteracting gas is specific to the harmonicconfinement quench considered here and is accidental.

(ii) Strongly interacting regime (γ 0 , 1/t 0 1): Here, the equation of state is that of an ideal Fermi gas so that the previous ideal gas results apply, and Eqs. ( 2)-( 4) are fulfilled with ν = 1. The breathing mode oscillations of the momentum distribution of a finite-temperature Tonks-Girardeau gas are discussed elsewhere [START_REF] Atas | [END_REF].

(iii) Quasicondensate regime (γ 0 , γ 3/2 0 t 0 1): In this case ν = 1/2, and the validity of the scaling solutions (2) can be demonstrated using the equation of state P = 1 2 gρ 2 . The latter can be derived from the quasicondensate chemical potential µ = gρ and the Gibbs-Duhem relation ρ = (∂P/∂µ) T . For a weak quench, 1, the solution to Eq. ( 4) oscillates at frequency ω B √ 3ω 1 and is nearly harmonic with an amplitude λ(t = π/ω B )-1 2 /3 [START_REF]For a small amplitude quench ( 1), Eq. (4) with ν = 1/2 be solved using the method of linearization [i.e., expanding λ(t) as λ(t) = δλ(t), with δλ(t) , yielding[END_REF]. According to Eq. ( 2), the density profile breathes self-similarly, maintaining its initial Thomas-Fermi parabolic shape, ρ 0 (x) = ρ 0 (0)(1x 2 /X 2 0 ) for x ≤ X 0 [ρ 0 (x) = 0 otherwise], with X 0 = 2gρ 0 (0)/mω 2 0 .

Finite-temperature effects are not seen in the dynamics of the density distribution ρ(x, t) [34] because in this regime the equation of state does not depend on the temperature. However, as we show below, such effects can be revealed in the dynamics of the momentum distribution.

2. Dynamics of the momentum distribution.-Let us consider a slice of the gas in the region [x, x + dx] of density ρ(x, t), velocity v(x, t), and entropy per particle s(x, t). In the laboratory frame its momentum distribution is n(ρ, s, kmv/ ), where n is the equilibrium momentum distribution of a homogeneous gas, which we normalize to ´dk n(ρ, s; k) = ρ. The total momentum distribution is then given by n(k, t) = ˆdx n(ρ, s; kmv(x, t)/ ).

There are two contributions to n(k, t): the hydrodynamic velocity field and the contribution of thermal velocities, which have different effects on the breathing mode oscillations. In order to see the sole effect of the hydrodynamic velocity field, let us first disregard the effect of the thermal velocities, taking n(ρ, s; kmv(x, t)/ ) = ρδ(kmv(x, t)/ ), where δ(k) is the Dirac delta function. If a scaling solution as in Eq. ( 2) exists, the hydrodynamic component of the momentum distribution evolves according to

n h (k, t) = m| λ| ρ 0 k m λ . (7)
For oscillatory λ(t), this implies that n h (k, t) collapses to a zero-width distribution twice per position-space density breathing cycle: when the width of the cloud in real space is both largest and smallest, both corresponding to λ = 0. Therefore, the oscillations of the hydrodynamic contribution to the momentum distribution will always display frequency doubling. Now consider the additional contribution of thermal velocities to n(k, t), which changes as each slice undergoes isentropic compression and decompression during the breathing cycle. Since one expects the momentum width to be a monotonically increasing function of the compression factor, the thermal momentum width of each slice [and hence of the overall momentum distribution n(k, t)] is expected to oscillate out-of-phase relative to the width of the real-space density profile, but at the same breathing frequency ω B .

The evolution of the overall momentum distribution n(k, t) results from the combination of the hydrodynamic and thermal parts. For a near-ideal gas at T > 0 this leads to a somewhat fortuitous cancellation of the hydrodynamic velocity field by the thermal component, and so the total momentum distribution always oscillates at ω B = 2ω 1 (see Ref. [31]) and never displays frequency doubling, consistent with the single-particle picture.

The situation is different, however, in the quasicondensate regime. The momentum distribution of a homogeneous quasicondensate of density ρ and temperature T , for wavelengths in the phononic regime (i.e., k / √ mgρ) is given by a

Lorentzian n(ρ, s; k) = (2ρl φ /π)/ 1 + (2l φ k) 2 [31]
. Substituting this Lorentzian into Eq. ( 6), we obtain the full momentum distribution of the trapped gas,

n(k, t) = 1 π ˆdx 2l φ (x, t)ρ(x, t) 1 + 4[l φ (x, t)] 2 [k -mv(x, t)/ ] 2 , (8)
where l φ (x, t) = 2 ρ(x, t)/mk B T (t). According to the scaling solutions ( 2) and ( 3) with ν = 1/2 (see also [31]), l φ (x, t) evolves as

l φ (x, t) = √ λ l (0) φ ρ0 (x/λ), (9)
where ρ0 (x) = 1x 2 /X 2 0 is the scaled initial density profile and l

(0) φ = 2 ρ 0 (0)/mk B T 0 = 2[ρ 0 (0)γ 2 0 t 0 ] -1 .
Combining the scaling solution for l φ (x, t) with that for ρ(x, t), and changing variables to u = x/λX 0 in Eq. ( 8), leads to the following final result

n(k, t) = B λˆ1 -1 du (1 -u 2 ) 2 1 + 4 λ(1 -u 2 ) 2 k -ω1 ω0 A λu 2 . ( 10
) Here, k = l (0) φ k, A = mω 0 X 0 l (0) φ / = √ 8/γ 3/2 0 t 0 , and B = 2ρ 0 (0)l (0)
φ X 0 /π is a normalization factor. In addition, we have introduced a dimensionless time τ ≡ ω 1 t, so that the dimensionless functions λ(τ ) ≡ λ(τ /ω 1 ) and λ = d λ/dτ , obtained from Eq. ( 4), depend only on the ratio ω 1 /ω 0 , or equivalently only on the quench strength = (ω 0 /ω 1 ) 2 -1. Thus, for a given , the evolution of n(k, t) is governed solely by the dimensionless parameter A, which itself depends only on the initial intensive parameters γ 0 and t 0 . Note that A 1 in the quasicondensate regime where γ 3/2 0 t 0 1 [25,28]. Using Eq. ( 10) for a given A and quench strength , we can now compute the evolution of the full momentum distribution and its half width at half maximum (HWHM); see Figs. 1(a)-(b). The HWHM can then be fitted with a sum of two sinusoidal functions: the fundamental mode oscillating at ω B ( √ 3ω 1 , for 1) and the first harmonic oscillating at 2ω B , with amplitudes c 1 and c 2 , respectively. Defining the weight of the fundamental mode as

K = c 2 1 /(c 2 1 + c 2 2
), we identify the frequency doubling phenomenon with K 1, whereas K 1 corresponds to the absence of doubling. The doubling crossover can, therefore, be defined as the value of A = A cr for which K = 1/2. As we show in [31], for small quench amplitudes one expects the frequency doubling to occur for A √ 1, while for A √ 1 the thermal effects dominate and the frequency doubling is absent; accordingly, A cr is expected to scale as A cr ∝ 1/ √ . Figure 1(c) shows the nonequilibrium phase diagram of the crossover from frequency doubling to no doubling and confirms that A cr , obtained using Eq. ( 10) and the fitting procedure described above, does indeed scale as ∝ 1/ √ .

In Fig. 1(d) we superimpose the conditions for observing frequency doubling on the equilibrium phase diagram of the 1D Bose gas. As we see, for small enough quench strengths, ) corresponds to the crossover between the two equlibrium regimes. The two lines with filled (red) and open (blue) circles, on the other hand, correspond to the frequency doubling crossover conditions for two different quenches, = 0.1 and = 0.02 (the respective data points in (c) are labelled in the same way). The (light red and light blue) shaded areas underneath these lines correspond to the conditions where the frequency doubling occurs.

Although the applicability of the hydrodynamic theory in this system might be questionable, our analytic results have been benchmarked against finite-temperature c-field simulations, which is an approach that is identical (in 1D) to the truncated Wigner approximation and whose validity for degenerate weakly interacting Bose gases is well established [19,20,[35][START_REF] Davis | [END_REF][37]. In this approach, the atomic field is treated as a classical field, whose evolution is governed by the timedependant Gross-Pitaevskii equation (GPE), with the initial stated being generated using a stochastic evolution that samples the Gibbs ensemble [START_REF]scaling solutions (2)-(3) remain applicable for arbitrary time[END_REF]. Qualitatively, the same behaviour as in Figs. 1(a)-(b) based on the hydrodynamic approach occurs in our c-field simulations [34]; quantitatively, the crossover from doubling to no-doubling is in line with the analytic predictions [see Fig. 1(c)]. Moreover, as we argue in Ref. [START_REF]scaling solutions (2)-(3) remain applicable for arbitrary time[END_REF], for sufficiently weak confinement (small ω 0 ), the cfield dynamics are governed by just two dimensionless param-eters, A and , as predicted from the hydrodynamic approach. Overall, the performance of the hydrodynamic theory-as validated by our c-field simulations-in modelling the harmonic confinement quench of a finite-temperature quasicondensate is remarkable. Moreover, even though the hydrodynamic results of Eq. ( 10) formally require A 1 to ensure the applicability of the quasicondensate regime, our comparison with c-field simulations shows that Eq. ( 10) continues to give accurate predictions even for moderate values of A 1: the predictions capture correctly the crossover values of A cr even for large to moderate quench strengths, for which A cr lies close to the equilibrium crossover to the ideal Bose gas regime.

In summary, we have developed a finite-temperature hydrodynamic approach for a harmonically trapped 1D Bose gas and applied it to the study of breathing mode oscillations in the quasicondensate regime. In addition to the standard scope of the hydrodynamic theory, which is to describe the evolution ) corresponds to the crossover between the two equilibrium regimes. The two lines with filled (red) and open (blue) circles, on the other hand, correspond to the frequency doubling crossover conditions for two different quenches, = 0.1 and = 0.02 (the respective data points in (c) are labelled in the same way). The (light red and light blue) shaded areas underneath these lines correspond to the conditions where the frequency doubling occurs. the crossover from doubling to no doubling lies well within the quasicondensate regime. We therefore conclude that this phenomenon is governed not by the crossover from the ideal Bose gas regime into the quasicondensate regime, but by the competition between the hydrodynamic velocity (which always displays doubling) and the broadening/narrowing of the thermal component of the gas due to adiabatic compression/decompression (which always oscillates at the fundamental frequency ω B ).

Although the applicability of the hydrodynamic theory in this system might be questionable, our analytic results have been benchmarked against finite-temperature c-field simulations, whose validity for degenerate weakly interacting Bose gases is well established [19,20,[START_REF] Davis | [END_REF][37][START_REF] Castin | Carusotto[END_REF]. In this approach, the Bose gas is approximated as a classical field, whose evolution is governed by the time-dependent Gross-Pitaevskii equation (GPE), with the initial state being sampled from the classical Gibbs ensemble for the given temperature and density [31]. Qualitatively, the same behaviour as in Figs. 1(a)-(b) based on the hydrodynamic approach occurs in our c-field simulations [35]; quantitatively, the crossover from doubling to nodoubling is in broad agreement with the analytic predictions [see Fig. 1(c)]. Moreover, as we argue in Ref. [31], for sufficiently weak confinement (small ω 0 ), the c-field dynamics are governed by just two dimensionless parameters, A and , as predicted from the hydrodynamic approach. Overall, the performance of the hydrodynamic theory-as validated by our c-field simulations-in modelling the harmonic confinement quench of a finite-temperature quasicondensate is remarkable. Moreover, even though the hydrodynamic results of Eq. (10) formally require A 1 to ensure the applicability of the quasicondensate regime, our comparison with c-field simulations shows that Eq. ( 10) continues to give accurate predictions even for moderate values of A 1.

In summary, we have developed a finite-temperature hydrodynamic approach for a harmonically trapped 1D Bose gas and applied it to the study of breathing mode oscillations in the quasicondensate regime. While the usual scope of the hydro-dynamic theory is to describe the evolution of the real-space density of a gas, our approach extends its utility to describe the evolution of its momentum distribution. The approach allowed us to discern the contribution of the hydrodynamic velocity field and that of thermal excitations in the oscillatory dynamics of the momentum distribution of the 1D quasicondensate, hence explaining the full mechanism behind the phenomenon of frequency doubling and the crossover to no doubling. The hydrodynamic predictions are in broad agreement with numerical simulations based on finite-temperature c-field simulations. Our approach can address not only the sudden quench scenario studied here, but also the dynamics under arbitrary driving of the trapping frequency ω(t), in which case the differential equation for the scaling parameter λ(t), Eq. ( 4), must be solved numerically. Future extensions of this work will concern the treatment of breathing mode oscillations in the strongly interacting regime [START_REF] Atas | [END_REF], and could also address collective behavior of 1D Bose gases in anharmonic traps, as well as of 2D and 3D quasicondensates in highly-elongated geometries. K .V. K. acknowledges stimulating discussions with D. M. Gangardt. I. B. acknowledges support by the Centre de Compétences Nanosciences Île-de-France; S. S. S. -by the Australian Research Council Centre of Excellence for Engineered Quantum Systems (grant No. CE110001013); M. J. D. and K. V. K. -by the ARC Discovery Project grants DP160103311 and DP140101763. Numerical simulations were performed using XMDS2 [39] on the University of Queensland School of Mathematics and Physics computer "Obelix", with thanks to I. Mortimer for computing support.

For a uniform ideal gas (either bosonic or fermionic) at temperature T and 1D density ρ = N/L, where N is the number of particles and L is length of the confinement box, the only two length scales are the mean interparticle separation ρ -1 and the thermal de Broglie wavelength Λ T = 2π 2 /(mk B T ); the corresponding energy scales are 2 ρ 2 /m and k B T . Using the thermodynamic definition of the 1D pressure, P = (∂U/∂L) s , where U is the internal energy, one can apply simple dimensional analysis to write down the equations of state for P and s:

P = k B T ρ F ( 2 ρ 2 /mk B T ), (S1)
s/k B = G(mk B T / 2 ρ 2 ). (S2)
Here F and G are dimensionless functions of the only dimensionless parameter-the ratio of the two energy scales. With this choice of expression for P , the classical ideal gas law in the high temperature limit is recovered with F 1. For a highly degenerate ideal Fermi gas, on the other hand, the equation of state

P = 2 π 2 ρ 3 /(3m) is recovered with F (π 2 /3)( 2 ρ 2 /mk B T ).
Now consider a confinement quench of the gas. Applying the general functional forms of P and s to small (locally uniform) slices of the gas, it can be shown by direct substitution that the scaling solutions, Eqs. ( 2) and (3) of the main text, satisfy Eqs. (1a) and (1c). For Eq. (1b), first note that Eq. (S1) together with the scaling solutions imply that P (x, t) = P 0 (x/λ)/λ 3 , where P 0 (x) ≡ P (x, 0). Since Eq. (1b) is assumed true at this initial time, then

∂ x P 0 = -ρ 0 (x)∂ x V (x, 0) = -mω 2 0 xρ 0 (x).
Together with the scaling solutions and Eq. ( 4) of the main text, this relation is sufficient to show that Eq. ( 1b) is true for all times.

B. Equilibrium momentum distribution of a uniform 1D

quasicondensate.

Here we outline the derivation of the Lorentzian shape of the equilibrium momentum distribution of a uniform 1D quasicondensate (for a more detailed derivation, we refer the reader to Refs. [1,2]). For a uniform and hence translationally invariant system, the momentum distribution is given by the Fourier transform of the first-order correlation function, G (1) (x) = Ψ(x) Ψ(0) , where Ψ(x) is the bosonic field operator. In the quasicondensate regime, corresponding to T √ γ 2 ρ 2 /(2mk B ) [3] (with ρ being the uniform density and γ ≡ mg/ 2 ρ), the first-order correlation function is dominated by the long-wavelength (low-energy) excitations whose Hamiltonian reduces (using the density-phase representation of the field operator, Ψ(x) = ρ + δ ρ(x)e i φ(x) ) to the Luttinger liquid form [2,4]

ĤL = ˆdx 2 ρ 2m (∂ x φ) 2 + g 2 (δ ρ) 2 . (S3)
Here, δ ρ(x) is the operator describing the density fluctuations, canonically conjugate to the phase operator φ(x).

In evaluating G (1) (x) we note that the density fluctuations are small in the quasicondensate regime and can be completely neglected as long as the relative distances of interest are much larger than the healing length ξ = / √ mgρ [1,2,5,6]. As the Luttinger liquid Hamiltonian (S3) is quadratic in φ, the correlation function G (1) (x) = ρ e i( φ(x)-φ(0)) can be expressed through the mean-square fluctuations of the phase via Wick's theorem:

G (1) (x) = ρe -1 2 [ φ(x)-φ(0)] 2 . (S4)
Denoting the Fourier component of φ(x) at wavevector k via φk , the corresponding term in the expectation value of the Hamiltonian is given by ´dk 2π Lρ 2 k 2 | φk | 2 /(2m), where L is the length of the uniform system. The relevant modes contributing to G (1) (x) are highly populated, such that a classical field picture is sufficient. The energy per quadratic degree of freedom is thus given by k B T /2 and therefore | φk | 2 = mk B T /(Lρ 2 k 2 ) [7]. Using this result one can then show, after a little algebra, that Eq. (S4) yields

G (1) (x) ρe -|x|/2l φ , (|x| ξ), ( S5
)
where l φ ≡ 2 ρ/(mk B T ). The Fourier transform of this exponentially decaying correlation function gives the Lorentzian momentum distribution, n(k) = (2ρl φ /π)/ 1 + (2l φ k) 2 , used in the main text.

C. Scaling solution for the temperature in the quasicondensate regime.

In order to calculate n(k, t) using the hydrodynamic scaling solutions, Eqs. ( 2) and ( 4) of the main text, we first need to determine how the temperature of the gas and hence the phase correlation length l φ = 2 ρ/(mk B T ) evolves during the breathing oscillations. To do this, we first note that the energy of the jth phonon mode in a quantization box of length L is given by E j = k j c, where k j = 2π L j is the phonon wave vector and c = (∂P/∂ρ)/m is the speed of sound. Using P ∝ ρ 2 , we find that c scales as c ∝ ρ 1/2 and therefore E j ∝ L -1 ρ 1/2 ∝ ρ 3/2 . Consider now an adiabatic compression/decompression cycle for a uniform slice of the gas confined to a box of length L. Such a compression does not change the mean occupation number of the mode j. The mode occupation is given by n j k B T /E j in the long-wavelength limit and scales as n j ∝ T /ρ 3/2 , whereas ρ ∝ λ -1 according to the scaling solution, Eq. ( 2) of the main text. Therefore, during the adiabatic breathing oscillations, the temperature T (t) evolves from the initial value T (0) ≡ T 0 to T (t) = T 0 /λ 3/2 , i.e., Eq. ( 3) of the main text with ν = 1/2.

D. Evolution of the momentum distribution for an ideal gas.

Applying the hydrodynamic approach and Eq. ( 7) of the main text to the ideal gas regime, we first note that the momentum distribution of a uniform ideal gas (normalized to ´dk n(k) = ρ) is given by n(ρ, s; k) = N (( 2 k 2 /2mµ)/k B T ), where N is a dimensionless function whose expression depends on the quantum statistics [8]. Since µ/k B T is a function of s (in the sense of a thermodynamic equation of state), which itself is a function of mk B T / 2 ρ 2 [see Eq. (S2)], one can assert that µ/(k B T ) = G(mk B T / 2 ρ 2 ), where G is a dimensionless function. Then, the scaling solutions ( 2) and (3) of the main text imply that µ(x, t

)/k B T (t) = G[mk B T (t)/ 2 ρ 2 (x, t)] = G[mk B T 0 / ρ 2 0 (x/λ)] ≡ G 0 (x/λ), or µ(x, t) = µ(x/λ, 0)/λ 2 = [µ 0 -1 2 mω 2 0 (x/λ) 2 ]/λ 2 , ( S6
) where µ 0 is the initial chemical potential in the trap center. Substituting n(ρ, s; k) along with this expression for µ(x, t) into Eq. ( 6) of the main text, and changing variables to x = αx-kλ λ/(mω 2 0 α), gives n(k, t) = n 0 (k/α)/α, (S7)

where n 0 (k) is the initial momentum distribution of the trapped gas and α 2 = (ω 2 0 + λ 2 λ2 )/(λω 0 ) 2 . Using Eq. ( 5) of the main text, we can explicitly write α as

α = [1 + cos 2 (ω 1 t)]/(1 + ), ( S8
)
which implies that the momentum distribution of a finitetemperature ideal gas in the hydrodynamic limit oscillates at ω B = 2ω 1 and never displays frequency doubling. One thus recovers the expected behaviour for an ideal gas, due to the position-momentum symmetry of the underlying harmonic oscillator Hamiltonian. The fact that this result is reproduced within the hydrodynamic approach is a result of a "fortuitous" exact cancelation of the effect of the hydrodynamic velocity field by the thermal component.

E. Scaling of the frequency doubling crossover Acr with the quench strength .

In this section we make qualitative arguments that provide an understanding of the dominant oscillation regimes in the dynamics of the momentum distribution of a quasicondensate, and derive an approximate scaling of A cr with . Let us first introduce typical momentum scales involved in the dynamics. For small-amplitude oscillations, corresponding to 1, the scaling parameter λ oscillates as λ(t)

1 + 3 -3 cos( √ 3ω 1 t) and therefore the magnitude of λ is of the order of ∼ / √ 3. This means that the characteristic hydrodynamic momentum, which can be estimated as kh ∼ mX 0 λ/ from Eq. ( 7) of the main text, and which can be rewritten as kh

∼ mX 0 ω 1 / ∼ (ω 1 /ω 0 )A /l (0) φ , is of the order of kh ∼ A /l (0) φ for 1.
Compared to this, the characteristic thermal momentum during the compression/decompression cycle oscillates above kth ∼ 1/l (0) φ with an amplitude variation of δ kth ∼ /l (0) φ kth . For A 1, the characteristic hydrodynamic momentum is much larger than both the characteristic thermal momentum and its variation, kh kth , δ kth . Then, as long as one is interested in momenta of the order of kh , the function n in Eq. ( 6) of the main text can be approximated by a δ-function, and therefore the breathing oscillations of the momentum distribution will be dominated by the hydrodynamic phenomenon of frequency doubling.

In the opposite regime of A 1, the characteristic hydrodynamic momentum is much smaller than the characteristic thermal momentum, kh kth , and therefore the above approximation, which neglects the effect of the width of n in Eq. ( 6) of the main text, breaks down. In this case, the contribution of the hydrodynamic momenta to n(k, t) can instead be estimated via a Taylor series of n as powers of k h (x, t) = mv(x, t)/ . In this series, the contribution of the first-order term to the integral vanishes because ∂ n/∂k is an even function of x, whereas k h (x, t) is odd; therefore, the hydrodynamic velocity field has no effect on n(k, t) in this order. The leading-order correction in n(k, t) thus comes from the second-order derivative term, proportional to k 2 h . To estimate this correction, let us consider the typical variations of the peak value of the momentum distribution n(0, t), which we denote via δn, induced solely by the hydrodynamic momenta. Using ∂ 2 n/∂k 2 | k=0 ∼ ρ 0 (0)(l (0) φ ) 3 in Eq. ( 6), one can find that δn ∼ -X 0 ρ 0 (0)l

(0) 3 φ k2 h , or δn ∼ -X 0 ρ 0 (0)l (0)
φ (A ) 2 , where we have used kh ∼ A /l

φ . Since the peak value n(0, t) is inversely proportional to the characteristic momentum width W , the hydrodynamic contribution to the typical change (δW ) of the momentum width fulfils δW/W -δn/n(0, 0).

Using W ∼ kth ∼ 1/l (0) φ and n(0, 0) ∼ ρ 0 (0)l (0) φ X 0 , we then find δW ∼ (A ) 2 /l (0)
φ . Comparing now this result with the typical variation of the thermal momentum width δ kth ∼ /l (0) φ , we can conclude that the hydrodynamic contribution δW ∼ (A ) 2 /l (0) φ will dominate the thermal contribution if A 1/ √ ; in this case, one would still observe the phenomenon of frequency doubling. If, on the other hand, A 1/ √ , the breathing oscillations will be dominated by the variations of thermal momenta and no frequency doubling will be observed. Accordingly, one can expect the crossover from doubling to no doubling to occur at A cr ∝ 1/ √ , with the proportionality factor to be found from the numerical results. From Fig. 1(c) of the main textwe see that, in the relevant region of A 5, A cr scales as 1/ √ as expected. Since A 1 in the quasicondensate regime, the frequency doubling crossover requires very small quench strengths .

II. NUMERICAL SIMULATIONS USING THE C-FIELD METHODOLOGY

A. The c-field method

The c-field (or classical field) method is a proven approach to studying the equilibrium properties and dynamics of degenerate Bose gases at finite temperature [9]. The crux of the technique is to treat the quantum Bose field ψ(x, t) as a classical field ψ C (x, t), thus ignoring the discrete nature of the particles that make up the field. The classical field approximation captures many features of weakly-interacting 1D Bose gases. For instance, for thermal equilibrium configurations it correctly describes the crossover from the ideal Bose gas regime to the quasicondensate regime [10].

The energy functional of a classical Bose field confined in a harmonic potential is

E({ψ C })= ˆdx E(x), (S9)
where

E(x) = 2 2m ∂ψ C (x) ∂x 2 + 1 2 mω 2 0 x 2 |ψ C (x)| 2 + g 2 |ψ C (x)| 4 -µ|ψ C (x)| 2 . (S10)
Configurations corresponding to thermal equilibrium are obtained from the Gibbs ensemble: the probability of a field configuration ψ C (x) is proportional to exp(-E({ψ C })/k B T ).

Here µ is the chemical potential that fixes the mean particle number.

A convenient method to sample configurations from the Gibbs ensemble is to integrate the projected stochastic Gross-Pitaevskii equation (SPGPE) for times long enough that the memory of the initial state is lost. The SPGPE is

dψ C (x, t) = P C - i L C + κ th (µ -L C )ψ C (x, t)dt + 2κ th T dW (x, t) , (S11)
where

L C ψ C = -2m ∂ 2 ∂x 2 + 1 2 mω 2 0 x 2 + g|ψ C | 2 ψ C , (S12)
dW (x, t) is uncorrelated complex white noise satisfying dW * (x, t)dW (x, t ) = δ(xx )dt, and P C {•} is the projector onto the computational basis. The value of the rate κ th has no consequence for the equilibrium configurations, and hence can be chosen for numerical convenience [11].

At this point we make some comments about the classical field model to present its limitations and provide some physical insight. For high energy modes of the classical field, the interaction energy can be neglected and the energy functional can be approximated by that for noninteracting particles. Given a mode of energy m , the classical field model predicts an energy of k B T , and thus a mean occupation number of n m = k B T /( mµ). Although this expression is a good approximation to the Bose-Einstein distribution for n m 1, it overestimates the population for n m 1, where the Maxwell-Boltzmann distribution is a more appropriate model. In 2D and 3D, this overestimation leads to an ultraviolet divergence of the field density |ψ C | 2 at a fixed temperature. To overcome this problem, one can introduce an energy cutoff: higher-energy modes are treated as an ideal Bose gas [12][13][14][15], while the classical field is restricted to the lowenergy modes. With this separation of the field into a classical and a quantum part (in the sense that the discreteness of particles is important), Eq. (S11) has a physical meaning: it represents the "real" time evolution of the field with κ th an effective collision rate that quantifies the thermal and diffusive damping experienced by particles in the classical field region due to the thermal reservoir (i.e. the quantum region) [9].

Although the inclusion of an energy cutoff is crucial in higher dimensions in order to prevent a divergence of the atomic density, its role is less crucial in 1D. In particular, the 1D classical field predictions for the atomic density do not diverge, even in absence of an energy cutoff, and are quantitatively correct for degenerate gases [16]. The results we present are correctly described solely by classical field theory without a cutoff: for a given peak linear density and temperature, the results do not depend on the cutoff once it is large enough. The cutoff merely determines the size of the basis used for the numerical calculations.

Our simulations were performed within the Hermite-Gauss basis, which is the single-particle eigenbasis of an harmonically-confined ideal gas, and therefore represents the natural, most computationally efficient basis (see Eq. ( 39) of [17] for an explicit expression of the SPGPE in the Hermite-Gauss basis). Since these eigenstates well-approximate the higher-energy, sparsely populated modes, simulation in the Hermite-Gauss basis imposes an energy cutoff that is in direct proportion to the basis size. As discussed above, the variation of this energy cutoff makes little difference to the resulting equilibrium states.

After generating the equilibrium ensemble, the simulations proceed by quenching the trapping frequency ω 0 → ω 1 and evolving the ensemble in time using the simpler projected Gross-Pitaevskii equation (PGPE) [which can be obtained by setting κ th = 0 in Eq. (S11)] [18,19]. This equation then conserves energy and number of particles. As pointed out above, this classical field approximation fails to correctly capture the behavior of high energy modes. In contrast to the equilibrium case, this can affect the dynamics, and in particular the damping rates, and the details can be sensitive to the cutoff [20,21]. Nevertheless, nonequilibrium properties have been studied with some success using this type of classical field approximation, including the collective oscillations of Bose gases [21,22]. In our case, the results of dynamics that we present here are found not to depend strongly on the choice of the energy cutoff, demonstrating that the role of the higher-energy states is relatively unimportant.

For consistency, we chose the same energy cutoff for both the PGPE evolution after the confinement quench and the SPGPE that generates the initial condition. However, the Hermite-Gauss modes that represent the single-particle eigenbasis depend upon the trapping frequency. This implies that the number of Hermite-Gauss modes used for the PGPE evolution, Mcut , is related to the number of modes in the SPGPE evolution M cut via

Mcut = ω 1 ω 0 (M cut + 1 2 ) -1 2 , (S13)
where x denotes the integer component of x.

Figure 1 illustrates PGPE evolution of the position density and momentum distributions after a confinement quench in three different parameter regimes. Although the density undergoes breathing oscillations at frequency ω B √ 3ω 1 in all three cases, the momentum distribution exhibits frequency doubling (top row), a crossover between quasi-doubling and no doubling (middle row), and no frequency doubling (bottom row). This is consistent with the breathing oscillations predicted by our finite-temperature hydrodynamic theory within the quasicondensate regime (see main text).

B. Details for comparison with finite-temperature hydrodynamic theory in the quasicondensate regime

The thermal equilibrium properties of a harmonically trapped Bose gas at temperature T 0 and peak density ρ 0 (0) can be parametrized by the three dimensionless quantities γ 0 = mg/[ 2 ρ 0 (0)], t 0 = 2 2 k B T 0 /(mg 2 ), and ω ≡ (l HO /ξ 0 ) = ω 0 /[gρ 0 (0)] 2 , where l HO = /(mω 0 ) and ξ 0 = / mgρ 0 (0) are the harmonic oscillator length scale and healing length, respectively. However, within the classical field approximation only two parameters are required since features on the order of mean interparticle separation 1/ρ 0 (0) are neglected. Specifically, if the classical field is scaled by ψ 0 = [mk 2 B T 2 0 /( 2 g)] 1/6 , the length scale by x 0 = ( 4 /(m 2 gk B T 0 )) 1/3 , and the trapping frequency by ω = (k B T 0 √ mg/ ) 2/3 / , then the thermal action E({ψ C })/k B T [see Eqs. (S9) and (S10)], which determines the grand canonical partition function, depends only on the dimensionless parameters η 0 = [ 2 /(mg 2 k 2 B T 2 0 ] 1/3 µ and ω 0 /ω. Equivalently, the gas can be described by the dimensionless parameters χ 0 = k B T /[ ρ 0 (0) gρ 0 (0)/m] and ω = ω/[gρ 0 (0)]; field correlation functions of order q have previously been shown to depend only on χ 0 and ω, provided they are scaled to ρ 0 (0) q and the lengths are scaled to 2 ρ 0 /(mk B T ) (cf. [10,23] which investigated the parameter dependence within the classical field approximation for a uniform Bose gas). For sufficiently weak trapping frequencies, ω → 0 and drops out of the problem. This occurs if the size of the atomic cloud is much larger than other microscopic correlation lengths of the gas, therefore implying that the local density approximation (LDA) is valid.

Consider now the post-quench dynamics investigated in this paper. They are a priori parametrized by the dimensionless parameters γ 0 , A = √ 8/χ 0 = √ 8/(γ 3/2 0 t 0 ), , and ω. Within the classical field approximation, if we rescale the PGPE [cf. Eqs. (S11) and (S12)] as done for the equilibrium case, we find that γ 0 drops out. Additionally, we chose parameters for our c-field simulations such that the dynamics only depend upon A and . That is, we required a sufficiently weak ω such that the size of the cloud was always much larger than the typical correlation length l (0) φ = 2 ρ 0 (0)/(mk B T 0 ) (recall that l (0) φ is the typical phase correlation length of the gas, itself larger or on the order of the density-density correlation length). Similarly to the equilibrium case, one then expects a dynamical LDA to be valid, ensuring that the parameter ω is irrelevant. This further fulfilled the high temperature condition required for the (S)PGPE, whilst still ensuring that the number of modes was numerically tractable.

In order to compute A cr for a given quench strength , we fixed t 0 and varied A. For each A, K was extracted by fitting B 1 exp(-b 1 t)[ √ K cos(νt) -√ 1 -K cos(2νt)] + B 2 exp(-b 2 t) to the HWHM of the momentum distribution at each time point t (here B 1 , B 1 , b 1 , b 2 , ν, and K are free fitting parameters); see the red curves of Fig. 1 for example fits. The rates b 1 and b 2 account for the damping present in the PGPE evolution, which is absent from our hydrodynamic theory. This gave a dataset (A, K); the point A cr where K = 1/2 was determined by fitting 1 2 [tanh[a(A -2/3 -A -2/3 cr )] + 1] to this dataset, with free fitting parameters a and A cr .

The equilibrium momentum distributions generated by the SPGPE differ slightly to those generated from Eq. (10) of the main text. This is not unexpected; various assumptions that go into Eq. ( 10) (such as a Thomas-Fermi density profile and the LDA) are relaxed in the SPGPE. A fairer comparison to our finite-temperature hydrodynamic theory is therefore obtained by fitting the equilibrium momentum distribution predicted by Eq. ( 10) of the main text to the c-field equilibrium momentum distribution (with A and ω 0 /[gρ 0 (0)] as free parameters). This shifts the dataset (A, K) → (A , K); the c-field n(k, t)/n(0, 0) values of A cr reported in Fig. 1 of the main text are computed with respect to this latter dataset.

4 n

 4 FIG. 1. (Color online). (a) Breathing mode oscillations and (b) the respective widths (HWHM) of the momentum distribution of a harmonically trapped 1D quasicondensate after a confinement quench as a function of the dimensionless time τ = ω1t. The three examples shown in (a) and (b) correspond, respectively, to: = 0.563, A = 104 (with t0 = 10 6 and ω0/[gρ0(0)] = 3.0 × 10 -3 is c-field simulations) -top row; = 0.101, A = Acr = 5.21 (t0 = 10 3 , ω0/[gρ0(0)] = 1.1 × 10 -2 ) -middle row; and = 0.0203, A = 3.95 -bottom row [34]. The grey triangles are the c-field data [30], with the error bars indicating 95% confidence interval. (c) Nonequilibrium phase diagram of the dynamical crossover from frequency doubling to no doubling in the A-parameter space; data points (circles) show the crossover values Acr for which the weights of the fundamental and the first harmonics are equal (K = 1/2); grey triangles are from c-field simulations. The dashed line is a fit in the region A > 5 with a power law Acr 1.58/ √ (see text), whereas the two dotted lines show the values of A corresponding to K = 0.05 and K = 0.95. (d) Frequency doubling conditions superimposed on the equilibrium phase diagram of the 1D Bose gas [25, 28], drawn in terms of the dimensionless temperature (t0) and interaction strength (γ0), and covering the quasicondensate and the neighbouring nearly ideal Bose gas regimes. The grey dashed line (t0 = γ -3/2 0

FIG. 1 .

 1 FIG. 1. (Color online). (a) Breathing mode oscillations and (b) the respective widths (HWHM) of the momentum distribution of a harmonically trapped 1D quasicondensate after a confinement quench as a function of the dimensionless time τ = ω1t. The three examples shown in (a) and (b) correspond, respectively, to: = 0.563, A = 104 (with t0 = 10 6 and ω0/[gρ0(0)] = 3.0 × 10 -3 is c-field simulations) -top row; = 0.101, A = Acr = 5.21 (t0 = 10 3 , ω0/[gρ0(0)] = 1.1 × 10 -2 ) -middle row; and = 0.0203, A = 3.95 -bottom row [35]. The grey triangles are the c-field data [31], with the error bars indicating 95% confidence interval. (c) Nonequilibrium phase diagram of the dynamical crossover from frequency doubling to no doubling in the A-parameter space; data points (circles) show the crossover values Acr for which the weights of the fundamental and the first harmonics are equal (K = 1/2); grey triangles are from c-field simulations. The dashed line is a fit in the region A > 5 with a power law Acr 1.58/ √ (see text), whereas the two dotted lines show the values of A corresponding to K = 0.05 and K = 0.95. (d) Frequency doubling conditions superimposed on the equilibrium phase diagram of the 1D Bose gas [25, 28], drawn in terms of the dimensionless temperature (t0) and interaction strength (γ0), and covering the quasicondensate and the neighbouring nearly ideal Bose gas regimes. The grey dashed line (t0 = γ -3/2 0

FIG. 1 .

 1 FIG. 1. (Color online). Results of PGPE simulations illustrating the breathing mode oscillations of a harmonically trapped, weakly-interacting 1D Bose gas after a confinement quench. (a) Density; (b) momentum distribution; (c) rms width ∆x of the density; and (d) half-width-at-halfmaximum ∆ kHWHM of the momentum distribution; all as a function of the dimensionless time τ = ω1t. The three examples correspond, respectively, to: = 0.563, A = 104, t0 = 10 6 and ω0/[gρ0(0)] = 3.0 × 10 -3 -top row; = 0.778, A ≈ Acr = 2.3, t0 = 10 3 , and ω0/[gρ0(0)] = 0.07) -middle row; and = 0.1, A = 1.15, t0 = 10 6 , and ω0/[gρ0(0)] = 0.1 -bottom row. In (d) the black points are c-field data, with the error bars indicating 95% confidence intervals, whilst the solid red curves are fits of the form B1 exp(-b1t)[ √ K cos(νt) -√ 1 -K cos(2νt)] + B2 exp(-b2t).

Supplemental Material

Finite-temperature hydrodynamics for one-dimensional Bose gases:

Breathing mode oscillations as a case study I. Bouchoule, 1 S. S. Szigeti, 2, 3 M. J. Davis, 2 and K. V. Kheruntsyan In this supplemental material we provide further details on the hydrodynamic scaling solutions, as well as a brief description and further results of our c-field simulations.