Mesoscopic impurities in generalized hydrodynamics

This study focuses on the impurities in integrable models from the viewpoint of generalized hydrodynamics (GHD). An impurity can be thought of as a boundary condition for the GHD equation, relating the states on the left and right sides. It was found that, in interacting models, it is not possible to disentangle the incoming and outgoing states, which means that it is not possible to think of scattering as a mapping that maps the incoming state to the outgoing state. A novel class of impurities, dubbed mesoscopic impurities, was then introduced, whose spatial size is mesoscopic (i.e. their size Lmicro≪Limp≪L is much larger than the microscopic length scale Lmicro but much smaller than the macroscopic scale L). Because of their large size, it is possible to describe mesoscopic impurities via GHD. This simplification allows one to study these impurities both analytically and numerically. These impurities show interesting non-perturbative scattering behavior, such as the nonuniqueness of the solutions and a nonanalytic dependence on the impurity strength. In models with one quasi-particle species and a scattering phase shift that depends only on the difference in momenta, the scattering can be described using an effective Hamiltonian. This Hamiltonian is dressed due to the interaction between the particles and satisfies a self-consistency fixed-point equation. In the example of the hard-rod model, it was demonstrated how this fixed-point equation can be used to find almost explicit solutions to the scattering problem by reducing it to a two-dimensional system of equations that can be solved numerically.


Introduction
Integrable models are an important set of fine-tuned (mostly one-dimensional) models in theoretical physics which can be solved exactly.This allows for the analytical study of the effect of interactions in these quantum or classical many-body systems in much greater detail than it would be possible in generic models.At first, it may seem that these fine-tuned models are not necessarily good models for actual physical systems.In particular, they show the unusual feature that they posses infinitely many (local) conserved quantities, while a generic system typically has only 3: particle number, momentum and energy.For instance, this manifests in the equilibrium the system approaches after a long time: it is characterized by a generalized Gibbs ensemble (GGE) [1,2], which takes into account all conserved quantities instead of the usual Gibbs ensembles.This failure of thermalization has been observed in experiments, such as the famous quantum-Newton cradle experiment [3] and others [4,5].However, in one dimension, also non-integrable systems often show effects linked to integrability (like long equilibration times, the most famous example here is the FPUT system [6]).The reason for that is that often these systems are close to an integrable model and therefore for short times they can be treated as integrable models (this also explains why experiments -which are never precisely fine-tuned -still can observe GGE's).Only after waiting for long times eventually integrability breaking effects will be observed.Therefore understanding the physics of integrable models provides a first step to understanding the physics of general many-body systems in one dimension.This idea has sparked huge interest in integrable systems, and in particular there has been plenty of work on how to perturb these systems with integrability breaking perturbations, see for instance [7][8][9][10][11].
In the last decade the progress on understanding integrable models was further fueled due to the establishment of a powerful theoretical toolbox: generalized hydrodynamics (GHD) [12,13].Generalized hydrodynamics is a framework to study large scale dynamics of integrable models (for a pedagogical introduction, see [14]).This provides access to the out-of-equilibrium dynamics of these models, which is particularly useful to compare with experiments [15][16][17][18].GHD is a coarse-grained theory which describes the state of a system locally only by its conversed quantities, all other degrees of freedom are averaged out.This provides an ideal starting ground for studying integrability breaking, since integrability breaking manifests through the breaking of conservation laws.In order to be able to account for an integrability breaking effect in GHD it has to be small in some sense.Small could for instance mean a) that the coupling to the integrability breaking terms is weak.A perturbative treatment typically gives rise to Boltzmann-type equations and has been used to describe various experimentally relevant effects, like the effect of particle losses [19,20], the effect of being only a quasi 1D system [21] and an integrability breaking background [22].In GHD it is also possible to go beyond weak perturbations and to consider certain strong perturbations, which b) vary slowly in space and time.This includes the experimentally important case of an external trapping potential [22][23][24] or also slowly varying interactions [25].
Another option to perturb an integrable model is c) to introduce impurities: Impurities can be potentially strong perturbations, but they will only affect the system in a small localized region.Therefore, the system can still be described by GHD away from the impurity, but at the position of the impurity there will be a non-trivial boundary condition relating the state to the left and the right of the impurity.Unlike other forms of integrability breaking, it has been observed they not lead to thermalization at late times [26][27][28].Impurities have been long studied in the context of integrable models both analytically [29][30][31][32][33][34][35] starting from the infamous Kondo impurity (see for instance the review [36]) and experimentally [37][38][39].
In this paper we study the out-of-equilibrium dynamics of integrable models in the prescence of an impurity from a generalized hydrodynamics perspective.The situation we have in mind is quite general: Given the initial state ρ(t = 0, x, p) at time t = 0 we want to predict the state ρ(t, x, p) at some later time t taking into account the impurity.First, we formalize and unify ideas already mentioned in the literature [28,40,41] to establish a coherent picture for the incorporation of impurities into GHD: The boundary conditions correspond to non-equilibrium stationary states (NESS) of the impurity.In our discussion we will find that the classification of impurities can be quite peculiar: Due to their interaction the reflected particles affect the effective velocities of the incoming particles, implying that it is not possible to think of the outgoing state as a function of the incoming state (unless the model is free) ‡.Because of that it is also not clear whether a boundary condition will always exist or whether it is unique (in fact we will give an example where they it is unique).Another problem for actual computations is that even if the NESS would be known, one would still need to map the states on the left and right onto a GGE, which requires the computation of expectation values of charge densities [40].Unless the model is non-interacting, this is an incredibly cumbersome task.Note that the problems discussed here are not problems connected with specific impurities, but rather problems connected to the interacting nature of the integrable model.
In this paper we would like to put forward the understanding of the peculiarities of the boundary conditions for impurities in interacting models.We already discussed that the effects mentioned above are not present in free models, which were studied extensively in the non-equilibrium context [27,28,[42][43][44].We also do not expect them to appear the perturbative treatment applicable to weak impurities [28,45,46], since effects like non-uniqueness are often non-perturbative in nature (and indeed we will observe this at an explicit example).They should appear in integrable impurities, but they are quite fine-tuned and analytical computations on them are involved (also note the recent work on how to incorporate them into GHD [41]).Instead, for our purposes, ‡ Note that if one would apply ordinary scattering theory to the impurity, this problem would not occur.There one can define the S-matrix which maps incoming states to outgoing states.The difference is that ordinary scattering theory corresponds to a infinite time limit t ≫ N , s.t.eventually all particles are non-interacting.In GHD however, we look at times t ∼ N , implying that particles are still interacting.
we are seeking a broad class of impurities, which a) exists for any integrable model, b) includes strong impurities, c) can be treated analytically and numerically and d) can easily be connected to the GHD language.In the second part of this paper we introduce a class of impurities which satisfies all of these requirements: mesoscopic impurities.
Mesoscopic impurities are impurities whose spatial scale is mesoscopic, i.e. they are very large (and slowly varying) impurities, but still much smaller than the system size.We choose them to be linear combinations of charge densities (note that this is a major difference from the impurity studied in [28], which is also mesoscopic in size).Due to their large spatial size one can describe such impurities using the framework of generalized hydrodynamics itself.This drastic simplification will allow us to study them both analytically and numerically in efficient manners.They are a universal description for large scale impurities built from charge-densities and can be used to model non-perturbative impurities in interacting integrable models.Furthermore, the nonequilibrium stationary states are directly available in the language of GHD.Beside this we also expect that it should be possible to implement them in actual experiments, for instance in cold atom setups.
Even more, following basic ideas from generalized hydrodynamics, one can view the regime of mesoscopic impurities also as a zero'th order term of an expansion in 1/L imp , where L imp is the spatial size of the impurity.In this sense our results can also be a starting point to systematically gain insights into smaller impurities by taking diffusion (or further higher order corrections) to GHD into account [47][48][49].
To keep the discussion simple we restrict ourselves to models in which φ(p, q) = φ(p − q) is a function of the difference only.This includes important classical and quantum mechanical systems, like hard rods and the Lieb-Liniger model.In these models we find an intuitive description of scattering at mesoscopic impurities: One can view the interacting system as non-interacting particles evolving in an effective Hamiltonian.The effective Hamiltonian is the single particle bare Hamiltonian plus a correction coming from the interaction of the particles.We show that this effective Hamiltonian satisfies a functional fixed point equation, which provides a convenient starting point for further analysis and explicit solutions.We use our analysis to establish the following general properties of scattering at mesoscopic impurities: • Scattering at mesoscopic impurities is invariant under local spatial rescalings of the impurity.
• Scattering identically vanishes if the strength of the impurity is below a certain cutoff.
• The trajectory of particles at the impurity are deterministic.This means that all particles with a certain momentum p are eighter transmitted or reflected.
• The solution to the scattering problem (i.e. the boundary condition) is not necessarily unique.
We further demonstrate how the fixed point equation can be used to obtain analytical solutions to the scattering problem at the example of hard rods scattering at a potential barrier, particles only approaching from one side.In the case where the length of the rods d > 0 is positive the solution to impurity problem is unique.However, this is not always the case as we demonstrate for hard rods with negative length d < 0 (negative length corresponds to a time-delay during scattering).We give an explicit example with two possible stable solutions.
We close the paper with a demonstration of our initial assumption of this paper, namely that replacing the impurity by a boundary condition of the GHD equation indeed gives the correct large scale evolution: For a simple impurity we simulate hard rods scattering at an impurity starting from some large scale initial state at time t = 0.At a later time t = T we compare the result to the result obtained from simulating the GHD with impurity boundary conditions and show that they coincide.For the simulation of the GHD equation and in particular the boundary condition we use an efficient algorithm outlined in Appendix A.
The paper is structured as follows: In Section 2 we discuss the relation between the impurity and the boundary condition in general.In Section 3 we then introduce mesoscopic impurities and derive the effective Hamiltonian and its fixed point equation in Section 4. In Section 5 we study the explicit example of hard rods scattering at an potential.In Appendix A we also describe an efficient numerical algorithm to solve the scattering at mesoscopic impurities.

General considerations
Given a system which consists of an (translation invariant) integrable model and an impurity we would like to understand the behavior of the system on large scales.First, let us specify what we mean by an impurity: An impurity is a (possibly strong) localized perturbation of the integrable model whose spatial size L imp is much smaller than the macroscopic scale L ≫ L imp .This is contrary to the case of an external potential, which is a perturbation whose spatial size L ext is comparable to the macroscopic scale L ∼ L ext .In this paper we will restrict ourselves to the case where the system far on the left and far on the right is equivalent.However, the same discussion and the tools can be extended to the case where the system on both sides is different.A simple example of that is the case of a potential V (x) which decays to different values V (∞) ̸ = V (−∞) as x → ±∞.A more complicated example would be systems with different interactions on both sides of the impurity.In such a situation the impurity acts as the connection between two different systems.Also note that our analysis also applies to boundaries (seen as a limiting case where one system is non-existent).
In this paper we want to study the effect of an impurity in a non-equilibrium setting.The problem is as follows: Given an initial state of the system at time t = 0, predict the state of the system at some later time.In interacting integrable models already the evolution without any impurity is complicated and in most cases cannot be computed explicitly.However, in many interesting cases and experiments (for instance in cold atom systems [15][16][17]) the system size L is much larger than the microscopic lengthscale of the particle, the system is observed on long timescales compared to the microscopic timescale and also the system contains many particles.In this case one can describe the out-of-equilibrium dynamics on large scales integrable models by generalized hydrodynamics [12,13].The (Euler scale) generalized hydrodynamics equation can be stated in two ways: In a conservation form where it describes the quasi-particle density ρ(t, x, p) or in a transport form where it describes the occupation function n(t, x, p) [12]: Both equations are linked via: where E(p) is the bare energy of a quasi-particle and the dressing operation is defined as the solution to: Here we defined the operator T, whose kernel T (p, q) = 1 2π φ(p, q) is given by the scattering phase shift, which depends on the model.One can interpret 1 2π φ(p, q) as the effective displacement which occurs when two particles with momenta p and q interact with each other (we will assume that φ(p, q) is symmetric).When φ(p, q) is positive the interaction corresponds to a time-delay or a backwards displacement of the trajectory of the particle.When φ(p, q) is negative the interaction speeds up the particle and thus corresponds to a forward displacement of the trajectories.Note that the above equations are written for an integrable model with a single particle species.It is also possible to formulate them for multiple particle species [12,13].
So far we discussed the system without an impurity.By definition the impurity is way smaller than the macroscopic scale and thus from the GHD viewpoint the impurity will shrink to a point.Throughout this paper we will place the impurity at x = 0.For x ̸ = 0 the system is still described by the GHD equation, but at x = 0 the impurity will lead to some boundary condition, which relates n L (p) = ρ(0 − , p) and n R (p) = ρ(0 + , p).We can now write the GHD equation (e.g. in transport form) including the impurity as: where M denotes the set of all physically allowed relations between n L and n R .Note that we choose to express the boundary conditions in terms of the occupation function n and not ρ (see discussion at the end of this section).
How can one characterize the set M? For that we need to study the scattering problem at the impurity.The idea is that since the impurity is small L imp ≪ L, the timescales for scattering at the impurity are also way smaller than the macroscopic timescales.Thus, at the scale of the impurity we study a long time limit where we can send t imp → ∞.Also since n L and n R change on the macroscopic time-scale (i.e.slowly), we can assume them as constant during scattering.This is the usual assumption on scattering at impurities in general and is, for instance, at the base of scattering theory in quantum mechanics.In the limit t imp → ∞ the solution converges to a stationary solution of the equation of motion.A prominent example of this is the Lippmann-Schwinger equation in quantum mechanics where the scattering state is an eigenstate of the Hamiltonian (and thus stationary in time) [50].We will also see this feature when we study mesoscopic impurities.
The scattering problem of the impurity therefore consists of finding a solution to the stationary equations of motion of the impurity system, so called non-equilibrium stationary states (NESS).This is well known, also in the context of GHD, and these NESS have been constructed and studied in special cases [28,40,42].
Given such a NESS, one can construct the boundary conditions in the following way: Evaluating a NESS far away from the impurity x imp → ±∞ it will become a solution to the unperturbed system.From the GHD viewpoint this state should correspond to a GGE, which is characterised by its expectation values of the charge densities.Therefore, one has to compute the expectation values of charge densities far away from the impurity.Given those expectation values one can in principle construct a GGE, which corresponds to a quasi-particle density ρ(p) or alternatively to an occupation function n(p) in GHD.The n(p) obtained by this procedure are the boundary conditions n L and n R .To summarize we look for a solution to the stationary impurity system, a NESS whose asymptotics for x imp → ±∞ are described by n R/L .The set of all those possible states then determines M.
Usually, one does not describe the boundary condition just as a set M, but one views the outgoing state as a function of the incoming state n out [n in ].For instance, in quantum mechanics one defines the S matrix which maps an incoming state to the corresponding outgoing state.This is a natural and physically intuitive description of the scattering at the impurity.
Unfortunately, viewing the outgoing state as a function of the incoming state is not possible in interacting models.The fundamental reason is that the effective velocity of particles gets affected by the presence of all other particles, incoming as well as outgoing.Therefore, given for instance the solution n L we can determine which particles are incoming and which are outgoing on the left side.But if we do not know the solution we do not know which particles are incoming a priori.In addition to that, even in cases where it is clear which particles are incoming, it is not clear whether the solution to the scattering problem is unique.In fact, we will demonstrate the nonuniqueness at an explicit example, see Section 5.3.To avoid potential confusion, this non-uniqueness refers to the non-uniqueness of boundary conditions.If one applies ordinary scattering theory to any impurity, it is always possible to define a S-matrix, which maps incoming states to outgoing states in a unique way.However, ordinary scattering theory is based on a long time limit and describes the incoming and outgoing states in terms of asymptotic states.These asymptotic states form only at very long times |t| ≫ N where the particles have separated so much that they do not interact anymore.In GHD on the other hand, we consider much shorter times t ∼ N : The states to the right and to the left of the impurity are finite density states, where particles are still interacting.In fact, this also shows that the S matrix is not the correct object to describe scattering on hydrodynamic timescales of interacting integrable models.
Due to the peculiarities outlined in the last paragraph, in general, we need to describe the boundary conditions by a set M. However, there are some restrictions on that set based on basic symmetries and conservation laws.For instance, if the impurity is PT symmetric then if (n L (p), n R (p)) ∈ M , also (n R (−p), n L (−p)) ∈ M .Unless the impurity is integrable as well, the impurity will break some, if not all, conservation laws.For a typical impurity we can expect that it will conserve the total particle number and the total energy.In this case we have the restrictions dp However, in principle also particle number or energy conservation violating impurities are possible.
We finish this general discussion with another subtlety which only appears in interacting models: In principle one can express the set M either in terms of ρ or in terms of n.Both are possible, however, we find that n is better suited for the following reason: Consider the unperturbed model and take as initial state two bumps: a bump on the left side with positive velocity and a bump on the right side with negative velocity.When the system evolves the two bumps will collide.During the collision n(p) will be transported along some complicated trajectory of the quasi-particle, but the value n(p) will not change.The quasi-particle density ρ = 1 2π 1 dr n, however, will change along the trajectory, since it is multiplied by 1 dr , which depends on the presence of other particles as well.The same effect appears during the scattering at the impurity: The reflected particles will alter the shape of the incoming ρ, but not of the incoming n.Therefore n is better suited to describe an incoming state.

Mesoscopic impurites
In general, microscopic impurities can become very complicated and break many conservation laws, which is why we cannot hope to find a general solution to their scattering problem.Usually one can only study very specific integrable impurities or try to study the problem perturbatively in the impurity strength.Even then, in order to connect with GHD, one has to compute the expectation values of all charge densities, which is incredibly cumbersome in interacting models [40].Therefore, most computations so far restricted to free models.As we discussed, in free models one does not encounter the delicate problems for interacting models outlined in the last section.In order to understand them it would be instructive to have a class of impurities, which can be treated analyt-ically for interacting models and in particular whose solutions can be easily connected to the GHD language.In the remainder of the paper we will introduce such a class of impurities, which we call mesoscopic impurities.This class of impurities exists for any integrable model, therefore allowing to model and study impurity effects quite generally.
The idea is to study impurities which are large in size (compared to microscopic), but smaller than the macroscopic length scale of the system, i.e.L imp ∼ L γ , where 1 2 < γ < 1 (we require γ > 1  2 in order to avoid diffusive effects).To be precise we consider impurities which are combinations of charge densities V = dx V n x L imp qn (x), where qn (x) are the charge densities of the integrable model.
We dub these impurities mesoscopic impurities, since their size is mesoscopic.On the length scale L imp we can still divide space in small fluid cells in which we can assume that local relaxation to a GGE has already happened.Since this assumption is the basic assumption of GHD, this in fact means we can fully describe the impurity using the GHD framework alone.This is a great simplification.Of course the results only hold for very large impurities.But still they provide a first insight into general features of impurities in GHD and also can be seen as a first approximation to an actual impurity.In order to study smaller impurities the next step would be to also take diffusive corrections into account, which would give access to effects of impurities on scales L imp ∼ L γ for . By adding more and more orders of the gradient expansion one can thus access smaller and smaller impurities.However, it is not clear whether the gradient expansion will give rise to a convergent series, so potentially it will fail to describe microscopic impurities.
In order to find an equation describing the mesoscopic impurity consider the general GHD equation for a space dependent bare energy function E(x, p) [23]: which is written in conservation form, but using the occupation function n.In our case is the size of the impurity as seen from the macroscopic system.We now want to study the above equation for small ℓ → 0. For that let us change the position variable to x → xℓ.This yields: and by taking the limit ℓ → 0 we find the stationary GHD equation: Note that this is a differential equation which requires boundary conditions at x → ±∞.These are given by the incoming/outgoing data n L (p) = n(−∞, p) and n R (p) = n(∞, p).This data has to be given externally.Again we have the problem that in order to determine which particles are incoming we need to know the complete state n(±∞, p).This means we cannot specify the boundary conditions directly, but rather have to assume the form of n(±∞, p) at infinity, then solve the stationary GHD equation using the incoming data and then check whether the outgoing data coincides with the incoming data.
We conclude that the impurity problem for a mesoscopic impurity reduces to the scattering problem of GHD at the impurity potential.Let us summarize some basic properties of solutions to the stationary GHD equation (8).

Particle number and energy conservation
We can establish particle number conservation by integrating (8) over p: which means that the total current dp v eff (x, p)ρ(x, p) is independent of x and thus the outgoing current is equal to the incoming current.
Similarly we can establish energy conservation by showing that the energy current vanishes: where the last line vanishes due to the well-known symmetry property of the dressing dp f ng dr = dp f dr ng (see for instance [14] equation ( 117)).

Rescaling invariance
There is also another immediate property following from the way we derived the stationary GHD equation: Consider a global rescaling of space V (x, p) → V (λx, p).
The solution to the stationary GHD equation in the rescaled potential is simply given by n(λx, p).This can either be checked explicitly or can simply be observed by noting that if ℓ in derivation ( 7) was a mesoscopic scale then ℓ/λ is a mesoscopic scale as well and thus the resulting GHD equation and its solution must be the same, only rescaled.This shows that the stationary GHD equation is invariant under global rescaling of space.Interestingly the same is true for local rescaling as well: Observation 1 Consider a smooth (or at least differentiable) map y : R → R s.t.y(x → ±∞) = ±∞.Then the stationary GHD equation with potential V (y(x), p) is solved by n(y(x), p).
This can be seen rather easily: We have Since the dressing operation does not depend on x we simply find their dressed expressions evaluated at y(x): We also find that 1 dr = 1 dr (y(x), p).Now inserting everything into the stationary GHD equation we immediately find: Therefore n(y(x), p) is a solution to the locally rescaled GHD equation.Note that y(x) is not required to be monotonically increasing.It is allowed to go back in space.
This rescaling property can be used to simplify problems.For instance, consider an impurity potential which depends only on space V (x) and is non-negative and has a single maximum of V = max x V (x) which is obtained at x 0 .Then we can start from the potential: Now define y(x) as: Then V tri (y(x)) = V (x) and thus we can obtain the solutions to all these potentials simply by solving one potential V tri .If V (x) only has no other local maxima than the global one this also produces the expected result.However, if V (x) has some local minima, y(x) is not monotonically increasing: We still find a solution to the stationary GHD equation, but the minima are filled with some quasi-particles, which are trapped in the minima (see Figure 1).We usually assume the impurity to be empty before scattering, but in principle one could also consider a filled impurity.

Impact of the sign of T (p, q)
In general, T (p, q) describes the effective displacement on the trajectories of two quasiparticles after they scattered.If it is negative it describes a sudden forward jump of the quasi-particles during scattering (like in the hard rods case T (p, q) = − d 2π ).Positive T (p, q) (like in the Lieb-Liniger model) physically corresponds to a time-delay during scattering, but can be thought off as a backward displacement of the particles (like in the flea gas algorithm [51]).
Imagine a particle moving 'up-hill' at a potential V (x) and consider the case where it scatters with a reflected particle going 'down-hill'.If T (p, q) is negative then the two particles will be displaced forward, meaning that the incoming particle gains potential energy ∆E = V ′ (x)2π|T (p, q)|, while the reflected particle looses this energy.Since this effectively lowers the height of the potential barrier for the incoming particle, we conclude that negative T (p, q) tends to help incoming particles to pass a potential barrier.
For T (p, q) > 0 the situation is reversed.Now both particles are displaced backwards, meaning the incoming particle looses potential energy, while the reflected one gains it.Therefore the potential barrier is effectively higher.
Note that this also has an impact on the uniqueness of solutions.For T (p, q) < 0 reflected particles will help other particles to pass the barrier, which leads to less reflected particles.Therefore, there is a competition between the number of reflected particles and their effect on the incoming particles, which should intuitively lead to a stable equilibrium.For T (p, q) > 0, however, both effects act in the same direction: If there are many reflected particles even more will get reflected.Alternatively, if only few particles are reflected more particles will be transmitted.This can lead to two solutions separated by an unstable equilibrium in between and thus the solution might not be unique.
We will see this phenomenon in the discussion of the hard rods (see Section 5).

Stationary GHD equation for T (p − q) and effective Hamiltonian
We already established some general properties of the solutions in the last section.Now we will specifically study the case where the scattering phase shift depends only on the difference of momenta T (p, q) = T (p − q), where great simplifications happen since the normal modes are known.

Stationary GHD equation in normal modes
The stationary GHD equation ( 8) is hard to solve because it is written in a continuity equation form ∂ x (An) = ∂ p (Bn).It would be more convenient to write the equation in transport form: where m(x, p) are called the normal modes.In general, it is not clear how to derive the normal modes, or whether they exist.However, in the case of a scattering phase shift which only depends on the differences of the momenta T (p, q) = T (p − q) the normal modes exists and are simply given by n(x, p) [23], similar to the GHD equation without external potential.This gives the stationary GHD equation in transport form: For other models where T (p, q) ̸ = T (p − q) this does not work.However, one could try to reparametrize p → f (p), which might lead to a reparameterized scattering phase shift depending only on the difference.Alternatively, if one is still able to find some other normal modes, the derivations done in this section will still be applicable with minor adaptions.
Let us briefly recap why the occupation function n are the normal modes of the GHD equation in the case T (p, q) = T (p − q).The stationary GHD equation in conservation form is given by: Let us write out the derivatives of the dressing explicitly: Here we used T (p, q) = T (p − q), which implies ∂ p , T = 0 to swap ∂ p and T. Using the definition of the dressing we can rewrite this as: In particular, we find that: From this we can easily see that if (21) holds (27) and thus also (23) will vanish.The advantage of ( 21) over ( 8) is that solutions to (21) can be described in terms of characteristics.Consider a particle moving along a solution of the characteristic ODE: Then it is easy to see that n(x(t), p(t)) is constant in time, i.e. n is constant along characteristics.This property is particularly useful for numerical simulations (see Appendix A.1.2).

Effective Hamiltonian
Writing the stationary GHD equation in transport form is already very useful.However, it becomes even more interesting if one notices an additional fact.As a byproduct of the derivation of the stationary GHD equation in normal modes we also showed that 27).This implies that there exists a function H(x, p) with the property: where ∇ = ∂ x ∂ p .We call this function H(x, p) the effective Hamiltonian for the following reason: The stationary GHD equation can be written as which is precisely the equation for a distribution of non-interacting particles evolving according to the Hamiltonian H.If there was another time derivative this equation would be the Lioville equation for that Hamiltonian.Since the time derivative is not there, this problem is now the scattering problem for non-interacting particles with Hamiltonian H(x, p).This idea is typical in integrable models and in GHD.One tries to rewrite the problem in terms of another effective problem for non-interacting particles where all quantities get dressed by the other particles (consider for instance the effective velocity v eff (x) which is not the bare velocity of the particles, but gets altered due to the presence of other particles).In that respect we can view the Hamiltonian H(x, p) as the effective (or dressed) Hamiltonian under which the particles evolve.The effective Hamiltonian H(x, p) will be given by the bare Hamiltonian E(x, p) plus some other terms which originate from the interactions between particles.
We want to note that H(x, p) is technically equal to E Dr defined in [52, supplementary material, equation ( 5)].However, there the authors only look at the p dependence of E Dr (p), which works in general for any T (p, q).Since we are also interested in the x dependence, the situation is more complicated and only gives (30) in the case T (p − q).
Let us now derive a formula for the effective Hamiltonian.We will start with the definition of the dressing operation: The curl of n∇H is given by: which vanishes due to the stationary GHD equation.Therefore, we can find a function N (x, p) with the property: and thus: Let us take this equation in x and integrate it from −∞ to x: By looking at x → −∞ we can fix H(−∞, p): where we are free to choose the constant C. We can interpret the effective Hamiltonian (35) as follows: The first part H(−∞, p) describes the evolution according to the GHD equation without impurity.The second term V (x, p) describes the bare contribution to the impurity.The third term: describes the contribution to the impurity coming from the interactions between particles.Note that while V (x, p) will vanish for large x, U (x, p) will vanish by definition for x → −∞, but in general does not need to vanish for x → ∞.In fact, this is important as U (x → ∞, p) = 0 implies that the state on the right side of the impurity is unaffected by the impurity, i.e. it signals full transmission.As we will establish in the next section, N (x, p) will not vanish for large p → ±∞ but rather approaches a constant.Therefore, if T (p) is an integrable function, i.e.
dp |T (p)| < ∞ this allows us to remove the second term in U by redefining H(−∞, p) and U (x, p) → Ũ (x, p) = TN (x, p).This gives the following simplified expression for the Hamiltonian: which agrees with (35), up to the arbitrary constant.However, if T (p) is not integrable (for instance in the hard rods case T (p) = −d) then this definition of U (x, p) would be ill-defined ( Ũ (x, p) = ∞) and therefore we will stick to (35) for the following general discussion.
Let us give two remarks about the general form of U (x, p): Second, note that U (x, p) is not a general function, but has to be in the image of T.An important example where this is useful is the hard rods case where U (x) = TN (x, p) = −d dq N (x, q) is a function of x only.The fact that U (x) depends only on x simplifies the problem so much that one can actually solve the impurity problem (see Section 5).
We will finish this section with a simple, yet physically very important observation: Observation 2 At a mesoscopic impurity all particles with momentum p are either transmitted or reflected with probability one.
This follows directly from the existence of the effective Hamiltonian H(x, p).The particles have to follow deterministic trajectories.Therefore all incoming particles with a specific momentum will follow the same trajectory and eventually leave the impurity at the same point.
Note that this is substantially different from scattering at microscopic impurities, where particles are reflected or transmitted with a certain probability (for instance recall scattering at an rectangular potential barrier from your undergraduate quantum mechanics course).We conclude that this is a limitation of the mesoscopic impurity model and that non-deterministic scattering is an effect which appears only if the impurity is small enough.On the other hand, the amount of non-deterministic scattering can be used to quantify how well an impurity can be described by a mesoscopic impurity.

The scattering problem in a general Hamiltonian
In the previous section we introduced the effective Hamiltonian: with:  So far this result not useful since we do not know N (x, p).For that we need to solve the scattering problem in the effective Hamiltonian H(x, p): We know that both n(x, p) and N (x, p) satisfy the following equation As we described earlier this equation implies N (x, p) is constant along characteristics: However, we also know that a basic property of Hamiltonian systems is that the Hamiltonian is conserved along trajectories.This means a trajectory follows a level set of the Hamiltonian H(x(t), p(t)) = const.Now we know that N and H are both constant along a trajectory which implies we can locally write N (x, p) = N (H(x, p)).
Locally here means that this is true for a domain Ω ⊂ R 2 which does not contain any critical points ∇H = 0 and does not contain more than one connected part of each level set.The reason why this is important is that globally the level sets of H(x, p) can have multiple disconnected components.If one starts on one of these components and follows the level set one never reaches the other disconnected components and thus there is no reason why N (x, p) should have the same value on these other components.Now let us assume we know all the incoming data n L (p) = n(−∞, p) and n R (p) = n(∞, p) on the left and the right side.As we discussed beforehand this typically requires the full knowledge about the state including the outgoing data as well.
For simplicity we will also restrict to the case where p dr (±∞, p) is a monotonically increasing function.This has the advantage that we can be sure that there is only one p which corresponds to zero velocity and on the left side all particles with higher momenta are incoming and all particles with lower momenta are outgoing (and similar on the right side).
If p dr (±∞, p) is not monotonically increasing then we might have several momenta corresponding to zero velocity which complicates the following reasoning.Nevertheless, by introducing more sets Σ, it will be possible to also include that case.
Starting from the incoming data on both sides we can first compute the Hamiltonians H L/R (p) on both sides (which we know up to a constant) and then we can define the two functions: where p L/R are the momenta corresponding to zero velocity on both sides and where we choose the incoming branch, i.e.
We will now describe how to construct N (x, p) from the incoming data in a general Hamiltonian.This construction is also illustrated in Figure 2 for the example Hamiltonian depicted in Figure 3.
Given a general Hamiltonian we can group its level sets into three type of sets: First, we define Σ L/R which are all the level sets that originate from incoming data on either the right or the left side.We think of this as the regions which are allowed for particles to enter.For instance, Σ L consist of level sets which start from the x = −∞ on the left side with positive velocity and then either extend all the way to either the right side (i.e. the particle is transmitted) or it bends back and extends to the left side again (i.e. the particle is reflected).Then we have a collection of sets Λ L/R i which are 'islands' in either Σ L/R where the energy H(x, p) is too high (or to low) for particles to enter.Λ L/R i consist of closed loops of level sets.Note that we associate to each of these Λ ) which corresponds to the energy of particles 'flowing' around Λ L/R i .The last type of sets are non-allowed regions similar to the Λ L/R i which lie between Σ L and Σ R .Again they consist of closed loops of level sets.However, unlike the Λ L/R i these islands will not give a contribution as we will define N (x, p) = 0 on them.
The two sets Σ L/R are separated by a curve γ(s) ∈ R 2 .This curve might not necessarily be unique, it is only required that Σ L lies on one side of the curve and Σ R on the other.We will parametrize this curve in such a way that γ x (s → ±∞) = ±∞.In many cases the curve γ(s) will be a function of x only in which case we write γ = (x, p 0 (x)).Note that p 0 (−∞) = γ p (−∞) has the interpretation of the largest momentum incoming from the left side which is reflected by the impurity and the smallest momentum coming from the right side which is transmitted by the impurity (and similarly for p 0 (∞) = γ p (∞) on the right side).We will fix the integration constants in such a ways that both N (γ(s)) = 0 and the Hamiltonian H(γ(s)) = 0 on this curve.
Once these sets have been characterized one can explicitly define N (x, p): The first two lines just copy the incoming data into the regions Σ L/R .For instance, on Σ L we know that N (x, p) = Ñ L (H(x, p) + C) + D, where C and D are constants we are allowed to choose.We set The constant is determined by its value on the boundary ∂Λ L/R i , which corresponds to a level set of the Hamiltonian with energy h L/R i .This complicated procedure allows us to compute N (x, p) given H(x, p).In turn we can now use this to compute U (x, p): Let us define Σ L/R (x) = p : (x, p) ∈ Σ L/R and Λ Then we can rewrite definition (37) as: Note that if the T (p) is not integrable, then the individual integrals might not converge since N (p → ∞) → const.In this case one has to combine the integrals into one integral over the sum of all three terms (which will then be finite since all constants at infinity cancel).We can also insert expression (39) for H(x, p): This is a functional fixed point problem for U (x, p) with fixed point functional G[U ].Note that this fixed point problem is quite complicated, in particular since the shape of the domains Σ L/R and Λ L/R i depend on U (x, p).In case we manage to find a solution to the fixed point equation we can write the solution to the stationary GHD equation (21) as follows: Remark: The fixed point problem also has contributions from the regions Λ L/R i , even though there n(x, p) = 0.This is because the defining relation ∇N = n∇H = 0 does not not imply that N (x, p) = 0, but only that N (x, p) = const in Λ L/R i .This is mathematically similar to the 'Aharonov-Bohm effect' in quantum mechanics [53].

Discussion of the fixed point equation
The fixed point equation ( 47) is quite complicated and can in general only be solved numerically.However, even without knowing the solution explicitly one can still establish some basic properties.

Existence of solution
It is physically clear that there should be some solution to the scattering problem given the incoming data.However, as we discussed earlier it is not clear how to specify the incoming data in interacting models and therefore it is also not clear whether a solution will always exits.In the following we will give a mathematical argument why we expect that the fixed point equation should always have solution for given incoming data N L/R .Let us rewrite this as a functional equation: First, observe that this functional changes continuously under smooth deformations of U (x, p).This is obvious as long as changing U does not create or destroy domains Λ This is particularly obvious if T (p) is integrable as the interaction induced term is then bounded.But also when T (p) is not integrable the interaction induced term will still approach a finite value as λ → ±∞.No matter the sign of λ as |λ| → ∞ we have that |U (x, p)+λ∆U (x, p)| will be very large around the support of ∆U and thus we have an island Λ L/R i formed around the support of ∆U .The value of N (x, p) on this island is constant and will approach N (∂supp(∆U )) which is independent of λ for large |λ|.Thus, as λ → ±∞, N (x, p) approaches a finite function and therefore the interaction induced term is bounded in λ, which implies (50).
We have now shown that for any ∆U we have dx dp F [U +λ∆U ](x, p)∆U (x, p) → ±∞.Let us take the formal limit where ∆U approaches a delta function which gives: Of course the Hamiltonian H(x, p) + λδ(• − x)δ(• − p) does not make proper sense which is why we choose to regularize them by considering bump functions.We now view U (x, p) as a vector ⃗ U and δ(• − x)δ(• − p) as a basis of this functional vector space.A finite-dimensional version of statement ( 51) is: for all k as λ → ∞.If (52) holds in one dimension and F is continuous then the intermediate value theorem shows that there must be a zero of F (U ) on the real line.In the finite-dimensional case there is a generalization of the intermediate value theorem, called Poincaré-Miranda theorem, which (up to some mathematical technicalities) establishes that (52) implies the existence of a zero of F [U ].This indicates that (51) should imply the existence of a solution F [U ] = 0 as well (however, this is not rigorous as results in finite dimension possibly do not carry over to the infinite-dimensional setting).
4.4.2.Discussion of uniqueness While the existence of solutions is physically expected, there is no reason why solutions should be unique.Non-uniqueness of solutions means that the solution will depend on the scattering history.Therefore studying the uniqueness of solutions gives physical insights.For that it is constructive to compute the Jacobian of F [U ], i.e. how F is perturbed by small perturbations of U (x, p) → U (x, p) + δU (x, p): Let us again consider a δU which is a bump function.The important observation here is that as long as δU does not affect any critical points of H(x, p) (i.e. it does not move them around or create/destroy them) the boundaries of the sets Σ L/R and Λ L/R i might change in a continuous fashion, but the topology of the sets does not change.In particular the energy levels h L/R i associated to Λ L/R i remain unchanged.For such a perturbation we therefore find: = δU (x, p) dq T (p − q)n r (H(−∞, q) + V (x, q) + U (x, q))δU (x, q) where we used (48) to identify the solution n(x, p).Again we let δU (x, p) → δ(x − x 0 )δU (p) approach a delta function in space which gives: where n(x 0 , p) is the density at x 0 .From this we can see that perturbations at x 0 only affect the result at x 0 but not at any other y ̸ = x 0 .This means that the fixed point equation decouples for different x (as long as there are no critical points at x).
Uniqueness of the fixed point equation is now linked to whether the Jacobian 1 − Tn(x 0 ) is invertible.We can identify this operator as the inverse dressing operator (4), which we know should be invertible if n represents a physical state.Therefore, as long as we stay inside the set of physical states if we find a solution at x 0 it will be a unique one.There can be multiple solutions but they either need to be separated by nonphysical states or by changes of the critical points.In particular quantum mechanical models (like the Lieb-Liniger model) where T (p) is integrable dp T (p) = 1 and in addition the occupation function is bounded n(x, p) < 1 (i.e.every quantum number can be occupied at most once) the Jacobian will always be invertible (in that case the operator norm Tn ∞ < 1 and it is a standard result that 1 − Tn is invertible).The above arguments only partially answer the question of uniqueness since there we excluded the case when δU affects critical points of the Hamiltonian.It only tells us that the Jacobian is invertible on that subspace.We also need to discuss what happens in case δU affects a critical point.This is hard to analyze in general since changing U at a critical point will affect the solution globally.For instance, when δU acts on a critical point on the boundary of a Λ L/R i it will change the energy h , which (since the boundary is a level set of H(x, p) with that energy h L/R i ) will move the whole boundary ∂Λ L/R i .However, the set of critical points of the Hamiltonian will usually be finite and determine the topology of the problem.If we know the positions of the critical points, the values p 0 (±∞) = γ p (±∞) and the energies h L/R i associated to the islands we can construct the sets Σ L/R and Λ L/R i .In a model where we can ensure that the dressing is always well defined, this finite amount of information is enough to construct a full unique solution U (x, p) from it.In practice this allows us to simplify the functional fixed point problem into a finite-dimensional fixed point problem, which is much easier to analyze.We will use this strategy later in the case of the hard rods model to show that for positive rod length the solution is unique, while for negative rod length we give a specific example where the solution is not unique, see Section 5.3.
Remark: The discussion of uniqueness here considers only the uniqueness of the fixed point problem.The construction of the fixed point problem requires the knowledge about which particles are incoming.As we already discussed one can only determine this if one knows the outgoing particles as well.Thus, even if the fixed point problem always has a unique solution there could be multiple configurations of incoming and outgoing particles.

Solution for weak potentials
Although the fixed point equation cannot be solved explicitly for a general model it is a non-perturbative method.A different, perhaps more standard, way of approaching the impurity problem would be to treat the impurity via perturbation theory in impurity strength λ.In general, the problem with perturbation theory is that it is often not clear whether the perturbative expansion indeed converges to the actual solution, or whether it misses non-perturbative effects.In this section we establish that for small impurity strength λ the scattering will vanish identically.This implies that a perturbative expansion in λ would yield 0 to all orders.We conclude that scattering at mesoscopic impurities is fully non-perturbative in nature.Consider any (bounded) potential V (x, p) and scale it λV (x, p) with λ → 0. Then there exists a λ 0 > 0 s.t. for all |λ| < λ 0 there exists a solution to (47) where the states on the left and right side coincide.This means that all particles are unaffected by the impurity and no particles are reflected.
Note that this is clear for an impurity for free particles: For instance, consider the impurity V (x, p) = V e − x 2 2 .Then all particles with momentum p > √ 2λ V will be transmitted, while all particles with momentum p < √ 2λ V will be reflected.If all particles have a momentum bounded from below |p| > p 0 , then for λ < λ 0 = p 2 0 2 V all particles are transmitted.
In interacting models the idea is similar: The interaction changes the Hamiltonian, but if λ is small enough this change is negligible and thus the free particle result applies as well.While this establishes that particles are not reflected, it does not explain why the transmitted particles are also unaffected.For instance, the momenta of the outgoing particles could be shifted.This can only be checked by explicitly constructing the solution, which we do in Appendix B.

The hard rod case
We would like to show how how the ideas from the last section can be applied to a specific model: the hard rods model.In this model we have φ(p − q) = −d and E 0 (p) = p 2 2 .The physical hard rods model is given by positive rod length d > 0, but we will also allow d < 0. In negative hard rods the d does not describe the rods size, but should rather be interpreted as the time delay for the scattering of two particles, which leads to an effective negative position shift −d.
Before we start let us note a speciality of the hard rods model.One can explicitly evaluate the dressing (4): and the effective velocity: which in particular means that the momentum that corresponds to zero velocity is d ⟨p⟩ := d dp pρ(x, p).Using p dr = p − d ⟨p⟩ = ∂ p H(±∞, p) we can explicitly compute the effective Hamiltonian outside the impurity: where C ± are the integration constants on both sides.Note that d ⟨p⟩ has the same value on both sides since the total current satisfies dp v eff (p)ρ(p) = dp pρ(p) = ⟨p⟩ and the total current is independent of x (see equation ( 9)).
Due to T (p − q) = − d 2π we find that U (x, p) = U (x) of ( 47) is a function of x only.Compared to the general case, where U (x, p) can depend on both x and p this already simplifies the problem.We will now argue that one can simplify the fixed point equation even further and finally arrive at a finite-dimensional fixed point problem: First, let us recall that U (x) satisfies the following fixed point equation: dq N r (H(−∞, q)) . ( As we discussed before (see Section 4.4.2),unless at x there is a critical point of H(x, p), then a change of U (x) will only affect the fixed point equation at x, but not at any other position.Therefore fixing one of those x, where is no critical point of H(x, p) we can rewrite the fixed point equation as: which is a function of one variable only.Furthermore by taking the derivative of F x we find for d > 0: dq n r (H(−∞, q) + V (x, q) + U ) ≥ 1 (62) and therefore F x (U ) is monotonically increasing.Here n L/R (h) = dN L/R (h) dh , which we identify with the solution (48), i.e. n(x, p) = n L/R (h(x, p)) on Σ L/R .This is implies uniqueness of a solution U = U (x).Given the location of the critical points we can therefore compute U (x) at all other positions x.The only missing pieces of information are locations of the critical points, d ⟨p⟩ and the integration constants C ± from (59).Note that this is a finite set of information.Therefore, it should be possible to reduce the fixed point equation to a finite-dimensional problem.This is an abstract statement, but we will make the ideas more explicit at an example.
Also note that the above discussion was for d > 0, for negative d there is no mathematical reason why (62) should be positive.Therefore, there could be potentially more than one solution.

Explicit example: Potential does not depend on p and incoming particles from the left
In the following we would like to explain this procedure on the most simple example: An impurity V (x) which only depends on position x.For simplicity let us also restrict to V (x) ≥ 0 which has one unique maximum V at x = 0 and V ′ (x) > 0 for x < 0 and V ′ (x) < 0 for x < 0. Note that due to the local rescaling invariance all of these impurities will give rise to the same scattering (see Section 3.2).To make the analysis even more simpler let us now look at the situation where there are no incoming particles from the right, i.e. n R (p) = 0.This example can easily be extended to the case where particles also come in from the right, but then the discussion of uniqueness becomes much more complicated.
Since there are only particles coming in from the left it is convenient to redefine the zero value of the Hamiltonian to be at the zero-velocity momentum d ⟨p⟩ at x = −∞.This means that the Hamiltonian is given by: where we combined W (x) = V (x) + U (x).To understand the simplicity of the Hamiltonian (63) is, the reader might find it helpful to define v = p − d ⟨p⟩ and write the Hamiltonian as: This is just the Hamiltonian of a classical particle moving in the effective potential W (x). Scattering at such a potential is very simple.Let us denote by W = max x W (x) the highest value of the effective potential.As we will observe later this maximum is at x = 0 which is the same position as the maximum of V (x).Then all incoming particles with kinetic energy smaller than W will get reflected and all other particles will get transmitted.
This can be expressed in the following way (note that we slightly redefines N (x, p) due to the redefinition of H(x, p)): Here N L and the boundary p 0 (x) are explicitly given by: Note that in the region p < p 0 (x) there are no particles, but still N (x, p) is non-zero since N (x, p) has to be a continuous function.
From the fixed point equation for U (x) we can write down a equation for W (x): The last term comes from the −N (−∞, q) term in (37).Note that since N (x, p) is a non-zero constant for p < p 0 (x) both integrals in (37) are formally infinite, however the divergent part in both integrals cancels.Equation ( 68) is particularly useful since V (x) only appears on the right hand side.
Let us take the derivative of F x (W ): where we defined n ).For d > 0 this derivative is always positive, for d < 0 it might become negative as well.However, note that for a physical state dFx(W ) dW coincides with 1  1 dr (x) , which has to be positive.From equation (69) we can compute the derivative of W (x) w.r.t.x: This equation has a nice interpretation: The slope of the effective potential (i.e. the force) is the slope of the bare potential, but modified due to the interaction with other particles.In particular for d > 0 the slope is decreased (the potential barrier is lowered, thus more particles will be transmitted), while for d < 0 the slope is increased (the potential barrier is higher, thus less particles will be transmitted).This is in line with the intuitive picture we had discussed before in Section 3.3.From equation (70) we can also infer that the highest point of the effective potential W (x) is at the highest point x = 0 of V (x) as well.
Let us evaluate (68) at x = 0: Since F 0 is continuous and F 0 (0) = 0 and F 0 ( W → ∞) → ∞ we know that there will exist at least one solution W > 0 for V > 0. Again let us compute the derivative: This is always positive for d > 0 and thus a unique solution exists.For d < 0 this is not clear: While dp n(x, p) < 2π |d| is bounded for physical reasons (1 dr is has to be positive) there is no bound on the second term n L ( W ) = n L (d ⟨p⟩ + √ 2 W ). It is not hard to construct a situation where there are multiple solutions W for some value V .
Recall, however, that this solution is only a solution to the fixed point equation, which is not the full solution to the scattering problem (see the discussion at the end of Section 2).In fact, even if there are multiple solutions to (71) these solutions will be physically different since they give rise to different d ⟨p⟩.

Determining d ⟨p⟩
So far we have discussed how we can find a solution to the impurity problem given the effective Hamiltonian at x → −∞, in particular on d ⟨p⟩.In fact, we do not know d ⟨p⟩ a priori, but instead we need to compute it in a self consistent fashion: We need to match the assumed d ⟨p⟩ with the d ⟨p⟩ = 1 2π 1 dr d dq qn(x → −∞, q) computed from the solution.
First note that any solution to the fixed point equation (71) can always be upgraded to a full solution by a shift of the momentum variable: In fact, if we define v = p − d ⟨p⟩ and parametrize the initial condition in terms of n v (v) = n(d ⟨p⟩+v) the parameter d ⟨p⟩ completely drops out of the equation.Then if we take any solution n v (x, p) of the new equation, compute its ⟨p⟩ = dv n v (v) and reintroduce the momentum p = d ⟨p⟩ + v we always find a solution n(x, p) stationary GHD equation.Now let us consider a scattering scenario where we know the distribtion of incoming particles as a function p.At this point we find that the problem is ill-defined.In order to specify which particles are incoming we need to know the momentum corresponding to zero velocity d ⟨p⟩.
Fortunately, in case of only incoming particles from the left there is a situation where we can uniquely specify the initial data.For d > 0 the trick is to observe that the more particles are reflected the smaller d ⟨p⟩ becomes.Thus d ⟨p⟩ is largest for full transmission.Therefore, if we consider a state n(p) which is identically zero for all momenta smaller than some cutoff p 0 and this cutoff is larger than p 0 > d ⟨p⟩ = d dppn(p) 1+d dpn(p) then no matter the amount of reflected particles all incoming particles will always have positive velocity.Note that for d < 0 we can always choose p 0 = 0 since d ⟨p⟩ ≤ 0.
In the present case of an impurity which does not depend on x and incoming particles only from the left, one can derive the following compact expression for the zero velocity momentum: Equation ( 73) together with (71) are a two-dimensional closed system of equations: By explicit computation one can check that the determinant of the Jacobian of this system is always positive if d > 0. We conclude that for d > 0 there always exists a unique solution to the scattering problem.

Hysteresis for negative hard rods d < 0
Contrary to hard rods with positive size in the case of negative rod length there can indeed exist multiple solutions to the two-dimensional system of equations.We would like to demonstrate that using an explicit example.We study the scattering of n L (p) = θ(p − 1 2 ) 1 2 e −50(p−1) 2 at the potential V (x) = V e − x 2 2 .We plot the values W of all solutions as function of V in Figure 4 (left) and compare it to a numerical In the depicted range there are three theoretical solutions W .The simulation follows the lower branch for increasing V and the upper branch for decreasing V (the middle branch is unstable).On the right side we picked V = 0.4 and show trajectories of particles at the impurity (left) and the effective potential W (x) as function of x (right).
simulation.In the range V ≈ [0.36, 0.41] there are three different solutions, only the outer two are stable.The numerical simulation was done using the algorithm described in Appendix A.1.2.In addition we adiabatically ramp V up and later down: after each iteration we increase/decrease V (∆ V = 0.001 per step).For increasing V the simulation follows the lower branch, while for decreasing V the simulation follows the upper branch.On the right side of Figure 4 we also show the trajectories of particles and the shape of the effective potential.For small V only few particles are reflected, but as 1 dr (x) = 1 1+d dqn(x,q) > 1 the effective potential is already larger than the bare one by (70).When we increase V over the threshold value V ≈ 4.1 the situation suddenly changes: Now a considerable fraction of the particles are reflected, which increases 1 dr (x) on the left side, which in turn increases W (x). Therefore even fewer particles can penetrate the impurity.This feedback loop only stops when most of the particles are reflected and only a negligible fraction of the particles is transmitted.If we now decrease V this situation is stable: The reflected particles contribute to a large 1 dr (x) and thus to a high potential W (x), which most particles cannot penetrate.Only when we decrease V below the threshold value V ≈ 3.6 too many particles are transmitted and the system jumps back to the original configuration with few reflected particles.This explains the hysteresis observed in the numerical simulation.We would like to finish this paper with a demonstration that the GHD equation including the boundary conditions obtained from the mesoscopic impurity indeed gives the correct prediction for the evolution of the quasi-particle density in the large-scale limit L → ∞.

Simulation of the GHD equation with impurity boundary condition
We study the evolution of hard rods starting from an initial state characterized by (in microscopic coordinates), where L imp = 1 10 L 3 4 is scaled with the macroscopic scale L. The relation between the microscopic and macroscopic coordinates is as follows: x micro = Lx macro and t micro = Lt macro .We choose the height of the impurity to be λ = 0.4, the hard rod length to be d = 0.3 and the simulation time to be T macro = 3, after which most of the particles have scattered.
The GHD simulation with impurity boundary condition are performed using the hybrid algorithm described in detail in Appendix A.2. Space, momentum and time are discretized with ∆x macro = 0.006, ∆p = 0.004 and ∆t macro = 0.001 (for further details refer to Appendix A.3).At these values the GHD simulation has well converged.
For the hard rods simulation for a given L we first distribute the hard rods randomly according to their initial distribution ρ(0, x macro , p) = 1 2π 1 dr (x macro )n(0, x macro , p) and then evolve them using simple molecular dynamics with ∆t micro = 0.01.Since we place the particles initially random, the resulting final state will also be random.Therefore, in order to gain sufficient statistics, this simulation is repeated 100 times for each L and averaged.
The results of both simulations are compared in Figure 5, where we plot a) the resulting ρ(t macro = 3, x macro , p) for both simulations as a density plot and b) the total density of particles in two regions R (reflected particles) and T (transmitted particles.Already in the density plot a) one can see that the quasi-particle density obtained via GHD with impurity boundary condition agrees quite well with the averaged density from the hard rods simulation.
The results can also be compared more quantitatively in b).The data points and their errorbars are the mean and standard deviation of the hard rods Monte Carlo results for several L. We can see that this agrees well with the GHD simulation (solid line) in both regions.
We conclude that the GHD equation with impurity boundary conditions indeed correctly describes the large scale evolution of an integrable model with a mesoscopic impurity.
Further numerical results are given in Appendix A.3: We simulate the hard rods GHD equation starting from the same initial state and the same impurity, but for different hard rod lengths d = −0.3,0, 0.3 and qualitatively discuss the difference of how they scatter at the impurity.

Conclusion
In this paper we studied mesoscopic impurities in the GHD framework.First, we discussed how to include impurities into GHD in general and found that they are described by boundary conditions.The boundary conditions correspond to nonequilibrium stationary states of the microscopic impurity model.A big complication in interacting models compared to non-interacting models is that reflected particles affect the incoming particles (in particular their effective velocities) and therefore one cannot specify the incoming particles without knowing the outgoing particles.This means the solution to the scattering problem is only given by a collection of possible solutions, but it is not directly given as a map from the incoming data to the outgoing data.
Since it is not known how to find non-equilibrium stationary states in general, we decided to study a specific class of impurities: Impurities on a mesoscopic scale (larger than the diffusive scale), which are linear combinations of charge densities.These can be described by the (Euler-scale) stationary GHD equation in an external potential.This provides a broad class of impurities, which are present in all integrable models and can be studied analytically also for strong impurities.
Additional simplification occurs when we restrict to models where the scattering phase shift is only a function of the difference of the momenta.Here the scattering problem can be interpreted as particles moving in a one-particle effective Hamiltonian H(x, p), given by the bare energy E(x, p) plus an interaction dependent correction.Given the incoming data we described how to solve the scattering problem at such a Hamiltonian H(x, p) in general.This allows us to write down a functional fixed point equation for the correction term in H(x, p).If the incoming data is known then a solution to the fixed point equation corresponds to a valid solution to the stationary GHD equation.However, as the outgoing data might affect the incoming data, in general a solution to the fixed point equation will not reduce to the assumed incoming state for x → ±∞.Instead, one would have to adapt the asymptotic data until it is matched by the incoming data of the solution to the fixed point problem.
Still, the fixed point equation provides a starting point to study mesoscopic impurities in general and even to find analytic solutions to the scattering problem.In general, mesoscopic impurities show the following features: They are invariand under local spatial rescaling in x, the scattering at them is deterministic and scattering vanishes for sufficiently weak impurities.Furthermore, it is possible to solve the scattering problem for mesoscopic impurities via an efficient algorithm.
We provided explicit solutions to an impurity in the hard rods model, given by a potential, which is independent of p.Here the difference between a positive (corresponding to a time delay) and a negative (corresponding to an instant jump) scattering phase shift φ(p) becomes apparent: Negative φ(p) tends to increase the transmission, while positive φ(p) decreases it.For the hard rods we also compared molecular simulations and GHD simulations with impurity and found agreement.This shows that the scattering can be indeed captured by a boundary condition.
Mesoscopic impurities are of course only an approximation for large impurities, in particular since we only work on the Euler scale.It would be interesting to go beyond that.In regard of the gradient expansion we can view this as a first order approximation in the impurity size 1/L imp .By adding more terms of the gradient expansion one can compute higher order corrections to our results which would allow the treatment of smaller mesoscopic impurities (for instance on the diffusive scale).However, it is not clear whether this expansion extends all the way to microscopic impurities, i.e. whether this expansion converges or is only an asymptotic expansion.
To study this we expect that it will be interesting to look at the transition between deterministic and non-deterministic scattering: We have established in this paper that scattering at (Euler scale) mesoscopic impurities is always deterministic: particles are always transmitted or reflected with probability one.This is of course not generally true for microscopic impurities.It would be interesting to see whether higher order gradient expansion [47][48][49] introduces non-deterministic scattering, otherwise microscopic impurities will be fundamentally different from mesoscopic ones.
But also on the level of mesoscopic impurities there are open questions left.We discussed that solutions to the scattering problem are not unique in general, but it is not clear whether this appears in all models and under which conditions.In particular, it would be interesting to have a more detailed look at the scattering problem in quantum mechanical models, for instance the Lieb-Liniger model.There the interaction between particles is momentum dependent and also φ(p) is integrable which is quite different from the hard rods model.Furthermore, we restricted to impurities which can be described by external potentials.Another possibility would be to have an impurity where the interaction changes locally.If the interaction changes sufficiently slowly, one can describe the system by GHD [25] and therefore it should be possible to study properties of such impurities using similar ideas as for external potentials.
It would also be interesting to go beyond these simple integrable models and also consider models where φ(p, q) ̸ = φ(p − q) is not a function of the difference of momenta only.In this case, the normal modes of the GHD equation with external potential are not known, so most of the analysis in this paper cannot be applied.It is not even clear whether normal modes will always exist in any model.If there are models where they do not exist, scattering could be quite different.For instance, since one cannot think about the system as particles moving along GHD characteristics anymore, scattering might be non-deterministic.Also the numerical algorithms outlined in Appendix A do not apply in that case and it would be interesting to extend them.
A completely different direction of research would be to implement mesoscopic impurities in actual experiments.We believe that this should be possible among others in cold atom experiments, which are well described by the Lieb-Liniger model.For instance, a potential barrier could be implemented similar to an external potential, which is already used in GHD experiments [3,15].A precise setup we have in mind would be to prepare the atoms in a 1D trap (similar to [15]).Then at time t = 0 the trap is removed and at the same time a sharply peaked potential is turned on which serves as the impurity.After releasing the atoms, the distribution of particles can be observed as function of time and then compared to theoretical computations.A feature that would be particularly interesting for experimental realization is the local rescaling property of mesoscopic impurities.This implies some kind of stability against perturbations: Only few specifications of the impurity determine the scattering behavior, the precise shape of the impurity is not important.the total number of characteristics are important, but also their distribution.In fact, it is important to make sure that there are sufficient characteristics in regions where the incoming densities n(±x 0 , p 0 ) are high.High densities have the biggest impact on the precision of the dressing, which in turn is important to be able to compute accurate characteristics.On the other hand it is also important to have at least some characteristics in low density regions.We found that the following way of distributing produces good results: First, we sample the majority of the p 0 from the probability measures which are proportional to the incoming densities.This naturally places many characteristics in high density regions, but places only few characteristics in the low density regions.In order to ensure sufficient coverage of those regions as well we place the remaining characteristics with a uniform spacing between them, i.e. p 0 = k∆p, where k ∈ N.
In our simulations we choose a small time-step ∆t = 0.001 to compute the characteristics.The dressed (∂ p E) dr and (∂ x E) dr we compute by partitioning the space into small boxes of size ∆x = 10∆t.Therefore, each characteristic will provide ∼ 10 points (x(t), p(t)) in each box.We collect all those points from all characteristics and sort them according to their momenta.At these points (x(t), p(t)) we now know the value of n(x(t), p(t)).In order to find an approximation for the full n(x, p) we (linearly) interpolate between those base points.model.In addition, using the following simple trick the algorithm will automatically deal with the possibility of multiple solutions: Instead of solving the impurity problem from scratch during each time step, one can start the iterative procedure from the solution obtained in the last time step.This way the state at the impurity will adiabatically follow the correct solution to the impurity problem over time.Furthermore, we find that it is sufficient to perform one iteration of the iteration scheme in each time-step.This gives an efficient hybrid algorithm which solves the GHD equation with impurity.For now this algorithm is first order, but it would be interesting to extend it to higher order methods as well [55].The density plots show the distribution of n(t, x macro , p) for t = 0, t = 1.5 and t = 3.For time t = 1.5 and t = 3 we also plot some trajectories of particles at the impurity, together with the current effective potential W (x) (inset).The red lines at p 0 = ± √ 2λ (solid) and p = 0 (dashed) are given for comparison with the non-interacting case.particles are reflected.Only towards the end of the simulation, at t = 3, the density decreases again, which decreases the effective potential and allows more particles to be transmitted.

Appendix B. Explicit construction of a solution for weak impurities
In this appendix we give the explicit construction of the solution for a sufficiently weak impurity discussed in Section 4.4.3,i.e. we want to show the following observation: Consider any (bounded) potential V (x, p) and scale it λV (x, p) with λ → 0. Then there exists a λ 0 > 0 s.t. for all |λ| < λ 0 there exists a solution to (47) where the states on the left and right side coincide.This means that all particles are unaffected by the impurity and no particles are reflected.
The computations below are a bit technical, the idea, however, is very simple.Since the impurity is weak the level sets of the effective Hamiltonian (i.e. the trajectories of particles) will be almost straight.Therefore no reflection occurs.This allows for the following construction: Let us call p L (−∞) the smallest incoming momentum from the left and p R (∞) the smallest incoming momentum from the right.Then there will be two curves p L (x) and p R (x) which separate the space into three distinct parts: In the region p > p L (x) we have particles coming from the left to the right and in the region p < p R (x) we have particles coming from the right.In the region p R (x) < p < p L (x) the density vanishes.Note that this construction only works if there are no reflected particles.
For convenience let us redefine the Hamiltonian by adding a constant, s.t.H(−∞, v p = 0) = 0. Also let us redefine: where in both definitions we use different branches of H −1 (−∞, •).Then we can write: N R (H(x, p)) p < p R (x) 0 p R (x) < p < p L (x).

(B.3)
We now have to find a solution to the adapted fixed point equation: dq T (p − q)N L (H(−∞, q) + λV (x, q) + U (x, q)) dq T (p − q)N R (H(−∞, q) + V (x, q) + U (x, q)) dq T (p − q)N L (H(−∞, q)) dq T (p − q)N R (H(−∞, q)) . (B.4) Note that p L/R implicitly depend on U (x, p).For λ = 0 a solution to this equation is simply given by U (x, p) = 0 and p L/R (x) = p L/R (−∞).We already established in Section 4.4.2 that the Jacobian is given by the inverse dressing operation.Because the dressing has to be defined for λ = 0, we know that the Jacobian is invertible at λ = 0. Therefore, there will be a small region |λ| < λ 1 in which a solution to the fixed point equation will exist (mathematically this follows from the implicit function theorem), which will also be close to the unperturbed solution |U | ≤ U 0 = O(λ).This establishes, for sufficiently small λ, the existence of the effective Hamiltonian H(x, p) describing the scattering situation (even though we do not know its explicit form).
However, this Hamiltonian was found on the assumption that there is no reflection.For consistency we now need to check that this Hamiltonian H(x, p) = H(−∞, p) + λV (x, p) + U (x, p) does not have any critical points in the regions p > p L (x) and p < p R (x) (the existence of such a critical point is necessary for reflection).Without loss of generality consider the region p > p L (x) (the other one can be treated equivalently).The boundary of the region is determined via: In case λ = 0 then p L (x) = p L (−∞) and ∂ p E 0 (p L (−∞)) ≥ C. By continuity, if λ is small (and therefore U is small) there will be a region |λ| < λ 2 in which |λ∂ p V (x, p) + ∂ p U (x, p)| < C and therefore H(x, p) does not have any critical points for |λ| < λ 2 .
The same is true for p < p R (x) and thus our ansatz of H(x, p) is consistent: For sufficiently small λ there cannot be any critical points and therefore no 'islands' Λ L/R i and also no reflection.
This finishes the construction of the solution to the scattering problem for sufficiently small λ.We now go on to show that for this solution the state on the left and on the right side are actually coincide (meaning that all particles just flow 'around' the impurity, without being affected by it).Note that this is not trivial: For instance, one could imagine that all particles are transmitted but their final momenta might be shifted compared to their initial momenta.In order to study this let us look again at equation (B.4).This fixed point equation is in fact independent for different x, which -in this case -holds globally since there are no critical points of the Hamiltonian in the regions where particles are.Now imagine starting at x = −∞ and observing how the fixed point U (x, p) varies as x runs from x = −∞ to x = ∞.U (x, p) change adiabatically as function of λV (x, p), but will stay in some finite region |U | ≤ U 0 , where U 0 depends on λ.Again, since the Jacobian is given by inverse dressing operation, and the dressing is invertible at V = 0, the Jacobian will be invertible for all impurity potentials V (x, p) in a neighbourhood around V = 0 as well (mathematically this follows from the inverse function theorem).In this neighbourhood the solution U (x, p) will be a unique function U (x, p) = U [λV (x, p)](x, p) of λV (x, p).If λ is small enough |λ| < λ 3 then λV (x, p) will stay inside this neighbourhood.Since for x → ∞ we have that V (x, p) → 0, far away from the impurity the solution is given by U (x → ∞, p) = U [0](x, p).By construction however, we know that U [0](x, p) = 0. Therefore, the Hamiltonians on the left and the right side will be equal, which in turn implies the same for the states on both sides.This concludes the derivation.

Figure 1 :
Figure 1: Schematic sketch of transformation of solution under local rescaling: a) Triangle potential and its stationary solution (sketch), b) Gaussian potential and solution obtained via rescaling of a), c) Potential with local minimum and solution obtained via rescaling of a).Note that the density increases again in the minimum.These particles there are trapped in the potential and only passively contribute to the scattering.A stationary state without trapped particles also exists, but may look very different.

Figure 2 :
Figure2: Contour plot of an effective Hamiltonian H(x, p) (for illustrative purpose only -Hamiltonian does not correspond to a specific system).The corresponding 3D plot is shown in Figure3.The contours coincide with the trajectory of quasi-particles.We denote by Σ L (Σ R ) all contours corresponding to incoming particles on the left (right).Both regions are separated by the curve γ(s) (red), which is the boundary between reflected and transmitted particles.The 'islands' (gray) are regions contours whose energy is too high or too low to be reached by the incoming particles.The density of particles there is zero.

.
Then it only shifts around the boundaries in a continuous fashion which gives a continuous change of the functional F [U ].Also in the other cases F is continuous.For instance, when the change in U creates a new domain Λ L/R i we tear a 'hole' into Σ L/R and then fill the interior with a constant value of N s.t. the function N is still continuous.This operation also only changes N (x, p) continuously and thus F [U ] is continuous.Now consider a bump function ∆U (x, p) with compact support and let us study how dx dp F [U + λ∆U ](x, p)∆U (x, p) behaves as λ → ±∞.In both cases it is easy to see that dx dp F [U + λ∆U ](x, p)∆U (x, p) → ±∞.

Observation 3
Consider a state n(p) s.t.n(p) is identically zero in a finite region around the zero velocity momentum, i.e. the density is only non-zero for |∂ p E 0 | > C.

Figure 5 :
Figure 5: Comparison of microscopic hard rods simulation (averaged over 100 samples) and the corresponding result obtained by simulation the GHD equation with impurity boundary conditions for d = 0.3: a) Comparison of the quasi-particle density ρ.For the hard rods simulation this was obtained by doing a 2D histogram with 100 bins in both dimensions.The red lines indicate the momenta p = ± √ 2λ whose kinetic energy corresponds to the height of the potential λ = 0.4.This is given for comparison with noninteracting particles d = 0 (as in Appendix A.3). b) Quasi-particle density ρ integrated over the regions R (reflected particles) and T (transmitted particles) indicated in a).The solid line is the GHD result and the datapoints with errorbars are the mean and standard deviation of the Monte Carlo result for different L.

Figure A1 :
Figure A1: Simulation of hard rods scattering at an impurity V (x) = λe − 1 2 x 2 with λ = 0.4, for d = 0 (non-interacting), d = 0.3 and d = −0.3.The density plots show the distribution of n(t, x macro , p) for t = 0, t = 1.5 and t = 3.For time t = 1.5 and t = 3 we also plot some trajectories of particles at the impurity, together with the current effective potential W (x) (inset).The red lines at p 0 = ±

Observation 4
Consider a state n(p) s.t.n(p) is identically zero in a finite region around the zero velocity momentum, i.e. the density is only non-zero for |∂ p E 0 | > C.