Introducing iFluid: a numerical framework for solving hydrodynamical equations in integrable models

We present an open-source Matlab framework, titled iFluid, for simulating the dynamics of integrable models using the theory of generalized hydrodynamics (GHD). The framework provides an intuitive interface, enabling users to define and solve problems in a few lines of code. Moreover, iFluid can be extended to encompass any integrable model, and the algorithms for solving the GHD equations can be fully customized. We demonstrate how to use iFluid by solving the dynamics of three distinct systems: (i) The quantum Newton’s cradle of the Lieb-Liniger model, (ii) a gradual field release in the XXZ-chain, and (iii) a partitioning protocol in the relativistic sinh-Gordon model. Copyright F. S. Møller and J. Schmiedmayer. This work is licensed under the Creative Commons Attribution 4.0 International License. Published by the SciPost Foundation. Received 09-01-2020 Accepted 10-03-2020 Published 13-03-2020 Check for updates doi:10.21468/SciPostPhys.8.3.041


Introduction
In recent decades great experimental advances in the field of ultracold atoms have enable the preparation and manipulation of one-dimensional many-body quantum system far from equilibrium [1][2][3][4][5][6][7][8][9]. Therefore, theoretical tools for understanding the complex dynamics of these systems have been highly sought after [10][11][12]. Some of these low dimensional systems exhibit integrability by abiding to an infinite set of local conservation laws [13]. For a long time, integrable models have been a theoretical playground, although several of these models have also been realized experimentally [1,6,7,9].
Recently the theory of generalized hydrodynamics (GHD) has emerged as a powerful framework for studying integrable models out of equilibrium [14,15]. In its most basic form, generalized hydrodynamics describes the flow of all the conserved charges of integrable models. Thus, an infinite set of advection equations emerge, which through the thermodynamic Bethe ansatz, can be formulated as a single Euler-scale equation for a quasiparticle distribution. Since the inception of GHD several applications have been added to the framework, such as calculations of entanglement spreading [16][17][18][19], correlation functions [20][21][22], diffusive corrections [23,24], and many others [25][26][27][28]. Recently it has also been demonstrated to capture the dynamics of a cold Bose gas trapped on an atom-chip [29]. An especially appealing feature of the GHD framework is how the main equations can be universally applied to a large set of integrable models including the Lieb-Liniger model [14,[30][31][32], XXZ chain [15,[33][34][35], classical [36] and relativistic [14] sinh-Gordon, and many more [37][38][39][40]. The theory has already proven its worth by providing exact predictions for the many-body dynamics in several cases [14,15,20]. Additionally, GHD appears to have great potential as a numerical tool, as the computational complexity of solving the many-body dynamics is entirely independent of the Hilbert space size. Despite this, only a couple of different numerical schemes have been implemented so far [27,28,34,38,41]. Thus, if GHD is to applied on larger scales, such as describing experimental observations, more powerful numerical methods must be developed.
Currently, no open-source code exists for solving GHD equations. The goal of iFluid (integrable-Fluid) is to provide a powerful and intuitive numerical framework for finding and propagating the root density distribution, which serves as the basic quantity for thermodynamic calculations in integrable systems [42]. Hence, iFluid supplies a platform for solving the core hydrodynamical equations on top of which user-specific applications can be built. The universality of the GHD equations enables a highly flexible code base, wherein any integrable model can be seamlessly integrated. iFluid already supports a couple of models (see Appendix A), and the implementation of a new model can be achieved in relatively few lines of code after extending core classes of the framework.
So far only little effort has been put into comparing different algorithms for solving the GHD equations. For very cold temperatures the underlying quasiparticle distribution resembles that of a Fermi sea, whose hard walls and edges complicates numerical solutions. However, such problems have already been studied for many years within the field of fluid dynamics [43]. Adopting these methods should greatly bolster the numerical capability of generalized hydrodynamics. Therefore, iFluid abstracts the algorithm for solving the main advection equation, whereby users are free to implement whatever algorithm is suitable for their specific problem. Again, iFluid already implements a couple of basic algorithms sufficient for solving most tasks.
The purpose of this paper is to introduce the user to the iFluid framework and demonstrate its applicability in various scenarios. Thus, the paper is organized as follows. In Section 2, we review the basic concepts of GHD serving as the core of the iFluid framework. In section 3, the core features of iFluid are discussed. In Section 4, the applicability of iFluid is demonstrated by solving three distinct problems: The quantum Newtons cradle in the Lieb-Liniger model, a gradual confinement release in the XXZ model, and a partitioning protocol in the relativistic sinh-Gordon model. Finally in Section 5, we conclude and give an give an outlook for the development of iFluid. Details of the exact numerical implementations are reported in appendices.

Review of generalized hydrodynamics (GHD)
Generalized hydrodynamics in its essence utilizes the quasiparticle formulation of the thermodynamics Bethe Ansatz to describe the flow of charges within an integrable system. Integrable systems abide to infinitely many local conservation laws [44], thus preventing a conventional hydrodynamical description which only captures conservation of particles, momentum, and energy. This infinite set of conserved charges imposes constraints on the dynamics of the system and inhibits thermalization. Hence, under the assumption of local thermal equilibrium, the systems relaxes to a generalized Gibbs ensemble (GGE) from which thermodynamic quantities can be derived [45]. Once it is at this stage, the system can be described via the quasi-particles based thermodynamic Bethe ansatz (TBA). Within TBA the eigenstates of the full set of local conserved charges are multiparticle states, with each particle labelled by a rapidity λ [13]. Under integrability all multiparticle scattering events factorize into two-body elastic scatterings. Thus, all interactions between the quasiparticles are captured by the twobody scattering matrix S(λ). In the thermodynamic limit the rapidity can be thought of as a continuous variable, while the coarse-grained root density ρ(λ) gives the density of particles within the interval [λ, λ + dλ) [42]. The root densities (like the GGE density matrix,ρ GGE ) fix the expectation value of the the local charges,Q j , such that [46] 1 where h j (λ) is called the one-particle eigenvalue of the chargeQ j , and L is the system length. Among the infinite set of conserved charges we find the particle numberQ 0 =N , the total momentumQ 1 =P, and the total energyQ 2 =Ê. Thus, the corresponding one-particle eigenvalues are h 0 (λ) = 1, h 1 (λ) = p(λ), and h 2 (λ) = ε(λ), where p(λ) and ε(λ) are the momentum and energy of a single quasiparticle respectively. In a similar fashion we can calculate the expectation values of the currents associated to the charges via where v eff (λ) is the velocity by which the quasiparticles move. Later we will see how exactly this velocity is computed.
Until recently, the thermodynamic Bethe ansatz was used only to describe the expectation values of a homogeneous, stationary state. Imagine however, a weakly inhomogeneous system, where physical properties vary on space-time scales much larger than the underlying time scales. In this case local equilibrium is established faster than the physical quantities can change, whereby the overall system remains in a quasi-stationary state. Thus, the system can be thought of as consisting of space-time fluid cells, each of which is described by a local GGE with minimal variations to those of neighbouring cells [14]. For a quasi-stationary state the quasiparticle description of the thermodynamic Bethe ansatz still applies, however, the root density is now weakly dependent on time and space, ρ = ρ(t, x, λ). The main feature of GHD is formalizing how the root density behaves under time evolution. Thus, the GHD details how the flow of an infinite set of charges of an integrable model is given by the semiclassical propagation of a phase-space density of a quasiparticle collection.
In practice, it is more convenient to express the hydrodynamical equations in terms of the filling function, ϑ(λ). The filling can be interpreted as the density of quasiparticles over the density of available states at a given rapidity, ρ s (λ), such that Note that we have omitted the spacial and temporal argument for the sake of lighter notation. The subscript dr in eq. (3) denotes the dressing of the quantity. Non-trivial interactions between the particles induces collective effects, which are captured by the dressing operation defined through the integral equation Due to the factorization of scatterings, the interparticle interactions are captured solely by the two-body scattering phase Θ(λ) = −i log S(λ), with S(λ) being the two-body scattering matrix. The interaction between particles also influences their equations of motion. In the non-interacting case the particles move with the group velocity, v(λ) = ∂ λ ε/∂ λ p. However, in the presence of interactions the bare quantities become dressed, whereby the particles move with an effective velocity [14,15] v eff (λ) = (∂ λ ε) dr Furthermore, space-time inhomogeneities in the parameters of the model induce force terms on the quasiparticles, which can change their rapidities. Once again the interparticle interactions collectively dress these force terms, whereby the quasiparticles experience an effective acceleration [27] where α is a coupling of the model (such as the interaction strength c in the Lieb-Liniger model [30]). Inhomogeneities of the couplings in space and time have their own associated force term given respectively by and Finally, the evolution of the phase-space quasiparticle density is captured with the single hydrodynamical equation [26,27] where the brackets explicitly indicates that the velocity and acceleration is dependent on the current state. Eq. (9) is a simple Eulerian fluid equation, which describes the flow of the infinite set of conserved charges through a single expression. It is the main equation of GHD along with its root density-based counterpart which identifies the root density as a conserved fluid density [14]. From a numerical perspective, eq. (9) is more convenient to work with than eq. (10).
While the thermodynamic states of the Lieb-Liniger model are characterized by a single root density, this is in general not the case. For instance, the XXZ chain supports bound states, which are captured by including multiple types of quasi-particles, each with their own corresponding root density, ρ k (λ). Thus, one must sum over the contribution from each quasi-particle type to obtain the charge densities as described in eq. (1). To simplify the notation we adopt the convention from [26] of writing the rapidity as a single spectral parameters λ = (λ, k), whereby the integrals above can be generalized to After accounting for multiple species of quasiparticles, all the equations above can be applied to any integrable model. In fact, the only model-specific parameters that enter the calculations is the scattering phase, Θ(λ), encoding the interactions of the quasiparticles and the oneparticle eigenvalues, h j (λ). These quantities can be obtained for a given model through the thermodynamic Bethe ansatz, and once plugged into the hydrodynamical equations the full framework of GHD can be applied to the problem. The equations above constitutes the core of generalized hydrodynamics. Under evolution detailed by eq. (9) the system is always in a quasi-stationary state, from which various physical quantities can be calculated. In addition to solving eq. (9), iFluid also computes the so called characteristics [41] encoding the trajectories of the quasiparticles. Thus, the characteristics can be used for computing the hydrodynamics spreading of entanglement [19] and correlations [20]. The characteristics U and W have the simple interpretation as the inverse space and rapidity trajectories of the quasi-particles respectively [19], yielding Thus, U(t, x, λ) is the position at time t = 0 of the quasi-particle λ, which at time t has the position x [26]. Interestingly, the characteristics follow the same hydrodynamical equation as the filling function with the initial conditions given by definition Hence, the characteristics can be propagated alongside the filling function with minimal numerical cost.

Core functionality of iFluid
iFluid implements all the equations listed in the section above along with several additional features. The main goal of the framework is to remain highly flexible while delivering fast performance. To achieve this goal, iFluid implements all its core methods in abstract classes, which must be extended by the user in order to supply the necessary model-specific functions required by the internal routines of iFluid. Several models and numerical solvers are already implemented in iFluid and will be demonstrated in Section 4. This section aims to introduce the iFluid base classes, while more in-depth information can be found in the documentation [47].
The hydrodynamical equations described in Section 2 are solved numerically by discretizing the integrals using appropriate quadratures. Thereby the linear integral equations are converted into linear matrix equations enabling very fast numerical calculation of the hydrodynamical quantities. One should note that iFluid employs a very strict convention for indices which is enforced through the iFluidTensor class. The discretized equations are found in Appendix B, while further information regarding the iFluidTensor is written in the documentation [47].

Implementing a model
A key feature of iFluid is its intuitive interface and extendibility. This is achieved through the abstract class iFluidCore, which implements all the equations of the thermodynamics Bethe ansatz from the previous section. Following the hydrodynamical principle the system is always in a quasi-stationary state, whereby all the methods of the iFluidCore can be applied for any given root density. However, in order to perform any specific calculations some explicit information regarding the model is required, namely the energy and momentum functions, ε(λ) and p(λ), and the two-body scattering phase, Θ(λ − λ ), along with their rapidity derivatives. Additionally, inhomogeneous systems require derivatives with respect to the couplings in order to compute eqs. 7 and 8. Hence, before any calculation can be undertaken one must extend the general iFluidCore class with a model-specific class myModel < iFluidCore, wherein the following abstract functions must be implemented 1 % Abstract methods which must be overloaded in model class 2 getBareEnergy(obj, t, x, rapid, type) 3 getBareMomentum(obj, t, x, rapid, type) 4 getEnergyRapidDeriv(obj, t, x, rapid, type) 5 getMomentumRapidDeriv(obj, t, x, rapid, type) 6 getScatteringRapidDeriv(obj, t, x, rapid1, rapid2, type1, type2) 7 getEnergyCoupDeriv(obj, coupIdx, t, x, rapid, type) 8 getMomentumCoupDeriv(obj, coupIdx, t, x, rapid, type) 9 getScatteringCoupDeriv(obj,coupIdx,t,x,rapid1,rapid2,type1,type2) A few things to note about the input parameters: t, x and rapid are the spatial, temporal and rapidity arguments respectively. They can be either scalars or vectors. The type argument specifies the index of the quasiparticle type (only relevant for TBAs with multiple quasiparticle species). This argument can be either a scalar or an array of scalars and should default to 1 for single-particle TBAs. Lastly, the coupIdx input is a scalar index specifying which coupling the derivative is taken with respect to. The couplings are passed to the model-specific class through an array upon initialization (more on this later), whereby the array-index of the given coupling must match the coupIdx.
Although this might seem like a lot of work at first glance, most of these functions can be written in a single line. Examples of this can be found in the documentation [47].

Solving the GHD equation
Having specified the model-specific functions, all the equations listed in Section 2 except the main GHD equation (9) can be solved. Solving eq. (9) is achieved through the abstract iFluidSolver class. As previously mentioned, iFluid abstracts the algorithm for solving eq. (9) in order to accommodate future advances in the GHD numerics. Algorithms already implemented in iFluid can be found in the documentation [47], while new algorithms can be added simply by extending the iFluidSolver class and implementing the abstract methods 1 % Abstract methods which must be overloaded in sub−class 2 step(obj, theta _ prev, u _ prev, w _ prev, t, dt) 3 initialize(obj, theta _ init, u _ init, w _ init, t _ array ) The method step() has the simple function of advancing the filling function by a single time step, d t, following eq. (9) Several different approaches exists for taking this step. The solvers already implemented in iFluid utilize the implicit solution of eq. (9) [27,34] where the functionsx(t , t) andλ(t , t) are given bỹ andλ The subscript of the effective velocity and acceleration denotes that the dressing is taken with respect to the filling function at time τ. Further, note that the functionsx(t, 0) andλ(t, 0) are in fact the characteristics U(t, x, λ) and W (t, x, λ) respectively. The step() function in iFluids FirstOrderSolver and SecondOrderSolver approximates solutions of eqs. (19) and (20) for a single time step at various orders.
Some algorithms require the filling function at only a single time to perform the step above, while others need the filling at multiple times in order to perform a more accurate step. For example, the class SecondORderSolver [27] employs a midpoint rule, whereby the midpoint filling function is stored within the class. However, in order to take the first step, the class needs to know the first midpoint, which is not given a priori. Thus, one must also implement the method initialize(), which prepares all the quantities necessary for starting the time evolution algorithm.
Once the abstracted methods above are implemented, one can solve eq.
where each filling is stored as an iFluidTensor. In order to implement the abstract methods listed above, one will need some of the hydrodynamical equations listed in Section 2. Thus, the iFluidSolver constructor takes an iFluidCore object as argument and stores it. Whenever the dressing of a quantity (or something similar) is needed, one simply calls the appropriate method from the stored iFluidCore object.

Solving problems with iFluid
The previous section demonstrated how to implement models and algorithms in iFluid. Once the specific model and algorithm is implemented the hydrodynamical calculations become almost trivial, as most problems can be formulated in only a couple of lines of code. In the following examples we solve three distinctly different hydrodynamical problems using the already implemented methods of iFluid. The example codes can all be found on the iFLuid git page [48] and can be run in a matter of minutes on a laptop.

Quantum Newton's cradle
The original experimental realization of the quantum Newtons cradle [1] beautifully demonstrated integrability in a one-dimensional Bose gas. Recently, generalized hydrodynamics was utilized in a numerical study of the experiment [28], where the numerical results were obtained using the flea gas algorithm first described in [38]. Here we simulate a similar scenario of a Bose gas oscillating in a harmonic confinement. In contrast to previous studies which used an initial Bragg pulse to imprint different momenta unto the system, we simply displace half of the cloud with regards to the center of the trap akin to lifting a bead in the classical cradle. The initially displaced cloud will oscillate back and forth in the harmonic trap for several periods, thus demonstrating the lack of thermalization in integrable models. Furthermore, by keeping part of the cloud in the center we clearly illustrate the effect of the interparticle interaction, as the central cloud will be distorted upon overlapping with oscillating one.
The one-dimensional Bose gas is described by the Lieb-Liniger model with the Hamiltonian [30,31] The interaction strength, c > 0, and the chemical potential, µ, are the two couplings of the model, while m is the particle mass and L the system length. By employing an inhomogeneous chemical potential we can describe external traps through the local density approximation. First we need to specify the problem at hand, namely the discretization grids, the couplings and the temperature. For this example we employ a rectangular quadrature for solving the integrals, whereby the quadrature weights are simply the grid spacings. Next, we have to specify the dynamical couplings used in the simulation. The couplings and their derivatives are declared as a cell array of anonymous functions with time and space arguments. The first row specifies the raw couplings, while the second and third row contains the temporal and spatial derivatives respectively. The class LiebLinigerModel requires the first column of the coupling array to be the chemical potential and the second column to be the interaction strength. The chemical potential is given by some central value, µ 0 , minus the harmonic confinement, while the interaction is simply set to unit: Note that all operations in the anonymous functions should be elementwise (signified by the dot-prefix). Furthermore, entries for derivatives equal to zero can be left empty for a boost in performance.
Having specified the problem, we now turn to calculating the initial state given by two clearly separated, identical clouds. To illustrate the interaction between the two clouds, any distortion of the central cloud should be caused by the interactions. Hence, the initial state of the central cloud should be stationary with respect to the harmonic trap. This can be achieved in several fashions: Here we simply create an initial "double well", by displacing a copy of the harmonic trap via a heaviside function. In this case, the initial couplings are very different from the dynamical ones. Therefore, the method calcThermalState() of the iFluidCore class takes an initial set of couplings as optional argument, whereby: Finally, we are ready to solve the GHD equation (9), using the SecondOrderSolver class [27].
Simply pass the model to the solver and run the simulation through the propagateTheta() method: The final output is a cell array, where the i'th entry is the filling function, ϑ(x, λ), at the time t_array(i). Note, each ϑ(x, λ) is stored as an iFluidTensor.
The 27 lines of code above is all there is needed for specifying and solving a typical problem in iFluid. According to the hydrodynamic principles, the system is in a quasi-stationary state at every point in time. Hence, once the filling function is computed, it can be used for any calculation within the thermodynamic Bethe ansatz.
We start out by illustrating the motion of the two clouds of Bose gases in the Newton's cradle by calculating the linear (atomic) density corresponding to the zeroth charge density in eq. (1). Given the filling function, we can also calculate the root density, ρ(t, x, λ), thus illustrating the motion of the quasiparticles. The interactions between the oscillating and stationary cloud transfers momentum between them. At the end of the evolution, the system moves as a single cloud exhibiting both a center-of-mass and breathing motion. Below: Snapshots of the root density at each period. Initially the two clouds have distinct root densities, which gradually merges into a single, binary distribution. 28 charge _ idx = 0 % linear density is 0th charge density 29 n _ t = LL.calcCharges(charge _ idx, theta _ t, t _ array) 30 rho _ t = LL.transform2rho(theta _ t, t _ array) The evolution of the linear density is plotted in Figure 1 along with the root density at selected times. In the root density picture we initially see two blobs well separated in space and centred on zero rapidity. The central blob is the stationary state of the harmonic trap and thus remains in place, while the offset blob is accelerated by the harmonic confinement. This causes the offset root density to encircle the central one, effectively resulting in an oscillating motion of the density. Every time the two clouds overlap interactions occur, effectively transferring a small portion of the oscillating clouds momentum to the central one. Their total interaction is primarily determined by two things: the interaction strength, c, and the amount of time at which the clouds overlap. By separating the clouds by only a small distance, the offset cloud will only accelerate to a low velocity before passing the stationary cloud. Thus, the overlap time becomes long leading to a large distortion of the root densities. Therefore, the two blobs partially merge after merely a few oscillations, creating a binary system orbiting the center while rotating around itself. In the density picture this produces a single cloud whose center of mass oscillates in the harmonic trap while the cloud itself exerts a breathing motion. Figure 1 clearly demonstrates the collective interactions in generalized hydrodynamics. In a non-interacting theory, the offset cloud would simply have encircled the center forever without any deformations of the root densities. Meanwhile, any non-integrable system would have rapidly thermalized producing a single, Gaussian cloud.
Next, we wish to calculate the characteristic, U and W of eqs. (12), and illustrate their interpretation. This is easier if the two blobs of quasiparticles stay separate, which can be achieved simply by starting them further apart, thus decreasing their effective interaction over Just like the filling, the characteristics are returned as cell arrays of iFluidTensor. According to eq. 12, the filling function after some evolution time can be inferred from initial state via the characteristics. To demonstrate this we interpolate the initial filling to the characteristics at some time t and compare with ϑ(t). Performing the interpolation is straightforward in Matlab:  36 37 % interpolate theta _ init to (U(t _ final) , W(t _ final)) 38 theta _ UW = interp2( x _ grid, rapid _ grid, plt(theta _ init), U(:), W(:) ) 39 theta _ UW = reshape( theta _ UW, length(rapid _ grid), length(x _ grid) ) The resulting filling function is plotted in Figure 2 along with the characteristics. Starting with the characteristics we observe a spiral pattern, which provides new information about the quasiparticle trajectories not accessible from the root densities themselves. Although the central blob is stationary with respect to the harmonic confinement, its quasiparticles still have a finite velocities causing them to move in an orbit around the center of the blob. In fact, in the case of a harmonic confinement and no interactions, all the quasiparticles move in closed orbit around (x, λ) = (0, 0). However, the interactions between the two clouds distort the quasiparticle trajectories resulting in the observed spiral patterns. Interpolating the initial filling to the characteristics does indeed reproduce the filling function at time t as seen in Figures 2. The small differences between the two fillings are due to inaccuracies in the interpolation and the finite number of grid points, making it very hard to resolve the fine structure of the characteristics. The interpretation of the characteristics can be further understood by plotting them as function of time. Going back to definition in Section 2, we recollect that the characteristics encode the positions and rapidities of the quasiparticles at an earlier time. Note that this is not the trajectory of the quasiparticle but rather the inverse trajectory. Figure 3 depicts the filling at different times along with characteristics of the quasiparticles in the center of the offset blob as a function of time. Starting with Figure 3.(a), the bullets mark represent the coordinates at time t − τ of the quasiparticles now located in the center of the offset blob. The characteristics depict a circular motion, as the two blobs have yet to overlap. However, as the two cloud pass through each other, the quasiparticle trajectories become distorted, as seen in Figure 3.(b). This becomes especially clear when looking at the point t = τ, which clearly has moved. Hence, due to the distortion of the trajectories, a different quasiparticle of the initial root density can now be found in the center of the blob. Every time the two clouds overlap the trajectories become increasingly distorted, which in the end produced the spiral patterns observed previously in Figure 2. Thus, the characteristics do indeed encode the original location of the quasiparticles, which can be used to infer correlations of entanglements of the system [19,20].

Charges and currents of XXZ model
The Heisenberg XXZ model is another integrable model already implemented in iFluid. It has the Hamiltonian [42] whereŜ σ j are the standard spin-1 2 operators, while the couplings are the magnetic field, B, and the interaction, ∆. The model differs from the Lieb-Liniger model in a number of different ways. In the case of ∆ ≥ 1 the quasiparticles are restricted to the first Brillouin zone in the rapidity. Furthermore, bound states within the chain can occur, which requires a TBA description consisting of multiple root densities, i.e. multiple types of quasiparticles otherwise known as strings. For B < 0 infinitely many root densities are needed in theory, however, in practice one can truncate this to a relatively small number, as each additional root density has diminishing effect.
In this example we examine a system initially confined by a strong magnetic field, which afterwards is gradually decreased. Setting up the problem in iFluid is very similar to the previous example. Since the expression for the coupling is a little longer in this case, we use Matlab's symbolic feature to take the derivative of the coupling for us. Furthermore, we solve the TBA integrals using a Legendre-Gauss quadrature obtained via the legzo() method [49].   In the final two lines of the code we calculate the expectation values of the zeroth and second charges and associated currents. We wish to calculate the kinetic energy, thus the energy without the contribution from the magnetic field. Hence, we simply set the field to zero before the calculation. Figure 4 shows the calculated quantities at select times. Starting with subfigure 4.(a) we see the evolution of the linear density. Initially the density profile consists of a single, smooth curve dictated by the parabolic, confining magnetic field. However, the temporal change in the coupling induces force terms on the quasiparticles, and each string experiences a different effective acceleration. Thus, the strings starts moving outwards at different velocities, whereby three distinct density profiles become visible after some time. These three profiles correspond Figure 4: Expectation values of charge densities and their associated currents as function of x for the XXZ model. The system is initialized in a parabolic, confining magnetic field, whereafter the field is gradually lowered. (a,c) Linear density and its current at different times. As the field is lowered the difference in velocity between the different quasiparticles becomes apparent. (b,d) Kinetic energy density and current at different times. The higher kinetic energy of first order quasiparticles hides the contribution from higher orders.
to the three root densities accounted for in the calculation. This further emphasises the point made earlier that each additional root density included has diminishing effect. The same threepart structure can be seen in the associated current, albeit to a much lesser degree. Meanwhile, the expectation value of the kinetic energy barely changes by taking the contribution from additional strings into account, as the string corresponding to no bound states has the largest kinetic energy.

Partitioning protocol in relativistic sinh-Gordon
In our final example we examine the relativistic sinh-Gordon model with the Hamiltonian [50,51] where the constant c is the speed of light, and : • : denotes normal ordering with respect to the ground state. The mass-parameter β is related to the physical renormalized mass m by m 2 = β 2 sin(απ) απ and α = c g 2 8π + c g 2 .
In the iFluid implementation of the model, the couplings are given by the renormalized interaction, α, and the parameter β. Additionally, iFluid adds a chemical potential to the Hamiltonian as a third coupling, such that users can control the initial linear density.
With this example we wish to illustrate how iFluid deals with extensive systems. Thus, we explore a well-known protocol in the GHD community [14,15], namely a partitioning protocol where two semi-infinite, homogeneous systems, or leads, of different temperature are joined together. Since iFluid uses finite-sized grids, we obviously cannot store an infinitely long system. However, we can toggle the option extrapFlag = true of the iFluidSolver to enable extrapolation of the filling functions upon propagation. Usually extrapolation is ill advised, but in this case the extrapolation works well since each lead is homogeneous.
There are several ways of realizing the two leads. One option is declaring a space-dependent temperature and then balancing the difference in density via the chemical potential. Achieving exactly equal density of the two leads can be tricky using this approach, however, the example here is merely a demonstration of using the model. In addition to performing the regular partitioning protocol we also gradually increase the interactions, such that α(t) = (1 + 0.5 tanh(2t))/(8π + 2), while β = 1. By now the process of setting up grids and couplings in iFluid should be known to the reader, so we merely show the part of the code unique to the problem: At the end of the calculation the linear density is calculated along with the expectation value of the vertex operator [51], Φ k (x) ≡: e kgφ(x) :, which are both plotted in Figure 5.
Starting with the linear density we confirm the density being homogeneous initially. However, the higher temperature of the right lead causes a larger number of quasiparticles to be initialized at higher rapidities, whereby they move at higher velocity. Thus, quasiparticles from the right travel into the left lead faster than the quasiparticles from the left can fill out the void left behind. Therefore, one observes a wave of increased density travelling left and another wave of decreased density travelling right.
The change in interaction, α, has barely any influence on the redistribution of density. However, it greatly influences the expectation values of the vertex operators, since the expression for 〈Φ k 〉(x) is directly dependent on α [51]. Thus, as the interaction increases, so does 〈Φ k 〉(x). Even the regions not yet reached by the density wave are affected, since α is changed globally.

Conclusion
We have demonstrated that iFluid enables the user to perform state of the art GHD calculations in only a few lines of code. Additionally we have shown that iFluid can be easily extended to encompass a large number of integrable models and numerical solvers. We hope that iFluid will be a help to both students and researchers, who wish to explore the numerical side of generalized hydrodynamics. Furthermore, the recent experimental evidence of GHD's ability to describe the dynamics of cold gas experiments [29] further increases the need for powerful numerical tools. Thus, iFluid represents a great step forward in making the theory more widely accessible, since no open-source software exist in the GHD community so far.
Aside from being easy to use, iFluid also offers great extendibility through its abstract classes. Thus, users can implement new models simply by extending the iFluidCore class and overloading a couple of methods. Similarly, new solvers of the GHD Euler-equation (9) can be added to the framework by extending the class iFluidSolver. Well-established algorithms from the field of fluid dynamics can thereby be seamlessly added and tested. The development of iFluid is an ongoing process, as more and more advances are being made in the theory of generalized hydrodynamics. By employing a modular layout, the framework aims to function as a fundamental platform on which further tools can be built upon. Applications for calculating the hydrodynamical spreading of correlations and entanglement seem especially promising.
The current version of iFluid is written is Matlab, however, plans are currently in the works to write the framework as either a Python package or a C++ library. While each language has its own advantages, the Matlab iteration of iFluid is easily accessible to most members of the GHD community while retaining decent performance. We also welcome anyone interested in contributing to the project to contact the authors through either email or the official iFluid repository on Github.

A Thermodynamic Bethe ansatz of implemented models
The thermodynamic Bethe ansatz is a textbook technique nowadays [42], which can be applied to a large range of integrable models. In this section we report the basic TBA quantities of the models highlighted in the main text and emphasise the details of the iFluid implementation.
The iFluid implementation of the Lieb-Liniger model is contained in the LiebLinigerModel class, which takes the chemical potential as the first coupling and the interaction strength as the second coupling. In addition to the standard TBA equations, the LiebLinigerModel class implements additional methods: 1 % Given external potential V _ ext, fit mu to get given number of atoms 2 fitAtomnumber(obj, T, V _ ext, Natoms, mu0 _ guess, setCouplingFlag) 3 4 % Calculate the n−body local correlator g _ n 5 calcLocalCorrelator(obj, n, theta, t _ array) The first method fitAtomnumber() fits the central chemical potential, µ 0 , where µ = µ 0 − V e x t (x), to achieve a thermal state whose root density integrates to a specified number of atoms. This is especially useful for experimental comparisons. The second method calcLocalCorrelator() computes the local n-body correlator through the approach detailed in [21].
Recently generalized hydrodynamics was shown to describe the dynamics of an experimentally realized one-dimensional Bose gas. Hence, the theory appears to be a powerful tool for simulating real systems. Therefore, iFluid implements a wrapper class LiebLinigerModel_SI for converting inputs in SI-units to the internal units of iFluid, which then calls the appropriate methods of the LiebLinigerModel class.
For more detailed description of these methods we refer the reader to the official iFluid documentation [47].

A.2 XXZ spin chain model
The XXZ spin chain model is a discrete integrable model of N sites with the Hamiltonian H = Here the standard spin-1 2 operators areŜ σ j . while B j denotes the magnetic field at site j and ∆ j the interaction. Although the model is discrete, it is treated exactly like the continuous models in the hydrodynamical description.
The iFluid implementation of the model, XXZchainModel, takes B as the first coupling, while the angle θ = arccosh ∆ serves as the second coupling.
Like the Lieb-Liniger model, the thermodynamic Bethe ansatz is determined by only a single root density. The TBA functions of the model read ε(λ) = m cosh λ − µ , p(λ) = m sinh λ , Θ(λ) = i log sinh λ − i sin(απ) sinh λ + i sin(απ) , where we have added a chemical potential, µ, to the energy function.
The iFluid implementation of the model, sinhGordonModel, takes α as the first coupling, β as the second coupling, and µ as the third coupling. In the addition to the standard TBA functions, the sinhGordonModel class also implements methods for calculating the expectation values of vertex operators using the approach detailed in [51].