Finite-time Landauer principle beyond weak coupling

Landauer's principle gives a fundamental limit to the thermodynamic cost of erasing information. Its saturation requires a reversible isothermal process, and hence infinite time. We develop a finite-time version of Landauer's principle for a bit encoded in the occupation of a single fermionic mode, which can be strongly coupled to a reservoir. By solving the exact non-equilibrium dynamics, we optimize erasure processes (taking both the fermion's energy and system-bath coupling as control parameters) in the slow driving regime through a geometric approach to thermodynamics. We find analytic expressions for the thermodynamic metric and geodesic equations, which can be solved numerically. Their solution yields optimal processes that allow us to characterize a finite-time correction to Landauer's bound, fully taking into account non-markovian and strong coupling effects.

The unattainability principle suggests that Landauer's bound cannot be saturated with finite resources, namely time and energy [30][31][32].In finite time, using tools from optimal transport theory [33][34][35] and thermodynamic geometry [36][37][38][39][40][41][42][43], optimal erasure protocols have been derived both for classical systems described by overdamped Langevin dynamics [44][45][46][47][48] and open quantum systems described by Lindblad master equations [35,[49][50][51][52][53][54].Such optimal protocols naturally lead to a finite-time correction to Landauer's bound in different physical set-ups, which has given rise to the term finite-time Landauer principle [45,48,52].For a slowly driven (quantum) two-level system weakly coupled to a thermal bath, the finite-time bound takes the simple form (see App. B.2 and Refs.[50,54]) where τ is the total time of the process and Γ is the thermalization rate.The finite-time correction is positive, in agreement with the second law of thermodynamics, and when Γτ → ∞ we recover the standard bound.We also note that the optimal protocol saturating the finite bound eq. ( 1) has been recently implemented in a semiconductor quantum dot [29].More general versions of eq. ( 1) have also been recently developed for Markovian systems driven at any speed [35,51,52].Despite this remarkable progress, previous works on the finite-time Landauer principle have focused in Markovian systems which, for quantum systems, can be guaranteed by a sufficiently weak interaction between system and bath.In the presence of strong coupling [55][56][57][58][59], we expect both new opportunities arising due to faster relaxation rates and non-Markovian dynamics [60][61][62][63][64][65][66][67], as well as challenges due to the presence of new sources of irreversibility [68][69][70][71][72][73][74][75].The goal of this work is to take a first step into this exciting regime by deriving the first order to a tight finite-time correction of Landauer's principle for a single fermion that can interact strongly with a reservoir, as described by the resonant-level model [76][77][78][79][80][81][82][83].Our main result is summarised in what follows:

Main result
Given a two-level system that can be strongly coupled to a thermal bath, we find that the finite-time version of Landauer's principle can be expressed as where a ≈ 2.57946, τ Pl = ℏ/k B T is the so-called Planckian time [84], and Γ is the average thermalization rate (see details below).This extends eq. ( 1) to strong system-bath couplings, with the transition between the two being characterized in Fig. 2. The finite-time correction in eq. ( 2) is of quantum-mechanical nature and independent of the coupling strength, hence prevailing even for arbitrarily strong system-bath coupling (roughly speaking, Γ → ∞ in eq. ( 1)).
The appearance of the Planckian time τ Pl = ℏ/k B T in eq. ( 2) is particularly interesting.This timescale encodes two fundamental constants in nature: Boltzmann's constant k B and Planck's constant ℏ.It arises in several contexts in manybody physics, including quantum transport and quantum chaos; see Ref. [84] for a review.In analogy with the "Planck time" in quantum gravity, it is associated with the shortest timescale of thermalization [84][85][86][87]; that is, the shortest time needed to redistribute energy between particles and reach thermal equilibrium.This gives an insightful context to our main result eq. ( 2): a fundamental finite-time quantum correction must appear to Landauer's bound due to a minimal time required for thermalization.This also suggests that the form of eq. ( 2) has a broader range of applicability, with the value a depending on the specific many-body thermalizing dynamics considered.

Framework
We consider a driven system S that can be put in contact with a thermal bath B, so that the total time-dependent Hamiltonian reads: ( Here Focusing on protocols where H int (0) = H int (τ ) = 0, we can naturally identify from the first law of thermodynamics W = Q + ∆E S with ∆E S = Tr[H S (τ )ρ(τ ) − H S (0)ρ(0))], the dissipated heat as the total energy absorbed by the bath [3].
Assuming that the initial state of SB is a thermal state: ρ(0) = e −βH(0) /Z(0) with Z(t) ≡ Tr[e −βH(t) ], eq. ( 4) can be re-expressed as [100]: where ∆F = k B T ln[Z(0)/Z(τ )] is the change of equilibrium free energy of SB, and the entropy production Σ can be expressed as: Σ = S ρ(τ ) e −βH(τ ) Z(τ ) . The entropy production Σ ≥ 0 accounts for the irreversible energetic contribution in finite-time processes, and depends on the particular driving path H(t) linking H(0) to H(τ ).Minimising Σ over all finite-time processes leads to thermodynamic protocols that minimize the work W . Furthermore, in an erasure process, ∆E S = 0 (see details below) therefore these protocols also minimize the dissipated heat Q.

Thermodynamic geometry
The framework of quantum thermodynamic geometry [42] allows us to minimize the entropy production Σ, and therefore the dissipated heat Q for erasure, for protocols that are slow compared to their relaxation time-scale.

Strongly coupled systems
Let us expand the Hamiltonian H(t) in eq. ( 3) as H(t) = j λ j (t)X j where {λ j (t)} are the externally controllable parameters and {X j } are the corresponding observables.In order to apply the geometric approach, we need to impose more structure on the possible evolutions U (t) generated by eq. ( 3).We require two basic ingredients: Requirement 1: Thermalization.In absence of driving, the conjugated observables X j thermalize.More precisely, for a frozen Hamiltonian H(t), we have where Ũt (s) ≡ e −iH(t)s , ⟨X j (t)⟩ eq = Tr[ω β (t)X j ], ω β (t) = e −βH(t) /Z(t), and β is implicitly defined by the initial energy of the total system.In the context considered here, namely purely unitary dynamics of SB, this condition is satisfied both by non-integrable systems satisfying the ETH hypothesis [101,102] and also for integrable systems typically appearing in open quantum systems [103][104][105][106].
Requirement 2: Slow driving, so that the system remains close to the instantaneous equilibrium state while being driven.This enables us to keep only leading terms when making a linearresponse expansion in the driving speed [107], which can be expressed as: The coefficients m ij , which depend on the point {λ i (t)}, can in principle be derived from the exact equations of motion.We should note that it is well known that optimal finite-time protocols feature jumps [108].However, these jumps disappear near the reversible limit (see App. E and [109]) and their contribution to the dissipated heat becomes negligible.Therefore this requirement becomes a natural assumption in the context of finding a first order correction to Landauer's bound.
Combining the expansion of eq. ( 8) with eq. ( 4) and eq.( 6), we obtain the standard expression for entropy production at leading order in the inverse of the driving speed [37,40,42]: where, in contrast to previous works, the metric m ij depends on the unitary dynamics of SB (this will later be solved for a specific model).
Because of the second law of thermodynamics, it follows that m ij can be expressed as a metric, i.e., a symmetric, positive-definite m ≥ 0 operator that depends smoothly on the point {λ i (t)}.
We can associate a length to a protocol by defining It is related to the entropy production via a Cauchy-Schwarz inequality [37,40,42]: where equality is satisfied by protocols with constant entropy production rate Furthermore, to minimize the entropy production of any (slow) protocol we have to find the shortest path between the desired initial and final value of the Hamiltonian's parameters.This corresponds to a geodesic path, with length L, which naturally defines a minimal entropy production We can find Σ min by solving the geodesic equation that is derived from the metric and computing its length [37,40,42].

Resonant-level model
Having explained the general ideas behind our work, we now focus on finite-time driving processes of a single fermionic mode coupled to a fermionic bath, which can e.g.describe a singleelectron quantum dot.The total Hamiltonian reads: where â † is the creation operator of the two-level system and b † k is the creation operator of a bath mode with frequency ω k , following the canonical anti-commutation relations: â} = 0; and finally λ k are the interaction weights which define the spectral density function of the bath The energy ε is the difference between the energy of the two-level system and the chemical potential of the bath 1 .We are assuming optimal control over the functions ε(t) and g(t) so that we can fully optimize the protocol and reach the fundamental limit for this system.While this level of control is, in principle, ambitious experimentally in regards to the coupling, it has been achieved in quantum dots [110] where the tunneling rate (i.e.interaction strength) can be modified by several orders of magnitude.We take the continuum limit and assume that the spectral density of the bath is a Lorentzian where Λ > 0 is a parameter characterizing its width.Exact and explicit solutions for the resonant-level model are known in the wide-band limit Λ → ∞ [76][77][78][79][80][81][82].This limit is commonly used to describe quantum systems in contact with fermionic macroscopic baths, e.g. in quantum dots or single-molecule junctions [111,112].In essence, it neglects the structure of the density of states in the bath and, as a consequence, a main limitation is that it fails to describe the shorttime dynamics [112].Nevertheless, this problem does not affect this study since we are interested in large times.We should further note that the energy of the system-bath interaction is proportional to Λ, and therefore is divergent in this limit.We will therefore take Λ to be finite but much larger than any other energy scale of the system.For our analysis to be valid we simply require dynamics much slower than Λ −1 [76].
The dynamics are solved via a quantum Langevin approach, which is detailed in App.A (see also [62,83]).Taking the initial state to have no correlations between S and B (ρ(0) = ρ S (0) ⊗ ρ B (0)) and the bath to be in a thermal state; we find the probability of occupation of the excited 1 The chemical ν potential of the bath is incorporated by subtracting νâ † â to the system's Hamiltonian.Since here HS = εa † â (with ε the energy of the system), we can simply redefine ε to be the difference between the system's energy and the chemical potential and set ν = 0 without loss of generality.level of the system p(t) = ⟨â † â⟩ and the systembath interaction energy v(t) = k λ k ⟨â † bk ⟩+h.c., where f β (ω) = (1 + e βω ) −1 is the Fermi-Dirac distribution and we defined the propagator with µ(t) := 1 2 g(t) 2 .From these expressions we can exactly compute the thermodynamic work eq.( 4), which reads: From the exact solutions for p(t) and v(t), in App.A we show that Requirement 1 is satisfied, and hence W = ∆F in the quasistatic limit.For slow but finite-time processes, we perform a slow driving expansion of eq. ( 14) and eq.( 15) (details in App.B) using that the thermalization rate of the system is Γ := 2 ℏτ τ 0 dt µ(t), so that the expansion can be performed in orders of 1/(τ Γ).We then obtain an expansion for W analogous to eq. ( 6) where the entropy production Σ is described by eq. ( 9) with ⃗ λ(t) = (ε(t), µ(t)) and the thermodynamic metric where This metric gives a geometrical description of slow thermodynamic protocols performed on the system.By solving the geodesic equations [113], we can find the geodesic length L and hence the minimal entropy production eq.( 11).

Special limits of the metric
Before attempting to solve the geodesic equations for the case of erasure, we now study the high  [109] with boundary conditions ε(0) = 0 and βε(τ ) = 100 to the lower bounds given by eq. ( 1), Van Vu et al. [52] and Zhen et al. [51].
and low temperature limits, as well as the limit of weak coupling, to gain further analytical insights on the form of optimal protocols and the associated entropy production.

High temperature limit (βε, βµ ≪ 1)
Since the terms of eq. ( 18) quickly decay at high frequencies we can perform the high temperature expansion f β (ω) = 1  2 − 1 4 βω + O(β3 ω 3 ) directly in the metric.At leading order, we find: This enables an analytical solution of the geodesic equations.
2 Usually, the initial condition for erasure would be ε(0) = ν for ν the chemical potential of the bath and ε the energy of the two-level system (so that the corresponding thermal state is the fully mixed state).But since here we defined ε to be the difference to the chemical potential we take ε(0) = 0 without loss of generality.

Zero temperature limit (βε or βµ → ∞)
In the limit of T = 0 3 we have , where Θ is the Heaviside step function.Therefore the metric becomes (cf.App.D) which coincides with the metric of an angle distance in the (ε, µ) space -hence the metric is singular.If we re-parameterize (ε, µ) as (r cos ϕ, r sin ϕ) we find 2 .Therefore any protocol that keeps φ(t) constant is a geodesic, leading to the minimal entropy production: with ϕ = arctan(µ/ε).Note that there are multiple (infinitely many) geodesics for any pair of boundary points.This fact prevents us from continuing the expansion to further orders in temperature.Nevertheless, this limit provides analytical insights on optimal protocols with βε or βµ ≫ 1.
In particular, we note that there is no need for a diverging coupling even when ε(τ ) → ∞ as, once µ has become large, eq.(25) shows that it is optimal to reduce the coupling while increasing the energy.Furthermore, eq. ( 25) shows that at zero temperature, while the reversible cost of the operation goes to zero, the dissipation remains strictly positive.This result is complementary to the findings of Ref. [13] which demonstrate a finitesize correction to Landauer's bound that does not disappear in the zero-temperature regime.

Weak coupling limit
Lastly, we take take the weak coupling limit to compare to previous erasure results that are obtained via Lindbladian dynamics, a common assumption in previous works on optimal thermodynamic control in the quantum regime [35,[49][50][51][52][53][54].In this limit the coupling is taken to be small and constant, therefore the metric becomes a scalar (cf.App.B.2): which matches with the results of [50,54] which were also obtained with thermodynamic geometry.Indeed, this metric can be obtained from the rate equation which can also be obtained by taking the weak coupling limit in the Heisenberg equations that define eq. ( 14).In this regime protocols that minimize dissipated heat at arbitrary speed were found by [109].Therefore we will compare the results one obtains in slow driving and the results of [51,52] to the exact minimization of [109].
We are interested in erasure processes, where ε(t) is driven from 2 ε(0) = 0 to ε(τ ) = ε * with ε * ≫ k B T in a time τ .Optimal finite-time protocols are those which minimize the work cost W = τ 0 dt ε(t)p(t), and hence the heat dissipated to the environment Q = W − ∆E.The results of [109] provide an exact solution to this problem, which is shown in Fig. 1.As it is well-known in finite time stochastic thermodynamics [108], jumps appear in the optimal solution.However, as we approach the quasistatic limit where τ Γ ≫ 1, the jumps progressively disappear.In App.E we prove why the jumps should also disappear in the long times limit at strong coupling.As detailed in App.B.2, and also discussed in previous references [50], the optimal driving solution in this limit has the simple analytical form leading to the work cost from where we can directly recover eq. ( 1) through the first law of thermodynamics (note that ∆E S ≈ 0).In Fig. 1 we notice that the exact solution of [109] agrees well with this analytical form in the slow driving limit.For completeness, we also show recent results of [51,52].These results apply more generally to any Markovian master equation (here we apply them to the particular case of eq. ( 27)), and one can see that they provide a bound to the exact numerical (and approximate analytical) solutions.

Optimized erasure
We now focus on erasure outside of any approximation, where we will optimize the driving over both the energy and coupling.In what follows, we focus on minimising Σ in an erasure process, which imposes specific boundary conditions to the geodesic equations.We assume that we have no prior knowledge of the system, therefore its initial state is ρ S (0) = 1/2.This translates in taking ε(0) = 0 so that it coincides with the thermal state of H S .For the qubit to be erased we want its final state to be ρ S (τ ) ≈ |0⟩⟨0| (i.e.p(τ ) ≈ 0).Since the driving is done slowly, p(t) is always close to its thermal expectation value.Therefore by choosing βε(τ ) → ∞ we ensure p(τ ) ≈ 04 .For the coupling, the boundary conditions are µ(0) = µ(τ ) = 0, because we want to think of this as an "erasure machine" that the qubit is "brought to" at the start and "retrieved from" at the end.Given this family of protocols, we recognise from eq. ( 6) that W = k B T (ln 2 + Σ), and similarly Q = k B T (ln 2 + Σ), thus justifying the minimisation of Σ as given in eq. ( 9).After the qubit has been decoupled (i.e. at t > τ ), we bring the Hamiltonian of the system back to its starting value (ε = 0) to close the cycle.Since p(τ ) ≈ 0, (left) A series of optimal protocols depicted in the (µ, ε) space.They all start with zero energy and coupling and end with finite energy and zero coupling.In the limit of large βε(τ ) they can be considered as erasure protocols.
(middle) The entropy production of the optimal erasure protocol as a function of the final energy, compared to the high temperature regime cost eq.( 23).(right) Comparison of the entropy production for a geodesic protocol in which one parameter is varied at a time (with µ being increased until µ * ) and the weak coupling approximation eq. ( 1); the minimal possible entropy production, τ Σ min /β = 2.57946 ± 1 • 10 −5 [ℏ], obtained when both parameters are changed simultaneously is also shown.
this step requires no work, and it can be done arbitrarily quickly.
The geodesic equations we obtain for this process are not solvable analytically.The integral of eq. ( 18) can be solved to give us an expression of the metric in terms of polygamma functions (cf.App.B) but it does not simplify the geodesic equations into an analytically solvable form.We therefore turn to numerical tools to obtain the optimal protocol and compute the the dissipated work.Though, in our case, we want to impose the aforementioned boundary conditions; this is known as a Boundary Value Problem (BVP), which is famously hard to solve numerically [114].Though we can use the fact that the high temperature limit is accurate at the start of an erasure protocol, therefore the initial conditions of the optimal protocol for erasure will match the initial conditions of eq. ( 21) and eq.( 22).This allows us to turn the BVP into an Initial Value Problem (cf.App.G) which is much simpler to solve.
In Fig. 2 we show optimal erasure protocols in the (µ, ε) space for different final values of βε.We can notice that the predictions of the high and low temperature limit are verified: at the start of the protocols the coupling is increased but once we reach a certain value there is no more need to increase it, regardless of the final value of βε we try to reach.Interestingly, the maximal value reached by βµ is larger than 1.This shows that reaching the strong coupling is needed to achieve optimal erasure, which is one of the main insights of our work.In the same figure we also show the value of τ Σ min /β for the optimal protocol as a function of the final energy.
We can see that for small values of βε(τ ) eq. ( 23) gives an accurate description of the work cost, but as we reach higher values it saturates around . This provides a finite-time correction to Landauer's principle in this set-up, thus leading to a generalisation of eq. ( 1): with a ≈ 2.57946 and τ Pl = βℏ.This is one of the main results of this work and can be seen as a generalization of eq. ( 1).As opposed to the results of [51] and [52], eq. ( 30) is only valid for large protocol times; yet, it has the advantage of taking into account strong coupling effects (including any possible variation of the coupling strength), having a much simpler form for the correction (which is independent of any chosen relaxation timescale), and we provide an explicit protocol to achieve it.By turning around eq. ( 30) one can highlight a quantum speed limit for erasure of a qubit, furthermore this speed limit is of the order of the Planckian time τ Pl = ℏ/k B T which is conjectured to be the fastest relaxation timescale for thermalization [84].In particular, one can see that eq. ( 30) bounds the speed of erasure by the order of τ Pl regardless of how large is the coupling strength used in the protocol.
Interestingly, we now argue that the form of the correction eq. ( 30) is in fact general of any Landauer erasure protocol with control on S and the SB coupling.Indeed, first note that the geodesic length L is dimensionless and can only depend on β and the boundary conditions as we optimize over µ and ε.In an erasure process, the boundary conditions read: ε(0) = 0, ε(τ ) → ∞, and µ i (0) = µ i (τ ) = 0 where i runs over all the possible control parameters on SB.But this implies that L is independent of them and hence of β.Therefore k B T Σ min will take the form of a constant, independent of any parameter of the system and bath, divided by τ .This is a crucial difference from eq. ( 1).This simple argument based on dimensional analysis thus shows that eq. ( 30) is rather general, with the value of a depending on the specific implementation (e.g. the ohmicity of the bath).It is important to highlight that the bound eq. ( 30) implies that, even when having access to arbitrary strong SB interactions (naively taking Γ → ∞ in eq. ( 1)), infinite time is still required for perfect erasure due to the quantum-mechanical correction derived here.
Finally, we analyze a scenario where the coupling is kept constant while ε(t) is driven, which is motivated both by experimental set-ups and for a comparison with the weakly interacting case.Therefore, we restrict to one-parameter protocols consisting of the three following steps: 1. while keeping ε at 0 we turn on the coupling to some value µ * ; 2. while keeping the coupling fixed we bring ε from 0 to some value ε * ≫ k B T ; 3. while keeping ε constant we turn off the coupling.Each step contributes positively to the entropy production, and its minimisation is discussed in App.F. In Fig. 2, we show Σ min for different values of µ * , ranging from the weak to the super-strong coupling regime.It can be appreciated how eq. ( 1) breaks down, and also how such one-parameter protocols become close to the fundamental limit eq. ( 30) for βµ * > 1.

Conclusions and outlook
Deriving finite-time corrections to the seminal Landauer bound is a challenging endeavour in stochastic and quantum thermodynamics [35,[44][45][46][47][48][49][50][51][52][53][54].Previous works have focused on markovian systems only, which in the quantum regime is obtained through the weak coupling limit (βg 2 → 0).However, should a general finite-time correction exist, it will require the presence of strong coupling at some point during the process as the dissipation generated in finite time is propor-tional to g −2 when g is small 5 .Motivated by this observation, we have developed new insights into the form of optimal protocols for erasure beyond the weak coupling limit.
We have focused on a bit encoded in the occupation of a single fermionic mode, which can be strongly coupled to a reservoir.We have derived analytically the thermodynamic metric, which governs the dissipation rate in the slow driving regime, and showed that it takes a simple form in the high and low temperature limits.From the general form of the metric we obtained the optimal erasure protocol, which requires increasing the coupling strength to g 2 ∼ k B T , which corresponds to a relaxation timescale of the order of the Planckian time τ Pl .The corresponding dissipation yields a finite-time correction to Landauer's bound for this setup, which is substantially lower than similar results in the weak coupling regime.Furthermore, by using the obtained bound as a quantum speed limit, this result adds further evidence to the conjecture [84] that τ Pl is fastest relaxation timescale many-body systems can achieve.
While our results were derived in a fermionic model there are some general insights that follow from our work.First there is a fundamental quantum correction that prevails, see eq. ( 30), which can be compared with eq. ( 1) derived in the weak coupling regime.While the specific value of a in eq. ( 30) will depend on the specific setup, it will never approach 0 (even for diverging system-bath coupling) due to the inherent cost of changing the interaction strength.Furthermore, to obtain these results we adapted the framework of thermodynamic geometry to system-bath unitary dynamics in which the coupling can be arbitrarily large or small.This is in contrast to recent claims of failure of this approach in closed quantum systems [115].Finally, as was argued before, our results make evident the need of strong coupling in a general finite-time correction to Landauer's principle.
This work opens exciting directions for the future.On the one hand, the level of experimental control required to implement such protocols is in principle possible in quantum dots [116][117][118][119], where the energy-level ε(t) and coupling g(t) can be independently controlled, even by several orders of magnitude [110].On the other hand, it would be interesting to characterise the dependence of a in the nature of the bath and the SB coupling (e.g. its spectral density), more generally to derive similar quantum-mechanical finitetime corrections that are independent of the specific implementation, and to gain further insights in the connection between Landauer erasure and the Planckian time.We consider a single fermionic mode coupled to a fermionic bath.Without loss of generality we can set the chemical potential to 0 and the ground state of the two-level system to 0. The Hamiltonians of the system and bath are where â † is the creation operator of the two-level system and b † k is the creation operator of a bath mode with frequency ω k .These ladder operators follow the canonical anticommutation relations.For the interaction between system and bath, the Hamiltonian is where the λ k are the interaction weights.
Will will consider that ε(t) and g(t) are the control parameters to then perform the erasure of information in the single mode of the system.We now proceed to solve the dynamics of the system and bath in the Heisenberg picture.For an operator Â in the Schrödinger picture, we denote by ÂH (t) the corresponding operator in the the Heisenberg picture.The evolution of ÂH (t) is defined by the Heisenberg equation of motion: where Ĥ(t) = ĤS (t) + Ĥint (t) + ĤB and ℏ = 1.Applying this equation to the ladder operators âH (t) and bk,H (t) we find the following system of n + 1 equations: By defining ûk (t) = e iω k t bk,H (t) we can see that eq. ( 36) becomes which is solved by ûk (t) = ûk (0) − iλ * k t 0 ds e iω k s g(s)â H (s). We therefore find bk, where we used that ûk (0) = bk,H (0) = bk .Therefore where we defined the noise operator ξ(t) = k e −iω k t λ k bk and χ(t) = k e iω k t |λ k | 2 .We can notice that χ(t) is the (symmetrized) noise correlation function: and its Fourier transform is the (unit-less) spectral density of the bath By inserting eq. ( 38) into equation eq. ( 35) we find an equation of motion for âH (t): In order to solve eq. ( 39) we will need to explicitly take the continuum limit so that our bath indeed becomes a bath.We can take its spectral density to be either a Lorentzian J(ω) = Λ 2 Λ 2 +ω 2 or a pass-band J(ω) = Θ(Λ − |ω|) (Θ is the Heaviside step function).We will need to assume that we are working in the wide-band approximation (Λ → ∞).More practically, we are assuming that the bath interaction is the same over the energies we are spanning with the system.This limit allows us to say that the noise correlation function is negligible for time differences larger than zero: In this limit eq. ( 39) becomes considerably simpler: Similarly to how we solved eq. ( 36), we define û(t) = exp t 0 z(s)ds âH (t) for z(t) := iε(t) + 1 2 g(t) 2 .Now we have which is solved by û(t) = û(0) − i t 0 ds g(s) exp [ s 0 z(r)dr] ξ(s).We therefore find the solution of the evolution of the ladder operator of the distinguished mode: where we defined the propagator G(t, s) = exp − t s z(r)dr .And from eq. ( 37) we find the solution for the bath modes:

A.2 Relevant observables
Since we are performing an erasure, we will assume the system starts in a factorized state and that the bath starts in a thermal state at inverse temperature β state with respect to its Hamiltonian: We are interested in computing the occupation probability of the excited level p(t) = ⟨â † â⟩ and the system-bath interaction potential v(t) = ⟨ V ⟩.From eq. ( 41) we have where we used the CAR to get Tr[â † ξ(s)] = 0 and drop the cross terms.Further using the CAR we simplify the remaining trace in the integral: where f β (ω) = (1 + e βω ) −1 is the Fermi-Dirac distribution.We can apply the continuum limit to eq. ( 44) by using the equality , which holds for any function h by definition of J.We can then apply the wideband limit by using lim Λ→∞ lim n→∞ J(ω) = 1.We find Applying eq. ( 45) to eq. ( 43) we get Using eq. ( 38), the CAR and eq.( 41), we have

So we have
A.3 Proof of requirement 1 We will now proceed to prove that, in absence of driving, p(t) and v(t) thermalize.We do so in two steps, we first simplify the expressions of eq. ( 46) and eq.( 47) for ε(t) = ε and g(t) = g and compute the infinite time limit.Then we compute the thermal expectation value of the corresponding observables and prove that the obtained expressions are the same.

A.3.1 Infinite time limit in absence of driving
By assuming that the driving parameters are kept constant the propagator becomes This allows us to compute the time integrals in eq. ( 46) and eq.( 47): As a side-note, it is interesting to note that the frequency integral of eq. ( 49) can be solved to give the following expression for the occupation probability where ψ (0) (z) = d dz ln Γ(z) is the digamma function (defined as the logarithmic derivative of the Gamma function) and B(x; a, b) = x 0 ds s a−1 (1 − s) b−1 is the incomplete beta function.This expression is useful for numerical implementations as it is faster to compute than the integral of eq.(49).
By taking the limit t → ∞ in eq. ( 49) and eq.( 50) we find Here we can notice that if we take the Laplace transform of the propagator we obtain which allows us to rewrite eq. ( 52) and eq.( 53) as

A.3.2 Thermal expectation value
We now compute the expectation value of â † â and V when the state is a Gibbs state Therefore we want to find p th := Tr[ω β â † â] and v th := Tr[ω β V ].Using the fact that the total Hamiltonian is quadratic, we can diagonalize it to rewrite it in the following way where ε k are eigen-energies and ĉk are fermionic ladder operators that follow the CAR: {ĉ † j , ĉk } = δ jk 1, {ĉ j , ĉk } = 0.They are related to the original ones in the following way where |k⟩ = ĉ † k |0⟩ are 1-particle eigenstates of the Hamiltonian with eigenvalue ε k .Inserting this in the expression for the thermal expectation of the probability of occupation we find From eq. ( 41) it is easy to see that we can write the propagator in the following way where we used the fact that the vacuum state does not evolve Û (t) |0⟩ = |0⟩ and that since we are performing no driving we have Û (t) = e −it k ε k ĉ † k ĉk .Note that the sum needs only to be over 1particle states as there is a scalar product with the 1-particle state â † |0⟩.By now defining φ(ω) := Considering eq. ( 55) it is clear that if φ(ω) = 1 π ℜ G(−iω) then we have proven p th = lim t→∞ p(t).Therefore we compute the Laplace transform of G(t, 0) using eq.( 62) where P. denotes the Cauchy principal value, Θ(t) is the Heaviside step function and we used that its Fourier transform is (in a distributional sense) dt e ist Θ(t) = πδ(s) + P. i s .Since φ(ω) is by definition a real function we can see that P. ∞ −∞ dω ′ φ(ω ′ ) ω−ω ′ is a real number.Therefore we can conclude φ(ω) = 1 π ℜ G(−iω) .Which concludes the proof of the thermalization of p(t).To prove the thermalization of v(t) we proceed in a similar fashion.We start by computing v th To proceed we have to define the following cross-propagators where we defined ψ j (ω) = k ⟨k| b † j |0⟩⟨0|â|k⟩ δ(ω − ε k ).By further defining ψ 0 (ω) := k λ * k ψ k (ω) and ψ(ω) = ψ 0 (ω) + ψ * 0 (ω) we can see that Therefore, similarly to the case of G(t, 0), we have and in particular πψ(ω) = ℜ K(−iω) (since ψ(ω) is real by definition).Hence, by eq. ( 56) and eq.( 68), the last step to prove that v(t) thermalizes is to check that ℜ K(−iω) = gℑ G(−iω) .To do so we start by computing the components of K(t): from eq. ( 41) we can see that and from eq. ( 42) Since the time time dependence is contained in the exponentials, it is straightforward to compute the Laplace transform Therefore we find which allows us to conclude ψ(ω) = gℑ G(−iω) .This concludes the proof of the thermalization of v(t).

B Slow driving expansion B.1 Deriving the thermodynamic metric
We are interested in performing an erasure protocol and minimizing the work cost of performing it.
The erasure protocol is one where ε(0) = 0, ε(τ ) ≫ β −1 and g(0) = g(τ ) = 0, for τ the total time of the protocol.The work cost of a protocol where we control ε and g is To get a correction to Landauer's bound for finite time protocols, and to work with more tractable expressions, we expand eq. ( 73) in the long times limit up to first order.To do that we first need to make some notation changes.First we make the time parameter in ε and g dimensionless, so that the protocol starts at "time" input parameter 0 and ends at "time" input parameter 1.So we have the following mappings: t → τ t, dt → τ dt and d dt → τ −1 d dt .Second we need to "extract" the evolution timescale of the system in order to make the slow driving expansion.From eq. ( 46), eq. ( 47) and the definition of the propagator (or more clearly form eq. ( 49) and eq.( 50)) it is quite clear that the relaxation timescale of the system, at any point of the evolution, is of the order (g(t) 2 ) −1 .Hence we are going take the average of the square of the coupling as normalizing factor, we therefore define (in normalized time) Γ := 1 0 dt g(t) 2 .We now define a normalized version of our control parameters: We can therefore write the expression for work cost in this new convention We can also rewrite the propagator and the expectation values of the observables where we were able to insert a phase in the time integral of eq. ( 77) because of the absolute value and rescaled ω by Γ.We can see that we have to expand in 1/τ Γ the same integral for both eq.( 77) and eq.( 78).To do that we do partial integration.First, we can notice that Therefore we can write where we can evaluate the first part as γ(s) ).We absorbed the G ω (t, 0) term in O(e −τ Γ ).For the second term we evaluate the derivative and continue integrating by parts Similarly as in eq. ( 79), we keep only the evaluation at s = t for the first term because the evaluation at s = 0 is of order O(e −τ Γ ).Whereas the remaining integral will also have to be evaluated by parts, and in doing so we will obtain another power of 1/τ Γ.But since we are only interested in the first order correction and the integral will only yield terms of order O(1/τ 2 Γ 2 ) we don't need to compute it.By combining eq. ( 79) and eq.( 80) we finally find The absolute value squared of eq. ( 81) is Combining this with eq. ( 77) we find the expansion of p(t) in the slow driving regime Whereas the imaginary part of eq. ( 81) is Therefore the slow driving expansion of v(t) is Therefore we can see that we can rewrite the work cost of the protocol as The leading order term is where we re-scaled ω by Γ.Here we were able to perform the integral independently of the function describing the control parameters.Therefore W (0) only depends on their initial and final value.More importantly, we can identify the instantaneous thermal expectation values of p(t) and v(t) (from eq. ( 52) and eq.( 53)) in the time integral of W (0) .This implies very directly that W (0) = ∆F .
We can notice that we can write W (1) as with ⃗ λ t = (ϵ(t), γ(t)) T and the metric for Since the leading order is independent of the path taken in parameter space, minimizing the work cost of erasure only implies minimizing W (1) , i.e., the entropy production k B T Σ.As we see from eq. ( 88) it is equivalent to finding the shortest path in a metric space described by the metric m( ⃗ λ).The length of this shortest path is known as thermodynamic length.In the main text eq. ( 18) and eq.( 19) represent the metric when the problem is rewritten in terms of the unit-full parameters.

B.2 Weak coupling limit
Previous works on optimization of finite-time Landauer erasure have focused on the Markovian regime [44][45][46][47][48][49][50][51][52][53][54], corresponding to the weak coupling limit.We analyze this regime in this section.First we assume that the coupling remains unchanged during the protocol, which means Γ = g 2 and γ = 1.We start by rewriting p(t) from eq. ( 83) in a more convenient manner under this first assumption: where we used eq.( 74) to go back to unit-full parameters and re-scaled ω by Γ. Integrating by parts the second integral we get 3 and dropped the O(1/τ 2 Γ 2 ) to make the notation lighter.We now take the weak coupling limit, but we do so while keeping the slow driving assumption: τ Γ ≫ 1.Using the results of [120] we get the following where the equalities are meant in a distributional sense.Therefore we find that the occupation probability in the weak coupling limit is This result coincides with applying a slow driving expansion to a simple exponential relaxation model with characteristic time Computing the work cost yields with ∆F = β −1 ln 1+e −βε(0) 1+e −βε (1) .We will now minimize the work cost of the erasure protocol, similar optimizations have been done before in [50,109,121].From variational calculus we know that the extremal function of the integral in eq. ( 94) will keep the integrand constant.So we can solve the variational problem as follows: with We therefore find and recover the result of eq. ( 1) from the main text:

B.3 Solving the integral of the thermodynamic metric and finding the symmetry
We will now try to find a more tractable version of the metric in eq. ( 89).First we notice that We can remark that m 0 coincides with a metric of an angle distance in the (ϵ, γ) space.To solve the integral of eq. ( 89) we will go in Fourier space.For a function h(ϵ) its Fourier transform h(ξ) has the defining property Even though f β (ω) is not an integrable function we can find its Fourier transform in a distributional sense For m 0 we find with m0 (0, γ) = π 2γ 1. Therefore we can rewrite eq. ( 89) as where we used eq.( 99), eq. ( 100) and the fact that ∞ −∞ dω e iωx = 2πδ(x).If we now insert eq. ( 101) and flip the sign in the second integral we find We will now be able to compute these integrals in terms of poly-gamma functions.The poly-gamma function of order m ≥ 0 is defined as ψ (m) (z) := d m+1 dz m+1 ln Γ(z).For m > 0 and ℜ[z] > 0 they have an integral representation: Using the fact that 1 sinh(x) = 2e −x 1−e −2x , a change of variable and eq.( 104) we find We can notice that the metric explicitly depends on Γ in such a way that it seems that the solution for the geodesic should depend on this scale factor.Though this dependence disappears if we reparameterize the problem in terms of its original unit-full parameters.We start by rewriting the work as W = ∆F + W (1) + O( where we redefined W (1) with the unit-full parameters ⃗ λ t = (ε(t), µ(t)) T (with µ(t) := 1 2 g(t) 2 = Γγ(t)): We remind the reader that z = µ + iε.This metric is the same as the one presented in the main text.
Despite not looking very approachable, eq. ( 108) is a much more tractable version of eq. ( 18) from the main text when it comes to numerical implementations (as the polygamma functions are computed much faster than integrals) and analytical studies of the geometric properties of thermodynamic protocols.
We can notice from eq. ( 107) and eq.( 108) is that there is a symmetry in the corrective term.If we perform the following transformation: for λ > 0; then W (1) remains unchanged.This symmetry allows us to conclude that the minimal value of W (1) to perform Landauer erasure will be of the form c/τ where c is a constant that does not depend on any physical quantity.

C High temperature limit
In order to find the minimal dissipation in the multi-variable case we need to numerically solve the equations of motion given by the exact metric, which (unsurprisingly) are very untractable analytically.But instead of solving an initial value problem (which we can always solve by numerical integration, in principle) we are trying to solve a boundary value problem.Generically, to solve a BVP numerically, the solver will try many IVPs until the wanted BVP is reached.But here we can notice that we can turn the BVP into an IVP by taking an analytical approximation of the problem around the point (ε, µ) = (0, 0).
As we have seen in eq. ( 109) there is an underlying symmetry in this problem, so a limit where ε and µ are infinitesimal is the same as a limit where β is infinitesimal but ε and µ finite.Formally we are requiring β|µ + iε| ≪ 1, which is a high-temperature limit.It is important to note that despite the fact that this approximation will yield some analytical results on how to optimize a protocol in the high temperature regime it will not give us a result that is relevant for Landauer erasure because to perform erasure we are assuming that we reach βε ≫ 1, which is a low temperature limit.
One way to obtain an analytical result in this framework is by going back to eq. ( 89) and apply the high-temperature expansion of the Fermi-Dirac distribution: f β (ω) = 1 2 − 1 4 βω + O(β 3 ω 3 ).Since m 0 (±∞, γ) = 0 the first term of the metric in this expansion is 0. But from the next order we find (in unit-full parameters) It might not be immediate why we also require βµ ≪ 1, but it becomes clear that it is required when we want to obtain the same result by applying the same expansion on eq. ( 108) (which would also allow us to get further orders).We will now compute and solve the equations of motion: where we assumed the Einstein tensorial notation and Γ i jk are the Christoffel symbols Here we have m il = 8µ β δ il and ∂ a m bc = − β 8µ 2 δ aµ δ bc .We therefore find which can be rewritten as We get the following differential equations for µ and ε: From the first equation we can see that d ε/ ε = dµ/µ, therefore ε = Cµ for some constant C (we can already see as a sanity check that ε never changes sign in an optimal protocol).The equation for when we consider that µ = g 2 /2 and μ = ġ2 + gg we can see that Taking into account that the boundary conditions for g are g(0) = g(1) = 0 we find that g(t) = A sin(kπt) for some constant A, k ∈ N * and C = 2kπ.By choosing ε(0) = 0 and ε(1) = ε * we have ε(t) = kπA 2 t 0 ds sin(kπs) 2 , therefore A 2 = 2ε * kπ .Thus the optimal protocol, portrayed in Fig. 3, is From eq. ( 116) we can compute the integrand of the dissipated work: We therefore find the dissipated work for the high temperature limit by inserting this in eq. ( 107) and taking k = 1: We can see that in this scenario the corrective term grows extensively with the final energy ε * , combining this with the fact that the exact metric goes to 0 faster than O(|z −1 |) we can presume that most of the dissipation in the exact protocol is caused by the part of the protocol that matches with the high temperature regime.

D Low temperature limit
We will now study the limit of T → 0 we have f β (ω) → f ∞ (ω) = Θ(−ω), where Θ is the Heaviside step function.Therefore the integral of eq. ( 89) becomes (in unit-full parameters) where we used eq.( 98) and the fact that m 0 (+∞, µ) = 0. Thus we have We can compute the integrand of W (1) to find By defining the coordinates ⃗ λ r,ϕ = (r, ϕ) T , such that ε = r cos ϕ and µ = r sin ϕ, we have ϕ From which we can deduce the metric in these new coordinates m (r,ϕ) Crucially, we can notice that this metric is singular: here changes in the coordinate r do not cause an increase in the work cost.Therefore we can parameterize geodesics at zero temperature as follows where r(t) is any function that satisfies the boundary conditions.By using eq.( 107) we can find the dissipated work for ∆ϕ = ϕ(1) − ϕ(0).

E Discontinuities in the protocol
As is mentioned in the main text, it is well known [108] that, at weak coupling, discontinuities appear at the beginning and at the end in optimal finite-time protocols.In the main text we gave numerical evidence to the fact that these jumps disappear as one approaches the quasistatic limit.Here we make an a posteriori argument as to why these jumps should also disappear for the system we studied in the strong coupling regime.
We can start by immediately discarding jumps in the coupling because these lead to a diverging work cost in the wideband limit.Then, if one makes a jump in the energy at the start of the protocol from 0 to ε * its work cost is where p(0) is the probability of occupation at t = 0 and is set to be 1/2.We can compare it to the work cost given by an optimal continuous protocol with the same boundary conditions: where a(βε * ) is a bounded function of βε * (c.f.figure in the main text).We note that by expanding ∆F around β = 0 we find Therefore at β > 0 and small enough one finds that ∆F < W jump .But one can also note that lim β→∞ ∆F = 0, and since ∆F is a monotonous function of β we conclude that ∆F < W jump for any temperature and any ε * > 0. At this point it is immediate that there exists τ large enough such that W cont < W jump .We finally consider jumps in the energy at the end of the protocol.These jumps would have to be performed once the coupling is very close to zero since it is at the end of the protocol.In this limit optimal protocols are found by [109].These optimal protocols feature jumps whose magnitude is controlled by a constant of integration K, in particular the magnitude of these jumps is O( √ K).This constant is defined as follows (in terms of our notation with unit-less time) From which it is clear that K goes to 0 in the limit of τ Γ ≫ 1 and therefore the jumps disappear.This behavior is also confirmed in Fig. 1 where these optimal protocols are shown.

F One-parameter case
Because of the metric we obtain in eq. ( 108) it is quite clear that, in the two-parameter case, we will not be able to solve analytically the geodesics for the full problem.This even prevents us form finding an analytical expression for the distance between two points in the parameter space, as it is the length of the shortest path (for which we have no expression).But by fixing one parameter to an arbitrary value and solving for the other we can use the fact that geodesics always have a conserved quantity along their path (the integrand: ⃗ λ T t m( ⃗ λ t ) ⃗ λ t ) to avoid solving the geodesic equation and finding an explicit formula for the length and geodesic.We point out the fact that the symmetry mentioned in App.B does not lead to a conserved quantity because β is a constant of the system instead of a function of time for which we are solving.
Here we will take the erasure protocol to be made of three parts, which will be optimized separately: 1. we turn on the coupling to some value µ * while keeping the energy at zero; 2. while keeping the coupling at µ * we increase the energy from zero to infinity; 3. we turn the coupling off.Incidentally, this type of protocols are more realistic for an experimental realization as often setups are not able to control optimally energy an coupling at the same time.And even if the control over the coupling is only to turn it on to some value and turn it off, step 2 will remain valid.Furthermore, previous studies done at weak coupling essentially are described by this type of protocol; but they are in a regime where step 1 and 3 can be neglected.Therefore we can compare the results of this section to those of the weak coupling limit.
We start by looking at step 3, in the limit of βε → ∞ we actually reach a scenario described by the T = 0 limit.Therefore the length of this step is described by eq. ( 123), for any finite value of µ * the angle span of this step is trivially 0. Therefore, up to first order, this step will not cause any extra dissipation, no matter how it is realized.
We notice that we can write the length of the first step as follows where we used the fact that, since the metric is not explicitly time-dependent, the sign of μ has to be always positive for this step.With eq. ( 108) we find the following expression for the length Now that we have an expression for L 1 we can recover an equation for µ(t).We can use the fact that the integrand of the time integral in eq. ( 130) is constant to obtain which gives an implicit definition of µ(t), or rather an explicit definition of its inverse t(µ).
By following the same procedure as in step 1 we can recover the length of step 2 These implicit definitions of ε(t) and µ(t) can be solved numerically, the results are shown in Fig. 4. The question remains about how to subdivide optimally the protocol times of step 1 (τ 1 ) and step 2 (τ 2 = τ − τ 1 ).The total excess work is given by W (1) = L 2 1 /τ 1 + L 2 2 /τ 2 , by taking the derivative and imposing it to be zero we find And therefore we find which is indeed what was to be expected, as L 1 + L 2 is the total length of the protocol.
We now discuss how we can obtain an exact version of eq. ( 1) of the main text, so that it applies also in the strong coupling regime.First we can notice that since L 1 > 0 and L 2 > 0 we have W (1) ≥ L 2 2 /τ .Then by considering that m εε (ε, µ * ) is a one-dimensional metric it has to be positive by definition.Therefore the integrand of eq. ( 133 and since the approximation becomes exact in the limit µ * → ∞ we have lim µ * →∞ L 2 = √ π/2.Therefore, In Fig. 5 we show the minimal value of W (1) for step 2 as a function of βµ * , and we compare it to eq. ( 1) of the main text and to eq. ( 140).We can see how, when the coupling becomes small, the exact curve agrees with eq.(1).133)) with the weak coupling approximation (eq.( 1) of the main text) and the lower bound of eq. ( 140).
G Numerical solution to the general case Figure 6: A series of optimal protocols depicted for multiple values of βε (1).They all start with zero energy and coupling and end with finite energy and zero coupling.In the limit of large βε(1) they can be considered as erasure protocols.Shown in the parameter space (left) and as a function of time (centre and right).
We now discuss how the numerical problem of finding the optimal erasure protocol was approached.
Having found the metric eq.( 108), all we have to do to find the optimal erasure protocol is to solve the geodesic equations with the Christoffel symbols defined as in eq.(112).Though the differential equations we get are quite untractable and cannot be solved analytically, we won't even write them here as they are very long and will not bring any insight.Therefore we will solve them numerically, and indeed eq. ( 141) is quite practical for numerical integration since the second derivative of the parameters can be easily isolated.
But at this point we can notice that close to t = 0, by continuity, we are satisfying the conditions for the high temperature approximation.And at t = 0 the approximation becomes exact.Therefore the initial conditions of an optimal erasure protocol in the general case must match with the initial conditions of the protocols that we previously studied in the high-temperature regime.From a numerical perspective it is much more preferable to solve an initial value problem instead of a boundary value problem.Therefore we used the numerical solver DOP853, implemented in the scipy library in python, to solve eq. ( 141) with the initial conditions given by eq.(116).
To be precise, we cannot start the integration from t = 0 as the metric is formally divergent at (ε, µ) = (0, 0), therefore we evaluate eq. ( 116) at an infinitesimal time and integrate from there.The specific value we choose for ε * sets the value of ε(1) that is reached in a monotonous way.When ε * is chosen small the protocol closes matches with those of eq. ( 116) (as long as ε( 1) is also small).Then for larger values of ε * we get more interesting behavior, as is shown in Fig. 6.
When one changes the value of k in eq. ( 116) the value of ε( 1) that is reached is different.But, as is shown in Fig. 7, by numerically searching values of ε * such that the same ε( 1) is reached for different values of k we find that the protocols end up being the same.
Finally we thought it might be interesting to compare the best one-parameter protocol to the geodesic erasure protocol we find numerically.And we can see from Fig. 8 that, despite seeming very different in the path taken in the parameter space, when we look at the functions of time they are actually quite similar.The apparent difference happens because the part of the protocol for βε ≫ 1 is done very quickly since the metric is vanishing in that region.

Figure 2 :
Figure2: (left) A series of optimal protocols depicted in the (µ, ε) space.They all start with zero energy and coupling and end with finite energy and zero coupling.In the limit of large βε(τ ) they can be considered as erasure protocols.(middle) The entropy production of the optimal erasure protocol as a function of the final energy, compared to the high temperature regime cost eq.(23).(right) Comparison of the entropy production for a geodesic protocol in which one parameter is varied at a time (with µ being increased until µ * ) and the weak coupling approximation eq.(1); the minimal possible entropy production, τ Σ min /β = 2.57946 ± 1 • 10 −5 [ℏ], obtained when both parameters are changed simultaneously is also shown.

Figure 3 :
Figure 3: Parametrization of µ(t) and ε(t) described by eq.(116) for multiple values of k.Shown in the parameter space (left) and as a function of time (centre and right).

Figure 5 :
Figure 5: Comparison of the excess work W(1) = k B T Σ for a slow erasure protocol at constant coupling in the exact description (eq.(133)) with the weak coupling approximation (eq.(1) of the main text) and the lower bound of eq.(140).

Figure 7 :
Figure 7: Comparing two optimal erasure protocol with βε(1) ≈ 21 for two different values of k.The same is found for other values of k.