Introduction

The 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 self-powered 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 medicine1,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 growth10. 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 performance11,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 physics15,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 ramifications18,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 eigenvectors22,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 invisibility24, hypersensitive sensing25,26, loss-induced transparency27, and EP lasers28,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 sources30,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 systems32,33,34,35,36,37.

Results and discussion

Mathematical formulation

The system consists of two thermal reservoirs at different temperatures TH > TC, 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 \({\mathcal{C}}\) in parameter space. The design is chosen in a way that the motion along the path \({\mathcal{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.

Fig. 1: Schematic representations of our thermal motor.
figure 1

a A schematic photonic circuit connected to two thermal baths at different temperatures, TH > TC, 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 L1(2) and capacitances C1(2), connected via a coupling capacitance Cc. The two thermal baths are represented via Thevenin equivalent transmission lines, with characteristic impedance R and grounded voltage sources VL(R), and they are attached to the circuit via capacitances Ce1(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.

Full size image

In typical circumstances, the MDF X = {X1, X2, ... , XM} 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

$${\mathcal{M}}\ddot{{\bf{X}}}(t)={{\Gamma }}\dot{{\bf{X}}}(t)+{{\bf{F}}}_{{\rm{av.}}}+{\boldsymbol{\xi }}(t),$$
(1)

where \({\mathcal{M}}\) is the generalized inertia tensor, ξ(t) is a fluctuating force, and Γ is the friction tensor that satisfies a fluctuation-dissipation 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 \({\bf{x}}=\left\langle {\bf{X}}\right\rangle\), where 〈 ⋅ 〉 indicates thermal averaging. Finally, the mean force Fav., drives the mechanical rotor diverting energy from the photonic thermal current to produce mechanical work. In the photonic framework (Fig. 1a), Fav. is analogous to the radiation pressure associated with the radiation inside the circuit. In the electromechanical framework of Fig. 1b, Fav. 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

$${{\bf{F}}}_{{\rm{av.}}}=-\hslash \left\langle {{{\Psi }}}^{\dagger }\frac{\partial {H}_{0}}{\partial {\bf{x}}}{{\Psi }}\right\rangle =-\hslash \mathop{\sum}\limits _{n,{n}^{\prime}}{\left(\frac{\partial {H}_{0}}{\partial {\bf{x}}}\right)}_{n{n}^{\prime}}\left\langle {\psi }_{n}{(t)}^{* }{\psi }_{{n}^{\prime}}(t)\right\rangle .$$
(2)

where \({{\Psi }}={\left({\psi }_{1},\ldots ,{\psi }_{N}\right)}^{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 H0(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

$${\rm{i}}\frac{{\mathrm{d}}{{\Psi }}(t)}{{\mathrm{d}}t} = {H}_{{\rm{eff.}}}{{\Psi }}(t)+{\rm{i}}{D}^{T}{\theta }^{(+)}(t),\\ {H}_{{\rm{eff.}}} = {H}_{0}({\bf{x}})-{\rm{i}}\frac{{D}^{T}D}{2},\\ {\theta }^{(-)} = -{\theta }^{(+)}+D{{\Psi }},$$
(3)

where the matrix D, with elements \({D}_{n,\alpha }=\sqrt{2{\gamma }_{\alpha }}{\delta }_{n,\alpha }\), 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)=\mathop{\int}\nolimits_{0}^{\infty }f(\omega ){{\rm{e}}}^{-{\rm{i}}\omega t}d\omega\)), the amplitudes \({\theta }_{\alpha }^{(+)}(\omega )\) satisfy the relation

$$\left\langle {[{\theta }_{{\alpha }^{\prime}}^{(+)}({\omega }^{\prime})]}^{* }{\theta }_{\alpha }^{(+)}(\omega )\right\rangle =\frac{1}{2{\rm{\pi }}}{\tilde{{{\Theta }}}}_{\alpha }(\omega )\delta (\omega -\omega ^{\prime} ){\delta }_{\alpha ,{\alpha }^{\prime}},$$
(4)

where \(\tilde{{{\Theta }}}(\omega )={{\Phi }}(\omega )\cdot {{\Theta }}(\omega )\), with Φ(ω) being a noise filter function and \({{{\Theta }}}_{\alpha }(\omega )={\left\{\exp \left[\hslash \omega /({k}_{B}{T}_{\alpha })\right]-1\right\}}^{-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 \({\mathcal{C}}\) as39,40,41,42,43,44,45 (see Supplementary Note 1)

$$W =\int_{0}^{\infty }\frac{{\mathrm{d}}\omega }{2\pi }\mathop{\sum}\limits _{\alpha }{\tilde{{{\Theta }}}}_{\alpha }(\omega ){P}_{\alpha }(\omega ),\\ {P}_{\alpha } =\frac{\hslash }{{\rm{i}}}{\oint }_{{\mathcal{C}}}{\mathrm{d}}{\bf{x}}{\left[{({S}^{{\bf{x}}})}^{\dagger }{\partial }_{{\bf{x}}}{S}^{{\bf{x}}}\right]}_{\alpha ,\alpha }=\hslash {\oint }_{{\mathcal{C}}}{R}^{{\bf{x}}}{\partial }_{{\bf{x}}}{\alpha }^{{\bf{x}}}{\mathrm{d}}{\bf{x}}.$$
(5)

where \({S}^{{\bf{x}}}(\omega )=-{I}_{{N}_{\alpha }}+{\rm{i}}D{G}^{{\bf{x}}}(\omega ){D}^{T}\) is the unitary instantaneous scattering matrix and \({G}^{{\bf{x}}}={[\omega {I}_{N}-{H}_{{\rm{eff.}}}(x)]}^{-1}\) is the Green’s function associated with the effective Hamiltonian Heff. (Im 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 \({\mathcal{C}}\), it is a geometric quantity46,47. It turns out that for the two-reservoir setup of Fig. 1, Pα can be written only in terms of the reflectance Rx 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 \({\mathcal{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

$${P}_{\alpha }=\iint_{ A}\frac{\partial {\omega }_{p}}{\partial {x}_{p}}\frac{\partial {\omega }_{q}}{\partial {x}_{q}}{{\mathcal{W}}}_{\alpha }{\mathrm{d}}{x}_{p}{\mathrm{d}}{x}_{q},$$
(6)

where we have used that \({\left(\frac{\partial {H}_{0}}{\partial {x}_{p}}\right)}_{n,m}=\frac{\partial {\omega }_{n}}{\partial {x}_{p}}{\delta }_{n,m}{\delta }_{n,p}\), and we have introduced the work density per unit area as

$${{\mathcal{W}}}_{\alpha }={\mathop{\mathrm{lim}}\limits_{A\to 0}}{P}_{\alpha }/A=4{\gamma }_{\alpha }\hslash {\mathcal{R}}e\left({G}_{p\alpha }^{* }{G}_{pq}{G}_{q\alpha }-{G}_{q\alpha }^{* }{G}_{qp}{G}_{p\alpha }\right),$$
(7)

with \(A={\iint }_{A}\frac{\partial {\omega }_{p}}{\partial {x}_{p}}\frac{\partial {\omega }_{q}}{\partial {x}_{q}}{\mathrm{d}}{x}_{p}{\mathrm{d}}{x}_{q}\to 0\) (see Supplementary Note 2).

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 the \({\partial }_{{\bf{x}}}\times {\left[{({S}^{{\bf{x}}})}^{\dagger }{\partial }_{{\bf{x}}}{S}^{{\bf{x}}}\right]}_{\alpha ,\alpha }\,\ne\, 0\), and (b) the closed path \({\mathcal{C}}\) must enclose a nonzero area in the parameter space {x1, x2}. A biproduct of the last condition is that variations of x1, x2 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 Heff. 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 = TH and Tα=2 = TC. The effective Hamiltonian of such a reduced system reads

$${H}_{{\rm{eff.}}}=\left(\begin{array}{ll}{\omega }_{1}-{\rm{i}}{\gamma }_{1}&\kappa \\ \kappa &{\omega }_{2}-{\rm{i}}{\gamma }_{2}\end{array}\right),$$
(8)

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 Heff. is \({\omega }_{\pm }={\omega }_{0}-{\rm{i}}{\gamma }_{0}\pm \frac{1}{2}\sqrt{{\left({{\Delta }}\omega -{\rm{i}}{{\Delta }}\gamma \right)}^{2}+4{\kappa }^{2}}\), where \({\omega }_{0}=\frac{{\omega }_{1}+{\omega }_{2}}{2}\), \({\gamma }_{0}=\frac{{\gamma }_{1}+{\gamma }_{2}}{2}\) and Δω = ω1 − ω2 and Δγ = γ1 − γ2 ≠ 0. The corresponding (non-normalized) eigenvectors are

$${u}_{1,2}={\left(2\kappa ,-{{\Delta }}\omega +{\rm{i}}{{\Delta }}\gamma \pm \sqrt{{\left({{\Delta }}\omega -{\rm{i}}{{\Delta }}\gamma \right)}^{2}+4{\kappa }^{2}}\right)}^{T}.$$
(9)

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 γ027. 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 \({{\mathcal{W}}}_{\alpha }\).

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 \({\mathcal{C}}\) associated with changes to the resonant frequencies being \({\omega }_{n}={\omega }_{0}\left[1+\delta \cos ({x}_{n}+{\phi }_{n})\right]-{(-1)}^{n}\epsilon\), 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-of-phase i.e., ϕ1 = π/2, ϕ2 = 0. For such a modulation scenario, the associated enclosed area in the parameter space \(\left({\omega }_{1}({x}_{1}),{\omega }_{2}({x}_{2})\right)/{\omega }_{0}\) is A = πδ2.

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 \({{\mathcal{W}}}_{\alpha }\) in terms of Green’s function Gx. In fact, for the 2 × 2 case, the calculations for Green’s function can be carried out explicitly for any perturbation, giving

$${G}^{{\bf{x}}} =\frac{1}{D}\left(\begin{array}{ll}\omega -({\omega }_{2}-{\rm{i}}{\gamma }_{2})&-\kappa \\ -\kappa &\omega -({\omega }_{1}-{\rm{i}}{\gamma }_{1})\end{array}\right)\\ \approx \frac{{\mathcal{A}}}{{{{\Delta }}}_{{\rm{EP}}}}+\frac{{\mathcal{B}}}{{{{\Delta }}}_{{\rm{EP}}}^{2}},$$
(10)

where D = [ω − (ω2 − iγ2)][ω − (ω1 − iγ1)] − κ2 and ΔEP = ω − ωEP. In the above expression, the generic perturbation p is hidden in the parameters that define Heff., 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 \({\mathcal{A}}\) and \({\mathcal{B}}\) are frequency-independent matrices (see Supplementary Note 3). It turns out that the functional dependence of \({{\mathcal{W}}}_{\alpha }\) 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 Heff. 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 \({{\mathcal{W}}}_{\alpha }\) acquires its maximum value (see below).

A more elaborated treatment can extend the above analysis of Gx, 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 Heff. do not span the Hilbert space. Instead, one has to complete the eigenvectors of Heff. into a basis by introducing the associated Jordan vectors48,49. Following this approach, we can recover the last expression of Gx in Eq. (10).

Substituting the expression for the Green’s function back in Eqs. (6), (7)) we get

$${{\mathcal{W}}}_{1} = 4\hslash {\gamma }_{1}{\gamma }_{2}{\mathcal{R}}e\left[\frac{-2{\rm{i}}{\kappa }^{2}{D}^{* }}{| D{| }^{4}}\right]\\ = \frac{-8\hslash {\gamma }_{1}{\gamma }_{2}{\kappa }^{2}\left(2{\gamma }_{0}{{{\Delta }}}_{0}\,+\,\epsilon {{\Delta }}\gamma \right)}{{\left\{{\left[{{{\Delta }}}_{0}^{2}\,-\,{\gamma }_{0}^{2}\,-\,{\epsilon }^{2}\,+\,c\right]}^{2}+{\left[2{\gamma }_{0}{{{\Delta }}}_{0}\,+\,\epsilon {{\Delta }}\gamma \right]}^{2}\right\}}^{2}},$$
(11)

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 \({\mathcal{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 Heff. 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

$${{\mathcal{W}}}_{1}=-\frac{16\hslash {\gamma }_{1}{\gamma }_{2}{\kappa }^{2}{\gamma }_{0}{{{\Delta }}}_{0}}{{\left\{{\left[{{{\Delta }}}_{0}^{2}+{\gamma }_{0}^{2}+c\right]}^{2}-4{\gamma }_{0}^{2}c\right\}}^{2}}.$$
(12)

In fact, by considering the EP condition c = 0, we are able to identify in the denominator of \({{\mathcal{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 \({{\mathcal{W}}}_{1}\) is nonmonotonic and antisymmetric with respect to the EP resonance frequency axis ω = ω0 (Δ0 = 0) for all κ-values. Furthermore, \({{\mathcal{W}}}_{1}({{{\Delta }}}_{0}=0)=0={{\mathcal{W}}}_{1}({{{\Delta }}}_{0}\to \pm \infty )\) while its extrema occur in the vicinity of EP (see the gray-filled circle) at \({{{\Delta }}}_{0}=\omega -{\omega }_{0}=\pm \sqrt{\frac{1}{7}}{\gamma }_{0}\), see Fig. 2a.

Fig. 2: Work in the proximity of an exceptional point (EP).
figure 2

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 Heff., 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 short-dashed 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 × 1012 rad s−1, the temperature of the hot (cold) bath is TH = 109 K (TC = 3 K), and the decay rates of the two modes are γ1 = 0.01 ω0 and γ2 = 0.002 ω0. In c–e, T = (TH + TC)/2 is the average temperature of the reservoirs, A is the area in the parameter space, and kB is the Boltzmann constant.

Full size image

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 Heff.. An example case is when the resonances of the two coupled modes are detuned by ϵ. In this case, the diagonal elements of Heff. take the form \({\omega }_{1,2}={\omega }_{0}\pm \epsilon +{\omega }_{0}\delta \cos (x+{\phi }_{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

$${{\mathcal{W}}}_{1}=-\frac{8\hslash {\gamma }_{1}{\gamma }_{2}}{| D{| }^{4}}{\left(\frac{{{\Delta }}\gamma }{2}\right)}^{2}\left[2{\gamma }_{0}{{{\Delta }}}_{0}+\epsilon {{\Delta }}\gamma \right],$$
(13)

where now the denominator takes the form

$$| D{| }^{4}={\left\{{\left[{{{\Delta }}}_{0}^{2}-{\gamma }_{0}^{2}-{\epsilon }^{2}\right]}^{2}+{\left[2{\gamma }_{0}{{{\Delta }}}_{0}+\epsilon {{\Delta }}\gamma \right]}^{2}\right\}}^{2},$$
(14)

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

$${{\mathcal{W}}}_{1}\approx -\frac{8\hslash {\gamma }_{1}{\gamma }_{2}}{| {{{\Delta }}}_{{\rm{EP}}}{| }^{8}}{\left(\frac{{{\Delta }}\gamma }{2}\right)}^{2}\Bigg[2{\gamma }_{0}{{{\Delta }}}_{0}\\ +\epsilon {{\Delta }}\gamma \left(1-\frac{16{{{\Delta }}}_{0}^{2}{\gamma }_{0}^{2}}{| {{{\Delta }}}_{{\rm{EP}}}{| }^{4}}\right)\Bigg],$$
(15)

where the term associated with the perturbation ϵ is an even function in Δ0. We conclude, therefore, that the work density \({{\mathcal{W}}}_{\alpha }\) 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 \({{\mathcal{W}}}_{\alpha }\) 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 \({{\mathcal{W}}}_{\alpha }\), weighted with the function \({\tilde{{{\Theta }}}}_{\alpha }(\omega )\), 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 \({{\mathcal{W}}}_{\alpha }\), 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 \({{\mathcal{W}}}_{\alpha }\) with \(\tilde{{{\Theta }}}(\omega )\) 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 \({{\mathcal{W}}}_{\alpha }\) 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\approx A{{{\Theta }}}_{1}({\omega }_{0})\int d\omega /(2\pi ){{\mathcal{W}}}_{1}\propto \epsilon\). The same argument applies also in the case of spectral filtering with Φ(ω) = H(ω − ω0) (see red dashed line). In both cases, the extreme work Wmax. 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.

Time-domain simulations and implementation using electromechanical systems

We validate the above proposal by performing time-domain (TD) simulations using COMSOL software50 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 × 104 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

$${W}_{{\rm{TD}}}=\int_{0}^{2\pi }\tau (x)dx=2\pi \cdot {{\Gamma }}\cdot {{{\Omega }}}_{{\rm{s}}},$$
(16)

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 Cc = 2κC0 and C0 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 \({\mathcal{C}}\). We find that in the long time limit, the MDF reaches a terminal angular velocity \({{\Omega }}(t\to \infty )\equiv {{{\Omega }}}_{{\rm{s}}}^{{\rm{TD}}}\) which can be used in Eq. (16) for the numerical evaluation of WTD. 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 COMSOL50 (see “Methods”).

Fig. 3: Langevin dynamics and work.
figure 3

a–c The dynamics of the angular velocity Ω(t) for some representative values of the coupling coefficient κ, whose terminal velocity \({{{\Omega }}}_{{\rm{s}}}^{{\rm{TD}}}\) determines the work delivered by the photonic circuit. The black dashed lines indicate the associated theoretical values of the saturation velocity Ωs, see Eq. (16). d Work performed by the two-resonator circuit setup versus the coupling parameter κ. The numerical evaluation for the work from time-domain (TD) simulations (dots with error bars) is based on the mean value of the terminal angular velocity \({{{\Omega }}}_{{\rm{s}}}^{{\rm{TD}}}\), see Eq. (16). The error bars indicate the range of values for the work per cycle due to fluctuations of Ω(t) after saturation, see “Methods”. The TD simulations match nicely with the theoretical predictions for the work (green line), given by Eq. (5) in terms of the S-matrix. The blue dashed line reports the work predicted by the coupled-mode theory (CMT) modeling, see Eqs. (4) and (8) and “Methods”. The vertical red dotted line indicates the position of the exceptional point (EP).

Full size image

In Fig. 3d, we report a summary of the extracted WTD 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 TH and TC (TH > TC) is bounded by the so-called Carnot efficiency ηC = 1 − TC/TH51,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 high53,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 power44. The latter can be used, e.g., for lifting a weight or charging a capacitor. The usable output power is

$${P}_{{\rm{out}}}=\frac{{{{\Omega }}}_{{\rm{s}}}}{2\pi }W-\frac{{{{\Omega }}}_{{\rm{s}}}}{2\pi }{{\Gamma }}\int_{0}^{\frac{2\pi }{{{{\Omega }}}_{{\rm{s}}}}}{\dot{x}}^{2}{\mathrm{d}}t\approx \frac{{{{\Omega }}}_{{\rm{s}}}}{2\pi }W-{{\Gamma }}{{{\Omega }}}_{{\rm{s}}}^{2},$$
(17)

where we have assumed that the MDF has large inertia, forcing the rotor to move with terminal velocity \(\dot{x}\approx {{{\Omega }}}_{{\rm{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 \({{{\Omega }}}_{{\rm{s}}}^{* }=W/(4\pi {{\Gamma }})\) leading to \({P}_{{\rm{out}}}^{* }={\left(\frac{W}{4\pi }\right)}^{2}\frac{1}{{{\Gamma }}}\), which is half of the total frictionless power \({\left(\frac{W}{4\pi }\right)}^{2}\frac{2}{{{\Gamma }}}\). For circuit parameters such that \({{{\Omega }}}_{{\rm{s}}}\gg {{{\Omega }}}_{{\rm{s}}}^{* }\), the motor dissipates most of the incident energy, while in the other limiting case, where \({{{\Omega }}}_{{\rm{s}}}\ll {{{\Omega }}}_{{\rm{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, Pout, that is extracted from the motor during one period of the cycle \(2\pi /{{{\Omega }}}_{{\rm{s}}}^{* }\) when it operates under optimal conditions (i.e., \({{{\Omega }}}_{{\rm{s}}}={{{\Omega }}}_{{\rm{s}}}^{* }\)), to the total input power Pin delivered to the photonic circuit. Specifically

$${\eta }^{* }\equiv \frac{{P}_{{\rm{out}}}^{* }}{{P}_{{\rm{in}}}},\quad {P}_{{\rm{in}}}\approx {\bar{I}}_{b}+{\bar{I}}_{p},$$
(18)

where for the evaluation of Pin we have also considered the fact that the slow variation of the photonic network’s parameters induces a pumping energy current \({\bar{I}}_{{\rm{p}}}\) in addition to the energy current \({\bar{I}}_{{\rm{b}}}\) due to the temperature bias59. Both currents above are measured at the hot reservoirs. Since typically \({\bar{I}}_{{\rm{p}}}\ll {\bar{I}}_{{\rm{b}}}\) we can omit the pumped current from the denominator while we can substitute in Eq. (18) the maximum usable power as \({P}_{{\rm{out}}}^{* } \sim {W}^{2}/{{\Gamma }}\). Therefore, \({\eta }^{* } \sim (\frac{W}{{\bar{I}}_{{\rm{b}}}})\frac{W}{{{\Gamma }}}\), which suggests that the maximum \({\eta }^{* }\le \frac{{\eta }_{{\rm{C}}}}{2}\) might be expected in the parameter domain where W acquires its maximum values, see Figs. 2 and 3.

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}_{{\rm{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 \({{{\Omega }}}_{{\rm{s}}}^{* }\) is small enough such that the Born–Oppenheimer approximation is valid. From Fig. 4, we see that both η* and \({P}_{{\rm{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 optimized46,53,58,60,61,62,63. Our proposal identifies, as an optimal domain for the design of cycle \({\mathcal{C}}\), the parameter space in the proximity of an EP.

Fig. 4: Maximum power and efficiency.
figure 4

Maximum power \({P}_{{\rm{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 \({{{\Omega }}}_{{\rm{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.

Full size image

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 \({\mathcal{W}}(\omega )\). In typical circumstances, \({\mathcal{W}}(\omega )\) is antisymmetric with respect to the position of ωEP, leading to a near-zero total work W. When, however, the spectral work density \({\mathcal{W}}(\omega )\) 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 space64,65,66,67. The possibility to extend these design schemes for the realization of optimal quantum motors39,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 Z0 = 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πf0, with f0 = 1 MHz, and ϵ = 0.0488 ⋅ ω0. The capacitors C1,2 are considered as a pair of conductive plates separated by a median air gap \({d}_{0}=\frac{{e}_{0}\cdot A}{{C}_{0}}\), where e is the vacuum permittivity, A = 100 cm2 is the area of the capacitor plates, and \({C}_{0}=\frac{1}{{Z}_{0}\cdot {\omega }_{0}}\) is a median capacitance. The upper plates of the capacitors are assumed to be attached to a wheel (the MDF) of radius r = d0/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\cdot \cos ({\phi }_{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 mr2 = 7.5810−15 kg m2, 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 Cc = 2κ ⋅ C0, where the coupling coefficient κ is a tunable parameter of the simulations. The left/right resonators are coupled to hot/cold baths via capacitors Ce1 = 0.1 ⋅ C0 and Ce2 = 0.03 ⋅ C0, 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}\cdot \mathop{\sum }\nolimits_{i = 1}^{200}\sin \left({\tilde{\omega }}_{i}\,t+{\varphi }_{i}\right)\), where V0 = 1 V is the amplitude of the noise, \({\tilde{\omega }}_{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 \({\tilde{\omega }}_{1}=2\pi \,0.85\) MHz with an upper limit of \({\tilde{\omega }}_{200}=2\pi \,1.1\) MHz.

The error bars in Fig. 3d reflect the fluctuations in the numerical evaluation of \({{{\Omega }}}_{{\rm{s}}}^{{\rm{TD}}}\) and are extracted from the temporal analysis of Ω(t) as \({{{\Omega }}}_{{\rm{min/max}}}^{{\rm{TD}}}={\rm{min/max}}\left({{\Omega }}(t\in [{t}_{1},{t}_{\max }])\right)\), where t1 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 setup69 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 scatterer59 where the voltage on the left (right) node, L (R), connecting the Thevenin equivalent transmission line to the capacitor Ce1 (Ce2) is driven by time- dependent voltage sources \({v}_{{\rm{L(R)}}}^{{\rm{s}}}(t)={\mathcal{R}}e\left({V}_{{\rm{L(R)}}}^{{\rm{s}}}{{\rm{e}}}^{{\rm{i}}\omega t}\right)\). Both the voltage and the current measured at the node L(R) can be split as a sum of two contributions

$$\begin{array}{lll}{v}_{{\rm{L(R)}}}(t,\omega )&=&{\mathcal{R}}e\left[({V}_{{\rm{L(R)}}}^{+}+{V}_{{\rm{L(R)}}}^{-}){{\rm{e}}}^{{\rm{i}}\omega t}\right],\\ {i}_{{\rm{L(R)}}}(t,\omega )&=&\frac{1}{{Z}_{0}}{\mathcal{R}}e\left[({V}_{{\rm{L(R)}}}^{+}-{V}_{{\rm{L(R)}}}^{-}){{\rm{e}}}^{{\rm{i}}\omega t}\right].\end{array}$$
(19)

Here \({V}_{{\rm{L}}}^{+}={V}_{{\rm{L}}}^{{\rm{s}}}/2\) (\({V}_{{\rm{R}}}^{-}={V}_{{\rm{R}}}^{{\rm{s}}}/2\)) is the complex voltage amplitude of the incoming wave from the left (right) Thevenin equivalent source, while \({V}_{{\rm{L}}}^{-}\) (\({V}_{{\rm{R}}}^{+}\)) is the voltage amplitude of the outgoing scattered wave. In the particular case where \({V}_{{\rm{R}}}^{-}=0\), \({V}_{{\rm{L}}}^{-}\) (\({V}_{{\rm{R}}}^{+}\)) is the amplitude of the reflected (transmitted) wave and, thus, the reflection (transmission) coefficient results \({r}_{{\rm{L}}}=\frac{{V}_{{\rm{L}}}^{-}}{{V}_{{\rm{L}}}^{+}}\) (\({t}_{{\rm{R,L}}}=\frac{{V}_{{\rm{R}}}^{+}}{{V}_{{\rm{L}}}^{+}}\)). Similarly, when \({V}_{{\rm{L}}}^{+}=0\), the reflection (transmission) coefficient is \({r}_{{\rm{R}}}=\frac{{V}_{{\rm{R}}}^{+}}{{V}_{{\rm{R}}}^{-}}\) (\({t}_{{\rm{L,R}}}=\frac{{V}_{{\rm{L}}}^{-}}{{V}_{{\rm{R}}}^{-}}\)).