Main

Two connected reservoirs exchanging particles and energy is a paradigmatic system that is key to understanding transport phenomena in diverse platforms of both fundamental and technological interest ranging from heat engines to superconducting qubits1 and even heavy-ions collisions2. Entropy and heat, both irreversibly produced and transported by the currents flowing between the reservoirs, are key quantities in superfluid and superconducting systems3,4. They help to reveal microscopic information in strongly interacting systems5,6 and more generally characterize far-from-equilibrium systems7. Yet in traditional condensed matter systems such as superconductors and superfluid helium, the entropy is not directly accessible and requires indirect methods to deduce it8,9,10.

In this work, we leverage the advantage of quantum gases of ultracold atoms as naturally closed systems well-isolated from their environments to study entropy transport and production in fermionic superfluid systems. Using the known equation of state11, we measure the particle number and total entropy in each of the two connected reservoirs as a function of evolution time, therefore directly obtaining the entropy current and production. In general, the nature of these currents depends on the coupling strength between the superfluids. On the one hand, two weakly coupled superfluids exhibiting the Josephson effect12,13 exchange an entropy-free supercurrent described by Landau’s hydrodynamic two-fluid model14,15,16,17,18. In quantum gases19 as well as superconductors20, this is accomplished with low-transparency tunnel junctions weakly biased in chemical potential or phase, while narrow channels are used to block viscous currents in superfluid helium21,22. On the other hand, superfluids strongly coupled by high-transparency channels23 can exhibit less intuitive behaviour since the supercurrent no longer dominates the normal current, making the system fundamentally non-equilibrium24,25. In particular, the signature of superfluidity in such systems is often large particle currents on the order of the superfluid gap which respond nonlinearly to chemical potential biases smaller than the gap. This is observed in ballistic junctions between superconductors20, superfluid He26 and quantum gases27,28,29. However, entropy transport in this setting has so far only been experimentally studied indirectly and at higher temperatures in the linear response regime where the superfluidity of the system is ambiguous30,31, leaving open the question of entropy transport between strongly coupled superfluids.

Here, we connect two superfluid unitary Fermi gases with a ballistic channel and measure the coupled transport of particles and entropy between them. We observe large subgap currents of both particles and entropy, indicating that the current is not a pure supercurrent and cannot be understood within a hydrodynamic two-fluid model. In particular, superfluidity counterintuitively enhances entropy transport in this system by enhancing particle current while maintaining a large entropy transported per particle. We also observe in our parameter regime that, while the system can always thermalize via the irreversible flow of this superfluid-enhanced normal current, thermalization via pure entropy diffusion is inhibited in one-dimensional (1D) channels, giving rise to a non-equilibrium steady state previously observed in the normal phase30. The observed nonlinear dynamics of particles and entropy are captured by a phenomenological model we develop whose only external constraints are the conservation of particles and energy and the second law of thermodynamics.

We begin the experiment by preparing a balanced mixture of the first- and third-lowest hyperfine ground states of 6Li at unitarity in an augmented harmonic trap shown in Fig. 1a (Methods and Supplementary Information section 3). To induce transport of atoms, energy and entropy between the left (L) and right (R) reservoirs, we prepare an initial state within the state space in Fig. 1b characterized by the conserved total atom number and energy, N = NL + NR and U = UL + UR, and the dynamical imbalances in the atom number ΔN = NL − NR and entropy ΔS = SL − SR. The imbalances in the extensive quantities induce biases in the chemical potential Δμ = μL − μR and temperature ΔT = TL − TR according to the reservoirs’ equations of state (EoS; Methods) which in turn drive currents of the extensive properties INN, ΔS) = − (1/2)dΔN/dt and ISN, ΔS) = − (1/2)dΔS/dt. Note that IS is an apparent current, not a conserved current like IN and IU, though we can place bounds on the conserved entropy current from the apparent current and entropy production rate \({I}_{S}^{{{\;{\rm{cons}}}}}\in [{I}_{S}-(1/2){{{\rm{d}}}}S/{{{\rm{d}}}}t,{I}_{S}+(1/2){{{\rm{d}}}}S/{{{\rm{d}}}}t]\) (ref. 29). These equations of motion, together with the initial state ΔN(0), ΔS(0), determine the path the system traces through state space ΔN(t), ΔS(t) as well as the speed with which it traces this path. The paths that we explore experimentally are shown as black lines overlayed on top of the entropy landscape S = SL(NL, UL) + SR(NR, UR) = S(N, U, ΔN, ΔS) computed from the EoS. The paths all exhibit a strictly positive entropy production rate dS/dt > 0, indicating that the transport is irreversible, until they reach either a non-equilibrium steady state (ΔN, ΔS ≠ 0) or equilibrium (ΔN = ΔS = 0) where S is maximized for the fixed N and U. The von Neumann entropy of a closed system does not increase in time under Hamiltonian evolution, though the measured thermodynamic entropy S can increase due to buildup of entanglement entropy shared between the two reservoirs32,33,34,35. We have verified that the system is nearly closed given the measured particle loss rate dN/dt/N < 0.01 s−1 and heating in equilibrium d(S/NkB)/dt < 0.02 s−1, limited by photon scattering of optical potentials, such that the entropy production observed during transport is due to the fundamental irreversibility of the transport.

Fig. 1: Irreversible particle, energy and entropy transport through a non-equilibrium channel connecting two superfluids.
figure 1

a, Slices along the x = 0 and z = 0 planes of the calculated local degeneracy T/TF during transport through the 1D channel, showing that the reservoirs (left L and right R) are in local equilibrium in the normal phase (red, T/TF > 0.167) over most of their volume and superfluid (blue, T/TF ≤ 0.167) at the contacts to the channel. The channel (green) is far from local equilibrium. Differences in the atom number N, energy U and entropy S between the reservoirs induce currents between them. b, Depending on the initial conditions (filled or open circles) and microscopically allowed processes, the reservoirs can exchange entropy advectively and diffusively, tracing out various paths through state space with the constraints that atom number and energy are conserved dN/dt = dU/dt = 0 and the entropy production is positive definite dS/dt ≥ 0. This evolution halts (stars) by reaching either a non-equilibrium steady state (ΔN, ΔS ≠ 0) or equilibrium (ΔN = ΔS = 0) where the total entropy S = SL + SR has a global maximum.

We measure the evolution of the system by repeatedly preparing the system in the same initial state, allowing transport for a time t, then taking absorption images of both spin states and extracting Ni, Si for both reservoirs i = L, R using standard thermometry techniques (Methods and Supplementary Information section 5B). Between the end of transport and imaging, we adiabatically ramp down the laser beams that define the channel to image the reservoirs in well-calibrated, half-harmonic traps. The beams do work on the reservoirs during this process and change Ui but Ni and Si remain constant. The cloud typically contains N = 270(30) × 103 atoms and S/N = 1.59(7)kB before transport, below the critical value of ~1.90kB for superfluidity in the transport trap (numbers in brackets represent statistical uncertainties). Figure 1a shows the local degeneracy T/TF in the x = 0 and z = 0 planes during transport calculated within the local density approximation36 using the three-dimensional (3D) equation of state11 to determine the local Fermi temperature TF(x, y, z) (Methods). The imbalance is illustrated using the representative values of νx = 12.4 kHz, Δμ = 75 nK × kB and ΔT = 0. Assuming local equilibrium, the most degenerate regions in the channel would be deeply superfluid due to the strong, attractive gate potential and reach T/TF ≈ 0.027, s ≈ 7.2 × 10−4kB and Δ/kBT ≈ 17 where s is the local entropy per particle and Δ is the superfluid gap. When νx 7 kHz (Methods and Supplementary Information Section 3) and the channel is two-dimensional (2D), normal regions appear at the edges of the channel that can also contribute to transport.

In a first experiment, we prepare an initial state ΔN(0), ΔS(0) ≠ 0 (filled circle in Fig. 1b) such that equilibrium is reached within 1 s. For the strongest confinement, ΔN(t), shown in Fig. 2a, clearly deviates from exponential relaxation and the corresponding IN is much larger than the value ~Δμ/h of a quantum point contact in the normal state, where h is Planck’s constant, indicating that the subgap current-bias characteristics are nonlinear (non-Ohmic) and the reservoirs are superfluid28,29,37. When reducing νx to cross over from a 1D to 2D channel, ΔN(t) relaxes faster (IN increases) and, although it is less pronounced, the nonlinearity persists. The dynamics for the two smallest values of νx are nearly identical, suggesting that there are additional resistances in series with the 1D region such as the viscosity of the bulk reservoirs or the interfaces between the 3D and 2D regions31 or between the normal and superfluid regions38.

Fig. 2: Observing advective entropy current from the non-exponential evolution of imbalances in particle number and entropy.
figure 2

a,b, Particles (a) and entropy (b) are transported between the reservoirs by currents that respond nonlinearly to biases in chemical potential and temperature, evidenced by the non-exponential decay. Both currents increase as the channel opens from 1D to 2D (νx decreasing). c, The net entropy increases very slightly while ΔS changes significantly, indicating that entropy is indeed transported by the particle flow. d, Showing IS = s*IN where the entropy advectively transported per particle s* is nearly independent of the channel geometry νx despite IN varying significantly. The observation that s* > 0 means the nonlinear particle current is not a superfluid current. Each data point is an average of three to five repetitions and error bars represent the standard deviation. Solid lines are fits of the phenomenological model.

Source data

Concurrently with ΔN(t), we observe non-exponential relaxation of ΔS(t) (Fig. 2b) which bears a remarkable resemblance to ΔN(t). We find that by plotting ΔS(t) against ΔN(t) in Fig. 2d, all paths collapse onto a single line. This demonstrates that the entropy current is directly proportional to the particle current IS = s*IN where the average entropy advectively transported per particle s* is nearly independent of νx even though IN itself varies significantly. Moreover, dS/dt (Fig. 2c) is barely resolvable and is significantly smaller than IS, meaning there is indeed a large conserved entropy current flowing between the reservoirs. The dependence of S/NkB on the confinement νx in Fig. 2c has a technical origin and does not affect the system during transport (Methods). The entropy transported per particle s* = 1.18(3)kB is near its value in the normal phase30 and is orders of magnitude larger than the local entropy per particle in the channel assuming local equilibrium s ≈ 7.2 × 10−4kB. Because superfluidity in the contacts enhances IN while only slightly suppressing s*, superfluidity increases IS.

The fact that s* > 0 directly shows that the large, non-Ohmic current between the two superfluids is itself not a pure supercurrent in the context of a two-fluid model39. The observation that the flow is resistive is insufficient alone to conclude that it is not superfluid as there are many mechanisms for resistance to arise in a pure supercurrent40,41. The observation that s*s suggests that the channel is far from equilibrium and hydrodynamics breaks down as is often the case in weak link geometries39 unlike previous assumptions31,38. We discuss in Supplementary Information section 1 how the degree to which hydrodynamics breaks down depends on the preparation of the system. The large entropy current suggests an irreversible conversion process from superfluid currents in the contacts to normal currents in the channel and back to superfluid, or the propagation of normal currents originating in the normal regions of the reservoirs through the superfluids while remaining normal. Moreover, the independence of s* from νx implies that this process is independent of the channel geometry. There is an analogy between this observation and the central result of Landauer–Büttiker theory that the conductance through a ballistic channel is also independent of the geometry and depends only on the channel’s transmission and the number of propagating modes.

In a second experiment, whose results are presented in Fig. 3, we prepare the system with a nearly pure entropy imbalance (ΔN(0) ≈ 0, ΔS(0) ≠ 0, open circle in Fig. 1b). The initial response of ΔN and ΔS from t = 0 to when ΔN reaches its maximum value is clearly non-exponential, resembling the advective dynamics in the first experiment with the same s*, while the dynamics that follow are much slower and consistent with exponential relaxation and therefore linear response. With decreasing νx, both dynamics become faster and the maximum values of ΔN and ΔS achieved at the turning point become smaller. For the largest value of νx, the initial response is still fast while the relaxation that follows is extremely slow and resembles a non-equilibrium steady state over experimentally accessible timescales: the relaxation time of this state is 8(2) s while it is reached from the initial state in only ~0.2 s. In this state, the non-vanishing imbalances ΔN and ΔS depend on the initial state as well as the path of the system through state space determined by s*, that is, the system is non-ergodic. This indicates that the non-equilibrium steady state previously observed at higher temperatures30 persists in the superfluid regime where the current with which it is reached is 6 times larger and non-Ohmic. Figure 3f shows the measured path in state space also illustrated in Fig. 1b. It demonstrates that the path is determined by the competition between the nonlinear and linear dynamics and varies with νx, in contrast to the first experiment (Fig. 2d).

Fig. 3: Competition between two modes of entropy transport: nonlinear advection and linear diffusion.
figure 3

a,b, A pure entropy imbalance induces currents of particles (a) and entropy (b) which are dominated by the nonlinear advective mode at early time and the linear diffusive mode at long times. For strong confinement (red) the system reaches a non-equilibrium steady state which persists beyond the experimentally accessible time, indicating the decoupling of advective and diffusive entropy transport. Reducing the confinement facilitates the diffusive mode, allowing equilibration and reducing the maximal response. c,d, Enlarged views of the initial dynamics of a,b for additional clarity. e, As in Fig. 2c, the increase in net entropy indicates the irreversibility of the transport process. A fit with our phenomenological model (solid lines) describes the data well over the explored parameter regime. f, The competition between the two modes, which changes with νx, determines the path traced through state space. Error bars (data points) indicate the standard deviation (average) of three to five repetitions.

Source data

In the following, we formulate a minimal phenomenological model to describe our observations which are not captured by the linear response approach that successfully describes this system in the normal state30,31 as it predicts purely exponential relaxation. We therefore turn to the formalism of generalized gradient dynamics42, a generalization of Onsager’s theory of irreversible processes (Methods and Supplementary Information section 2). While it does not provide a microscopic theory, this formalism can describe general, irreversible, non-equilibrium processes and provides a convenient way to impose macroscopic constraints such as the second law of thermodynamics and conservation laws for the particle number and energy. Within this framework, we make the Ansatz

$$\begin{array}{rcl}{I}_{N}&=&{I}_{{{{\rm{exc}}}}}\tanh \left(\frac{\Delta \mu +{\alpha }_{\rm{c}}\Delta T}{\sigma }\right)\\ {I}_{S}&=&{\alpha }_{\rm{c}}{I}_{N}+{G}_{\rm{T}}\Delta T/T\end{array}$$
(1)

which produce entropy via the irreversible flow dS/dt = (INΔμ + ISΔT)/T. The non-trivial result that αc appears in both IN and IS is a generalization of Onsager’s reciprocal relations to nonlinear response and is a consequence of the irreversibility of these currents. The system exhibits two modes of entropy transport: a nonlinear advective mode \({I}_{S}^{\,\rm{a}}={\alpha }_{\rm{c}}{I}_{N}\) characterized by the excess current Iexc, Seebeck coefficient αc and nonlinearity σ, wherein each transported particle carries entropy s* = αc on average, and a linear diffusive mode \({I}_{S}^{\,\rm{d}}={G}_{\rm{T}}\Delta T/T\) characterized by the thermal conductance GT which enables entropy transport without net particle transport according to Fourier’s law. The linear model is reproduced in the limit of large σ with conductance G = Iexc/σ. In a Fermi liquid, these two modes are related by the Wiedemann–Franz law where the Lorenz number L = GT/TG has the universal value \({\uppi }^{2}{k}_{\rm {B}}^{2}/3\). The nonlinearity implies the breakdown of the Wiedemann–Franz law since the advective and diffusive modes are no longer linked30. The excess current Iexc is the particle current with the nonlinearity saturated, as in superconducting weak links40.

Figure 4 shows the parameters of the model as functions of νx extracted from the fits shown as curves in Figs. 2 and 3. Panel a shows Iexc normalized by the fermionic superfluid gap Δ/h along with the number of occupied transverse modes at equilibrium nm (Methods). Filled (open) circles were extracted from the first (second) experiment. Iexc follows nmΔ/h, increasing as νx decreases until νx ≈ 4 kHz where it plateaus, likely due to additional resistances in series with the 1D region. The fitted Iexc is apparently reduced in the second experiment relative to the first because the initial current is suppressed by the diffusive mode, making it more difficult to fit. It is intriguing that, consistent with previous studies28,29, the equilibrium superfluid gap Δ within the local density approximation is still the relevant scale for the current, despite the evidence that the channel region is far from equilibrium.

Fig. 4: Characterizing the linear and nonlinear transport modes by fitting the phenomenological model across the 1D–2D crossover.
figure 4

Fitted parameters from the first (second) experiment corresponding to Fig. 2 (Fig. 3) are plotted with filled (open) markers here. a, The excess current normalized to the superfluid gap. b, The thermal and spin conductances normalized to the non-interacting values \({G}_{\rm{T}}^{0}=2{\uppi }^{2}{k}_{\rm {B}}^{2}T/3h\), \({G}_{\sigma }^{0}=2/h\) for a single transverse mode (see main text). For reference, the number of occupied transverse modes in the channel nm is also shown in a and b. c, The Seebeck coefficient (advectively transported entropy per particle) extracted directly from the fits to the first experiment. d, The slope of the path through state space in the first experiment and in the advective (boxes) and diffusive (diamonds) limits of the second experiment. Error bars represent standard errors from the fits.

Source data

Figure 4b shows GT and the spin conductance Gσ separately measured in the same system by preparing a pure spin imbalance (Methods) normalized to their values for the single-mode non-interacting ballistic quantum point contact \({G}_{\rm{T}}^{0}=2{\uppi }^{2}{k}_{\rm{B}}^{2}T/3h,\,{G}_{\sigma }^{0}=2/h\). Both conductances are suppressed relative to the non-interacting values and increase monotonically with decreasing νx, possibly due to the appearance of non-degenerate transverse modes at the edges of the channel (Methods and Supplementary Information section 3). The non-equilibrium steady state arises from the fact that GT → 0 in the 1D limit. The relative increase of GT with decreasing νx is larger than that of Gσ, suggesting that more types of excitations can contribute to diffusive entropy transport than spin transport, for example, both collective phonon and quasiparticle excitations can contribute to GT (ref. 43) while only quasiparticle excitations can contribute to Gσ (ref. 44).

Figure 4c shows the fitted Seebeck coefficient αc, while Fig. 4d shows the slope of the path through state space dΔS/dΔN during the advective and diffusive dynamics. The fitted αc and dΔS/dΔN match for the purely advective transport in the first experiment, showing that αc is remarkably insensitive to νx, while dΔS/dΔN more clearly shows how the two modes compete in the second experiment to determine the net response of the system. Figure 4d shows that, while both modes are generally present in the system’s dynamics, their relative prevalence depends on νx as well as the initial state: the initial state in the first experiment was carefully chosen to allow the system to relax to equilibrium via the advective mode alone by preparing ΔS(0) = αcΔN(0) while the initial state in the second was chosen to contain both modes.

In summary, we have observed that the conceptually simple system of two superfluids connected by a ballistic channel exhibits the highly non-intuitive and currently unexplained effect that the presence of superfluidity increases the rate of irreversible entropy transport between them via nonlinear advection. This contrasts with the more familiar case of superfluid and superconducting tunnel junctions where the reversible, entropy-free Josephson current dominates. The entropy advectively carried per particle is nearly independent of the channel’s geometry, while the timescales of advective and diffusive transport depends strongly thereon, raising the question of the microscopic origin of the observed entropy transported per particle s* 1kB. Our phenomenological model that captures these observations, in particular the identification of advective \({I}_{S}^{\,\mathrm{a}}\propto {I}_{N}\) and diffusive \({I}_{S}^{\,\mathrm{d}}\propto \Delta T\) modes along with the sigmoidal shape of INμ + αcΔT), may help guide future microscopic theories of the system. While extensive research has been conducted on the entropy producing effects of topological excitations of the superfluid order parameter8,12,21,39,41,45,46,47,48,49, less attention has been given to their influence on entropy transport and the possible pair-breaking processes they can induce. Early studies of superconductors found that mobile vortices can advectively transport entropy by carrying pockets of normal fluid50 with them24,51,52,53. More generally, entropy-carrying topological excitations, which give rise to a finite chemical potential bias according to the Josephson relation Δμ = hNv/dt (refs. 39), where Nv is the number of vortices, can result from a complex spatial structure of the order parameter49. Alternatively, it is possible that an extension of microscopic theories of multiple Andreev reflection37,54, which reproduce the finding that the excess current scales linearly with the number of channels and the gap28,29, may explain our observations. Clearly, a proper microscopic theory of this system is a challenge for the future. A complete understanding of the particle and entropy transport in superfluid systems is essential for both fundamental and technological purposes.

Methods

Transport configuration

The atoms are trapped magnetically along y and optically by a red-detuned beam along x and z, with confinement frequencies νtrap,x = 171(1) Hz, νtrap,y = 28.31(2) Hz and νtrap,z = 164(1) Hz. A pair of repulsive TEM01-like beams propagating along x and z, which we call the lightsheet (LS) and wire respectively, intersect at the centre of the trapped cloud and separate it into two reservoirs connected by a channel. The transverse confinement frequencies at the centre are set to νz = 9.42(6) kHz (kBT/hνz = 0.21) and νx = 0.61…12.4(2) kHz (kBT/hνx = 0.16…3.3) with the powers of the beams. An attractive Gaussian beam propagating along z with a similar size as the LS beam acts as a gate potential in the channel. The peak gate potential is \({V}_{{{{\rm{gate}}}}}^{\;0}=-2.17(1)\,\upmu {{{\rm{K}}}}\times {k}_{\rm{B}}\). We use a wall beam, which is thin along y and wide along x, during preparation and imaging to completely block transport with a barrier height \({V}_{{{{\rm{wall}}}}}^{\;0}\) much larger than μ and kBT. The repulsive LS, wire and wall are generated using blue-detuned 532 nm light while the attractive gate is created with red-detuned 766.7 nm light.

The effective potential energy landscape along y at x = z = 0 (Supplementary Information section 3) is approximately Veff(y, nx, nz) ≈ hνx(y)(nx + 1/2) + hνz(y)(nz + 1/2), where νx and νz vary along y due to the beams’ profile and nx, nz are the quantum numbers of the harmonic potential in x and z directions. The number of occupied transverse modes nm (Fig. 4) is calculated via the Fermi–Dirac occupation with local chemical potential set by Veff(y, nx, nz),

$${n}_{\rm{m}}=\mathop{\sum }\limits_{{n}_{x},{n}_{z}=0}^{\infty }\mathop{\min }\limits_{y}\frac{1}{1+\exp \left\{[{V}_{{{{\rm{eff}}}}}(\;y,{n}_{x},{n}_{z})-\mu ]/{k}_{\rm{B}}T\right\}}$$
(2)

where the minimum occupation of each mode is used to account for modes that are not always occupied throughout the channel (Supplementary Information section 3).

The complete potential energy landscape V(r), where r = (x, y, z), was used to produce Fig. 1a via the local density approximation for the density n(r) = n[μ − V(r), T] (refs. 11,36) that determines the local Fermi temperature \({k}_{\rm{B}}{T}_{\rm{F}}({{{\bf{r}}}})={\hslash }^{2}{[3{\uppi }^{2}n({{{\bf{r}}}})]}^{2/3}/2m\), where m is the atomic mass. The superfluid gap Δ assuming local equilibrium is estimated using the calculation in a homogeneous system Δ(μc, T) (ref. 55) where \({\mu }_{\rm{c}}=\mathop{\max }\limits_{{{{\bf{r}}}}}[\mu -V({{{\bf{r}}}})]\) is the maximum local chemical potential in the system. The crossover between 1D and 2D regimes (νx ≈ 7 kHz) of the channel is estimated by comparing the local degeneracy along x (at y = z = 0) to the superfluid transition. In the 2D limit, non-superfluid modes can pass through the edges of the channel while in the 1D limit the degeneracy across the channel is below the superfluid transition (Supplementary Information section 3).

Transport preparation

To prepare imbalances ΔS(0), ΔN(0) we ramp up the channel beams to separate the two reservoirs followed by forced optical evaporation. Using a magnetic field gradient along y, we shift the centre of the magnetic trap with respect to the channel beams before separation to prepare ΔN(0). By shifting the trap centre during evaporation, we can compress one reservoir and decompress the other, thereby changing their evaporation efficiencies and inducing a controllable ΔS(0). See Supplementary Information Section 4 for more details. To measure the spin conductance Gσ (Fig. 4b), we prepare a ‘magnetization’ imbalance ΔM = ΔN − ΔN. To do this, we ramp down the magnetic field before separating the reservoirs at 52 G where the spins’ magnetic moments are different and modulate a magnetic gradient along y until the two spins oscillate out of phase. We then separate the two reservoirs and ramp back the magnetic field.

Imaging and thermometry

Between the end of transport and the start of imaging, we ramp down the channel beams while keeping the wall on. At the end of each run, we obtain the column density \({n}_{i\sigma }^{{{{\rm{col}}}}}(\;y,z)\) of both reservoirs i = L, R and both spin states σ = ,  (first and third-lowest states in the ground state manifold) from two absorption images taken in quick succession in situ. We fit the degeneracy qiσ = μiσ/kBTiσ and temperature Tiσ of both reservoirs for each spin state using the EoS of the harmonically trapped gas11. However, we use the fitted temperature from the first image () for both spins since the density distribution in the second image is slightly perturbed by the first imaging pulse. The thermometry is calibrated using the critical S/N of the condensation phase transition on the BEC side of the Feshbach resonance. See Supplementary Information section 5B for more details.

Generalized gradient dynamics

To ensure that the phenomenological nonlinear model satisfies basic properties such as the second law of thermodynamics, it is formulated it in terms of a dissipation potential Ξ (ref. 42). The thermodynamic fluxes are defined as derivatives of the dissipation potential Ξ with respect to the forces IN = TΞ/∂Δμ and IS = TΞ/∂ΔT. In this formalism, Onsager reciprocity and the conservation of particles and energy are fulfilled. Our model (equation (1)) is the result of the following dissipation potential, which is constructed based on the experimental observation IS = s*IN and that IN follows a sigmoidal function of Δμ (ref. 28),

$$\varXi =\frac{\sigma {I}_{{{{\rm{exc}}}}}}{T}\log \left[\cosh \left(\frac{\Delta \mu +{s}^{* }\Delta T}{\sigma }\right)\right]+\frac{{G}_{\rm{T}}}{2}{\left(\frac{\Delta T}{T}\right)}^{2}$$
(3)

where the first part describes the advective and the second part the diffusive transport mode. See Supplementary Information section 2 for more details.

Reservoir thermodynamics

To formulate equations of motion in state space (ΔN, ΔS), we relate Δμ, ΔT to ΔN, ΔS in terms of thermodynamic response functions

$$\left(\begin{array}{c}\Delta N\\ \Delta S\end{array}\right)\approx \frac{\kappa }{2}\left(\begin{array}{cc}1&{\alpha }_{\rm{r}}\\ {\alpha }_{\rm{r}}&{\ell }_{\rm{r}}+{\alpha }_{\rm{r}}^{2}\end{array}\right)\left(\begin{array}{c}\Delta \mu \\ \Delta T\end{array}\right)$$
(4)

where κ is the compressibility, αr is the dilatation coefficient and r is the ‘Lorenz number’ of the reservoirs56 (Supplementary Information section 5A). To obtain the spin conductance, we use ΔM = (χ/2)Δb, where Δb = (Δμ − Δμ)/2. The spin susceptibility χ ≈ 0.32κ, following the computed EoS of a polarized unitary Fermi gas57.

The potential landscape V(r) during the transport experiment deviates from simple harmonic potential due to the confinement beams as well as the anharmonicity of the optical dipole trap. We estimate from numeric simulations based on our knowledge of the V(r) that T and κ agree within 1% to those determined from absorption imaging in near-harmonic traps while μ is 24% higher during transport. However, αr and r are more sensitive to the trap potential and can deviate by a factor of 3. We therefore fit these response coefficients in our model. See Supplementary Information section 5A for more details.

The total entropy S, being a state variable (contour plot in Fig. 1b), depends on the imbalances ΔN(t) and ΔS(t) but not the currents IN and IS. With the linearized reservoir response, the entropy produced by equilibration is given by

$$S(t)\approx {S}_{{{{\rm{eq}}}}}-\frac{\Delta {N}^{2}(t)}{2T\kappa }-\frac{{[\Delta S(t)-{\alpha }_{\rm{r}}\Delta N(t)]}^{2}}{2T{\ell }_{\rm{r}}\kappa }$$
(5)

where Seq is the maximum entropy at equilibrium given fixed total N and U. The increase in S/NkB with decreasing νx in Figs. 2c and 3c is caused by switching on the wall beam to block transport. For lower νx, there are more atoms in the channel to be perturbed by this process.

Fitting procedure

We fit the phenomenological model (equation (1)) with linear reservoir responses (equation (4)) to each dataset—the set of different transport times at fixed νx—independently for both the first and second experiment. We do this by solving the initial value problem for ΔN(t) and ΔS(t) given the parameters αc, Iexc, σ and GT along with the reservoir response functions κ, αr and r. From these solutions, we also compute the total entropy S(t) as a function of time (equation (5)). We fit ΔN(t), ΔS(t) and S(t) simultaneously to the data using a least-squares fit. For the first experiment, only Iexc, σ, αc, ΔS(0) and Seq are free parameters. Other parameters are fixed to their theoretical values for better fit stability. For the second experiment, we fit σ, ΔS(0), Seq, GT, αr, r and an offset in ΔS to account for drifts in alignment. αc is fixed to the averaged value obtained in the first experiment. See Supplementary Information section 6 for more details. The slopes shown in Fig. 4d are obtained by simple linear fits in the state space ΔS versus ΔN (Figs. 2d and 3f). The advective and diffusive modes in the second experiment are separated in time at the maximum ΔN(t) (Fig. 2a).