Enhanced energy harvesting near exceptional points in systems with (pseudo-)PT-symmetry

Exceptional point degeneracies, occurring in non-Hermitian systems, have challenged many well established concepts and led to the development of remarkable technologies. Here, we propose a family of autonomous motors whose operational principle relies on exceptional points via the opportune implementation of a (pseudo-)PT-symmetry and its spontaneous or explicit violation. These motors demonstrate a parameter domain of coexisting high efficiency and maximum work. In the photonic framework, they can be propelled by thermal radiation from the ambient thermal reservoirs and utilized as autonomous self-powered microrobots, or as micro-pumps for microfluidics in biological environments. The same designs can be also implemented with electromechanical elements for harvesting ambient mechanical (e.g., vibrational) noise for powering a variety of auxiliary systems. We expect that our proposal will contribute to the research agenda of energy harvesting by introducing concepts from mathematical and non-Hermitian wave physics. Optimising energy harvesting is relevant to a variety of micro and nanotechnologies with applications ranging from self powered nanodevices to microfluidic pumps in biological environments. Here, the authors show that an opportunely designed opto(/electro)-mechanical circuit can simultaneously demonstrate high efficiency and power production when operating in the vicinity of an exceptional point degeneracy implemented via a pseudo-PT symmetry.

T he objective of energy harvesting is to supply power to a variety of systems ranging from structural health monitoring systems to the billions of autonomous wireless sensors associated with the Internet of everything and to selfpowered micro-/nanorobots and micropumps for microfluidics in biological environments. In particular, the manipulation of microscopic objects via currents, for example, has become an indispensable tool in many disciplines of science and technology, revolutionizing a variety of applications in areas as diverse as microengineering and microrobotics, to biology and medicine [1][2][3][4][5][6][7][8][9][10] . Recently, for example, the use of a photon-driven micropump was reported to control the flow of fluid in the vicinity of a neuron, thus affecting its growth 10 . Depending on the application, the source of these currents varies from thermal radiation and vibrations to electrical and chemical energy extracted in biological processes. On the fundamental level, such applications require the development of design principles that will allow us to realize powerful and efficient engines that operate between two reservoirs at different effective temperatures (or chemical potentials) and produce useful work with maximum efficiency. The situation is even more complex when one abandons the convenience of macroscopic frameworks (where recipes from traditional thermodynamics are available) and delves into the challenges of modern nanodevices, where wave interferences and thermal fluctuations dominate their performance [11][12][13][14] . In this realm, the management of wave interferences via the implementation of appropriate symmetries in composite structures can be proven crucial for the design of high-performance engines.
Along these lines, the development of engineering schemes that manipulate dynamical symmetries in order to enhance wave-matter interaction has been extremely fruitful over the last few years. A prominent example is parity-time (PT) symmetric wave physics [15][16][17] , which has initiated a paradigm shift and triggered increasing attention to non-Hermitian wave transport. During the last 10 years, this activity has led to a number of surprises with immediate technological ramifications [18][19][20] . Many of these concepts are directly related to the notion (and subsequent manipulation) of non-Hermitian spectral degeneracies known as exceptional points (EP) 20,21 . As opposed to traditional (diabolic) degeneracies occurring in the spectrum of Hermitian systems, the EP degeneracy is characterized by coalescence of both eigenfrequencies and the corresponding eigenvectors 22,23 . This collapse of the eigenbasis has many important consequences in the transport properties of a system and can lead to a plethora of fascinating phenomena like unidirectional invisibility 24 , hypersensitive sensing 25,26 , loss-induced transparency 27 , and EP lasers 28,29 . At the heart of many of these phenomena is the fact that, in the proximity of an EP degeneracy, the wave-matter interaction is enhanced. Recently, however, there was raising concern about the efficiency of EP design schemes under the influence of ambient noise sources 30,31 . This viewpoint is consistent with the general physics/engineering guidelines that consider noise as an anathema-an undesirable, but unavoidable, contribution of nature that tends to mitigate physical processes. Here, we will advocate for exactly the opposite viewpoint, aiming to change this narrative and highlight the coexistence of noise, with an opportunely designed EP, as an alternative spectral engineering paradigm for efficient energy harvesting from ambient noise sources.
Here, we show that the near-field thermal radiation, emitted from a hot reservoir toward a cold reservoir, can be harvested by an optomechanical circuit as a nonconservative (wind-)force. Under its influence, a (slow) mechanical degree of freedom (MDF) undergoes a closed path periodic motion. We show that for a long-but finite-driving period of the MDF, these circuits act as autonomous radiative motors. Their extracted power and efficiency are maximized when they are designed to operate in a domain of their parameter space that is in the vicinity of an EP degeneracy. The origin of this enhancement is traced to the coalescence of the corresponding eigenmodes, which maximizes the spectral power density in the vicinity of the EP. Such EP appears in the spectrum of the effective non-Hermitian Hamiltonian that describes the coupling of these motors with the thermal reservoirs. We have engineered such EPs by enforcing a (pseudo)-PT symmetry via an opportunely arranged coupling with the reservoirs, which leads to a differential loss between modes of the system. In this case, the effective non-Hermitian Hamiltonian is PT-symmetric invariant, after performing a renormalization with respect to the average losses. We have exploited the influence of the EP in the motor performance enhancement, by appropriately violating (spontaneously or explicitly) the designed (pseudo)-PT symmetry and by manipulating the thermal emissivity of the attached reservoirs. Our predictions can guide the design toward optimal operational conditions of autonomous motors. Their applicability extends beyond the photonic framework to other platforms like electromechanical circuits for harvesting mechanical (e.g., vibrational) ambient noise for the power supply of a variety of auxiliary systems [32][33][34][35][36][37] .

Results and discussion
Mathematical formulation. The system consists of two thermal reservoirs at different temperatures T H > T C , which are brought in contact via a circuit. The latter is coupled to a MDF from which we extract work. For simplicity, we will assume that the circuit incorporates two single-mode resonators whose frequencies ω n (where n = 1, 2) are modulated by the motion of the MDF. The temperature gradient between the two reservoirs produces a thermal current through the circuit that, in turn, exerts a force to the MDF engaging them in slow periodic motion along a given closed trajectory C in parameter space. The design is chosen in a way that the motion along the path C creates out-of-phase variations in ω 1,2 , thus leading to a violation of spatiotemporal symmetries of the structure. One possible implementation of the above setup is in the photonic framework Fig. 1a, while a parallel proposal in the electromechanical framework is shown in Fig. 1b. Below, we will mainly use the photonic terminology associated with Fig. 1a, while we will also have the electromechanical scenario of Fig. 1b in mind.
In typical circumstances, the MDF X = {X 1 , X 2 , ... , X M } describes a change in position or angle of the mechanical element due to a respective force or torque. For concreteness of our presentation, we will assume that M = 2. The coordinate X abides by the Langevin equation where M is the generalized inertia tensor, ξ(t) is a fluctuating force, and Γ is the friction tensor that satisfies a fluctuationdissipation relation. In our analysis below, we will assume that the fluctuating force ξ can be neglected due to the large inertia of the MDF. Consequently, we can approximate the dynamics of X by its mean value x ¼ X h i, where 〈 ⋅ 〉 indicates thermal averaging. Finally, the mean force F av. , drives the mechanical rotor diverting energy from the photonic thermal current to produce mechanical work. In the photonic framework (Fig. 1a), F av. is analogous to the radiation pressure associated with the radiation inside the circuit. In the electromechanical framework of Fig. 1b, F av. is associated with a torque acting on the mechanical rotor.
The interaction between the mechanical part and the radiation is obtained from the variation of the energy inside the photonic circuit due to displacement x of the MDF. Specifically, the thermal averaged force is where Ψ ¼ ψ 1 ; ; ψ N À Á T , ψ n is the field amplitude at the n − th resonator of the circuit, (ℏω n |ψ n | 2 represents the energy density in the n − th mode/resonator, and H 0 (x) is the effective Hamiltonian of the circuit that provides a description of the dynamics of the radiation field in the single-mode resonators. The dynamics of the open system (circuit coupled with reservoirs) is described in terms of a temporal coupled-mode theory (CMT) 38 i dΨðtÞ dt ¼ H eff : ΨðtÞ þ iD T θ ðþÞ ðtÞ; where the matrix D, with elements D n;α ¼ ffiffiffiffiffiffiffi 2γ α p δ n;α , describes the coupling of the circuit with the thermal baths. The thermal excitations from (toward) the α-th reservoir are given by the incoming (outgoing) complex fields θ (+) [θ (−) ]. In the frequency domain (using the Fourier transform f ðtÞ ¼ whereΘðωÞ ¼ ΦðωÞ Á ΘðωÞ, with Φ(ω) being a noise filter function and Θ α ðωÞ ¼ exp _ω=ðk B T α Þ ½ À1 f g À1 is the Bose-Einstein statistics describing the mean number of photons that are emitted from reservoir α with frequency ω. Finally, T α is the temperature of the α-th reservoir. Work density in the adiabatic limit. We assume that the dynamics of the MDF Eq. (1) occurs on timescales much larger than the ones associated with the field dynamics Eq. (4). Under this assumption, we can invoke the Born-Oppenheimer approximation and obtain the work performed by the motor along the path C as [39][40][41][42][43][44][45] (see Supplementary Note 1) where S x ðωÞ ¼ ÀI N α þ iDG x ðωÞD T is the unitary instantaneous scattering matrix and G x ¼ ½ωI N À H eff : ðxÞ À1 is the Green's function associated with the effective Hamiltonian H eff. (I m is the m × m identity matrix). In Eq. (5), the kernel P α (ω) indicates the spectral response of the system at a frequency ω. Since P α (ω) only involves a parametric integral along the path C, it is a geometric quantity 46,47 . It turns out that for the two-reservoir setup of Fig. 1, P α can be written only in terms of the reflectance R x and the corresponding reflection phase α x (see the right part of Eq. (5)). As a matter of convention, a positive W in Eq. (5) indicates that the dynamics of x follows the positive direction of the path C.
An analytically useful expression of P α is achieved by substituting in Eq. (5), the scattering matrix in terms of Green's function. We get where we have used that ∂H 0 ∂x p n;m ¼ ∂ω n ∂x p δ n;m δ n;p , and we have introduced the work density per unit area as . Direct inspection of Eq. (5) allows us to establish the following two conditions for the implementation of our proposal as a motor: (a) the force has to be nonconservative, which means that ≠ 0, and (b) the closed path C must enclose a nonzero area in the parameter space {x 1 , x 2 }. A biproduct of the last condition is that variations of x 1 , x 2 with a phase difference 0 or π cannot produce work.
Engineering EP degeneracies. In the frequency range near an EP degeneracy, the resolvent of the effective Hamiltonian H eff. can be approximated by a 2 × 2 subspace involving only the resonant modes associated with the EP. We therefore consider a minimal model consisting of two coupled modes with resonant frequencies ω 1 , ω 2 . Alternatively, one can consider, as a concrete example, the setup of Fig. 1a. The system consists of two single-mode resonators coupled asymmetrically to two reservoirs at temperatures T α=1 = T H and T α=2 = T C . The effective Hamiltonian of such a reduced system reads where κ describes the coupling between the two modes and γ 1 , γ 2 are the (asymmetric) decay rates of the two modes due to their coupling with the two reservoirs. The spectrum of H eff. is T H > T C , that is able to divert part of the thermal radiative energy into useful work in the form of the motion of a mechanical degree of freedom (MDF) described by a rotor. b An electromechanical motor consisting of two LC resonators, with grounded inductances L 1(2) and capacitances C 1(2) , connected via a coupling capacitance C c . The two thermal baths are represented via Thevenin equivalent transmission lines, with characteristic impedance R and grounded voltage sources V L(R) , and they are attached to the circuit via capacitances C e1(e2) . The capacitance plates of the LC resonators are coupled to pistons whose motion is out-of-phase with one another, thus exercising a torque to a rotor. When the system operates in the vicinity of exceptional points (EPs) and violates specific spatiotemporal symmetries (or the baths are subjected to spectral filtering), the motor operates at optimal performance.
(non-normalized) eigenvectors are It is easy to show that when Δω = Δω EP = 0 and κ = κ EP = Δγ/2 the system supports an EP degeneracy with ω + = ω − = ω EP = ω 0 − iγ 0 . In fact, under the condition Δω = Δω EP , the Hamiltonian Eq. (8) respects a pseudo-parity-time (PT) symmetry that reveals itself after renormalizing the losses with respect to their mean value γ 0 27 . Below, we will be discussing in detail two distinct scenarios involving perturbations around the EP that violate this pseudo-PT symmetry either spontaneously or explicitly. We will show that each of these cases affects in a different manner the characteristic features of the work density W α .
Work density in the presence of an EP. We analyze the extracted work density of the motor when the center of the modulation cycle is in the proximity of an EP. To this end, we consider a modulation cycle C associated with changes to the resonant frequencies being ω n ¼ ω 0 1 þ δ cosðx n þ ϕ n Þ Â Ã À ðÀ1Þ n ϵ, where ϵ describes a resonance detuning that displaces the unmodulated system Eq. (8) from the EP by violating explicitly its pseudo-PT symmetry. In order to satisfy the criteria for nonzero work, we have assumed that the two resonances are modulated out-ofphase i.e., ϕ 1 = π/2, ϕ 2 = 0. For such a modulation scenario, the associated enclosed area in the parameter space Next, we assume a generic perturbation p, which displaces the center of the modulation cycle with respect to EP. Using Eq. (7), we have evaluated the work density W α in terms of Green's function G x . In fact, for the 2 × 2 case, the calculations for Green's function can be carried out explicitly for any perturbation, giving In the above expression, the generic perturbation p is hidden in the parameters that define H eff. , e.g., in the frequencies ω 1,2 = ω 1,2 (p) and/or the coupling κ = κ(p) between the two resonant modes. When p → 0, the Green's function can be approximated with the last expression, where A and B are frequencyindependent matrices (see Supplementary Note 3). It turns out that the functional dependence of W α on ω, in the vicinity of EP, is affected by the presence of the square Lorentzian term on the last part of Eq. (10). This unique spectral feature is a consequence of the degeneracy of the eigenvectors of H eff. at the EP. Furthermore, a squared Lorentzian lineshape implies a narrower emission/absorption peak and greater resonant enhancement in comparison with a nondegenerate resonance at the same complex frequency. We will show that the competition between the two terms appearing at the right equality of Eq. (10) determines the conditions under which W α acquires its maximum value (see below). A more elaborated treatment can extend the above analysis of G x , in order to include any number of modes, by using a degenerate perturbation theory that takes into consideration the singular nature of EPs. In this case, the standard modal decomposition of the Green's function is not applicable since the biorthogonal eigenvectors of H eff. do not span the Hilbert space. Instead, one has to complete the eigenvectors of H eff. into a basis by introducing the associated Jordan vectors 48,49 . Following this approach, we can recover the last expression of G x in Eq. (10).
Substituting the expression for the Green's function back in Eqs. (6), (7)) we get where Δ 0 = ω − ω 0 , and for the evaluation of the contour integral in Eq. (6), we have explicitly written ω 1,2 in terms of the parameters ω 0 and ϵ that define the position of the path C. The constant c = (Δγ/2) 2 − κ 2 and/or the detuning ϵ indicate the degree of deviation from EP. Let us exploit further Eq. (11) by considering two specific examples corresponding to perturbations that preserve/violate the pseudo-PT symmetry of the effective unmodulated Hamiltonian H eff . In the first case, we displace the system away from the EP by varying the coupling κ ≠ κ EP while keeping ϵ = 0. We find that the work density takes the form In fact, by considering the EP condition c = 0, we are able to identify in the denominator of W 1 above, the signature of the square-Lorentzian anomaly associated with the collapse of the eigenvector basis. Equation (12) allows us to conclude that W 1 is nonmonotonic and antisymmetric with respect to the EP resonance frequency axis ω = ω 0 (Δ 0 = 0) for all κ-values.
while its extrema occur in the vicinity of EP (see the gray-filled circle) at The situation is different when we choose to perturb the system away from the EP using a parameter that enforces an explicit pseudo-PT-symmetry violation of the unmodulated Hamiltonian H eff. . An example case is when the resonances of the two coupled modes are detuned by ϵ. In this case, the diagonal elements of H eff. take the form ω 1;2 ¼ ω 0 ± ϵ þ ω 0 δ cosðx þ ϕ 1;2 Þ. Furthermore, the work density does not have a definite symmetry with respect to (ω − ω 0 ). To be concrete, we consider the particular case κ = κ EP = Δγ/2 for which the work density is where now the denominator takes the form demonstrating the traces of the square-Lorentzian anomaly. The latter is better appreciated in the limit of ϵ = 0 (EP condition). For ϵ ≪ ω 0 , we can further expand up to leading order in ϵ the denominator and get where the term associated with the perturbation ϵ is an even function in Δ 0 . We conclude, therefore, that the work density W α loses the parity as soon as ϵ is turned on, see also Fig. 2b. Below we will be discussing the consequences of such an effect in the power extraction of the autonomous motor.
Work in the presence of an EP. We are now ready to exploit the properties of W α for the design of autonomous motors with optimal performance. To this end, we remind that the extracted work W is essentially the frequency integral of W α , weighted with the functionΘ α ðωÞ, see Eq. (5). Let us first discuss the family of perturbations that preserve the pseudo-PT symmetry of the unmodulated effective Hamiltonian. In this case, the antisymmetric form of the work density W α , with respect to the ω 0 − axis, results in a near-zero total work, see Fig. 2c. The slight deviation from zero (toward positive W > 0) is due to the fact that Eq. (5) involves a product of W α withΘðωÞ which slightly desymmetrizes the integrand toward smaller frequencies (see the blue solid line). We can revert the situation by introducing a spectral filtering function Φ(ω), which enhances the unbalanced contribution of positive and negative work densities in the integral of Eq. (5). The resulting extracted work, for the example case of a filter function Φ(ω) = H(ω − ω 0 ), is reported in Fig. 2c with a black dashed line (H(x) is the Heaviside function). Our results indicate that such a spectral filtering approach can lead to an increase in W, which is higher by two orders of magnitude with respect to the unfiltered case. The same data indicate that the maximum work occurs in the vicinity of the EP where W α acquires its maximum value (gray vertical line) and where the desymmetrization strategy via spectral filtering is more impactful.
An alternative way to induce an asymmetric integrand in Eq. (5) is by deviating the system away from the EP via a perturbation that will explicitly violate the pseudo-PT symmetry of the unmodulated effective Hamiltonian. In the previous section, we have identified one such perturbation being the frequency detuning ϵ between the two resonators. In this case, the work density itself becomes asymmetric (see Fig. 2b), leading to a frequency integral Eq. (5), which is different from zero. In fact, the maximum W occurring in the proximity of κ EP , is again two orders of magnitude enhanced in comparison to the ϵ = 0 case, see the blue dashed line in Fig. 2c.
The enhancement of the extracted work W via engineered perturbations that violate the pseudo-PT symmetry of the motor is better appreciated in Fig. 2d. Here, we report the extracted work W (for fixed κ = κ EP ) for both spectrally unfiltered/filtered noise versus the perturbation ϵ. For the unfiltered case (red solid line), we find that in the vicinity of the EP, the total work is proportional to ϵ, a relation that is a direct consequence of the expansion Eq. (15) for the work density. Specifically, assuming for simplicity that Θ 1 (ω) ≈ Θ 1 (ω 0 ), the integration over ω leads to the conclusion that W % AΘ 1 ðω 0 Þ R dω=ð2πÞW 1 / ϵ. The same argument applies also in the case of spectral filtering with Φ (ω) = H(ω − ω 0 ) (see red dashed line). In both cases, the extreme work W max. occurs at perturbation strengths ϵ max. in the vicinity of the EP, where the linear approximation Eq. (15) breaks down. An additional conclusion that we extract from the above analysis is that the spectral filtering method (combined with perturbations that violate the pseudo-PT symmetry) leads to a slight (twofold) increase of the extracted work as compared to the unfiltered case (see red solid line).
A panorama of the extracted work W versus ϵ and κ is shown in Fig. 2e. Here, we report only the unfiltered case, i.e., Φ(ω) = 1. The data demonstrate nicely that the extreme value of the extracted work occurs in the vicinity of (ϵ, κ) = (0, κ EP ) where the EP is located. The case of spectral filtering with a function Φ(ω) (e.g., Φ(ω) = H(ω − ω 0 )) shows the same qualitative features (with the only difference that W is flat in the negative ϵ semiplane due to the specific filter function) and therefore is not reported here. . Work density (color scale) as a function of the frequency ω of the incident radiation and the coupling constant κ (a) and resonance detuning ϵ (b). Blue and green solid lines are the eigenfrequencies of H eff. , which demonstrate an EP degeneracy (gray dot) at κ = κ EP = 0.5 Δγ (Δγ is the differential loss) in a (where ϵ = 0) and at ϵ = 0 in b (where κ = κ EP ). c The work W for ϵ = 0 (ϵ = 0.006ω 0 ) is indicated with a blue solid (dashed) line. We also report the work in the case where we have introduced spectrally engineered (SE) reservoirs (black shortdashed line) via the spectral filtering function Φ(ω) = H(ω − ω 0 ) (H(x) is the Heaviside function). d Work with SE reservoirs for ϵ = 0 (red solid line) and ϵ ≈ 0.006 ω 0 (red dashed line). In (c, d), the vertical dashed gray lines indicate the position of the EP at κ = κ EP and ϵ = 0, respectively. e We report the total work W as a function of κ and ϵ. The EP is indicated with a gray dot (the black dashed line indicates the value of κ EP ). In the vicinity of EP, the work becomes extreme (positive or negative). In these calculations, the average resonance frequency is ω 0 = 200 × 10 12 rad s −1 , the temperature of the hot (cold) bath is T H = 10 9 K (T C = 3 K), and the decay rates of the two modes are γ 1 = 0.01 ω 0 and γ 2 = 0.002 ω 0 . In c-e, T = (T H + T C )/2 is the average temperature of the reservoirs, A is the area in the parameter space, and k B is the Boltzmann constant.
Time-domain simulations and implementation using electromechanical systems. We validate the above proposal by performing time-domain (TD) simulations using COMSOL software 50 with a realistic electromechanical system, see Fig. 1b. See "Methods" for a detailed description of the circuit setup.
In the simulations, the wheel is given an initial angular velocity Ω 0 = 2.5 × 10 4 rad/s. Its angular velocity is monitored as a function of time until it saturates at a particular value, Ω s . From here, we evaluate the work per cycle via the relation where τ is the torque produced by the capacitor plates on the wheel, x is the angular displacement, and Γ is the friction coefficient. The subindex TD indicates that the evaluated work is extracted from our time-domain simulations. In Fig. 3a-c, we show the transient dynamics of the angular velocity Ω(t) for three typical coupling constants κ, where κ is such that C c = 2κC 0 and C 0 is a median capacitance (see Methods). Notice that in some cases (e.g., Fig. 3a), the angular velocity Ω(t) acquires negative values indicating that the wheel rotates opposite to the direction of the closed path C. We find that in the long time limit, the MDF reaches a terminal angular velocity Ωðt ! 1Þ Ω TD s which can be used in Eq. (16) for the numerical evaluation of W TD . In each of the subfigures 3a-c, we also indicate the theoretical values of the saturation velocity Ω s (see dashed black line). The latter has been extracted via Eq. (16), where the work W on the left-hand side has been calculated using Eq. (5). For the theoretical evaluation of W, we have extracted the elements of the instantaneous S-matrix of the circuit using a frequency-domain analysis of COMSOL 50 (see "Methods").
In Fig. 3d, we report a summary of the extracted W TD versus the coupling constant κ. In the same figure, we also plot the theoretical predictions for the work W (green line) that have been derived using Eq. (5) with instantaneous scattering matrix elements given by the COMSOL frequency analysis of the electromechanical system. Finally, in the same figure, we are presenting the predictions of the CMT modeling of Eqs. (4) and (8). In the latter case, the various parameters (coupling, resonance frequencies, and linewidths) of the CMT model have been extracted from the transmission spectrum of the electronic circuit (see "Methods"). The nice agreement between CMT and TD simulations confirms the validity of our CMT modeling and establishes the influence of the EP protocols in extracting maximum work from thermal autonomous motors.
Efficiency. In the framework of thermal engines, the question of maximum efficiency has been addressed by the pioneering work of Carnot that pointed out that the efficiency of a thermal engine that performs a cycle between two reservoirs with temperatures T H and T C (T H > T C ) is bounded by the so-called Carnot efficiency η C = 1 − T C /T H 51,52 . Of course, this thermodynamic bound is of limited practical importance since the corresponding heat engine must work reversibly, and thus its output power is zero. A more practical direction is to identify conditions under which the power of irreversible thermal engines, working under finite-time Carnot cycles, is optimized while their efficiency is still high [53][54][55][56][57][58] .
For the setup shown in Fig. 1b, the temperature gradient between the two thermal reservoirs induces a thermal current that goes through the motor. Part of the associated input power is dissipated due to friction, resulting in a reduction in the amount of usable output power 44 . The latter can be used, e.g., for lifting a weight or charging a capacitor. The usable output power is where we have assumed that the MDF has large inertia, forcing the rotor to move with terminal velocity _ x % Ω s . The optimal terminal angular velocity that maximizes the usable work is dictated by the parameters of the setup and can be found from Eq. (17) to be Ω Ã s ¼ W=ð4πΓÞ leading to P Ã out ¼ W 4π À Á 2 1 Γ , which is half of the total frictionless power W 4π À Á 2 2 Γ . For circuit parameters such that Ω s ) Ω Ã s , the motor dissipates most of the incident energy, while in the other limiting case, where Ω s ( Ω Ã s , the friction can be neglected, but the device does not generate much power. In both limits, the usable output power is nearly zero. It is, therefore, useful to quantify the performance of an autonomous motor by introducing its efficiency η. The latter is defined as the ratio of the net usable average output power, P out , that is extracted from the motor during one period of the cycle 2π=Ω Ã s when it operates under optimal conditions (i.e., Ω s ¼ Ω Ã s ), to the total input power P in delivered to the photonic circuit. Specifically where for the evaluation of P in we have also considered the fact that the slow variation of the photonic network's parameters induces a pumping energy current I p in addition to the energy current I b due to the temperature bias 59  An efficient way to test the above expectations of the performance of our EP-influenced motor is by simultaneous evaluation of its efficiency, Eq. (18), together with the corresponding power P Ã out . These quantities are plotted in Fig. 4 as a function of the parameters κ, ϵ associated with the coupling and the resonance detuning between the two LC resonators of the electromechanical system of the previous section. For these calculations, we have used the CMT modeling with parameters that reproduce the results of the direct TD simulations of COMSOL for the electromechanical motor (see Fig. 3d). Furthermore, we have ensured that the angular frequency Ω Ã s is small enough such that the Born-Oppenheimer approximation is valid. From Fig. 4, we see that both η * and P Ã out acquire their maximum values in the vicinity of the EP-albeit at slightly different (κ, ϵ)-parameter values. This is because of a natural trade-off between efficiency and extracted power, which has triggered a number of recent studies to identify conditions where this trade-off is optimized 46,53,58,[60][61][62][63] . Our proposal identifies, as an optimal domain for the design of cycle C, the parameter space in the proximity of an EP.
Conclusion. We have theoretically proposed and numerically demonstrated a large enhancement of the performance of thermal motors when they are operating in a parametric domain that is in the proximity of an EP degeneracy. The latter appears in the spectrum of the effective non-Hermitian Hamiltonian that describes the open circuit and is achieved via an opportunely arranged (differential) coupling of the isolated circuit with the ambient baths. In the proximity of the EP, the eigenvector basis collapses (eigenvector degeneracy), leading to an enhanced spectral work density WðωÞ. In typical circumstances, WðωÞ is antisymmetric with respect to the position of ω EP , leading to a near-zero total work W. When, however, the spectral work density WðωÞ is desymmetrized, the total extracted power and the motor efficiency can acquire their maximum values in the domain of the parameter space that is in the vicinity of the EP. We have shown that this desymmetrization can occur either via an explicit PT-symmetry violation of the unperturbed system or via a spontaneous symmetry where one needs to supplement it with additional spectral filtering of the radiation of the bath.
Our results pave the way toward the development of a generation of optimal thermal motors that utilize engineered non-Hermitian spectral degeneracies. The proposed scheme can find applications for on-chip photonics (e.g., self-powered microrobots or micropumps in microfluidics) and electromechanical systems for harvesting ambient noise for powering a variety of auxiliary systems. It will be interesting to extend our study of motor efficiency to cases where the closed path in the parameter space is in the proximity of an EP degeneracy of higher order. It is plausible that the higher-order divergence of the resolvent will lead to further enhancement of the total work. Similar questions emerge in the case where there are more than one EP, hybrid EPs, and anisotropic EPs, in the proximity of the closed path in the parameter space [64][65][66][67] . The possibility to extend these design schemes for the realization of optimal quantum motors 39,68 is also another promising direction. These, and other, questions will be addressed in a separate publication.

Methods
Circuit setup and time-domain simulations. The setup consists of a pair of capacitively coupled resonators with impedance Z 0 = 70 Ohm tuned at different frequencies, ω 1,2 = ω 0 ± ϵ, which enforces violation of the pseudo-PT symmetry of the unmodulated system. In our simulations, we have considered that ω 0 = 2πf 0 , with f 0 = 1 MHz, and ϵ = 0.0488 ⋅ ω 0 . The capacitors C 1,2 are considered as a pair of conductive plates separated by a median air gap d 0 ¼ e 0 ÁA C 0 , where e is the vacuum permittivity, A = 100 cm 2 is the area of the capacitor plates, and C 0 ¼ 1 Z 0 Áω 0 is a median capacitance. The upper plates of the capacitors are assumed to be attached to a wheel (the MDF) of radius r = d 0 /10 in a way that during the wheel rotation with angular velocity Ω, the plates will undergo a motion described by the displacements d 1;2 ¼ d 0 þ r Á cosðϕ 1;2 Þ with ϕ 1 = Ωt and ϕ 2 = ϕ 1 + π/2. The wheel is assumed to have mass m = 1 g, moment of inertia I = 0.5 mr 2 = 7.5810 −15 kg m 2 , and experiences friction with the ambient medium with friction coefficient Γ = 2.5 ⋅ 10 −13 N m s rad −1 . The coupling capacitance between the two LC resonators is C c = 2κ ⋅ C 0 , where the coupling coefficient κ is a tunable parameter of the simulations. The left/right resonators are coupled to hot/cold baths via capacitors C e1 = 0.1 ⋅ C 0 and C e2 = 0.03 ⋅ C 0 , respectively, which yields the following value for the critical coupling κ EP = 0.001625 (red dotted line on Fig. 3d).
To enhance further the extracted work from the MDF, we have introduced, in addition to the detuning ϵ, spectral filtering of the thermal baths. Specifically, the hot bath is producing a noise signal consisting of 200 spectrally uniformly distributed harmonics vðtÞ ¼ V 0 Á ∑ 200 i¼1 sinω i t þ φ i À Á , where V 0 = 1 V is the amplitude of the noise,ω i is the frequency of each noise harmonic, and φ i is a random phase shift. The lower frequency of the noise considered in the simulations isω 1 ¼ 2π 0:85 MHz with an upper limit ofω 200 ¼ 2π 1:1 MHz.
The error bars in Fig. 3d reflect the fluctuations in the numerical evaluation of Ω TD s and are extracted from the temporal analysis of Ω(t) as Ω TD min=max ¼ min=max Ωðt 2 ½t 1 ; t max Þ ð Þ , where t 1 is the time during which Ω(t) reaches the theoretical value of Ω s for the first time for a given value of κ; t max ¼ 0:13 s is maximum time used in a TD analysis.
Finally, we point out that the CMT Hamiltonian's parameters that best fit the scattering spectrum of the circuit setup (see next section) are γ 1 ≈ 3.6 × 10 −3 , γ 2 ≈ 3.2 × 10 −4 , ω 1 ≈ 1 − κ, and ω 2 ≈ 0.94 − κ, in units of ω 0 = 2π × 1 MHz. We have used a detailed CMT modeling associated with the circuit setup 69 for both the calculation of work in Fig. 3d and for our extensive simulations of Fig. 4.
Scattering parameters for the electrical circuit. We consider the circuit setup in Fig. 1b as a scatterer 59 where the voltage on the left (right) node, L (R), connecting the Thevenin equivalent transmission line to the capacitor C e1 (C e2 ) is driven by Here is the complex voltage amplitude of the incoming wave from the left (right) Thevenin equivalent source, while V À L (V þ R ) is the voltage amplitude of the outgoing scattered wave. In the particular case where V À R ¼ 0, V À L (V þ R ) is the amplitude of the reflected (transmitted) wave and, thus, the reflection (transmission) coefficient results . Fig. 4 Maximum power and efficiency. Maximum power P Ã out (z axis) and efficiency η * normalized with respect to the maximum efficiency η C /2 (color scale) at operational conditions corresponding to optimal terminal angular velocity Ω Ã s . These quantities are plotted as a function of the coupling κ and resonance detuning ϵ for a fixed temperature gradient. A perturbation on κ respects the pseudo-PT-symmetric nature of the unmodulated system, while a perturbation on ϵ violates this symmetry. In these extensive simulations, we have used the coupled mode theory modeling with parameters associated with the circuit setup: decay rates γ 1 ≈ 3.6 × 10 −3 , and γ 2 ≈ 3.2 × 10 −4 , resonance frequencies of each resonator ω 1 ≈ 1 − κ, and ω 2 ≈ 0.94 − κ, all measured in units of the natural frequency ω 0 = 2π × 1 MHz.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.