Coexistence of stable and unstable population dynamics in a nonlinear non-Hermitian mechanical dimer

Non-Hermitian two-site ``dimers'' serve as minimal models in which to explore the interplay of gain and loss in dynamical systems. In this paper, we experimentally and theoretically investigate the dynamics of non-Hermitian dimer models with non-reciprocal hoppings between the two sites. We investigate two types of non-Hermitian couplings; one is when asymmetric hoppings are externally introduced, and the other is when the non-reciprocal hoppings depend on the population imbalance between the two sites, thus introducing the non-Hermiticity in a dynamical manner. We engineer the models in our synthetic mechanical set-up comprised of two classical harmonic oscillators coupled by measurement-based feedback. For fixed non-reciprocal hoppings, we observe that, when the strength of these hoppings is increased, there is an expected transition from a $\mathcal{PT}$-symmetric regime, where oscillations in the population are stable and bounded, to a $\mathcal{PT}$-broken regime, where the oscillations are unstable and the population grows/decays exponentially. However, when the non-Hermiticity is dynamically introduced, we also find a third intermediate regime in which these two behaviors coexist, meaning that we can tune from stable to unstable population dynamics by simply changing the initial phase difference between the two sites. As we explain, this behavior can be understood by theoretically exploring the emergent fixed points of a related dimer model in which the non-reciprocal hoppings depends on the normalized population imbalance. Our study opens the way for the future exploration of non-Hermitian dynamics and exotic lattice models in synthetic mechanical networks.

Within this field, significant effort has been devoted to the study of so-called "dimer" models, in which two sites are coupled together in a (non-Hermitian) Hamiltonian that can be represented generally as a 2 × 2 matrix. Such models are useful as they serve as minimal systems in which to understand (often analytically) the effects of non-Hermitian terms on dynamical behavior [42][43][44][45]. Some of the interesting dynamical effects that have been investigated in nonlinear dimer systems include sensitivity to input power [7,46], non-reciprocal dynamics and directed transport [47,48], and confinement in phase space [49]. Many of the previous works focused on the dimer models where the non-Hermiticity is introduced and controlled through on-site loss and gain terms (see, e.g., [5-7, 42-45, 50-52] and references therein), that is, through the diagonal components in the Hamiltonian.
In this paper, we instead explore dynamics of dimer models in which the non-Hermitian effects are introduced by making the hopping between the two sites asymmetric [53]. Such systems are inspired by the Hatano-Nelson (HN) model [54][55][56], which is a one-dimensional lattice model with asymmetric hoppings between the sites exhibiting the so-called non-Hermitian skin effect [21,28,57]. Such asymmetric, or non-reciprocal, hoppings are more difficult to experimentally implement than on-site gain and loss. However, there has been significant recent progress in realizing such asymmetric hoppings using setups based on optical systems [58][59][60][61], electrical circuits [62][63][64] and synthetic mechanical metamaterials [20,65], opening up a possibility to experimentally study such non-Hermitian models.
In this paper, we engineer non-Hermitian dimer models by taking advantage of the flexibility of a mechanical set-up consisting of two coupled harmonic oscillators with measurement-based feedback [65]. As we have previously demonstrated, this approach can be used to simulate near-arbitrary mean-field lattice Hamiltonians, with controllable on-site gain and loss, non-reciprocal couplings and (exotic) synthetic nonlinearities amongst other effects [20,65]. In addition to its tunability, a key advantage of this set-up is that it provides full access to the dynamics, allowing us to observe the evolution of the system in real time.
Here, we exploit this mechanical set-up to investigate the dynamics of HN dimer models. We first discuss the linear HN dimer in which the asymmetric hopping is externally fixed [53]. This linear model is mathematically arXiv:2302.03572v3 [cond-mat.mes-hall] 1 Aug 2023 equivalent to the aforementioned dimer models with onsite gain and loss, with the result that the trajectories of dynamics in the PT -symmetric regime show closed stable orbits, while the behavior in the PT -broken regime is unstable, with the population exponentially exploding or decaying. Building on the understanding of the linear model, we then explore a nonlinear version of the Hatano-Nelson dimer in which the asymmetric coupling is induced by population imbalance between the two sites. For this model, we again find both a stable regime at low coupling strengths, in which population oscillations are bounded, analogous to the PT -symmetric regime described above, and an unstable regime at high coupling strengths, in which the population grows/decays, similar to the PT -broken regime. However, we also find a new regime at intermediate coupling strengths, where the two types of behavior coexist, allowing us to tune the population dynamics from a stable oscillation to an unstable divergence by simply tuning the initial phase difference between the oscillators. As we discuss, this behavior can be understood by studying the fixed points of a variant of the nonlinear HN dimer model, which we call the instantaneous HN dimer model, in which the non-reciprocal coupling depends on the normalized population difference between the two sites as introduced below. Our work lays the foundation for exploring the dynamics of more exotic lattice Hamiltonians with non-Hermiticity and mean-field interactions.
The paper is structured as follows: In Sec. II we introduce and define the models we explore. We introduce the linear HN dimer model, the nonlinear HN dimer model, and the instantaneous HN dimer model that helps elucidate the dynamics of the original nonlinear HN dimer model. In Sec. III we describe the experimental approach and set-up, and the details and parameters chosen. In Sec. IV, we analyze the linear HN dimer model. We first analytically study the dynamical behavior of the model, and then we compare the results with numerically and experimentally obtained dynamics. In Sec. V, we analyze the dynamical behavior of the instantaneous HN dimer model, and we discuss the emergence of multiple fixed points and the structure of the transition between the weakly and strongly interacting regimes. In Sec. VI, we finally study the nonlinear HN dimer model. We give qualitative explanations of the phase diagram using results obtained in previous sections, and we discuss the coexistence of different phases in the dynamics from both numerical and experimental approaches. Finally, in Sec. VII, we draw conclusions and discuss the outlook for this work.

II. MODELS
We experimentally realize two types of non-Hermitian dimer models. The first is the linear Hatano-Nelson (HN) dimer model, in which Hermiticity is broken by externally-tuneable non-reciprocal couplings between the a) 1. a) Sketch of the Hatano-Nelson dimer as described in Eq. 1 where δJ is the hopping asymmetry parameter. b) Sketch of the dimer with population-dependent hopping asymmetry as described in Eq. 2, where z describes the population imbalance between the two sites and g is a control parameter.
two sites. The second is the nonlinear HN dimer model, in which the non-reciprocal couplings are instead induced by the population imbalance between the two sites, and hence evolve dynamically, depending on the interparticle interaction strength. To obtain an analytical understanding of the nonlinear HN dimer model, we also theoretically introduce a variant of the nonlinear HN dimer model in which the non-Hermiticity only depends on the normalized population imbalance between the two sites; such a model describes the dynamics of the nonlinear HN dimer model for a short period of time, and thus we call it an instantaneous HN dimer model. We shall now introduce the linear, nonlinear, and instantaneous HN dimer model in turn.

A. Linear Hatano-Nelson Dimer Model
The linear HN dimer model is a model in which two sites are coupled by non-reciprocal hopping amplitudes [28,53,66,67]. By writing the complex-valued wave-function of the two sites as α = (α 1 , α 2 ) T , the dynamics of the linear HN dimer is described by where ω is an overall energy offset and ∆ determines the on-site energy difference between the two sites (with ℏ = 1), as shown in Fig. 1 (a). Without loss of generality, from now on we will consider J, ∆ and δJ to be nonnegative real values. The coupling between the sites is split into the reciprocal part of the hopping amplitude −J and the non-reciprocal part ±δJ, the latter being responsible for breaking Hermiticity. This model is the two-site version of the famous Hatano-Nelson model for a one-dimensional chain [54][55][56] in which all the nearestneighbor inter-site couplings take the form as in Eq. 1.
As mentioned above, this model is mathematically equivalent to the more commonly studied non-Hermitian dimer model with reciprocal coupling and on-site gain and loss, c.f. e.g. Refs. [42][43][44][45]. To see this, the Hamil-tonian (i.e. the two-by-two matrix in Eq. 1) can be expressed with Pauli matrices as: H = −∆ σ z −J σ x −iδJ σ y , up to an overall energy shift. A suitable unitary transformation can bring this into the form H = √ ∆ 2 + J 2 σ x − iδJσ z , which describes a reciprocal hopping of strength √ ∆ 2 + J 2 and on-site imaginary terms of magnitude δJ with opposite signs corresponding to the gain and loss; the known behavior of this model can therefore be used to infer that of the linear HN model. We will instead analyze the linear HN model directly in Sec. IV in order to lay the groundwork for the rest of the paper.

B. Nonlinear Hatano-Nelson Dimer Model
In the linear HN dimer model, the non-Hermiticity was included via a constant off-diagonal contribution ±δJ.
In the nonlinear HN dimer model, the non-Hermiticity is introduced dynamically via the population imbalance between the two sites. The nonlinear HN dimer model is defined by the following equation of motion: in which the non-reciprocal part of the hopping amplitude is now set by ±gz, where g is the interaction strength and z ≡ |α 1 | 2 − |α 2 | 2 represents the population imbalance between the two sites. A sketch of this model is made in Fig. 1 (b). As z depends on the complex wavefunction, this is no longer a linear model and so does not have linear eigenstates and eigenenergies similar to the linear HN dimer. Understanding the rich behavior of this model is a goal of the current paper. Because of the nonlinearity, it is not possible to obtain a good analytical understanding of the model. In order to approximately understand the dynamics of this model for a short period of time, we introduce the following instantaneous HN dimer model.

C. Instantaneous Hatano-Nelson Dimer model
The nonlinearity in the nonlinear HN dimer model was in the term ±gz in the off-diagonal terms. We can rewrite this term as ±gn(z/n), where n is the norm of the wavefunction. In non-Hermitian models, the norm n generally depends on time, and so do gn and z/n. However, as we discuss in more detail in Sec. V, the model which is obtained by replacing gn in the nonlinear HN dimer model byḡ, allows a description of the dynamics without complication arising from the time-dependence of n ifḡ is taken to be a constant value. Physically, we can view the model as describing the dynamics of the nonlinear HN dimer model at a time close to t 0 if we set Sketch of the experimental apparatus used for implementing the described models. We perform real-time measurements of the effective position X k and momentum P k variables of our oscillators (labeled by index k), which are used to compute real-time feedback forces that are applied magnetically through a voltage-controlled current. The form of the feedback forces used to implement a desired Hamiltonian H is given simply by the relationship F k ∝ −∂H/∂X k . g = gn(t 0 ). We thus call the model the instantaneous HN dimer model. Before moving on to describe the dynamical features of the linear, nonlinear, and instantaneous HN dimer models, we briefly explain how the models can be experimentally realized in our system of mechanical oscillators.

III. EXPERIMENTAL SET-UP
Our experimental setup is described in detail in Ref. [65], and it consists of two almost identical oscillators, whose acceleration a(t) and its numerical derivative, the jerk j(t) = da(t)/dt, are measured in real time. These measurements can be used to produce a real-time feedback force on each oscillator that allows us to synthetically couple the two oscillators. This procedure enables us to effectively generate a wide variety of two-level tight-binding Hamiltonians. Within the rotating-wave approximation, the mean-field Schrödinger equations of motion for a generic tight-binding Hamiltonian can be mapped onto Newton's equations for the harmonic oscillators. Hence, we can engineer a combination of self and cross feedback [65]: the former is responsible for local on-site terms (e.g. site-dependent potential energy shifts, site-dependent gain or loss, and on-site nonlinear interactions), while the latter introduces off-site terms that allow energy to hop from site to site (e.g., complex hopping, non-reciprocal coupling, and even densitydependent hopping).
Experimentally, feedback forces are implemented by using real-time voltage output signals to control the currents through gradient solenoids surrounding each oscillator. These control signals result in magnetic forces on the oscillators, each of which features an embedded dipole magnet. The real-time voltage measurement signals re-lating to acceleration and jerk are normalized to a common scale. Given the simple harmonic nature of our oscillators, these normalized measurement signals X k and P k , where k is the oscillator index, serve as direct proxies for the oscillators' positions and momenta, x k and p k . As depicted in Fig. 2, a desired Hamiltonian H is implemented by feeding back on the oscillators with forces that are true to Hamilton's equation (dp k /dt = −∂H/∂x k ). Hence, feedback forces are of the form F k ∝ −∂H/∂X k . Naturally, site energy shifts are implemented by forces of the form F k ∝ X k , while gain and loss terms are implemented through feedback of the form F k ∝ P k , and so on [65]. Owing to our co-normalization of the X k and P k variables, the relative magnitude of all the linear terms in our experimentally implemented Hamiltonians is defined simply by the ratio of the applied feedback coefficients. The absolute calibration of our feedback forces (i.e., how the control voltage signals generated by our measurement-and-feedback system relate to the actual mHz-scale terms of the implemented Hamiltonian) is performed by investigating the frequency shift of the individual oscillators based on self-feedback forces (F k ∝ X k ), as detailed in Ref. [65].
To simulate the dynamics as in Eq. 1, we introduce a combination of self-feedback (to cancel natural loss terms and to shift the site energies) and linear cross-feedback (to introduce hopping between sites). To capture nonlinear terms, such as the population-dependent hopping contributions of Eq. 2, we introduce feedback forces of the form where X k,i and P k,i are the initial values taken by the effective position and momentum variables. Through this normalization, these nonlinear terms are governed by the same absolute calibration as the linear terms.
In addition to allowing for implementation of the desired linear and nonlinear HN dimer models, our control over these applied forces enables us to set the initial state of the oscillators for each experiment. Starting with oscillators nominally at rest, we sinusoidally drive the two oscillators at their common resonance frequency of ∼ 3.05 Hz for several seconds. By controlling the relative strength and phase of these two sinusoidal "initialization" drives, we control the initial amplitudes and phases of the two oscillators (or correspondingly, the initial complexvalued wave-function α).

IV. LINEAR HN DIMER MODEL
With the experimental setup we just described, we study the dynamical behavior of the linear HN dimer model. To understand the dynamics, we employ two descriptions, one is in the "phase space" of the population imbalance and the phase difference between the two sites, and the other is the Bloch sphere representation. Although these dynamics can already be obtained either by looking at the eigenstates of the Hamiltonian or through the mathematical equivalence between the linear HN dimer model and the on-site gain-loss model mentioned in Sec. II, the description in terms of the phase space and the Bloch sphere is worth discussing here as it gives us an intuitive picture for understanding, and will also lay the groundwork for later Sections to understand the dynamics of the nonlinear and instantaneous HN dimer models.

A. Dynamical Equations for the linear HN dimer
We rewrite here the equation of motion for the linear HN dimer model, Eq. 1, In standard Hermitian quantum mechanics, the overall normalization of the wavefunction is not an observable, and we can assume the wavefunction to be normalized with the overall phase left as a gauge degree of freedom. However, in non-Hermitian quantum mechanics, time evolution is non-unitary in general and the change of norm of the wavefunction over time plays an important role, describing phenomena such as decay and lasing. In our classical experimental setup, all the information of the wavefunction, α 1 and α 2 , can in principle be measured. Following the practice of non-Hermitian quantum mechanics, we consider the norm of the wavefunction to be an observable, but we choose to not consider the overall phase to be significant. In describing the dynamics of two-site (two-level) systems, it is both convenient and conventional to recast the dynamical equations either in terms of "phase space" dynamics [68] or in terms of the Bloch vector [43,44]. As we shall use both pictures interchangeably, we now briefly introduce each in turn.
In the phase-space picture, we look at the dynamics in the space of phase difference between two sites, φ ≡ arg(α 1 /α 2 ) ∈ (−π, π], and the normalized population imbalance between two sitesz ≡ (|α 1 | 2 − |α 2 | 2 )/n, where n ≡ |α 1 | 2 + |α 2 | 2 is the total population. For two-site Hermitian models, the dynamics is fully characterized in the space of {z, φ}. In non-Hermitian models, however, the overall norm of the wavefunction can change and have significance, so the full dynamics is characterized in the space of {n,z, φ}. Equations of motion for {n,z, φ} can be obtained from Eq. 1 in a straightforward manner, and they arė We see that the time evolutions of the variablesz and φ are closed by themselves and do not depend on the total population n. As a result, we can consider the dynamics in the restricted 2D phase space {z, φ}; the dynamics of n is separately determined from the information of {z, φ} [43,44]. Most importantly, when considering the dynamics of {z, φ}, we do not need to worry about the fact that the total population n can change in time, and thus we can analyze the dynamics in {z, φ} in a way analogous to the Hermitian dimer models. Secondly, we can alternatively visualize the dynamics by means of the normalized Bloch vector, s = (s x , s y , s z ) T , whose components are defined as: subject to the condition |s| 2 = 1/4, implying that the head of the Bloch vector always lies on the surface of a sphere of radius 1/2 called the Bloch sphere [43,44]. As the total population n also changes, the full dynamics is described by the four variables {n, s x , s y , s z }. It can then be shown that [53]: where again the time dependence of the Bloch vector (s x , s y , s z ) does not depend on the total population n, and thus, upon considering the dynamics of the Bloch vector (s x , s y , s z ), we do not need to worry about the time dependence of n. We note that since a vector (s x , s y , s z ) lies on the surface of a two-dimensional sphere, the dimension of the space in which the variables {n, s x , s y , s z } move is three-dimensional, agreeing with the three-dimensional description in the phase space {n,z, φ}. Naturally, the phase-space and the Blochvector equations are fully equivalent, as can be shown by noting that φ = arctan(s y /s x ) andz = 2s z .

B. Dynamics of the linear HN dimer
Linear non-Hermitian dimer models exhibit a PTsymmetry breaking transition, in which the eigenstates of the Hamiltonian coaelesce at an exceptional point [1,2,69] as the strength of the non-Hermitian terms is increased. The energy eigenvalues, E ± , of the Hamiltonian are [2,69]: The energies are purely real when J 2 + ∆ 2 > δJ 2 , corresponding to a weakly non-Hermitian regime called the PT -symmetric regime. When J 2 + ∆ 2 = δJ 2 the two eigenvalues become degenerate, at which point the eigenvectors also coalesce, yielding the exceptional point in the parameter space. When J 2 + ∆ 2 < δJ 2 , the energies acquire a nonzero imaginary part, which is called the PTbroken regime. The appearance of the imaginary part of eigenenergies indicates that the norm of the wavefunction will not be conserved over time, and the population exponentially grows or decays [28,53,66,67,70]. In the following, we first review the PT -symmetric region J 2 + ∆ 2 > δJ 2 where the dynamics is described by oscillation around two fixed points. We then discuss the PT -broken regime J 2 + ∆ 2 < δJ 2 where the fixed points turn into a source and sink of dynamics.

PT -symmetric Regime
In the PT -symmetric regime, J 2 + ∆ 2 > δJ 2 , where the eigenvalues of the Hamiltonian are real, the dynamics consist of Rabi oscillations between the two sites, which can be seen as closed orbits in phase-space and on the Bloch sphere [43,44]. The population can also be biased toward one of the two sites; the emergence of population imbalance comes from non-Hermiticity, and its mechanism is different from the self-trapping known in interacting Hermitian Josephson dimers.
The dynamics can be understood through the fixed points of motion, which are obtained by setting the time derivative of the variables to zero. In the non-Hermitian models that we analyze in this paper, we find it convenient and useful to look for fixed points of variables other than n, leaving the possibility for n to change in time.
In what follows, when we refer to fixed points, this refers to fixed points in the parameter space of either {z, φ} in the phase space description or {s x , s y , s z } in the Bloch sphere description without including n.
In terms of the Bloch-sphere description, the fixed points are thus obtained by settingṡ x =ṡ y =ṡ z = 0 in Eqs. 12, 13 and 14. Looking at the equation forṡ z = 0, one sees that the fixed point should satisfy either s y = 0 or s z = J/2δJ. The latter solution implies s z > 1/2 in the PT -symmetric region, J 2 + ∆ 2 > δJ 2 , which is not compatible with the condition |s| 2 = 1/4 and thus not a valid solution. Therefore, the only fixed points in the PT -symmetric region satisfy s y = 0. Solving the other equationsṡ x =ṡ y = 0, we obtain two fixed points: where Ω ≡ √ J 2 + ∆ 2 − δJ 2 > 0 in the PT -symmetric regime.
Equivalently, the fixed points can be expressed in terms of the phase-space variables by settingż =φ = 0. We find two fixed points, In the first fixed point, we should choose φ = 0 when δJ < J and φ = π when δJ > J. This discontinuous change of φ may look odd, but this is due to the fact that when J = δJ, the first fixed point becomesz = 1, which corresponds to the north pole in the Bloch sphere description. As the fixed point crosses the north pole, the phase angle φ changes discontinuously from 0 to π. Note that there is no discontinuity in the Bloch sphere description.
The dynamics around a fixed point can be understood by looking at the eigenvalues of the Jacobian matrix, J , of the fixed point [71][72][73][74], which in the phase space {z, φ}, is If the real part of the two eigenvalues are both zero, the fixed point acts as a center of motion around which the system oscillates. If the real part of the eigenvalues are both positive, the fixed point is called an unstable fixed point, and the dynamics flows away from the point. If the real part of the eigenvalues are both negative, the fixed point is called a stable fixed point, and the dynamics sinks into the point. If one of the eigenvalues has a positive real part and the other eigenvalue has a negative real part, the fixed point behaves as a saddle point of dynamics. The fixed points of Eq. 17 behave as centers for the dynamics [44,45]. This oscillatory behavior is analogous to the Rabi oscillation in coherently coupled two-level quantum systems. Such oscillatory dynamics is confirmed also experimentally, as seen in Fig. 3. Depending on the value of δJ, which is the strength of the non-Hermiticitiy, the value ofz at the fixed points varies, and it can even reach the maximum value ofz = 1. This imbalance of population is a two-site version of the non-Hermitian skin effect known in the extended Hatano-Nelson model with edges, where all the eigenstates are known to be localized on one edge [54][55][56].
These dynamical features are both numerically and experimentally confirmed, as described in Fig. 3. The left panel describes the Hermitian limit (δJ = 0) of the linear HN dimer model, whereas the central panel is in the PTsymmetric regime and the right panel is the PT -broken regime. Fixed points which serve as centers of dynamics are indicated by blue triangles, whereas the source and sink are indicated by a green diamond and an orange square, respectively. The roles of fixed points found in the dynamics obtained from numerical and experimental means are in accord with what we found above.
If we insert the fixed points into the equation for n, we seeṅ = 0, indicating that the total population does not change in time on the fixed points. This behavior is expected because the eigenvalues of the Hamiltonian in the PT -symmetric region are both real. However, if we look at the dynamics around the fixed points, the total population n changes in time. In the left panel of Fig. 4, we plot the numerically-calculated dynamics within a region of the 3D phase-space of {n,z, φ}. In the 3D phase-space, the fixed point in the phase space {z, φ} corresponds to a line perpendicular to thez − φ plane. The variation of the population during dynamical evolution even in the PT -symmetric region reflects that the eigenstates are not orthogonal, and their overlap can lead to a change of population during dynamics.

PT -broken Regime
Approaching the PT -breaking transition from the PTsymmetric region corresponds to taking Ω → 0, and the fixed points of the PT -symmetric region in Eq.(16) merge in this limit, giving rise to an exceptional point [44,45].
Beyond the PT -breaking transition, J 2 +∆ 2 < δJ 2 , we enter the PT -broken regime where the eigenvalues of the Hamiltonian are no longer real. The fixed point solution with s y = 0 from the PT -symmetric region is no longer a valid solution in the PT -broken regime. Instead, there are two fixed points at Equivalently, in the phase-space description, the fixed points are: The fixed point at φ < 0 is a stable fixed point acting as a sink of dynamics, whereas the fixed point at φ > 0 is an unstable fixed point acting as a source [44,45]. In the linear HN dimer, fixed points are merely the eigenstates of the 2-by-2 Hamiltonian; the eigenstate corresponding to the stable fixed point has the eigenvalue whose imaginary part is larger than the other eigenstate. We can physically understand the source/sink nature of these fixed point by noting that, starting from a state that is a superposition of the two eigenstates, the weight of the eigenstate corresponding to the sink will grow because of the larger imaginary part. We note that the total population n at the fixed points is no longer a constant, but rather grows in time. In this sense, these fixed points are not fixed points in the full dynamics taking n into account. As we noted earlier in this paper, we keep this terminology that the fixed points are points where variables other than n are kept constant.

V. INSTANTANEOUS HN DIMER MODEL
Now we turn to models with nonlinearities. Unlike the linear HN dimer model where the fixed points correspond to the eigenstates of the 2-by-2 Hamiltonian, obtaining the fixed points of the dynamics of models with nonlin-earity is generally nontrivial. Although the model we experimentally implement is the nonlinear HN model, this model does not allow for a simple description in terms of fixed points. However, as we shall explain, we can obtain the fixed points of the instantaneous HN model, which also contains nonlinearity and thus provides a good approximation to understand the experimentally realized nonlinear HN model. In this section, we give a theoretical description of the dynamics of the instantaneous HN dimer model, which we will use in later Sections to understand the dynamics of the nonlinear HN dimer model.
We rewrite here the equation of motion for the instantaneous HN dimer model Eq. (3): where we used the notationz = z/n we introduced in the previous section. The equations for the phase-space variables {n,z, φ} arė while the equations for the Bloch sphere variables are: n = −8ḡn s y s z (25) s x = 2∆s y + 8ḡ s x s y s z (26) s y = −2∆s x + 2(J −ḡ)s z + 8ḡ s 2 y s z (27) s z = −2Js y + 8ḡ s y s 2 z .
Just like the linear HN dimer model, time derivatives of the variablesz, φ, s x , s y , and s z have no n dependence, which implies that we can consider the dynamics in the phase space {z, φ} and on the Bloch sphere {s x , s y , s z } independent from how n depends on time.
Nonetheless, the instantaneous HN dimer model is not a linear model; indeed the dynamics hasz dependence in the 2-by-2 nonlinear Hamiltonian of Eq. (3). The situation should be contrasted with what happens for the nonlinear HN dimer model, whose equations of motions are Eqs. (40)(41)(42)(43)(44)(45)(46) in the next section. The nonlinear HN dimer model has an explicit n dependence in the equations of motion forz, φ, s x , s y , and s z , and a closed description in terms of either {z, φ} or {s x , s y , s z } is not allowed, making the analysis more difficult. The dynamics of the instantaneous HN dimer model we explore in this section serves as the basis to understand the dynamics of the nonlinear HN dimer model in the next section.

A. Dynamics of the instantaneous HN dimer model
As before, we first look for the fixed points of motion. Fixed points, as before, are points determined byż =φ = 0 for the phase-space dynamics andṡ x =ṡ y =ṡ z = 0 for the Bloch sphere dynamics, without imposingṅ = 0. The fixed points we discuss are, therefore, fixed points in the restricted space in which we do not look at the time evolution of n.
From the conditionṡ z = 0, we obtain s y = 0 or s 2 z = J/(4ḡ). Let us first examine the case s y = 0. From the other two equations,ṡ x =ṡ y = 0, we obtain the following two fixed point solutions A striking feature of the instantaneous HN dimer model is that these fixed points are always valid fixed points irrespective of the values of J, δJ, and ∆, which is to be contrasted with the linear HN dimer model where the fixed point s y = 0 was valid only in the PT -symmetric region. These fixed points in the phase-space description are The sign of the fixed points ofz depends on sign(J −ḡ); this is because when J =ḡ, the fixed point reaches the north pole in the Bloch sphere description and thus φ changes between 0 and π around J =ḡ. Note that there is no discontinuous change in the Bloch sphere representation. Now we look for the fixed point s y ̸ = 0. Combining s 2 z = J/(4ḡ) withṡ x =ṡ y = 0, we obtain the following additional four fixed points (31) where the signs of the first and the second terms can be chosen independently, and we defined These additional fixed points with s y ̸ = 0 are valid fixed points only whenḡ − g T > 0. Thus, the instantaneous model behaves qualitatively differently depending on whetherḡ < g T orḡ > g T , where the former has only two fixed points but the latter has six fixed points. We will refer to the caseḡ < g T as the weak interaction regime andḡ > g T as the strong interaction regime. These fixed points with s y ̸ = 0 in the phase-space description have the following expressions: where the sign ofz and the sign inside arccos should be chosen to be opposite, but the sign in front of arccos is independent of the other signs. We now consider the weak interaction regimeḡ < g T and the strong interaction regimeḡ > g T in turn to inspect the nature of fixed points and dynamics around them.

Weak interaction regime
In the weak interaction regimeḡ < g T , the instantaneous HN dimer model has two fixed points given by Eq. (29). These fixed points satisfyṅ = 0, indicating that the total number does not change in time. Thus, the weak interaction regime is a direct analog of the PTsymmetric regime in the linear HN dimer model.
We find that the eigenvalues of the Jacobian of the two fixed points in the phase-space description are λ ± = ±2 J(ḡ − g T ) for both fixed points. Sinceḡ < g T in the weak interaction regime, λ ± are complex conjugate pairs, indicating that these fixed points serve as the center of oscillation in the phase space of {z, φ}, analogous to the fixed points in the PT -symmetric regime in the linear HN dimer model. In Fig. 5, we plot numerically obtained dynamics of the weak interaction regime. All orbits are closed as expected.

Strong interaction regime
The strong interaction regimeḡ > g T of the instantaneous HN dimer model behaves differently from the PT -broken regime of the linear HN dimer model. First of all, the two fixed points present in the weak interaction regimeḡ < g T remain as fixed points in the strong interaction regime. The eigenvalues of the Jacobian of these two fixed points are λ ± = ±2 J(ḡ − g T ), which are now both real and have opposite signs; this indicates that the fixed points now act as saddle points of dynamics in the {z, φ} phase space.
In addition to the two fixed points inherited from the weak interaction regime, there are four additional fixed points in the strong interaction regime. Calculating the eigenvalues of the Jacobian of these four fixed points, we can group them into two categories.
Among the four additional fixed points, the eigenvalues of the Jacobian of the following two fixed points are in both cases which are both positive in the strong interaction regimē g > g T . Therefore, these two fixed points are unstable fixed points which act as sources of dynamics.
On the other hand, the eigenvalues of the Jacobian of the other two fixed points are both which are both negative, indicating that these two fixed points are stable fixed points that act as sinks of dynamics. This behavior is confirmed by the numerically obtained dynamics as plotted in Fig. 6.
Calculating the time dependence of the total population at these four fixed points, we finḋ where the positive (negative) sign corresponds to the two stable (unstable) fixed points. This indicates that the total population of the stable fixed points exponentially increases, whereas that of the unstable fixed points exponentially decreases; such behavior is consistent with these stable points acting as sinks and sources of dynamics. For all four of these fixed points, the population imbalance is |z| = J/ḡ, which indicates that |z| → 0 as g → ∞. The non-Hermitian localization on one of the two sites becomes smaller and smaller as the interaction g becomes larger, which is in stark contrast to the selftrapping phenomenon known in the two-site nonlinear Josephson model where the localization becomes stronger as the interaction becomes larger.

Transition between the weak and strong regimes
We have just seen that, as the interactionḡ increases and crosses g T , the number of fixed points changes from two to six. We now examine this transition.
We first note that such an increase of the number of fixed points beyond two is not possible with any linear dimer model, in which the fixed points are determined by Plots in the phase space and the Bloch sphere are for different initial conditions (marked with open markers); the evolution of the population over time for each initial condition is plotted in a-I and b-I using the corresponding marker. The saddle (blue triangles), stable (orange squares) and unstable (green diamonds) fixed points are highlighted. Note the trajectories converge towards the stable fixed points, after being slightly bent by the presence of saddle points. As expected, the population diverges quickly as a stable point is approached.
(at most) two eigenstates of the two-by-two Hamiltonian. Therefore, the appearance of six fixed points in the strong regime is an intrinsically nonlinear phenomenon.
Approaching the transition point from the weak regime, each of the two fixed points split into three as one crosses the transition pointḡ = g T . During this process, a fixed point on the weak interaction side, which is a center of dynamics, turns into three fixed points which are a saddle point, a stable fixed point, and an unstable fixed point. This process is consistent with the Poincaré-Hopf index theorem as we shall explain. On a two-dimensional parameter space, such as the phase space {z, φ} and the surface of the Bloch sphere, the tangent vectors of the dynamics define a vector field. Such a vector field can have singularities, corresponding to the fixed points. For each of these singularities, a topological index called the Poincaré index can be defined, which assigns the value of −1 for saddle points and +1 for centers, stable, and unstable fixed points. The index theorem states that the sum of the Poincaré indices on the two-dimensional parameter space should be equal to the Euler characteristics of the parameter space [44,75]. The Euler characteristics of our parameter space, which is a two-dimensional sphere as evident from the Bloch sphere description, is +2. In the weak interaction regime, we have two centers as singularities, and thus the sum of the Poincaré indices is +2, which is equal to the Euler characteristics as ex-pected. As one crosses the transition point, a center, which has the Poincaré index of +1, turns into a saddle point, a stable point, and an unstable point, whose Poincaré indices are −1, +1, and +1, respectively, conserving the sum of the Poincaré indices. Thus the strong interaction regime also satisfies the index theorem.

VI. NONLINEAR HN DIMER MODEL
With the understanding of the instantaneous HN dimer model, we can now understand the dynamics of the nonlinear HN dimer model, which is the model experimentally implemented.
We rewrite the equations of motion for the nonlinear HN dimer model Eq. (2), The dynamical equations in the phase-space arė and in terms of the Bloch sphere variables, the dynamics obeyṅ = −8gn 2 s y s z , (43) s x = 2∆s y + 8gn s x s y s z , We see that the time dependence ofz, φ and s x , s y , s z depend explicitly on n, implying that we can no longer use the two-dimensional phase space {z, φ} and the surface of the Bloch sphere {s x , s y , s z } as a parameter space within which the equations of motion are closed. Instead, we should consider the time dependence of n together with the other variables to understand the dynamics. This explicit n dependence in the equations of motion makes it difficult to analytically approach the dynamics of the nonlinear HN dimer model when the interaction is strong. However, as we shall see, the theory of the instantaneous HN dimer model can provide a good qualitative understanding of what happens in the nonlinear HN dimer model.

A. Dynamics of the nonlinear HN dimer model
We first look for fixed points of the dynamics. Since replacingḡ in the instantaneous HN dimer model by gn recovers the equations of motion of the nonlinear HN dimer model, the two fixed points in the weak interaction regime and the six fixed points in the strong interaction regime of the instantaneous HN dimer model still satisfẏ z =φ = 0 andṡ x =ṡ y =ṡ z = 0. However, these points are not truly fixed anymore in the parameter space because n can depend on time, and the time dependence of n itself affects the position of the fixed points in the parameter space.
We first note that the two fixed points satisfying s y = 0 are still fixed in the parameter space even in the nonlinear HN dimer model because they obeyṅ = 0, namely n is time independent. The additional four fixed points are not fixed anymore in the parameter space because n changes in time. The transition between the weak and strong interaction regimes is clear in the case of the instantaneous HN dimer model, given byḡ = g T . On the other hand, in the case of the nonlinear HN dimer model, there appears also a distinctive intermediate interaction regime between the weak and the strong regimes as we shall explain now.
In the nonlinear HN dimer model, one measure of the interaction strength is gn 0 , where n 0 is the total population at the initial time. Replacingḡ by gn 0 , the two fixed points with s y = 0 of the instantaneous HN dimer model are still the fixed points in the nonlinear HN dimer model. If gn 0 is small enough compared to g T , the dynamics we found for the weak interaction regime of the instantaneous HN dimer model applies also to the nonlinear HN dimer model. However, one should remember that, although the total population n does not change in time exactly at the fixed points, the total population does change during the periodic Rabi oscillation around the fixed points. This implies that even though the system initially satisfies gn 0 < g T , the total population changes and at some time t we may enter the regime with gn(t) > g T at which the dynamics should be compared to the strong interaction regime of the instantaneous HN dimer model.
We refer to the regime in which gn(t) < g T holds for any time t ≥ 0 as the weak interaction regime of the nonlinear HN dimer model. In the weak interaction regime, all the dynamics are described by orbital motion around the two fixed points which serve as the centers of dynamics.
We refer to the regime in which both gn(t) < g T and gn(t) > g T happen at some time t ≥ 0 during the evolution as the intermediate interaction regime. In this intermediate interaction regime, which is a unique feature of the nonlinear HN dimer model, there simultaneously exist two types of orbits: one is a closed orbit similar to the weak interaction regime, and the other is a diverging orbit which is reminiscent of the strong interaction regime of the instantaneous HN dimer model.
We finally refer to the regime in which gn(t) > g T holds for any time t ≥ 0 as the strong interaction regime, where the dynamics diverges as in the strong interaction regime of the instantaneous HN dimer model. A crucial difference between the strong interaction regime of the instantaneous HN dimer model and the nonlinear dimer model is that, in the latter, the source and the sink of the dynamics are no longer fixed points in the space {z, φ} or {s x , s y , s z } because of the change of the total population n.
We now examine these three regimes experimentally, and we compare them to numerical simulations. Figure 7 shows the dynamics of the nonlinear HN dimer model, experimentally and numerically obtained, for different values of gn 0 . Firstly, panels with gn 0 = 0 and gn 0 = 0.4 J correspond to the weak interaction regime. We observe Rabi-like oscillations, similar to the ones found in the linear and instantaneous HN dimer models. Secondly, the panel with gn 0 = 1.1J corresponds to the strong interaction regime, where the population of the trajectories diverges over time towards either positive or negative population imbalances, depending on the initial conditions. Finally, in the intermediate interaction regime at gn 0 = 0.7J the two behaviors co-exist, i.e., stable orbits and unstable trajectories are possible, depending on the initial conditions used, i.e. the initial phase difference between the two sites. This co-existence of both regimes is the unique feature of the nonlinear HN dimer model.

B. Experimental results and Numerical simulations
We can further see the unique feature of the nonlinear HN dimer model if we examine the dynamics on the normalized Bloch sphere, as shown in the central panel of Fig. 8. We observe that the trajectories bend in arcs on the sphere, before eventually converging towards points on the equator. This behavior arises because instantaneously the trajectory is attracted towards the stable fixed points of the instantaneous HN dimer model for that particular value of gn →ḡ. However, as n keeps growing, the coordinates of the fixed points in the instantaneous HN dimer model also keep changing (see, e.g., Fig. 6), so that the trajectories appear to effectively "chase down" the stable points by following these arcs. Indeed, in the limit that n → ∞, the coordinates on the Bloch sphere of the stable fixed points in the instantaneous HN dimer model become: corresponding to the points of convergence on the equator of the normalized Bloch sphere in the experimental model, e.g., as can be seen in the central panel of Fig. 8. Physically, this corresponds to the modes changing from being localized primarily on one of the two sites at small g (and hence small n), to being an equal superposition of the two sites asḡ → ∞, as the non-reciprocal coupling term dominates over all other terms in the Hamiltonian.
Of particular note in Fig. 8, two trajectories starting close to each other (black points and red points in the top panels) actually behave very differently: one describes a closed loop around the center (full blue triangle) and its population remains finite, although oscillating, whereas the other about half-way through the loop, diverts and ends up in the basin of attraction of the stable points (full orange squares), where the population diverges exponentially. As in the strong interaction case, note that the trajectories are bent in arcs to "follow" the always changing position of the fixed point in the corresponding instantaneous model.
Another helpful way to visualize this physics is shown in Fig. 9, where we plot arrows corresponding to the tangent vectors of dynamics on each point of the three- corresponding markers and colors) shows an evident decrease of population in the neighborhood of the unstable points (green diamonds). As the system evolves, it is then attracted towards the stable points (orange squares) where n grows exponentially.
dimensional phase space. The panels show the change in behavior of the fixed points on slices of fixed n: as long as for some time t 1 , gn(t 1 ) < g T , the former are centers (left panel), whereas as soon as for some later time t 2 the condition gn(t 2 ) > g T is met, the latter behave like saddles (see the middle and right panels).
In the intermediate regime, in the case of initial conditions close to φ = 0, the system is pushed by the equations in a closed orbit around the center (top left panel in Fig. 9) and the population oscillates, as shown in Fig. 5. If the initial conditions are not sufficiently close to the fixed point of the blue up-triangle in the left panel, while moving around the center, the trajectory ends up close enough to the basin of repulsion of unstable points (green diamonds) or to the basin of attraction of the stable points (orange squares).

VII. CONCLUSIONS AND OUTLOOK
In this paper, we have investigated both theoretically and experimentally the dynamics of two-site models with hopping asymmetry. In a linear model where the hopping asymmetry is externally fixed, we experimentally observed the transition from PT -symmetric to PT -broken regimes as the hopping asymmetry is increased. While all the orbits are closed in the PT -symmetric regime, the population diverges in the PT -broken regime due to the non-Hermiticity from the hopping asymmetry. In a nonlinear model where the hopping asymmetry is dynamically induced by population imbalance between the two sites, we experimentally observed three different regimes in behavior, depending on the initial coupling strength.
In the weak and strong regimes, we observe stable population oscillations and exponential growth/decay of the population, respectively, similar to the behavior in the linear model described above. However, in the inter-mediate regime, we observe a coexistence of these dynamics, meaning that we can tune from stable oscillations to divergent behavior by simply varying the initial phase-difference between the two sites. As we explain, all three different regimes can be understood by studying the emergent fixed points of a closely-related nonlinear model in which the non-reciprocal hopping depends on the normalized population imbalance between the two sites.
In the future, this work will pave the way towards the further exploration of non-Hermitian dynamics in more exotic systems. As demonstrated here and in our previous works [20], this mechanical platform can be used to simulate a wide-variety of lattice models, which would not be easy to realize in other systems. Going further, it will be interesting to explore, for example, the addition of other mean-field nonlinear effects [76] as well as extensions to larger systems, e.g., such as non-Hermitian three-site trimer models [77][78][79] or large lattices with many sites [80], where the interplay of gain and loss with artificial gauge fields and topological phenomena can also be explored [81,82]. and RGF\R1\180071. This work was also supported by the BRIDGE Seed Fund for collaboration between the University of Birmingham and the University of Illinois at Urbana-Champaign.