Active matter under control: Insights from response theory

Active constituents burn fuel to sustain individual motion, giving rise to collective effects that are not seen in systems at thermal equilibrium, such as phase separation with purely repulsive interactions. There is a great potential in harnessing the striking phenomenology of active matter to build novel controllable and responsive materials that surpass passive ones. Yet, we currently lack a systematic roadmap to predict the protocols driving active systems between different states in a way that is thermodynamically optimal. Equilibrium thermodynamics is an inadequate foundation to this end, due to the dissipation rate arising from the constant fuel consumption in active matter. Here, we derive and implement a versatile framework for the thermodynamic control of active matter. Combining recent developments in stochastic thermodynamics and nonequilibrium response theory, our approach shows how to find the optimal control for either continuous- or discrete-state active systems operating arbitrarily far from equilibrium. Our results open the door to designing novel active materials which are not only built to stabilize specific nonequilibrium collective states, but are also optimized to switch between different states at minimum dissipation.

Whilst the collective effects of active matter have been studied extensively, how to efficiently control active matter has only recently started to receive growing attention.Indeed, experiments have now demonstrated the ability to manipulate active nematics with a magnetic field [19], trigger spatial phase separation in bacteria [20], guide rotational patterns in magnetic rotors [21], and drive phase transitions in living cells [22].Moreover, the interplay of external control and activity has shown to induce interesting nonlinear behaviours, such as negative mobility, where optimal protocols could result in novel collective functions [23].Thus, the control of active matter opens unprecedented perspectives on developing novel nonequilibrium materials which selectively change collective states in response to perturbations.To this end, a first step is to establish generic principles for the reliable and optimal control of active systems.
Previous studies have explored the optimal control of active systems through the lens of minimal models.Such works have considered wet one-body control [24][25][26], one-body navigation strategies [25,[27][28][29][30], and more complex many-body scenarios and field theories [31][32][33][34][35].In these works, optimality typically involves penalizing protocols that steer away from the target state and/or rewards those finishing in the least time, without any constraint on how much energy is dissipated.Therefore, how to optimize control under thermodynamic constraints, such as minimizing the dissipated heat, is still largely unexplored for active systems.
For passive systems, recent advances in stochastic thermodynamics have led to a versatile framework for control at minimum dissipation [36][37][38][39][40][41][42].In this context, the operator controls the parameters of a given potential energy.The corresponding dissipation decreases with protocol time [Fig.1], so that the slowest protocol is always optimal, as expected from the thermodynamics of pas-FIG.1. Illustration of the heat dissipated for a control protocol with external perturbation in finite time.In passive systems, the protocol achieving the least dissipation is always the slowest one, as expected from standard thermodynamics.In active systems, the optimal protocol has instead a finite duration t * p which achieves the best trade-off between the dissipation stemming from the external perturbation, and that coming from internal activity.arXiv:2305.11078v2[cond-mat.stat-mech]15 Feb 2024 sive systems.For passive systems, the control framework amounts to geometric optimization [36,42]: optimal trajectories are geodesics on curved manifolds of a thermodynamic metric.The success of this geometric approach [38,39,41,43,44] fosters the hope that the underlying framework could also inform control beyond inherently passive systems.How such a framework can be adapted to active matter remains to be explored.
Underpinning the geometric approach for optimal control is linear response theory.It allows one to predict the effect of a weak perturbation in terms of correlations in the unperturbed dynamics.For passive systems, the relation between response and correlation functions follows the celebrated fluctuation-dissipation theorem [45], which relies on the time-reversal symmetry (TRS) of the dynamics.Interestingly, such a relation can be straightforwardly extended when TRS is broken [46,47].Indeed, it has already been used to study the response of various active systems [48][49][50][51][52].These recent developments establish a clear roadmap to extend the framework of thermodynamic control from passive to active systems.
In active systems, the energy consumption of individual constituents results in a constant rate of heat dissipation.Then, in contrast with the control of passive systems, dissipation now stems from two sources: (i) the perturbation by an external operator, and (ii) the internal energy consumption of particles.For a sufficiently large protocol duration, where the dissipation from the perturbation is small, the contribution from internal activity is ever-growing.For a sufficiently small protocol duration, the fast protocol incurs a large dissipation from the large perturbation, as with passive systems.As a result, slow protocols are no longer optimal for the active case, and, now, a finite protocol duration minimizes the dissipation [Fig.1].Therefore, the control of active systems involves not only finding optimal trajectories at fixed duration, but also finding which duration achieves the best trade-off between internal and external dissipation.An important challenge is then to rationalize how the optimal duration depends on the interplay between the external perturbation and the internal drive.
In this paper, inspired by previous efforts on passive systems, we herein utilize stochastic thermodynamics and response theory to systematically derive an optimal thermodynamic control framework for active systems.In so doing we build upon recent developments in and applications of linear response theory for active systems [48][49][50][51][52] and nonlinear response theory for driven thermal systems [53][54][55][56][57][58].First, we present, in Sec.II, our theoretical framework for the control of both continuous and discrete-state active systems.Second, we show, in Sec.III, how our framework can be deployed to inform the control of some specific model systems.We consider two cases: (i) an active Ornstein-Uhlenbeck particle confined in a harmonic trap with varying stiffness, and (ii) an assembly of active Brownian particles with purely repulsive interactions whose size is varied.In both cases, we obtain the optimal protocol for the corresponding con-trol parameter, and we discuss shared implications that the presence of self-propulsion has on deriving the optimal protocols.Overall, our results provide a systematic framework to address the control of active systems, and illustrate the main differences with respect to the control of passive systems.

II. OPTIMAL CONTROL
We propose a systematic recipe for optimizing the control of active systems.It consists in predicting how an external operator should vary the parameters of the potential energy to drive the system between two states at minimal cost.We consider the cost function to be the heat dissipated by the system to a thermal reservoir.Assuming weak and slow driving, we decompose the heat into specific correlations and averages, whose dependence on the control parameter determines the optimal protocol.Interestingly, the heat always features a minimum at the protocol time which achieves the best trade-off between internal and external dissipation [Fig.1].

A. Thermodynamic cost function
We consider a system of active particles immersed in a thermal bath at temperature T , where each particle i = {1, . . ., N } has an independent self-propulsion velocity v i that does not depend on position nor on any details of the protocol.We describe the motion of particles by an overdamped Langevin equation given as where r i is the position of particle i, µ is the mobility, ϕ(α, {r}) is the total potential energy that depends on the control parameter α and the set of all particle positions {r}, D = µT is the diffusion coefficient (with Boltzmann constant k B = 1), and η η η i is a Gaussian white noise with zero mean and unit variance.The potential energy can account for both particle interactions and external fields.In all what follows, the protocol consists in varying α(t) from its initial value α 0 = α(t = 0) to its final value α 1 = α(t = t p ) within the duration t p .
For our control problem, we choose the cost function to be the average heat dissipated into the thermal bath (held at temperature T ) [59][60][61][62][63], defined as where, unless stated otherwise, we perform a sum over repeated indices throughout this paper.Hereafter, the dot product (•) is interpreted within Stratonovich convention, for which standard rules of differential calculus carry over to stochastic variables [64].Substituting the dynamics from Eq. (1) into Eq.(2), and using the chain rule ṙi where ⟨X⟩ s denotes the steady-state average of observable X, ϕ 0 and ϕ 1 are the initial and final potential energies respectively, and We have assumed that the system is in a steady state at t = 0, yet it need not be the case at t = t p .Equation (3) can be regarded as the extension of the first law of thermodynamics (namely, the conservation of energy) for active systems [3,65].In addition to the potential change ⟨ϕ 0 ⟩ s − ⟨ϕ 1 ⟩ and the work rate α ⟨∂ α ϕ⟩, that are also present for passive systems, the heat features the dissipation rate J.In the passive limit (v i = 0 for all i) J vanishes, and one recovers the first law of thermodynamics in its standard form.The self-propulsion term v i /µ in Eq. ( 1) describes the force on particle i that results from the conversion of energy into directed motion.In active systems such a conversion involves multiple degrees of freedom, deliberately discarded in our framework, that typically provide additional contributions to the heat.Other models of active matter have been proposed which resolve these underlying nonequilibrium processes [66][67][68][69][70]. Interestingly, for the case where the chemical reactions between fuel molecules are always transduced into the motion of the active particle, any extra dissipation due to the presence of chemical reactions is only a background contribution, independent of the potential ϕ.Therefore, all the relevant contributions to the heat arising from ϕ are already accounted for in Eq. (3).
In the absence of external control ( α = 0), the heat reduces to the last term in Eq. ( 3), which scales linearly with time: ⟨Q⟩ = ⟨J⟩ s t p .Such a scaling illustrates that an active system, even when at rest, dissipates heat at a constant rate in order to sustain the self-propulsion of the particles.For a quasi-static protocol, namely for protocol durations much greater than the slowest relaxation timescale of the system (t p ≫ τ max ), α will change more slowly than any relaxation timescale of the system.Thus, the system goes through a series of steady states, so that all averages in Eq. ( 3) are now steady state averages: At large times, the heat scales linearly with protocol duration [Fig.1].
We are interested in describing corrections to the quasistatic heat that permit the prediction of optimal control protocols.To this end, we focus on cases where α varies weakly and slowly throughout the protocol.Thus, we can then express the average of a given observable ⟨X⟩ in terms of the linear (R 1 ) and second-order (R 2 ) response functions as with the time ordering t ≥ t ′ ≥ t ′′ , and where ∆α t,t ′ ≡ α(t) − α(t ′ ) is the external perturbation to the control parameter.The response functions will be explicitly defined later in Eq. (18).We then expand the perturbation as where ∆t = t − t ′ is the time increment.The expansions in Eqs.(6)(7) are inspired by a previous work on the control of passive systems [37], although, here, we expand to higher order in both ∆α and ∆t.In all what follows, we assume that these expansions are valid at all times within the protocol duration.It amounts to neglecting any abrupt change in the trajectory of α, and also regarding the averages in Eq. ( 3) as smooth functions of time t.Within this setting, we show in Appendix A that the heat can be cast as where α0 = α(t = 0) and α1 = α(t = t p ), and we have introduced the boundary term and the Lagrangian The functions {V, Φ, Σ, Λ, Ψ} are given by The decomposition of heat in Eqs.(8)(9)(10)(11) is one of the central results of this paper.By measuring response functions and averages at different values of α, one can systematically construct the dependence of the Lagrangian L in terms of α and α, and deduce the expression of the corresponding heat for any protocol.Indeed, experimentally measured linear and nonlinear response functions [71] can be injected into Eqs.(8)(9)(10)(11) to explicitly determine the functional dependence of the heat on the control parameter.The optimal protocol readily follows from L by straightforward optimization [Sec.II B].Importantly, Eqs.(8)(9)(10)(11) hold for any potential ϕ and selfpropulsion v i , which highlights the versatility of our approach.
In the absence of self-propulsion (v i = 0), the Lagrangian reduces to L = α2 Ψ(α), in agreement with previous results on controlling passive systems [37].As detailed in Sec.II B, the difference in L between active and passive systems drastically affects the corresponding dependence of ⟨Q⟩ in terms of t p .The heat monotonically decreases in the passive case, whereas it features a global minimum in the active case [Fig.1].Therefore, although our approach only accounts for leading-order corrections to the quasistatic regime, it is actually sufficient to capture the main qualitative change between optimizing either passive or active matter.

B. Optimal protocol
Our aim is to demonstrate that the heat decomposition in Eqs.(8)(9)(10)(11) entails a series of generic properties, both for the optimal protocol α(t) and for the corresponding optimal heat ⟨Q⟩.Interestingly, we derive such properties without actually specifying the explicit dependence of the functions {V, Φ, Σ, Λ, Ψ} on α [Eq.(11)].For an arbitrary potential ϕ and self-propulsion v i , we then predict how ⟨Q⟩ scales with the protocol time t p , and we provide an approximate estimation of the value t * p minimizing ⟨Q⟩.Our derivations essentially rely on drawing an analogy between the Lagrangian L in Eq. (10), and that of a harmonic oscillator with position-dependent mass [72].In that respect, V (α) and α2 (Ψ − Σ ′ )(α) respectively stand for the potential and kinetic energies, where the effective mass Ψ − Σ ′ here depends on the control parameter α.
In what follows, we assume that V (α) is bounded and (Ψ − Σ ′ )(α) ≥ 0 for all α between α 0 and α 1 , which ensures the validity of our analogy.The optimal trajectory then obeys the following Euler-Lagrange equation: As in standard Hamiltonian mechanics, one can also display the optimal trajectory in terms of a first integral of motion.Multiplying both sides of Eq. ( 12) by α, and integrating with respect to time t, we get The term E is analogous to the total energy in Hamiltonian mechanics.It is constant throughout the protocol, and it is constrained by E > max α V .
We deduce from Eq. ( 13) a relation between the protocol trajectory α(t) and the protocol speed α(t): which shows that the optimal trajectory obeys a firstorder ordinary differential equation.This differential equation can be solved by separation of variables, and it is fully determined by {E, V, Ψ, Σ}.In practice, the constant of motion E implicitly depends on the initial and final values {α 0 , α 1 }, and on the protocol duration t p .Indeed, separating variables in Eq. ( 14), and integrating throughout the protocol, we get Therefore, combining Eqs.(14)(15), the solution of the optimal trajectory α(t) follows directly.Although its explicit expression cannot be given analytically for generic potential ϕ and self-propulsion v i , numerical integration readily yields α(t) for given expressions of {V, Ψ, Σ}, and for arbitrary {α 0 , α 1 } and t p .Interestingly, we demonstrate in Appendix B that Eqs.(14)(15) are actually sufficient to predict how ⟨Q⟩ scales with protocol duration t p : where the expression of K(α 0 , α 1 ) follows from that of {Φ, Ψ, Σ} [Eq.(B7)].These scalings correspond to regimes where the dissipation is dominated either by the external driving of the control parameter α (namely for t p ≪ t * p ), or by the internal self-propulsion v i (namely for t p ≫ t * p ).The former reproduces the scaling expected for passive systems [37].The latter agrees with the quasistatic prediction [Eq.( 5)], where here the optimal protocol achieves the minimum value of ⟨J⟩ s , as expected.
In between the regimes of large and small t p , the heat reaches a minimum at t p = t * p [Fig.1].The protocol duration t * p achieves the best compromise between two regimes of high ⟨Q⟩.By extending the scalings in Eq. ( 16) to regimes where t p ≈ t * p , we approximate t * p as the value of t p where these scalings have same ⟨Q⟩: Equation (17) illustrates that t * p sets a trade-off between (i) the heat due to external drive, involving the system relaxation through response functions in the definition of K(α 0 , α 1 ) (equivalently in the definition of {Φ, Ψ, Σ}, see Eq. ( 11)), and (ii) the heat due to internal selfpropulsion, involving the steady state average ⟨J⟩ s .A qualitatively similar trade-off was reported in a different context for systems combining autonomous cycles and external drives [73].
In the passive limit (v i = 0), we deduce from Eq. ( 17) that t * p diverges.The heat now reaches its minimum asymptotically at large t p [Fig.1], and the scaling corresponding to t p ≪ t * p in Eq. ( 16) now extends to all t p .Indeed, as expected from the thermodynamics of passive systems, the protocol with smallest dissipation is always the slowest, in which case the heat reduces to its quasistatic value.
Interestingly, when controlling active systems, the boundary term B in Eq. ( 9) involves the protocol speeds α at initial and final times, respectively α0 and α1 .Within our framework, optimizing the Lagragian term L already leads to fixing α0 and α1 [Eqs.(14)(15)] for a given choice of {α 0 , α 1 } and t p .To optimize simultaneously B and L, one could potentially consider abrupt changes in the protocol trajectory α(t) at initial and final times [74].Yet, such discontinuous jumps would not be consistent with the expansions in Eqs.(6)(7).Therefore, in what follows, we focus on the optimal trajectory given by Eqs.(14)(15).

C. Response functions
We now demonstrate that the response functions in the heat decomposition of Eqs.(8)(9)(10)(11) can be expressed in terms of some specific correlations functions.Indeed, evaluating response functions generally requires measuring how the system is affected by weak perturbation, which can prove difficult both in numerical simulations and in experiments.Instead, our aim is to show that the response is actually already encoded in spontaneous fluctuations of the unperturbed dynamics.
In passive systems, the fluctuation-dissipation theorem enforces generic relations between linear response and correlation functions [45], and similar types of relations also exist for the second-order response [55,56,75].Interestingly, recent works have shown that responsecorrelation relations can be extended beyond passive systems [49], and in particular to active matter [48,[50][51][52].
Here, we explicitly derive such relations for the active dynamics in Eq. (1).
From the definition of the response functions in Eq. ( 6) for a generic observable X, we express R 1 and R 2 as Evaluating the response functions then amounts to predicting how the average ⟨X(t)⟩ is perturbed when varying the control parameter α.To this end, we express this average in terms of a path integral as where P η and P v denote the probabilites to observe a given trajectory for the noise terms η i and v i , respectively.We have used that these noises are independent.
Using that P η is Gaussian, integrating Eq. ( 19) with respect to position trajectories {r i } yields where Note that, although the dynamics is given in Stratonovich convention throughout the paper, S does not feature the term (µ/2) ´dt∇ i • f i usually given in the Onsager-Machlup action.This is because the path integral in Eq. ( 20) is written with respect to noise trajectories, instead of position trajectories.Substituting the path integral of Eq. ( 20) into Eq.( 18), we deduce where we have used that D[η, v] is independent of α, since noise trajectories are not affected by the control parameter.Then, predicting how S varies with α is a route to expressing R 1 and R 2 in terms of correlation functions.The perturbation α → α − ∆α affects the action as S → S + ∆S, where As detailed in Appendix C 1, combining Eqs.(22)(23) The connected correlation ⟨⟨X(t)Y (t ′ )⟩⟩ ≡ ⟨⟨X t Y t ′ ⟩⟩ is defined for arbitrary observables X and Y as In the absence of self-propulsion (v i = 0), we recover the celebrated fluctuation-dissipation theorem: Indeed, this theorem is readily obtained from Eq. ( 24) by considering the time-symmetrized response R 1 (X; t, t ′ ) − R 1 (X; t ′ , t), and using the time-reversal symmetry of equilibrium correlations [49] (see Appendix C 1).Similarly, we show in Appendix C 1 that the second-order response can be written as where the connected correlation ⟩⟩⟩ is defined for arbitrary observables X, Y , and Z as The relations in Eqs.(24)(25)(26)(27), coupled with Eqs.(8)(9)(10)(11), are the main relations of this paper.Indeed, they provide explicit connections between the response functions R 1 and R 2 , defined for an arbitrary observable X [Eq.( 6)], and correlation functions in the unperturbed dynamics.Similar linear and nonlinear response relations have been previously derived in the context of general Markov and non-equilibrium systems, with specific application to the Langevin equation and Ising spin variables [53,54].Our response relations Eqs.(24-27) provide a route to evaluating response functions resulting from external control without actually perturbing the system.Overall, we are now in a position to deploy these relations to the specific observables which define the response functions of interest in Eq. (11).

D. Discrete-state dynamics
In this section, we show that our generic framework for controlling the potential in Eq. ( 1), corresponding to a continuous-state dynamics, can be extended to nonequilibrium systems with discrete states.To this end, we consider a generic dynamics that is governed by the following master equation: where p i (t) is the probability for the system to be in the state i at time t, and K ij is the transition rate to go from state j to i.The control parameter α determines the expression of K ij , analogously to how α controls the potential ϕ in Eq. ( 1).The average of a given observable X and its time derivative are expressed as where X i is the projection of X over the state i.As for the continuous-state case, the system is in contact with a heat bath at temperature T , and subject to active driving.In practice, we consider transition rates with Arrhenius form, ensuring that the dynamics is thermodynamically consistent [77,78].Such a transition rate is given by where c ij = c ji is the rate amplitude, ϕ i is the potential energy of state i, and ε ij = −ε ji is the energy exchange associated with the active driving.
The heat dissipated into the bath while varying α from α 0 to α 1 , within a time duration t p , reads [79,80] Substituting the definition of the transition rates from Eq. ( 30) into Eq.( 31), we get where we have used the definition of averages given in Eq. ( 29), and we have introduced the active dissipation rate: The first law of thermodynamics in Eq. ( 32) takes a similar form as for the continuous-state case in Eq. ( 3).Therefore, within the same set of assumptions as in Sec.II A, the decomposition of heat as the sum of boundary and Lagrangian term [Eqs.(8)(9)(10)(11)] remains valid for the discrete-state case, as with all the results for the optimal protocol in Sec.II B. Again, to connect response functions to correlation functions, as for the case of continuous-state dynamics [Sec.II C], we rely on a path-integral approach.First, we define the path weight of a given trajectory ω([0, t p ]) which goes through a series of discrete states within the protocol time t p [60,81,82]: with the normalization factor N .The action A reads where we have introduced the reference transition rate Kij that does not depend on α.The stochastic density ρ j (t) is 1 if the trajectory ω([0, t]) goes through the state j at time t, and 0 otherwise: where δ ω,j is the Kronecker delta function.The stochas-tic transition rate n ij is defined as where γ (λ) ij denotes the series of times where the system jumps from state j to i, each one of them labelled by the index λ.The variation A → A + ∆A due to the perturbation α → α − ∆α then follows as (38) Using the same procedure as that yielding the responses in the continuous-state case, we show in Appendix C 2 that the linear response R 1 in the discrete case can be written as and the second-order response R 2 follows as where, for a given observable Y , we have introduced the notations Thus, Eqs.(39,40) provide explicit relations between response and correlation functions, which mirror those of the continuous-state case [Eqs.(24,26)] (also see [53,54]).Substituting these relations into the expressions in Eq. ( 11) then allows one to determine the boundary and Lagrangian terms of the dissipated heat.

III. FROM ONE-TO MANY-BODY CONTROL
We now apply our systematic recipe to the control of continuous-state systems in specific examples.First, in Sec.III A, we analytically derive the optimal protocol for an active particle subject to a harmonic potential whose stiffness varies.Second, in Sec.III B, we address the case of an assembly of repulsive active Brownian particles (RABPs) with a controllable size, that undergoes motility-induced phase separation (MIPS) at high packing fraction.

A. Active particle in a harmonic trap
Having established the main framework, we now test it on a toy model consisting of a single active particle in a one-dimensional harmonic trap [Fig.2(a)].The motion of the particle is governed by where ϕ(α, r) = αr 2 /2 is the trap potential, with α being the stiffness.For the self-propulsion, we choose the simplest Ornstein-Uhlenbeck process with ⟨v⟩ = 0 and ⟨v(t)v(s)⟩ = τ −1 D 1 e −|t−s|/τ , where τ is the persistence time and D 1 is the amplitude of the active noise.The heat [Eq.( 3)] associated with varying α from α(0) = α 0 to α(t p ) = α 1 reads where r 0 and r 1 are respectively the initial and final positions of the particle.Since the dynamics in Eq. ( 42) is linear, the time evolution of the moments ⟨r 2 ⟩ and ⟨rv⟩ is closed, and can be obtained using Itô's lemma [64] as In steady state, the left-hand side of Eq. ( 44) is zero, resulting in Therefore, the heat given in Eq. ( 43) can be numerically measured for a given protocol α(t) simply by simulating the dynamics of the moments in Eq. ( 44).
To evaluate the decomposition of heat [Eqs.(8)(9)(10)(11)], we use the continuous-state response functions in Eqs.(24) and (26).We show in Appendix D 1 how to obtain the following expressions (46) In the passive limit (D 1 = 0), the expressions reduce to Ψ = D/2α 3 µ 2 , Σ = V = 0, and Λ = αµΦ = D/(2αµ), in agreement with [37].The protocols for changing the stiffness α of the trap can be readily obtained following the procedure detailed in Sec.II B. First, we numerically solve the elliptic integral equation in Eq. ( 15) which connects the constant E with the protocol time t p .As E approaches max α V (namely for increasing t p ), one needs increasingly high precision to numerically solve Eq. (15).To this end, we utilize the arbitrary-precision numerical library Arb [83], that provides a state-of-the-art numerical integration procedure based upon adaptive bisection and adaptive Gaussian quadrature [84].Second, we obtain the protocol trajectory α(t) by numerically solving Eq. ( 14), for given protocol time t p and boundary conditions {α 0 , α 1 }.We report the results of optimal protocols for different t p in Figs.2(b-c).Very fast (t p ≪ t * p ) and very slow (t p ≫ t * p ) protocols collapse onto two separate master curves.This is in contrast with the thermodynamic control of passive systems, where there is only a single protocol master curve valid for all protocol durations [42].
For the optimal protocols in Fig. 2(c), we then compare two separate evaluations of heat.We either (i) simulate the dynamics in Eq. ( 44) with the optimal protocol, and measure the corresponding heat [Eq.( 43)], or (ii) directly substitute the optimal protocol in the expression of heat given in Eq. ( 8), which relies on response theory.Comparing these two results allows us to delineate the regime where the response theory is indeed valid.We report in Fig. 2(d) that the two evaluations of heat match very well not only at large t p , where the heat scales linearly with t p , but also in the regime where the heat is minimum (t p ≃ t * p ).At times smaller than t * p , a discrepancy arises, showing that the assumption of slow driving, underpinning the response framework, breaks down in this regime.
For a general potential, we expect that the protocol time at minimum heat t * p increases as the activity gets weaker (|v i | → 0).For the harmonic case treated here, we evaluate analytically the approximated estimation of FIG. 3. Optimal protocol duration t * p as a function of the active noise amplitude D1, as predicted by the scaling relation Eq. ( 17) in conjunction with Eq. (46).Parameters: t * p , which stems from matching the asymptotic behaviors of the heat, by substituting Eq. ( 46) into Eqs.( 17) and (B7).This estimation of t * p captures the divergence t * p ∼ 1/ √ D 1 at small activity (D 1 /D ≪ 1).Interestingly, it also predicts that t * p should plateau to a finite value at large activity (D 1 /D ≫ 1), as shown in Fig. 3.

B. Repulsive active Brownian particles (RABPs)
Having successfully tested our framework on a rather simple model, we now apply it to a many-body active system which features richer physics.Specifically, we consider a system of N RABPs in two spatial dimensions that can exhibit MIPS [Fig 4] [18].We take the dynamics in Eq. ( 1) with the potential energy given by where the pair potential U ij depends on the inter-particle distance r ij = |r i − r j |.To impose repulsive interactions, we take and U ij = 0 otherwise.The control parameter α here embodies the range of the repulsive interaction, which is akin to the excluded-volume radius.Therefore, changing α at a fixed number of particles amounts to changing the packing fraction, which can lead to a transition between homogeneous and phase-separated states [ Fig 4].Moreover, we define the self-propulsion of particle i as where v 0 is the constant magnitude of self-propulsion speed, and ζ i is a Gaussian white noise with zero mean and unit variance.The persistence time τ , which controls the angular noise, determines the self-propulsion correlations as ⟨v in (t)v im (0)⟩ = δ ij δ nm v 2 0 e −|t|/τ , where n and m refer to spatial coordinates [85].
Due to the difficulty in deriving analytical expressions for the decomposition of heat [Eqs.(8)(9)(10)(11)], we resort to numerically evaluating the corresponding averages and correlations given in Appendix D 2. We perform simulations at different particle sizes {α} in the homogeneous case (2.1 × 10 3 < αv/D < 2.4 × 10 3 ), and in the case of phase separation (3.0 × 10 3 < αv/D < 3.3 × 10 3 ).After obtaining data for various α values, we fit averages and correlations to deduce their continuous dependence in terms of α.We then use these fitting functions to obtain the control protocols through solving Eq. ( 15), where again we distinguish protocols for the homogeneous and MIPS cases.In practice, we consider N = 2000 for RABPs in a square box of side length L = 200, with periodic boundary conditions to mimic bulk conditions.The parameters are v = 1, µ = 1, D = 10 −3 , τ = 10 2 , and ε = 10 2 T .Considering that lengthscales and timescales are respectively taken in units of micrometers and seconds, the parameter values are chosen to be relevant to actual microswimmers [85], such as moving bacteria.We integrate the equation of motions [Eqs.( 1) and ( 48)] using a standard Euler discrete-update rule.
The resulting control protocols are shown in Figs.5(ab).In the homogeneous case, the control protocols appear to lie on a single master curve, as for the control of passive systems [42].In contrast, for the MIPS case, we obtain various protocol curves for different protocol durations t p , which collapse into distinct master curves either for very fast and very slow protocols.Interestingly, plotting the control speed α against the value of the control parameter α reveals that the protocol curves do not actually collapse onto a single master curve in the homogeneous case [Fig.5(c)], although the change in control speed is clearly smaller than that of the MIPS case [Fig.5(d)].Since the existence of a single master curve is a signature of passive systems [42], we speculate that RAPBs in the homogeneous state have macroscopic signatures that resemble a thermal system, at least as far as the heat decomposition [Eqs.(8)(9)(10)(11)] is concerned.Finally, we com-  pute the heat as predicted from the response framework [Eq.(8)].The optimal protocol duration t * p /τ ≈ 1 × 10 3 in the homogeneous case [Fig.5(e)] is approximately 200 times larger than in the MIPS case, where t * p /τ ≈ 2 × 10 5 [Fig.5(f)].This difference illustrates that the response to changes in particle size is slower in the MIPS case, where the cluster size adapts to the packing fraction, than in the homogeneous case.

IV. DISCUSSION
We have derived a thermodynamic framework for the optimal control of active systems that operate arbitrarily far from thermal equilibrium.The systematic nature of our approach relies on recent advances in stochastic thermodynamics and response theory, as inspired by previous works on the thermodynamic control of passive systems [36,37].Applications of our framework have revealed new insights into the control of active systems that are in direct contrast to the passive case.Specifically, we have found that: (i) for non-zero activity, expanding the dissipation to linear order in the protocol duration requires knowledge of the second-order response, whereas linear response is sufficient to obtain the same order in the passive case [37]; (ii) there is a cost (dissipation) as-sociated with the self-driving meaning that the optimal duration is not the longest one, at variance with passive systems [37,43]; (iii) we find that each protocol duration gives rise to a unique protocol curve which eventually collapse onto either a fast or slow master curve, in contrast to the strict many-to-one mapping as seen in passive control [42,86].
As a demonstration of the generality of our approach, we have built our framework to describe both continuous and discrete-state active systems.This opens the door to future work that may quantitatively compare optimal control scenarios between these different categories of active matter.Interestingly, different definitions of irreversibility and of dissipation have been proposed for various models of active matter [65,68,69,87,88].Indeed, the measure of irreversibility provides a legitimate quantification of dissipation only if the underlying dynamics satisfies specific thermodynamic constraints.Our results on discrete-state dynamics can be regarded as a motivation to develop active discrete-state models that satisfy these constraints thermodynamically.Moreover, our approach could also be adapted to field theories of active matter [89][90][91], some of which can be embedded within linear irreversible thermodynamics [92].Furthermore, although we have focused on optimizing heat, our framework can be straightforwardly adapted to other cost functions, for instance the extracted work (which has been optimized in the literature of passive systems [37-39, 42, 74, 86, 93-97]).Note that a recent study has considered optimizing the work in active matter by directly applying the framework of passive systems [98].
Our framework relies on the assumption of weak and slow driving of the control parameter.It leads to an explicit decomposition of heat which can serve as a seed, or a testing ground, for machine learning approaches beyond the smooth-driving assumptions [34,99].Interestingly, one could extend our framework to take into account corrections from higher-order responses, and thus to capture faster driving, within a systematic approach.Moreover, although we have focused here on controlling a single parameter for simplicity, it would be interesting to consider an arbitrary number of control variables in future works.Lastly, encouraged by the experimental implementations of the control framework in passive systems [40,100] and recent progress in measuring the response of living systems [101][102][103][104], our framework can be readily deployed by experimentalists to explore the control of experimental active systems.Indeed, even without knowing the nonequilibrium response-correlation relation in Sec.II C, our decomposition of heat in Eqs.(8)(9)(10)(11) already delineates which perturbation to apply and which response to measure.We anticipate that the potential success of applications of our framework will heavily depend on the quality of measured response functions.
Overall, our work paves the way towards designing active materials able to optimally switch between collective states.Indeed, while most studies of active matter focus on establishing phase diagrams, by associating con-trol parameter values with collective states [1,2], how to optimally induce transitions has remained largely unexplored.While our framework sets the stage to this end, there remains a series of open challenges.First, the assumption of slow driving precludes crossing a critical transition, since the relaxation of the system is always slower than any perturbation at criticality [105].Considering non-critical transitions instead, the corresponding order parameter is typically discontinuous, and so are the averages and correlations in the heat decomposition.Accordingly, one expects that the optimal protocol is no longer smooth when crossing phase boundaries, which challenges the assumption of weak driving.Therefore, it remains to be explored how the assumptions underlying the response framework need to be adapted when switching between collective states in active matter.
where we have assumed α 0 < α 1 .Therefore, it follows that the optimal trajectory α(t) follows a master curve when scaling t with t p , as for the control of passive systems [37].Moreover, the condition E ≫ max (B7) The term K accounts for contributions from both the boundary and Lagrangian terms.
Second, we deduce from Eqs. (14-15) that the protocol duration t p is large whenever E ≈ max α V .Therefore, for slow protocols (t p ≫ t * p ), the constant of motion E depends only weakly on the protocol duration t p .It follows from Eq. ( 14) that the protocol speed α(t) depends weakly on t p , and so does the boundary term B [Eq. ( 9 In this Appendix, we show how to explicitly relate response and correlation functions using a path probability representation of the dynamics.We consider two cases: (i) continuous-state dynamics [Eq.( 1)], and (ii) discretestate dynamics [Eq.( 28)].

Continuous-state dynamics
The linear response function R 1 , defined in Eq. ( 18), can be written in terms of the dynamic action S [Eq.(23) , (C3) where we have used Eq. ( 25).From Eq. ( 23), we get where we have used the chain rule ∂ α φ = ṙi •∇ i ∂ α ϕ in the unperturbed dynamics (namely, at constant α).Finally, substituting Eq. (C4) into Eq.(C3) yields the expression for the linear response function C5) in agreement with in Eq. (24).Symmetrizing the linear response function results in the following R 1 (X; t, t ′ ) − R 1 (X; t ′ , t) In the passive case (v = 0), the force acting on particles reduces to f i = −∇ i ϕ.
and the normalization factor is given by 1/N = ´e−S D[η, v].The term S is the standard Onsager-Machlup action [76]:

FIG. 5 .
FIG. 5.Optimal thermodynamic control of a many-body system consisting of repulsive active Brownian particles in the homogeneous (HOM) (left) and phase-separated (MIPS) (right) cases.(a-b) Derived protocols for driving the particle size α as a function of time.(c-d) Control speed α against control parameter α. (e-f) Scaled heat against protocol duration, computed using protocols shown in (a-b).