# Electron-hole asymmetry of quantum collective excitations in high-$T_{c}$
copper oxides

Maciej Fidrysiak<EMAIL_ADDRESS>Institute of Theoretical Physics,
Jagiellonian University, ul. Łojasiewicza 11, 30-348 Kraków, Poland

###### Abstract

We carry out a systematic study of collective spin- and charge excitations for
the canonical single-band Hubbard, $t$-$J$-$U$, and $t$-$J$ models of high-
temperature copper-oxide superconductors, both on electron- and hole-doped
side of the phase diagram. Recently developed variational wave function
approach, combined with the expansion in inverse number of fermionic flavors,
is employed. All three models exhibit a substantial electron-hole asymmetry of
magnetic excitations, with a robust paramagnon emerging for hole-doping, in
agreement with available resonant inelastic $x$-ray scattering data for the
cuprates. The $t$-$J$ model yields additional high-energy peak in the magnetic
spectrum that is not unambiguously identified in spectroscopy. For all
considered Hamiltonians, the dynamical charge susceptibility contains a
coherent mode for both hole- and electron doping, with overall bandwidth
renormalization controlled by the on-site Coulomb repulsion. Away from the
strong-coupling limit, the antiferromagnetic ordering tendency is more
pronounced on electron-doped side of the phase diagram.

##  I Introduction

High temperature superconductivity in layered copper oxides may be induced by
either hole (h) or electron (e) doping of the parent antiferromagnetic (AF)
insulating state Spałek  _et al._ (2022). Nonetheless, the equilibrium phase
diagrams of e- and h-doped cuprates exhibit substantial asymmetry, whose one
of the most evident manifestations is broader regime of the AF-phase
appearance on the e-doped side Armitage  _et al._ (2010). Thanks to
developments in spectroscopic techniques, particularly resonant inelastic
$x$-ray scattering (RIXS), this static picture has been recently supplemented
with a detailed account of non-equilibrium properties, including collective
excitation spectra Le Tacon  _et al._ (2011); Dean  _et al._ (2013a, b, c); Le
Tacon  _et al._ (2013); Ishii  _et al._ (2014); Lee  _et al._ (2014); Dean
_et al._ (2014); Guarise  _et al._ (2014); Jia  _et al._ (2014); Wakimoto  _et
al._ (2015); Minola  _et al._ (2015); Peng  _et al._ (2015); Ellis  _et al._
(2015); Huang  _et al._ (2016); Gretarsson  _et al._ (2016); Minola  _et al._
(2017); Ivashko  _et al._ (2017); Meyers  _et al._ (2017); Chaix  _et al._
(2018); Robarts  _et al._ (2019a); Fumagalli  _et al._ (2019); Peng  _et al._
(2018); Robarts  _et al._ (2021); Ishii  _et al._ (2017); Nag  _et al._
(2020); Singh  _et al._ (2022); Hepting  _et al._ (2018); Ishii  _et al._
(2019); Lin  _et al._ (2020). The low-energy dynamics of parent AF state is
dominated by spin-waves of localized spins, evolving with doping into
broadened magnetic excitations of the paramagnetic state (paramagnons).
Counter-intuitively, the paramagnons appear to be _more robust_ on the h-doped
side, where AF phase is rapidly suppressed away from half-filling. RIXS has
been also utilized to probe collective charge (plasmon) excitations that are
observed in both hole- Ishii  _et al._ (2017); Nag  _et al._ (2020); Singh
_et al._ (2022) and electron-doped Hepting  _et al._ (2018); Ishii  _et al._
(2019); Lin  _et al._ (2020) cuprates. Since magnetic and charge fluctuations
have been considered as candidates for the pairing mechanism in high-$T_{c}$
superconductors (see, e.g., Chubukov  _et al._ (2003); Scalapino (2012); Le
Tacon  _et al._ (2011); Yu  _et al._ (2020)), identification of a microscopic
origin of their distinct degree of e-h asymmetry is now in demand.

Parent compounds of high-$T_{c}$ copper-oxides are not canonical Mott
insulating systems, since they may be regarded as charge transfer insulators
within the Zaanen-Sawatzky-Allen scheme Zaanen  _et al._ (1985). In effect,
the e-h asymmetry manifests itself directly at the atomic-orbital level.
Nominally, the doped electrons position on the copper $3d_{x^{2}-y^{2}}$
states, while holes locate predominantly on $2p_{\sigma}$ oxygen orbitals.
Yet, there is a systematic variation of $d$-$p$-electron redistribution
between different cuprate families Rybicki  _et al._ (2016). Whereas minimally
the $d$-$p$ model of $\mathrm{CuO_{2}}$ plane is needed to address those
orbital-selective effects Zegrodnik  _et al._ (2019, 2020, 2021); Li  _et al._
(2021), doping dependence of the principal experimental equilibrium properties
and their e-h asymmetry is reproduced already within the effective single-band
picture Spałek  _et al._ (2017); Zegrodnik and Spałek (2017); Fidrysiak  _et
al._ (2018); Wysokiński and Fabrizio (2017a). The latter model is also
appropriate for studies of the correlated state dynamics following interaction
quench Schiró and Fabrizio (2010); Wysokiński and Fabrizio (2017b). The
single-band Hamiltonians remain thus a common starting point for describing
the physics of high-$T_{c}$ cuprates and related materials.

Here we carry out a study of the e-h asymmetry of collective magnetic and
charge excitations for the canonical Hubbard, $t$-$J$-$U$, and $t$-$J$ models
of high-$T_{c}$ copper-oxide superconductors, which allows us to single out
the microscopic parameters controlling relevant aspects of their many-particle
dynamics. The recently introduced concept of effective exchange interaction
Fidrysiak  _et al._ (2021) is invoked to relate all three Hamiltonians in
proper limits so that their collective excitations can be directly compared
and analyzed on equal footing. In order to account for the effects of strong
electronic correlations, we supplement Variational Wave Function (VWF)
approach with the expansion in inverse number of fermionic flavors
($1/\mathcal{N}_{f}$), which results in VWF+$1/\mathcal{N}_{f}$ scheme
Fidrysiak and Spałek (2021a). The latter has been demonstrated to be effective
in semi-quantitative analysis of observed spin- and charge dynamics of h-doped
cuprates, and to compare favorably with other computational techniques,
including determinant quantum Monte-Carlo Fidrysiak and Spałek (2020, 2021b);
Fidrysiak  _et al._ (2021); Spałek  _et al._ (2022).

For the cases of Hubbard and $t$-$J$-$U$ models up to moderately large value
of on-site interaction, we find that intense dispersive peak in magnetic
spectrum persists along the anti-nodal ($\Gamma$-$X$) Brillouin-zone direction
down to heavily h-overdoped regime, whereas the paramagnons along the nodal
($\Gamma$-$M$) line are rapidly suppressed. On the other hand, spin
excitations are shifted to larger energies with e-doping and become less
coherent away from half-filling. These findings indicate a substantial e-h
asymmetry at the dynamical level and are in agreement with available RIXS data
for the cuprates. Moreover, we find that the magnetic spectra calculated in
the strong-coupling limit ($t$-$J$ model) differ significantly from those
obtained for the Hubbard and $t$-$J$-$U$ Hamiltonians, and exhibit
characteristics that are not observed experimentally. This might suggest that
the $t$-$J$ model overestimates local electronic correlations, as has been
also noted within former theoretical survey of equilibrium and single-particle
quantities Spałek  _et al._ (2022).

For completeness, we also address the doping evolution of collective charge
excitations that are demonstrated to persist on both h- and e-doped sides of
the phase diagram. The dynamical charge response unambiguously separates into
the incoherent continuum part and a coherent charge mode. With increasing on-
site Coulomb repulsion, $U$, the collective mode energy undergoes a systematic
downward renormalization. Thus, the electronic band-narrowing effects are
instrumental for the charge sector, whereas magnetic part remains primarily
sensitive to effective exchange interaction, $J_{\mathrm{eff}}$. Remarkably,
within the $t$-$J$-$U$ model, the parameters $U$ and $J_{\mathrm{eff}}$ remain
independent and may be simultaneously tuned to match experiment. This is not
possible within the Hubbard and $t$-$J$ models. Such flexibility is needed to
account for non-trivial superexchange pathways via oxygen in a charge-transfer
insulator that alter the correspondence between kinetic exchange and $U$
Spałek  _et al._ (2017).

Finally, we carry out analysis of the leading ordering instabilities of the
paramagnetic state both against spin- and charge fluctuations. For the Hubbard
and $t$-$J$-$U$ models, we find that the static magnetic susceptibility at the
$M$ Brillouin-zone point is enhanced for e-doping relative to the h-doped
side, implying stronger tendency towards antiferromagnetism. This reflects the
experimentally observed e-h asymmetry of the equilibrium AF order.

##  II Model and method

We employ a generic square-lattice $t$-$J$-$U$ model of the copper-oxygen
plane, given by the Hamiltonian

$\displaystyle\hat{\mathcal{H}}=\sum_{i\neq
j}t_{ij}\hat{a}^{\dagger}_{i\sigma}\hat{a}_{j\sigma}+J\sum_{\langle
i,j\rangle}\hat{\mathbf{S}}_{i}\cdot\hat{\mathbf{S}}_{j}+U\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$
(1)

where $\hat{a}_{i\sigma}$ ($\hat{a}^{\dagger}_{i\sigma}$) are annihilation
(creation) operators of spin-$\sigma$ electrons on site $i$,
$\hat{n}_{i\sigma}\equiv\hat{a}^{\dagger}_{i\sigma}\hat{a}_{i\sigma}$, and
$\hat{\mathbf{S}}_{i}\equiv(\hat{S}^{x}_{i},\hat{S}^{y}_{i},\hat{S}^{z}_{i})$
denotes local spin operator. Out of the hopping integrals, $t_{ij}$, we retain
only those connecting nearest- and next-nearest neighbors,
$t\equiv-0.35\,\mathrm{eV}$ and $t^{\prime}\equiv 0.25|t|$, respectively. It
should be remarked that, within the $t$-$J$-$U$ model, the degree of e-h
asymmetry is controlled exclusively by the magnitude of next-nearest neighbor
hopping, $t^{\prime}$. Up to a trivial energy and chemical potential shift,
the remaining terms remain unchanged after application of the e-h
transformation in the usual form

$\displaystyle a^{\dagger}_{i\sigma}\rightarrow
a_{i\sigma}\cdot(-1)^{i},a_{i\sigma}\rightarrow
a_{i\sigma}^{\dagger}\cdot(-1)^{i},$ (2)

with $(-1)^{i}$ being phase factor alternating between neighboring sites. The
electronic interactions are governed by the magnitude of on-site Coulomb
repulsion, $U$, and by nearest-neighbor AF exchange integral, $J$.

The $t$-$J$-$U$ Hamiltonian (1) may be regarded as a generalization of the
Hubbard- and $t$-$J$ models, and encompasses both of them as particular cases.
The Hubbard model is obtained for $U>0$ and $J=0$, whereas the $t$-$J$ model
for $U=\infty$ and $J>0$. In order to analyze all three Hamiltonians on equal
footing, we refer to the concept of _effective exchange interaction_
$J_{\mathrm{eff}}\equiv J+\frac{4t^{2}}{U}$ that combines direct Heisenberg-
type interaction with second-order kinetic exchange Fidrysiak  _et al._
(2021). Hereafter its value is set to $J_{\mathrm{eff}}\equiv\frac{2}{3}|t|$
which has been used previously to study collective dynamics in the
paramagnetic state of h-doped cuprates Fidrysiak and Spałek (2021b). The
selected representative parameter values are: $U=6|t|$, $J=0$ (Hubbard model),
$U=12|t|$, $J=\frac{1}{3}|t|$ ($t$-$J$-$U$ model), and $U=\infty$,
$J=\frac{2}{3}|t|$ ($t$-$J$ model).

Figure 1: Calculated VWF+$1/\mathcal{N}_{f}$ collective excitation spectra for
the Hubbard model ($t=-0.35\,\mathrm{eV}$, $t^{\prime}=0.25|t|$, $U=6|t|$,
$J=0$) along the high-symmetry $X$-$\Gamma$-$M$ Brillouin-zone contour.
Imaginary parts of respective dynamical susceptibilities are represented as a
color map, with blue- and white colors corresponding to low- and high
intensity, respectively. The color scale is used consistently, so the
intensities can be compared. Panels (a)-(d) and (e)-(h) show the magnetic
excitations for the hole- an electron-doped system, respectively. Panels
(i)-(p) illustrate the charge response, calculated for the same model
parameters. The h- and e-doping levels ($\delta_{h}$ and $\delta_{e}$),
detailed inside the plot, are measured relative to half-filling.

We employ VWF+$1/\mathcal{N}_{f}$ approach in the local-diagrammatic variant
($\mathrm{LD}_{f}+1/\mathcal{N}_{f}$ in the notation of Ref. Fidrysiak and
Spałek (2021a)) that allows us to study collective excitations around the
strongly-correlated ground state. The saddle-point reference point is
constructed based on the variational wave function

$\displaystyle|\Psi_{\mathrm{var}}\rangle\equiv\frac{\hat{P}|\Psi_{0}\rangle}{{||\hat{P}|\Psi_{0}\rangle||}},$
(3)

where $|\Psi_{0}\rangle$ represents a Slater determinant, and $\hat{P}$ is an
operator introducing correlations into the trial state (the so-called
correlator). The denominator in Eq. (3) is needed for normalization, since
$\hat{P}$ is not unitary in general. We take
$\hat{P}\equiv\prod_{i}\hat{P}_{i}$ with

$\displaystyle\hat{P}_{i}\equiv\lambda^{0}_{i}|0\rangle_{i}{}_{i}\langle
0|+\sum_{\sigma\sigma^{\prime}}\lambda^{\sigma\sigma^{\prime}}_{i}|\sigma\rangle_{i}{}_{i}\langle\sigma^{\prime}|+\lambda^{d}_{i}|d\rangle_{i}{}_{i}\langle
d|,$ (4)

where states $|0\rangle_{i}$, $|{\uparrow}\rangle_{i}$,
$|{\downarrow}\rangle_{i}$,
$|d\rangle_{i}\equiv|{\uparrow\downarrow}\rangle_{i}$ span the local Hilbert
space for the lattice site $i$, and
$\sigma,\sigma^{\prime}=\uparrow,\downarrow$ enumerate spin configurations. At
zero temperature, both the variational parameters
$\boldsymbol{\lambda}=\bigcup_{i}\\{\lambda^{0}_{i},\lambda^{\uparrow\uparrow}_{i},\lambda^{\uparrow\downarrow}_{i},\lambda^{\downarrow\uparrow}_{i},\lambda^{\downarrow\downarrow}_{i},\lambda^{d}_{i}\\}$
and wave function $|\Psi_{0}\rangle$ are determined by optimization of the
energy functional
$E_{\mathrm{var}}(|\Psi_{0}\rangle,\boldsymbol{\lambda})\equiv\langle\Psi_{\mathrm{var}}|\hat{\mathcal{H}}|\Psi_{\mathrm{var}}\rangle$,
subjected to the requirement of fixed total particle number and other
constraints. Here we actually employ a thermal generalization of this
variational scheme based on free energy functional Fidrysiak and Spałek
(2021a), with temperature set to $k_{B}T=0.42|t|$ in order to stay clear of
density-wave orders. Stability of the paramagnetic state for all considered
parameter sets is demonstrated in Appendix A, where we also identify leading
ordering tendencies. It should be remarked that the transformation (2)
preserves the variational space spanned by the parameters of the correlator
(4). The present methodology is thus unbiased and suitable for the analysis of
intrinsic e-h asymmetry of collective excitations, as detailed in Appendix B.

In essence, VWF+$1/\mathcal{N}_{f}$ approach relies on promoting all the
variational parameters to imaginary-time dynamical quantities, extending the
number of fermionic families from one to $\mathcal{N}_{f}$, and subsequently
expanding the variationally determined free energy functional in the powers of
$1/\mathcal{N}_{f}$. We emphasize that, even in the employed here
large-$\mathcal{N}_{f}$ approximation, the residual interactions between the
Landau quasiparticles are incorporated. The latter are instrumental for a
proper description of collective paramagnon and plasmon excitations in the
strong-coupling regime.

Finally, we note that the strong-coupling limit ($U=\infty$, $J>0$) is
described within our formalism on the same footing as the finite-$U$ case. The
exclusion of doubly occupied sites at the dynamical level is achieved by
imposing the local constraint $\lambda^{d}_{i}(\tau)\equiv 0$, where $\tau$
denotes imaginary time. Since e-doped side of the phase diagram is
inaccessible within the low-energy sector of the $t$-$J$-model, in this regime
we actually carry out the analysis for the h-doped system with inverted sign
of next-nearest hopping, $t^{\prime}$, as dictated by e-h transformation (2).

We calculate the (imaginary-time) dynamical susceptibilities in the form
$\chi_{s}(\tau,\mathbf{k})\equiv\langle\mathcal{T}_{\tau}\hat{{S}}_{\mathbf{k}}^{z}(\tau)\hat{{S}}_{-\mathbf{k}}^{z}\rangle$
and
$\chi_{c}(\tau,\mathbf{k})\equiv\langle\mathcal{T}_{\tau}\hat{n}_{\mathbf{k}}(\tau)\hat{n}_{-\mathbf{k}}\rangle$.
The results have been obtained for a $300\times 300$ square lattice with
periodic boundary conditions, and analytic continuation of the Fourier
transformed susceptibilities, $\chi_{s/c}(i\omega_{n},\mathbf{k})$, is
performed as $i\omega_{n}\rightarrow\omega+0.02|t|$, where $\omega_{n}$ are
bosonic Matsubara frequencies.

##  III Results

Figure 2: Summary of calculated VWF+$1/\mathcal{N}_{f}$ collective excitation
spectra for the $t$-$J$-$U$ model ($t=-0.35\,\mathrm{eV}$,
$t^{\prime}=0.25|t|$, $U=12|t|$, $J=\frac{1}{3}|t|$). The plot is arranged in
the same way as Fig. 1 and the color scale is consistent among panels.

###  III.1 Hubbard model

We start the discussion of e-h asymmetry of collective excitations with the
Hubbard model case ($U=6|t|$, $J=0$). In Fig. 1, the calculated
VWF+$1/\mathcal{N}_{f}$ dynamical spin [panels (a)-(d)] and charge [panels
(e)-(h)] susceptibilities for hole (electron) doping levels
$\delta_{h}(\delta_{e})=0.1,0.2,0.3,0.4$ are displayed along the high symmetry
$X$-$\Gamma$-$M$ Brillouin-zone contour. The color map represents
susceptibility magnitude with blue- and white colors mapping to low- and high-
intensities, respectively. Relevant model parameters are given inside the
panels.

For h-doping (a)-(e), an intense and disperse signal emerges throughout entire
phase diagram along the anti-nodal ($\Gamma$-$X$) direction. This peak is
interpreted as robust paramagnon in high-$T_{c}$ literature. On the other
hand, the spectrum remains highly incoherent along the nodal ($\Gamma$-$M$)
line. Such a directional anisotropy of spin excitations is in agreement with
RIXS data collected for hole-doped copper oxides Robarts  _et al._ (2019b).

We now proceed to comparison of the h- and e-side of the phase diagram. As
follows form Fig. 1(e)-(h), a clear resonant peak appears for e-doping only
along the $\Gamma$-$X$ direction and for $\delta_{e}=0.1$ [panel (e)]. Yet, it
is less intense than the corresponding feature on the h-side [panel (a)]. With
increasing electron concentration, the magnetic spectral weight shifts to
larger energies and rapidly becomes incoherent. Such a hardening effect has
been reported experimentally for e-doped cuprates Lee  _et al._ (2014); Ishii
_et al._ (2014). The substantial e-h asymmetry of spin excitations is thus
reproduced already within the single-band Hubbard model.

For completeness, in Fig. 1(i)-(p), we also compare the h- and e-doping
evolution of charge excitations. The spectrum unambiguously separates into the
incoherent continuum part and sharp disperse collective mode. The latter
remains gapless at the $\Gamma$ point, since the long-range Coulomb
interactions are not included in the Hamiltonian (1). Contrary to the spin
part, the charge modes emerge for both e- and h-doping, which is also backed
by RIXS experiments Ishii  _et al._ (2017); Nag  _et al._ (2020); Singh  _et
al._ (2022); Hepting  _et al._ (2018); Ishii  _et al._ (2019); Lin  _et al._
(2020).

Figure 3: Summary of calculated VWF+$1/\mathcal{N}_{f}$ collective excitation
spectra for the strong-coupling ($t$-$J$ model) limit ($t=-0.35\,\mathrm{eV}$,
$t^{\prime}=0.25|t|$, $U=\infty$, $J=\frac{2}{3}|t|$). The plot is arranged in
the same way as Fig. 1 and the color scale is consistent among panels.

###  III.2 $t$-$J$-$U$ model

We now proceed to analysis of the $t$-$J$-$U$ model ($U=12|t|$ and
$J=\frac{|t|}{3}$). In Fig. 2 we display the calculated imaginary parts of
dynamical spin- and charge susceptibilities for h- and e-doping. The results
are organized in the same manner as those in Fig. 1. The dispersive magnetic
peak appears on the h-doped side along the $\Gamma$-$X$ line, whereas
$\Gamma$-$M$ excitations are highly incoherent [panels (a)-(d)]. Also, the
magnetic mode rapidly dissolves in the continuum on the e-doped side for both
analyzed high-symmetry Brillouin-zone directions [panels (e)-(h)]. All those
features of the calculated spectra are qualitatively consistent with those
obtained for the Hubbard model, which points towards essentially universal
manifestation of the e-h asymmetry in magnetic dynamics.

The advantage of the $t$-$J$-$U$-model over either Hubbard- or $t$-$J$-model
limits is that the scale of local correlations and kinetic exchange
interactions are independently controlled by microscopic Hamiltonian
parameters. Due to the same adopted value of $J_{\mathrm{eff}}$, the
$t$-$J$-$U$-model dispersion of the coherent magnetic peak closely follows
that of the Hubbard model [cf. Fig. 1(a)-(d)], even though the on-site
interaction $U$ is larger by the factor of two. The physical consequence of
increased on-site interaction is substantial narrowing of the continuum part
of the spectrum and reduction of the charge-mode bandwidth, particularly close
to half-filling.

###  III.3 $t$-$J$ model as strong-coupling limit

Finally, we address the strong-coupling limit by setting $U=\infty$ and
$J=\frac{2}{3}|t|$. Figure 3 summarizes the calculated imaginary parts of the
spin- and charge dynamical susceptibilities. A dispersive paramagnon-like peak
with similar characteristics to those of Hubbard and $t$-$J$-$U$ models is
visible only for the $\Gamma$-$X$ direction and h-doping [panels (a)-(d)].
Interestingly, the magnetic spectra on the e-doped side [panels (e)-(h)] are
more robust and intense than those obtained for finite-$U$ models.

The principal qualitative feature of the $t$-$J$ model, absent for weak- and
intermediate coupling, is the emergence of sharp peaks at the magnetic
continuum threshold, both along the $\Gamma$-$X$ and $\Gamma$-$M$ directions.
Those are not unambiguously seen in spectroscopy, which suggests that the
$t$-$J$ model might overestimate electronic correlations in the context of
high-$T_{c}$ cuprates. The latter finding is also independently supported by
former studies, showing that the $t$-$J$-$U$ model yields better overall
agreement with experiment for those materials Spałek  _et al._ (2017).

The charge-excitation spectra, displayed in Fig. 3(i)-(p), are renormalized
with respect to Hubbard and $t$-$J$-$U$ models, yet no qualitatively new
behavior is observed.

Figure 4: Energy profiles of the imaginary part of the dynamical spin
susceptibility, $\chi^{\prime\prime}_{s}(\omega,\mathbf{k})$, obtained for the
Hubbard, $t$-$J$-$U$, and $t$-$J$ models. The data are plotted for
representative wave vectors, both on the h- and e-doping sides of the phase
diagram. The parameters and models are indicated inside the panels (wave
vectors are listed in the units $\frac{2\pi}{a}$ with $a$ being lattice
spacing).

###  III.4 Structure of the paramagnon peaks

The relation between Hubbard, $t$-$J$-$U$, and $t$-$J$ models may be further
characterized by comparing the energy profiles of the collective magnetic
excitation spectra. This is carried out in Fig. 4, where we display the
imaginary part of dynamical spin susceptibility, plotted for all three models
as a function of energy. Several representative doping levels and wave vectors
are selected, with top and bottom panels corresponding to the anti-nodal and
nodal directions, respectively (parameters are detailed inside the plot).

As noted above, paramagnon-like peaks emerge only for h-doping along the
$\Gamma$-$X$ line [panels (a), (c), and (e)]. However, their lineshapes vary
substantially depending on the on-site interaction strength, $U$. The magnetic
peak contains a cusp for the $t$-$J$ model up to $\delta_{h}\lesssim 0.2$
[blue and orange curves in Fig. 4(e)] and $t$-$J$-$U$ model up to
$\delta_{h}\lesssim 0.1$ [blue line in Fig. 4(c)]. This non-analytic feature
is reminiscent of a Van Hove singularity in two-particle density of states,
and thus _is not_ of resonant origin. Disentanglement of the resonant
(paramagnon) part of the spectrum from non-resonant (continuum) excitations
poses thus a major challenge in the strong-coupling regime. Those aspects will
be addressed in detail elsewhere Fidrysiak and Spałek .

##  IV Discussion

Figure 5: Calculated static spin- and charge susceptibilities for the Hubbard,
$t$-$J$-$U$, and $t$-$J$ models, both on the h- and e-doped side of the phase
diagram. The relevant parameters are detailed inside the panels. The obtained
susceptibilities remain positive along the high-symmetry $X$-$\Gamma$-$M$
contour which signals stability of the high-temperature paramagnetic state
against fluctuations.

In Sec. III we presented the VWF+$1/\mathcal{N}_{f}$ characteristics of the
collective spin- and charge excitations for the Hubbard, $t$-$J$-$U$, and
$t$-$J$ models of high-$T_{c}$ copper oxides, both on h- and e-doped side of
the phase diagram. The microscopic parameters of those three Hamiltonians have
been related to each other with the help of the effective exchange concept
Fidrysiak  _et al._ (2021) so that Figs. 1-3 can be directly compared.

The main finding of the present analysis is that robust paramagnon peak
emerges in the magnetic spectrum only along the anti-nodal $\Gamma$-$X$
direction for h-doping, whereas e-doped systems exhibit hardening of spin
excitations and less coherent dynamics. On the other hand, coherent charge
mode is found on both sides of the phase diagram. Those results are consistent
among all models analyzed, and remain in agreement with RIXS experiments for
multiple copper-oxide families.

We also recall that the only term explicitly breaking the e-h symmetry in the
Hamiltonian (1) is the next-nearest-neighbor hopping $\propto t^{\prime}$
which controls single-particle excitations and is necessary match the
high-$T_{c}$ copper oxide fermiology. Indirectly, $t^{\prime}$ impacts also
the structure of continuum particle-hole excitations that contribute to the
magnetic response. Thus, the obtained universal e-h asymmetry of the
paramagnon-like spectrum for wide range of interaction parameters suggests
that this particular effect is governed predominantly by the single-particle
kinematics rather than by multi-particle correlations.

We have also identified several principal differences between the spectra of
the Hubbard, $t$-$J$-$U$, and $t$-$J$ models. As demonstrated in Sec. III,
continuum spin excitations undergo a substantial renormalization with
increasing value of the on-site interaction, $U$. Moreover, in the strong
coupling ($t$-$J$ model) limit, an additional sharp peak in the magnetic
response appears at the continuum threshold. There are also more subtle
differences in the spectral lineshapes. We observe that, on approaching the
strong-coupling limit, the magnetic peaks evolve into the cusp-like structures
characterized by discontinuous first derivative of magnetic intensity. This,
in turn, points towards a non-negligible admixture of incoherent continuum
excitations to magnetic signal.

We now comment on apparently opposite effect of the e-h asymmetry on the
collective spin excitations and static magnetic order. In Appendix A we have
presented the VWF+$1/\mathcal{N}_{f}$ static susceptibilities. The magnetic
response attains maximum at the antiferromagnetic $M$ point for both h- and
e-doping, and is systematically larger on the e-side for the Hubbard and
$t$-$J$-$U$ models. This signals an enhanced tendency towards AF ordering for
e-doping, as seen in experiment Armitage  _et al._ (2010) (interestingly, this
hierarchy is reversed in the strong coupling limit). Our present analysis of
non-equilibrium quantities may be reconciled with the above result by noting
that robust magnetic excitations emerge only along the $\Gamma$-$X$ direction,
thus AF order and the paramagnons are governed by distinct regions of the
Brillouin zone.

Finally, we note that electronic correlations beyond mean-field theory must be
incorporated in order to accurately describe the magnetic dynamics, even for
relatively small values of $U$. This has been demonstrated for the Hubbard-
model case by comparison of the VWF+$1/\mathcal{N}_{f}$ and random-phase-
approximation dynamical spin susceptibilities, and subsequently relating them
both to experiment and quantum Monte-Carlo simulations Fidrysiak and Spałek
(2020, 2021a). The present analysis suggests that both those correlation-
related effects and non-resonant (incoherent) contributions are needed to
consistently describe the magnetic spectra throughout the high-$T_{c}$ phase
diagram.

## Acknowldgments

I thank Prof. J. Spałek for discussion. This work was supported by Grants
MINIATURA 5 No. DEC-2021/05/X/ST3/00666 and OPUS No. UMO-2021/41/B/ST3/04070
from Narodowe Centrum Nauki.

##  Appendix A Stability of the paramagnetic state

Here we analyze the stability of high-temperature paramagnetic state against
both spin- and charge fluctuations for the models considered in the main text.
In Fig. 5 we display calculated $X$-$\Gamma$-$M$ scans of spin [panels
(a)-(f)] and charge [panels (g)-(k)] static susceptibilities,
$\chi_{s}(\omega=0,\mathbf{k})$ and $\chi_{c}(\omega=0,\mathbf{k})$,
respectively. The doping levels and models are indicated inside the panels,
whereas the remaining parameters are the same as those taken in Sec. III. Both
spin and charge static susceptibilities remain positive in entire doping
regime, showing that the paramagnetic ground state is stable against
fluctuations on both h- and e-doped sides of the phase diagram.

There are several additional remarks that might be made, based on the analysis
of Fig. 5. First, for the Hubbard and $t$-$J$-$U$ models, the tendency towards
antiferromagnetic ordering is stronger for e-doping, as evidenced by a larger
value of static spin susceptibility close to the $M$ point, cf. panels
(a)-(d). This is qualitatively consistent with experimental asymmetry of
equilibrium magnetic phase diagram of the cuprates, with AF state being more
robust on the e-doped side. Interestingly, the $t$-$J$ model yields opposite
result [panels (e)-(f)], which favors finite-$U$ Hamiltonians for quantitative
studies of high-$T_{c}$ systems. Second, for the $t$-$J$-$U$ model, the
largest spin susceptibility magnitude is obtained for intermediate doping
levels so that the ordering tendency is shifted away from half-filling.

On the other hand, charge susceptibility remains fairly featureless for all
models and doping levels exceeding $0.1$. The enhancement of static charge
response around the $\Gamma$ point should be regarded as an unphysical
feature, since the Hamiltonian (1) does not include long-range Coulomb
interactions. The algebraic tail of the Coulomb potential is singular at
$\mathbf{k}\rightarrow 0$, which is known to suppress the magnetic response in
the vicinity of the Brillouin-zone center Spałek  _et al._ (2022).

##  Appendix B Electron-hole symmetry of the variational space

Here we show that the correlator (4) defines a variational space that is e-h
symmetric. As a consequence, the obtained asymmetry of the collective
excitations is an intrinsic feature of the $t$-$J$-$U$ model and is not
induced by the particular choice of trial wave function.

We proceed by rewriting the correlator terms in second-quantization language:

$\displaystyle|0\rangle_{i}{}_{i}\langle 0|$
$\displaystyle=(1-\hat{n}_{i\uparrow})(1-\hat{n}_{i\downarrow}),$ (5)
$\displaystyle|\sigma\rangle_{i}{}_{i}\langle\sigma|$
$\displaystyle=\hat{n}_{i\sigma}(1-\hat{n}_{i\bar{\sigma}}),$ (6)
$\displaystyle|\sigma\rangle_{i}{}_{i}\langle\bar{\sigma}|$
$\displaystyle=\hat{a}^{\dagger}_{i\sigma}\hat{a}_{i\bar{\sigma}},$ (7)
$\displaystyle|d\rangle_{i}{}_{i}\langle d|$
$\displaystyle=\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$ (8)

where $\bar{\sigma}=\downarrow(\uparrow)$ for $\sigma=\uparrow(\downarrow)$.
By applying the prescription (2) and noting that $\hat{n}_{i\sigma}\rightarrow
1-\hat{n}_{i\sigma}$, one arrives at

$\displaystyle\hat{P}_{i}\rightarrow\lambda^{0}_{i}|d\rangle_{i}{}_{i}\langle
d|+\sum_{\sigma\sigma^{\prime}}\sigma\sigma^{\prime}\lambda^{\sigma\sigma^{\prime}}_{i}|\bar{\sigma}\rangle_{i}{}_{i}\langle\bar{\sigma}^{\prime}|+\lambda^{d}_{i}|0\rangle_{i}{}_{i}\langle
0|.$ (9)

The transformed correlator (9) may be now brought to its original form (4) by
rearranging the variational parameters as
$\lambda^{0}_{i}\rightarrow\lambda_{i}^{d}$,
$\lambda^{\sigma\sigma^{\prime}}_{i}\rightarrow\sigma\sigma^{\prime}\cdot\lambda_{i}^{\bar{\sigma}\bar{\sigma}^{\prime}}$,
and $\lambda^{d}_{i}\rightarrow\lambda_{i}^{0}$, which concludes the
reasoning. Note that the transformation rules for the single-particle part of
the correlator, $\lambda^{\sigma\sigma^{\prime}}_{i}$, could be simplified by
supplementing the usual e-h transformation with spin-flip operation, and
including the spin-dependent phase factor in Eq. (2).

## References

  * Spałek  _et al._ (2022) J. Spałek, M. Fidrysiak, M. Zegrodnik, and A. Biborski, “Superconductivity in high-$T_{c}$ and related strongly correlated systems from variational perspective: Beyond mean field theory,” Phys. Rep. 959, 1–117 (2022).
  * Armitage  _et al._ (2010) N. P. Armitage, P. Fournier, and R. L. Greene, “Progress and perspectives on electron-doped cuprates,” Rev. Mod. Phys. 82, 2421–2487 (2010).
  * Le Tacon  _et al._ (2011) M. Le Tacon, G. Ghiringhelli, J. Chaloupka, M. Moretti Sala, V. Hinkov, M. W. Haverkort, M. Minola, M. Bakr, K. J. Zhou, S. Blanco-Canosa, and et al., “Intense paramagnon excitations in a large family of high-temperature superconductors,” Nat. Phys. 7, 725 (2011).
  * Dean  _et al._ (2013a) M. P. M. Dean, A. J. A. James, R. S. Springell, X. Liu, C. Monney, K. J. Zhou, R. M. Konik, J. S. Wen, Z. J. Xu, G. D. Gu, V. N. Strocov, T. Schmitt, and J. P. Hill, “High-Energy Magnetic Excitations in the Cuprate Superconductor $\mathrm{Bi_{2}Sr_{2}CaCu_{2}O_{8+\delta}}$: Towards a Unified Description of Its Electronic and Magnetic Degrees of Freedom,” Phys. Rev. Lett. 110, 147001 (2013a).
  * Dean  _et al._ (2013b) M. P. M. Dean, G. Dellea, R. S. Springell, F. Yakhou-Harris, K. Kummer, N. B. Brookes, X. Liu, Y-J. Sun, J. Strle, T. Schmitt, and et al., “Persistence of magnetic excitations in $\mathrm{La_{2-\mathit{x}}Sr_{\mathit{x}}CuO_{4}}$ from the undoped insulator to the heavily overdoped non-superconducting metal,” Nat. Mater. 12, 1019 (2013b).
  * Dean  _et al._ (2013c) M. P. M. Dean, G. Dellea, M. Minola, S. B. Wilkins, R. M. Konik, G. D. Gu, M. Le Tacon, N. B. Brookes, F. Yakhou-Harris, K. Kummer, J. P. Hill, L. Braicovich, and G. Ghiringhelli, “Magnetic excitations in stripe-ordered La1.875Ba0.125CuO4 studied using resonant inelastic $x$-ray scattering,” Phys. Rev. B 88, 020403 (2013c).
  * Le Tacon  _et al._ (2013) M. Le Tacon, M. Minola, D. C. Peets, M. Moretti Sala, S. Blanco-Canosa, V. Hinkov, R. Liang, D. A. Bonn, W. N. Hardy, C. T. Lin, T. Schmitt, L. Braicovich, G. Ghiringhelli, and B. Keimer, “Dispersive spin excitations in highly overdoped cuprates revealed by resonant inelastic x-ray scattering,” Phys. Rev. B 88, 020501 (2013).
  * Ishii  _et al._ (2014) K. Ishii, M. Fujita, T. Sasaki, M. Minola, G. Dellea, C. Mazzoli, K. Kummer, G. Ghiringhelli, L. Braicovich, T. Tohyama, and et al., “High-energy spin and charge excitations in electron-doped copper oxide superconductors,” Nat. Commun. 5, 3714 (2014).
  * Lee  _et al._ (2014) W. S. Lee, J. J. Lee, E. A. Nowadnick, S. Gerber, W. Tabiś, S. W. Huang, V. N. Strocov, E. M. Motoyama, G. Yu, B. Moritz, and et al., “Asymmetry of collective excitations in electron- and hole-doped cuprate superconductors,” Nat. Phys. 10, 883 (2014).
  * Dean  _et al._ (2014) M. P. M. Dean, A. J. A. James, A. C. Walters, V. Bisogni, I. Jarrige, M. Hücker, E. Giannini, M. Fujita, J. Pelliciari, Y. B. Huang, R. M. Konik, T. Schmitt, and J. P. Hill, “Itinerant effects and enhanced magnetic interactions in Bi-based multilayer cuprates,” Phys. Rev. B 90, 220506 (2014).
  * Guarise  _et al._ (2014) M. Guarise, B. Dalla Piazza, H. Berger, E. Giannini, T. Schmitt, H. M. Rønnow, G. A. Sawatzky, J. van den Brink, D. Altenfeld, I. Eremin, and et al., “Anisotropic softening of magnetic excitations along the nodal direction in superconducting cuprates,” Nat. Commun. 5, 5760 (2014).
  * Jia  _et al._ (2014) C. J. Jia, E. A. Nowadnick, K. Wohlfeld, Y. F. Kung, C.-C. Chen, S. Johnston, T. Tohyama, B. Moritz, and T. P. Devereaux, “Persistent spin excitations in doped antiferromagnets revealed by resonant inelastic light scattering,” Nat. Commun. 5 (2014), 10.1038/ncomms4314.
  * Wakimoto  _et al._ (2015) S. Wakimoto, K. Ishii, H. Kimura, M. Fujita, G. Dellea, K. Kummer, L. Braicovich, G. Ghiringhelli, L. M. Debeer-Schmitt, and G. E. Granroth, “High-energy magnetic excitations in overdoped $\mathrm{La_{2-\mathit{x}}Sr_{\mathit{x}}CuO_{4}}$ studied by neutron and resonant inelastic $x$-ray scattering,” Phys. Rev. B 91, 184513 (2015).
  * Minola  _et al._ (2015) M. Minola, G. Dellea, H. Gretarsson, Y. Y. Peng, Y. Lu, J. Porras, T. Loew, F. Yakhou, N. B. Brookes, Y. B. Huang, J. Pelliciari, T. Schmitt, G. Ghiringhelli, B. Keimer, L. Braicovich, and M. Le Tacon, “Collective Nature of Spin Excitations in Superconducting Cuprates Probed by Resonant Inelastic X-Ray Scattering,” Phys. Rev. Lett. 114, 217003 (2015).
  * Peng  _et al._ (2015) Y. Y. Peng, M. Hashimoto, M. Moretti Sala, A. Amorese, N. B. Brookes, G. Dellea, W.-S. Lee, M. Minola, T. Schmitt, Y. Yoshida, K.-J. Zhou, H. Eisaki, T. P. Devereaux, Z.-X. Shen, L. Braicovich, and G. Ghiringhelli, “Magnetic excitations and phonons simultaneously studied by resonant inelastic x-ray scattering in optimally doped ${\mathrm{Bi}}_{1.5}{\mathrm{Pb}}_{0.55}{\mathrm{Sr}}_{1.6}{\mathrm{La}}_{0.4}{\mathrm{CuO}}_{6+\delta}$,” Phys. Rev. B 92, 064517 (2015).
  * Ellis  _et al._ (2015) D. S. Ellis, Y.-B. Huang, P. Olalde-Velasco, M. Dantz, J. Pelliciari, G. Drachuck, R. Ofer, G. Bazalitsky, J. Berger, T. Schmitt, and A. Keren, “Correlation of the superconducting critical temperature with spin and orbital excitations in $({\mathrm{Ca}}_{x}{\mathrm{La}}_{1-x})({\mathrm{Ba}}_{1.75-x}\text{La}$ ${}_{0.25+x}){\mathrm{Cu}}_{3}{\mathrm{O}}_{y}$ as measured by resonant inelastic x-ray scattering,” Phys. Rev. B 92, 104507 (2015).
  * Huang  _et al._ (2016) H. Y. Huang, C. J. Jia, Z. Y. Chen, K. Wohlfeld, B. Moritz, T. P. Devereaux, W. B. Wu, J. Okamoto, W. S. Lee, M. Hashimoto, and et al., “Raman and fluorescence characteristics of resonant inelastic X-ray scattering from doped superconducting cuprates,” Sci. Rep. 6, 19657 (2016).
  * Gretarsson  _et al._ (2016) H. Gretarsson, N. H. Sung, J. Porras, J. Bertinshaw, C. Dietl, Jan A. N. Bruin, A. F. Bangura, Y. K. Kim, R. Dinnebier, Jungho Kim, and et al., “Persistent Paramagnons Deep in the Metallic Phase of $\mathrm{Sr_{2-\mathit{x}}La_{\mathit{x}}IrO_{4}}$,” Phys. Rev. Lett. 117, 107001 (2016).
  * Minola  _et al._ (2017) M. Minola, Y. Lu, Y. Y. Peng, G. Dellea, H. Gretarsson, M. W. Haverkort, Y. Ding, X. Sun, X. J. Zhou, D. C. Peets, L. Chauviere, P. Dosanjh, D. A. Bonn, R. Liang, A. Damascelli, M. Dantz, X. Lu, T. Schmitt, L. Braicovich, G. Ghiringhelli, B. Keimer, and M. Le Tacon, “Crossover from Collective to Incoherent Spin Excitations in Superconducting Cuprates Probed by Detuned Resonant Inelastic X-Ray Scattering,” Phys. Rev. Lett. 119, 097001 (2017).
  * Ivashko  _et al._ (2017) O. Ivashko, N. E. Shaik, X. Lu, C. G. Fatuzzo, M. Dantz, P. G. Freeman, D. E. McNally, D. Destraz, N. B. Christensen, T. Kurosawa, and et al., “Damped spin excitations in a doped cuprate superconductor with orbital hybridization,” Phys. Rev. B 95, 214508 (2017).
  * Meyers  _et al._ (2017) D. Meyers, H. Miao, A. C. Walters, V. Bisogni, R. S. Springell, M. d’ Astuto, M. Dantz, J. Pelliciari, H. Y. Huang, J. Okamoto, and et al., “Doping dependence of the magnetic excitations in $\mathrm{La_{2-\mathit{x}}Sr_{\mathit{x}}CuO_{4}}$,” Phys. Rev. B 95, 075139 (2017).
  * Chaix  _et al._ (2018) L. Chaix, E. W. Huang, S. Gerber, X. Lu, C. Jia, Y. Huang, D. E. McNally, Y. Wang, F. H. Vernay, A. Keren, and et al., “Resonant inelastic $x$-ray scattering studies of magnons band bimagnons in the lightly doped cuprate $\mathrm{La_{2-\mathit{x}}Sr_{\mathit{x}}CuO_{4}}$,” Phys. Rev. B 97, 155144 (2018).
  * Robarts  _et al._ (2019a) H. C. Robarts, M. Barthélemy, K. Kummer, M. García-Fernández, J. Li, A. Nag, A. C. Walters, K. J. Zhou, and S. M. Hayden, “Anisotropic damping and wave vector dependent susceptibility of the spin fluctuations in ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$ studied by resonant inelastic $x$-ray scattering,” Phys. Rev. B 100, 214510 (2019a).
  * Fumagalli  _et al._ (2019) R. Fumagalli, L. Braicovich, M. Minola, Y. Y. Peng, K. Kummer, D. Betto, M. Rossi, E. Lefrançois, C. Morawe, M. Salluzzo, and et al., “Polarization-resolved Cu $L_{3}$-edge resonant inelastic $x$-ray scattering of orbital and spin excitations in $\mathrm{NdBa_{2}Cu_{3}O_{7-\delta}}$,” Phys. Rev. B 99, 134517 (2019).
  * Peng  _et al._ (2018) Y. Y. Peng, E. W. Huang, R. Fumagalli, M. Minola, Y. Wang, X. Sun, Y. Ding, K. Kummer, X. J. Zhou, N. B. Brookes, and et al., “Dispersion, damping, and intensity of spin excitations in the monolayer $\mathrm{(Bi,Pb)_{2}(Sr,La)_{2}CuO_{6+\delta}}$ cuprate superconductor family,” Phys. Rev. B 98, 144507 (2018).
  * Robarts  _et al._ (2021) H. C. Robarts, M. García-Fernández, J. Li, A. Nag, A. C. Walters, N. E. Headings, S. M. Hayden, and K.-J. Zhou, “Dynamical spin susceptibility in ${\mathrm{La}}_{2}{\mathrm{CuO}}_{4}$ studied by resonant inelastic x-ray scattering,” Phys. Rev. B 103, 224427 (2021).
  * Ishii  _et al._ (2017) K. Ishii, T. Tohyama, S. Asano, K. Sato, M. Fujita, S. Wakimoto, K. Tustsui, S. Sota, J. Miyawaki, H. Niwa, and et al., “Observation of momentum-dependent charge excitations in hole-doped cuprates using resonant inelastic $x$-ray scattering at the oxygen $K$ edge,” Phys. Rev. B 96, 115148 (2017).
  * Nag  _et al._ (2020) A. Nag, M. Zhu, M. Bejas, J. Li, H. C. Robarts, H. Yamase, A. N. Petsch, D. Song, H. Eisaki, A. C. Walters, M. García-Fernández, A. Greco, S. M. Hayden, and K.-J. Zhou, “Detection of Acoustic Plasmons in Hole-Doped Lanthanum and Bismuth Cuprate Superconductors Using Resonant Inelastic X-Ray Scattering,” Phys. Rev. Lett. 125, 257002 (2020).
  * Singh  _et al._ (2022) A. Singh, H. Y. Huang, Christopher Lane, J. H. Li, J. Okamoto, S. Komiya, Robert S. Markiewicz, Arun Bansil, T. K. Lee, A. Fujimori, C. T. Chen, and D. J. Huang, “Acoustic plasmons and conducting carriers in hole-doped cuprate superconductors,” Phys. Rev. B 105, 235105 (2022).
  * Hepting  _et al._ (2018) M. Hepting, L. Chaix, E. W. Huang, R. Fumagalli, Y. Y. Peng, B. Moritz, K. Kummer, N. B. Brookes, W. C. Lee, M. Hashimoto, and et al., “Three-dimensional collective charge excitations in electron-doped copper oxide superconductors,” Nature 563, 374 (2018).
  * Ishii  _et al._ (2019) K. Ishii, M. Kurooka, Y. Shimizu, M. Fujita, K. Yamada, and J. Mizuki, “Charge Excitations in $\mathrm{Nd_{2-\mathit{x}}Ce_{\mathit{x}}CuO_{4}}$ Observed with Resonant Inelastic $X$-ray Scattering: Comparison of Cu $K$-edge with Cu $L_{3}$-edge,” J. Phys. Soc. Japan 88, 075001 (2019).
  * Lin  _et al._ (2020) J. Lin, J. Yuan, K. Jin, Z. Yin, Gang Li, K.-J. Zhou, X. Lu, M. Dantz, T. Schmitt, H. Ding, H. Guo, M. P. M. Dean, and X. Liu, “Doping evolution of the charge excitations and electron correlations in electron-doped superconducting $\mathrm{La_{2-\mathit{x}}Ce_{\mathit{x}}CuO_{4}}$,” npj Quant. Mater. 5, 4 (2020).
  * Chubukov  _et al._ (2003) A. V. Chubukov, D. Pines, and J. Schmalian, “A Spin Fluctuation Model for $d$-Wave Superconductivity,” in  _The Physics of Superconductors_, Vol. 1, edited by K. H. Bennemann (Springer Berlin Heidelberg, 2003) pp. 495–590.
  * Scalapino (2012) D. J. Scalapino, “A common thread: The pairing interaction for unconventional superconductors,” Rev. Mod. Phys. 84, 1383 (2012).
  * Yu  _et al._ (2020) B. Yu, W. Tabiś, I. Biało, F. Yakhou, N. B. Brookes, Z. Anderson, Y. Tang, G. Yu, and M. Greven, “Unusual Dynamic Charge Correlations in Simple-Tetragonal ${\mathrm{HgBa}}_{2}{\mathrm{CuO}}_{4+\delta}$,” Phys. Rev. X 10, 021059 (2020).
  * Zaanen  _et al._ (1985) J. Zaanen, G. A. Sawatzky, and J. W. Allen, “Band gaps and electronic structure of transition-metal compounds,” Phys. Rev. Lett. 55, 418–421 (1985).
  * Rybicki  _et al._ (2016) D. Rybicki, M. Jurkutat, S. Reichardt, C. Kapusta, and J. Haase, “Perspective on the phase diagram of cuprate high-temperature superconductors,” Nature Communications 7, 11413 (2016).
  * Zegrodnik  _et al._ (2019) M. Zegrodnik, A. Biborski, M. Fidrysiak, and J. Spałek, “Superconductivity in the three-band model of cuprates: Variational wave function study and relation to the single-band case,” Phys. Rev. B 99, 104511 (2019).
  * Zegrodnik  _et al._ (2020) M. Zegrodnik, A. Biborski, and J. Spałłek, “Superconductivity and intra-unit-cell electronic nematic phase in the three-band model of cuprates,” Eur. Phys. J. B 93, 183 (2020).
  * Zegrodnik  _et al._ (2021) M. Zegrodnik, A. Biborski, M. Fidrysiak, and J. Spałek, “Superconductivity in the three-band model of cuprates: nodal direction characteristics and influence of intersite interactions,” J. Phys.: Condens. Matter 33, 415601 (2021).
  * Li  _et al._ (2021) S. Li, A. Nocera, U. Kumar, and S. Johnston, “Particle-hole asymmetry in the dynamical spin and charge responses of corner-shared 1d cuprates,” Communications Physics 4, 217 (2021).
  * Spałek  _et al._ (2017) J. Spałek, M. Zegrodnik, and J. Kaczmarczyk, “Universal properties of high-temperature superconductors from real-space pairing: $t$-$J$-$U$ model and its quantitative comparison with experiment,” Phys. Rev. B 95, 024506 (2017).
  * Zegrodnik and Spałek (2017) M. Zegrodnik and J. Spałek, “Universal properties of high-temperature superconductors from real-space pairing: Role of correlated hopping and intersite Coulomb interaction within the $t$-$J$-$U$ model,” Phys. Rev. B 96, 054511 (2017).
  * Fidrysiak  _et al._ (2018) M. Fidrysiak, M. Zegrodnik, and J. Spałek, “Realistic estimates of superconducting properties for the cuprates: reciprocal-space diagrammatic expansion combined with variational approach,” J. Phys.: Condens. Matter 30, 475602 (2018).
  * Wysokiński and Fabrizio (2017a) M. M. Wysokiński and M. Fabrizio, “Mott physics beyond the Brinkman-Rice scenario,” Phys. Rev. B 95, 161106 (2017a).
  * Schiró and Fabrizio (2010) M. Schiró and M. Fabrizio, “Time-Dependent Mean Field Theory for Quench Dynamics in Correlated Electron Systems,” Phys. Rev. Lett. 105, 076401 (2010).
  * Wysokiński and Fabrizio (2017b) M. M. Wysokiński and M. Fabrizio, “Interplay of charge and spin dynamics after an interaction quench in the Hubbard model,” Phys. Rev. B 96, 201115 (2017b).
  * Fidrysiak  _et al._ (2021) M. Fidrysiak, D. Goc-Jagło, and J. Spałek, “Collective spin and charge excitations in the $t$-$J$-$U$ model of high-$T_{c}$ cuprates,” J. Magn. Magn. Mater. 539, 168395 (2021).
  * Fidrysiak and Spałek (2021a) M. Fidrysiak and J. Spałek, “Universal collective modes from strong electronic correlations: Modified $1/{\mathcal{N}}_{f}$ theory with application to high-${T}_{c}$ cuprates,” Phys. Rev. B 103, 165111 (2021a).
  * Fidrysiak and Spałek (2020) M. Fidrysiak and J. Spałek, “Robust spin and charge excitations throughout the high-${T}_{c}$ cuprate phase diagram from incipient Mottness,” Phys. Rev. B 102, 014505 (2020).
  * Fidrysiak and Spałek (2021b) M. Fidrysiak and J. Spałek, “Unified theory of spin and charge excitations in high-${T}_{c}$ cuprate superconductors: A quantitative comparison with experiment and interpretation,” Phys. Rev. B 104, L020510 (2021b).
  * Robarts  _et al._ (2019b) H. C. Robarts, M. Barthélemy, K. Kummer, M. García-Fernández, J. Li, A. Nag, A. C. Walters, K. J. Zhou, and S. M. Hayden, “Anisotropic damping and wave vector dependent susceptibility of the spin fluctuations in ${\mathrm{La}}_{2-x}{\mathrm{Sr}}_{x}{\mathrm{CuO}}_{4}$ studied by resonant inelastic x-ray scattering,” Phys. Rev. B 100, 214510 (2019b).
  * (53) M. Fidrysiak and J. Spałek, unpublished.


# Composition Dependent Instabilities in Mixtures With Many Components

Filipe C. Thewes  Institut für Theoretische Physik, Georg-August-Universität
Göttingen, 37077 Göttingen, Germany  Matthias Krüger  Institut für
Theoretische Physik, Georg-August-Universität Göttingen, 37077 Göttingen,
Germany  Peter Sollich  Institut für Theoretische Physik, Georg-August-
Universität Göttingen, 37077 Göttingen, Germany  King’s College London,
Department of Mathematics, Strand, London WC2R 2LS, U.K.

###### Abstract

Understanding the phase behavior of mixtures with many components is important
in many contexts, including as a key step toward a physics-based description
of intracellular compartmentalization. Here, we study the instabilities of a
mixture model where the second virial coefficients are taken as random
Gaussian variables. Using tools from free probability theory we obtain the
exact spinodal curve and the nature of instabilities for a mixture with an
arbitrary composition, thus lifting the assumption of uniform mixture
component densities pervading previous studies. We show that, by controlling
the volume fraction of only a few components, one can systematically change
the nature of the spinodal instability and achieve demixing for realistic
scenarios by a strong composition imbalance amplification. This results from a
non-trivial interplay of entropic effects due to non-uniform composition and
complexity in the interactions. Our approach allows for the inclusion of any
finite number of structured interactions, leading to a competition between
different forms of demixing as density is varied.

Phase separation is an important phenomenon and is especially rich in mixtures
with many components. In particular, biological mixtures such as the cytoplasm
show a complex phase behavior believed to be a key driver in the formation of
nucleoli and other intracellular structures Alberti (2017); Lafontaine  _et
al._ (2020); Sear (2007); Berry  _et al._ (2018); Choi  _et al._ (2020).
Typically, these systems demix into liquid droplets with different
compositions, where each phase is enriched in a number of components and
depleted in others Berry  _et al._ (2018); Keating (2012).

Within a mean-field approach, previous studies have explored the rich behavior
of complex mixtures by investigating the number of phases formed for
appropriately tuned Jacobs (2021), evolutionarily optimized Zwicker and Laan
(2022) and random Jacobs and Frenkel (2013, 2017) interactions. Within the
latter approach Sear and Cuesta (2003) dynamical properties have also been
investigated, including early time instabilities using random matrix theory
Sear and Cuesta (2003) and direct numerical simulations Shrinivas and Brenner
(2021) of continuum field theories (model B).

Although these studies have been able to illuminate some of the behavior of
complex mixtures, a key restrictive assumption common to all of them is
uniform composition, meaning all species are present in the system in equal
amounts. However, biological mixtures rarely satisfy this condition Sear
(2005); Gasic  _et al._ (2021); Neidhardt and Umbarger (1987), and the size of
intracellular structures is in fact heavily dependent on the composition of
the cytoplasmic pool Weber and Brangwynne (2015); Goehring and Hyman (2012).
In addition, experimental protocols often rely on controlling the composition
of the mixture to study its phase behavior Gasic  _et al._ (2021); Riback  _et
al._ (2020); Zhang  _et al._ (2015); Fritsch  _et al._ (2021).

With this in mind, our aim in this letter is to open the door to exploring the
full composition-dependent complexity of multi-component mixtures. To achieve
this, we build on the model of Sear and Cuesta Sear and Cuesta (2003) but lift
the drastic simplification of uniform composition. We show that systematically
changing the number density of only a few components enables one to control
the nature of the instabilities and, as a consequence, the phases that can be
formed.

Model and general results: Following Sear and Cuesta (2003) we study a mixture
of $M$ different components, labeled by Greek letters, with interactions
described by the second virial coefficients $\epsilon_{\alpha\gamma}$. The
mean-field free energy density $f$ is given by Weber  _et al._ (2019)

$f=\frac{1}{2}\sum_{\alpha,\gamma=1}^{M}\rho_{\alpha}\epsilon_{\alpha\gamma}\rho_{\gamma}+T\sum_{\alpha=1}^{M}\rho_{\alpha}\ln\rho_{\alpha}+T\rho_{0}\ln\rho_{0}$
(1)

where $\rho_{\alpha}=N_{\alpha}/V$ is the number density of species $\alpha$,
$T$ is temperature and we use $k_{\rm B}=1$. The last term in (1) is the
entropic contribution of an implicit solvent, interacting only via volume
exclusion. We define the total density as $\rho=\sum_{\alpha}\rho_{\alpha}$,
the average density per component as $\bar{\rho}=\rho/M$, and work with units
such that $\rho_{0}=1-\rho\geq 0$.

From (1) we obtain the $M\times M$ Hessian matrix
$H_{\alpha\gamma}=\partial^{2}f/(\partial\rho_{\alpha}\partial\rho_{\gamma})$
as
$\bm{H}=(T/\rho_{0})\bm{uu}^{\mathsf{T}}+(T/\bar{\rho})\textrm{diag}\left(1/y_{\alpha}\right)+\bm{\epsilon}$
with $\bm{u}=(1,1\dots,1)^{\mathsf{T}}$ the constant vector and
$\textrm{diag}(1/y_{\alpha})$ a diagonal matrix with entries determined by the
relative densities $y_{\alpha}=\rho_{\alpha}/\bar{\rho}$; by definition the
latter have average $M^{-1}\sum_{\alpha}y_{\alpha}=1$. Thermodynamic stability
requires all eigenvalues of $\bm{H}$ to be non-negative. Otherwise, i.e. if
the lowest eigenvalue $\lambda_{\min}$ is negative, the system is unstable to
phase separation by spinodal decomposition.

The phase diagram in the ($\rho$,$T$)-plane for a fixed composition
$\\{y_{\alpha}\\}$ splits into stable and unstable regions, separated by a
spinodal line determined by the condition $\lambda_{\min}=0$ Sollich (2002).
The nature of the spinodal instability is determined by the eigenvector
$\bm{v}$ corresponding to $\lambda_{\min}$. We will be interested in
instabilities of condensation type ($\bm{v}\sim\bm{u}$), where the densities
of all species change by similar amounts, and of demixing type
($\bm{v}^{\mathsf{T}}\bm{u}\approx 0$), where some species are enhanced while
others are depleted. We will show that the demixing case can be further split
into delocalized or random, where all components of $\bm{v}$ are of similar
order, and localized where a few species have much larger entries in $\bm{v}$
and thus dominate the demixing.

Following Ref. Sear and Cuesta (2003), we model the second virial coefficients
$\bm{\epsilon}_{\alpha\gamma}$ as Gaussian random variables of mean $-b$ and
variance $s^{2}$, drawn independently except for the symmetry constraint
$\epsilon_{\alpha\gamma}=\epsilon_{\gamma\alpha}$. The Hessian matrix then
reads

$\displaystyle\begin{split}\bm{H}&=\bm{R}_{1}+\bm{D}+s\bm{\eta}\qquad\mbox{with}\\\
\bm{R}_{1}&=\left(-b+\frac{T}{\rho_{0}}\right)\bm{uu}^{\mathsf{T}},\quad\bm{D}=\frac{T}{\bar{\rho}}\textrm{diag}\left(\frac{1}{y_{\alpha}}\right)\end{split}$
(2)

and $\bm{\eta}$ a Wigner matrix with entries of zero mean and unit variance
Akemann  _et al._ (2011). Understanding thermodynamic instabilities then
requires us to obtain the eigenvalue distribution or spectrum of $\bm{H}$, and
specifically its lower edge $\lambda_{\min}$; we focus throughout on the
interesting multi-component limit $M\gg 1$.

The $s\bm{\eta}$ term in $\bm{H}$ produces a continuous spectrum of
eigenvalues, and this extends to $\bm{D}+s\bm{\eta}$ 111Only if the entries of
$\bm{D}$ are relatively close to each other. Otherwise, the spectrum may split
intro multiple bulk pieces.. The first term $\bm{R}_{1}$ can be viewed as a
rank one perturbation; due to the so-called interlacing property of
eigenvalues Hwang (2004), the spectrum of $\bm{H}$ for large $M$ is then
either the same as the spectrum of $\bm{D}+s\bm{\eta}$, or $\bm{R}_{1}$ may
give rise to a single outlier Benaych-Georges and Nadakuditi (2011), which
separates from the continuous bulk spectrum of eigenvalues. We therefore have
two regimes: if an outlier exists to the left of the bulk, then it is the
lowest eigenvalue $\lambda_{\min}$. Otherwise the lowest eigenvalue is given
by the lower edge of the bulk itself.

Free probability Voiculescu  _et al._ (1992); Akemann  _et al._ (2011) is a
powerful tool to obtain the statistics of eigenvalues and eigenvectors of
large random matrices, provided they obey the so-called freeness criteria. A
key insight is that freeness generically holds between $\bm{D}$, $s\bm{\eta}$
and $\bm{R}_{1}$ (see Supplemental Material at for discussion). We can thus
use free probability to analyze the spectrum of the scaled Hessian $\bm{H}/M$;
adopting also the scaling $s=M^{1/2}\tilde{s}$ Sear and Cuesta (2003) ensures
that all matrices involved have eigenvalues of $O(1)$. We find (see
Supplemental Material at for this and subsequent derivations) for the spinodal
equation, which determines where $\lambda_{\min}=0$,

$\psi(-\rho\tilde{s}^{2}z/T)=\tilde{s}^{2}z^{2}\,.$ (3)

Here $\psi(x)=\langle xy_{\alpha}/(1-xy_{\alpha})\rangle$ with the angular
brackets denoting an average over the distribution $p(y_{\alpha})$ of
$y_{\alpha}$, which specifies the mixture composition. Equation (3) is to be
solved for $T$ as a function of the total density $\rho$. The difference
between the two regimes discussed above, i.e. outlier and bulk, lies in the
way $z=z(\rho,T)$ is determined as we explain next.

First, in the outlier regime one has $z=\theta^{-1}$, where
$\theta=T/\rho_{0}-b$ is the non-zero eigenvalue of $\bm{R}_{1}/M$; notice
that $\theta$ has to be negative to give rise to an outlier to the left of the
bulk. We can also determine the overlap between the (normalized) instability
vector $\bm{v}$ and the normalized uniform vector
$\hat{\bm{u}}=\bm{u}/\sqrt{M}$ as

$|\bm{v}^{\mathsf{T}}\hat{\bm{u}}|^{2}=\max\left\\{-\frac{\tilde{s}^{2}}{\theta^{2}}\left[1+\frac{T^{2}}{\rho^{2}\tilde{s}^{2}G^{\prime}(-\rho\tilde{s}^{2}\theta^{-1}/T)}\right],0\right\\}$
(4)

where $G(x)=\langle(x-1/y_{\alpha})^{-1}\rangle$ is the so-called resolvent,
the prime denotes derivative and all quantities, i.e. $\theta$, $T$ for a
given $\rho$, are evaluated on the spinodal curve. Following Ref. Sear and
Cuesta (2003) we shall refer to instabilities with non-vanishing overlap
between $\bm{v}$ and $\hat{\bm{u}}$ as condensation (C), and as demixing
otherwise. Equation (4) thus provides information on the nature of the
instability and is one of our key results.

Turning next to the bulk regime, we have already written Eq. (4) in a form
that applies also there. The first argument of the $\max$ function in (4) is
then negative and we have demixing behavior, with the lowest eigenvalue of the
bulk $\bm{D}+s\bm{\eta}$ determining thermodynamic stability. For the
restricted case of uniform composition, $\bm{D}=(T/\bar{\rho})\bm{I}$ is
proportional to the identity matrix and just shifts the spectrum of
$s\bm{\eta}$ by $T/\bar{\rho}$; the entire demixing regime is then described
by a linear spinodal line $T\propto\bar{\rho}$ Sear and Cuesta (2003). For
non-uniform composition, on the other hand, the spread of eigenvalues in
$\bm{D}$ can dominate at high enough $T$ Bouchbinder  _et al._ (2021) as
illustrated in Fig. 1. The edge of the bulk of the spectrum is then determined
by $\bm{D}$ and in the spinodal condition (3) one has $z=-T/(\tilde{s}^{2}\rho
y_{\max})$, where $y_{\max}=\max_{\alpha}y_{\alpha}$. The corresponding
eigenvector of $\bm{D}$ only has a single non-zero entry, and we find as a
result that the instability direction $\bm{v}$ becomes concentrated on a few
species. For low enough $T$, on the other hand, the $s\bm{\eta}$ term
dominates in $\bm{D}+s\bm{\eta}$ and $z$ has a larger value maximizing
$F(g)=-\tilde{s}^{2}g+(1/g)\psi(-\rho\tilde{s}^{2}g/T)$, i.e. $z=g^{*}$ with
$g^{*}$ determined from $F^{\prime}(g^{*})=0$. The threshold temperature
$T^{*}$ separating these two cases (see Fig. 1) is the one where
$F^{\prime}(-T/\tilde{s}^{2}\rho y_{\max})=0$, and is given explicitly by

$T^{*}=\tilde{s}\rho y_{\max}\sqrt{\left\langle
y_{\alpha}^{2}/(y_{\max}-y_{\alpha})^{2}\right\rangle}.$ (5)

Figure 1: Exemplary realization of eigenvalue distributions of $\bm{D}$ and
$s\bm{\eta}$. The sum $\bm{D}+s\bm{\eta}$ becomes dominated by the interaction
complexity for $T<T^{*}$, and by the composition $p(y_{\alpha})$ for
$T>T^{*}$. We used for $p(y_{\alpha})$ a Beta distribution. For visualization
purposes the spectra have been shifted to have matching means.

Finally, the transition between condensation (outlier) and demixing (bulk)
regimes occurs when the two solutions meet each other, namely when
$g^{*}=\theta^{-1}$. Along the spinodal curve this condition defines a
threshold density $\rho^{*}$ such that, for $\rho<\rho^{*}$, the spinodal is
condensation-like and otherwise of demixing type.

Taken together, our results give a complete characterization of the spinodal
line for any composition $p(y_{\alpha})$ via equation (3), and of the nature
of the spinodal instability via (4). In addition to the condensation-demixing
transition Sear and Cuesta (2003) Eq. (5) reveals the novel possibility of
composition-driven demixing. We show below that this is not a simple entropic
effect, but arises instead from an interplay between the composition and the
complexity in the interactions.

Example 1: Uniform composition. The case $p(y_{\alpha})=\delta(y_{\alpha}-1)$,
where all mixture components have the same density, has been studied in Sear
and Cuesta (2003); Moran and Bouchaud (2019). We revisit it briefly in order
to illustrate the formalism developed so far. Starting from Eq. (3),
straightforward algebra yields $T/\rho+z^{-1}+z\tilde{s}^{2}=0$. In the
condensation regime, $z=\theta^{-1}$ with $\theta=T/(1-\rho)-b$, giving a
quadratic 222Sear and Cuesta Sear and Cuesta (2003) omitted the explicit
entropic solvent contribution, resulting in the spinodal curve always being
linear in $\rho$. Including the solvent entropy as done in other recent
studies Jacobs and Frenkel (2017); Jacobs (2021); Shrinivas and Brenner (2021)
makes $\theta$ dependent on $T$ and $\rho$ and gives a nonlinear spinodal
curve in the condensation regime. equation for the spinodal line $T(\rho)$. In
the demixing regime, $T^{*}$ from Eq. (5) diverges so that thermodynamic
stability is always governed by the interaction complexity $s\bm{\eta}$. One
then finds from the condition $F^{\prime}(g^{*})=0$ that
$\tilde{s}g^{*}=1-T/(\rho\tilde{s})$, and inserting $z=g^{*}$ into the
spinodal equation yields $T=2\tilde{s}\rho$.

To understand the nature of the instabilities we refer to Eq. (4). In the
uniform distribution case we have
$G^{\prime}(-\rho\tilde{s}^{2}\theta^{-1}/T)=-(\tilde{s}^{2}\theta^{-1}+T/\rho)^{-2}T^{2}/\rho^{2}$
and the term in brackets equals $-\theta$ by the spinodal equation. This
yields
$|\bm{v}^{\mathsf{T}}\hat{\bm{u}}|^{2}=\max\\{1-\tilde{s}^{2}/\theta^{2},0\\}$
so that we are in the condensation regime as long as
$-\theta=b-T/(1-\rho)>\tilde{s}$. This is always the case at low densities,
where $T$ and hence $T/(1-\rho)$ vanishes along the spinodal, provided that
$b>\tilde{s}$. Therefore, at low densities
$|\bm{v}^{\mathsf{T}}\hat{\bm{u}}|=O(1)$ with densities of all species
changing by similar amounts at the spinodal instability. As the total density
is increased, $-\theta$ decreases and can approach $\tilde{s}$; $\bm{v}$ and
$\hat{\bm{u}}$ then become orthogonal, resulting in random demixing (RD), with
a delocalized instability vector $\bm{v}$ Benaych-Georges and Nadakuditi
(2011).

Example 2: One dominant species. Next, we investigate the case of one single
dominant species ($\alpha=1$) with relative concentration $y_{1}>1$, while all
other species have $y_{2}=(M-y_{1})/(M-1)<1$ (see Supplemental Material at for
$M$-dependent effects). This example will show how tuning the density of a few
species can change the nature of the spinodal instability at high densities.

As a consequence of the single distinct entry in $\bm{D}$, two possibilities
exist for the demixing regime, depending on whether the lowest eigenvalue of
$\bm{D}+s\bm{\eta}$ is controlled by $s\bm{\eta}$ or by the distinct entry in
$\bm{D}$. Solving the spinodal equation (3) and computing the nature of the
instability (See Supplemental Material at ) yields three different regimes
depending on $y_{1}$. At low densities, the spinodal is dominated by the
average interaction $-b$ and by entropic effects, yielding condensation
behavior. Increasing $\rho$ results in a transition to demixing. We find
explicitly for the instability direction in the demixing regime
$|\bm{v}^{\mathsf{T}}\hat{\bm{u}}|=0$ and

$|\bm{v}^{\mathsf{T}}\bm{e}_{1}|^{2}=\max\left\\{\frac{y_{1}-2}{y_{1}-1},0\right\\}$
(6)

where $\bm{e}_{1}=(1,0,0,\dots,0)^{\mathsf{T}}$ is the direction of the
dominant species. The $O(1)$ overlap between $\bm{v}$ and $\bm{e}_{1}$
demonstrates that, whenever $y_{2}>2$, we have composition-driven demixing
(CD) controlled by the dominant species. If, on the other hand, $y_{2}<2$ the
instability is controlled by $s$ and the mixture will undergo random demixing.
The transition from C to CD happens at
$\rho^{*}=1-y_{1}\tilde{s}/(b\sqrt{y_{1}-1}+\tilde{s})$ and the CD-spinodal
for $\rho>\rho^{*}$ follows $T=y_{1}\tilde{s}\rho/\sqrt{y_{1}-1}$.

Figure 2: One dominant species case. Spinodal lines for different
$y_{1}=\rho_{1}/\bar{\rho}$ at fixed $M=100$, $s=1$ and $b=1$ (lines: theory,
symbols: average over $50$ numerical realizations of Hessian matrix). For
$y_{1}>2$ the demixing is controlled by the dominant species. Insets:
Projection of instability direction onto the constant vector $\hat{\bm{u}}$
and the dominant species direction $\bm{e}_{1}$, showing the transitions from
condensation (C) to random (RD) and composition-driven demixing (CD),
respectively. We keep the $M$-dependence of $\theta$ to account for finite-
size effects when comparing to numerics.

In Fig. 2 we compare the predictions for the spinodal curves for different
$y_{1}$ and the corresponding instability direction; the comparison to results
from numerical realizations of the Hessian matrix shows excellent agreement.
For $y_{1}>2$ the instability vector is strongly concentrated on the dominant
species at high densities, with $v_{1}/v_{2}=O(\sqrt{M})$. What is striking in
this CD region is the fact that the share of the dominant species in the
instability direction is much larger than expected from entropic
considerations, which would predict $v_{1}/v_{2}\sim y_{1}/y_{2}=O(1)$. This
strong composition imbalance amplification is our key insight into
instabilities in complex mixtures. It results from an interplay of entropic
effects and complexity in the interactions ($\tilde{s}>0$); indeed the CD
regime would be absent in the limit $\tilde{s}\to 0$ (where $\rho^{*}\to 1$).

Figure 3: (a) Beta distribution $p(y)$ with $y_{\min}=0$, $r=0$ and
$t=y_{\max}-2$. The inset shows the complement of the cumulative distribution
function CDF(y), highlighting the different tail behaviors. (b) Spinodal line
at fixed $M=600$, $b=1.0$, $s=1.5$ for same distribution parameters
(corresponding colors) as in (a). Inset: inverse participation ratio (IPR),
showing the localization transition to CD on the right ($y_{\max}=10$). (c)
Spinodal lines for Beta distribution with $y_{\max}=5$ and systematic
interactions, for different values of
$\delta=\tilde{s}/\sqrt{\textrm{Var}_{\sigma}+\tilde{s}^{2}}$. Inset:
projection of instability direction onto the two principal interaction
directions, i.e. the uniform vector $\hat{\bm{u}}$ resulting from the volume
exclusion and the remaining systematic interaction
$\tilde{\bm{\sigma}}=\bm{\sigma}-\langle\sigma\rangle\bm{u}$; for small enough
$\delta$ a new regime with instability direction dominated by systematic
interactions appears.

Example 3: Beta distribution. We now turn our attention to a continuous
distribution $p(y_{\alpha})$ and show how its shape affects instabilities. We
focus on the Beta distribution defined by

$p(y)=Z_{r,t}^{-1}(y-y_{\min})^{r}(y_{\max}-y)^{t},\quad
y\in[y_{\min},y_{\max}]$ (7)

with $Z_{r,t}=(y_{\max}-y_{\min})^{r+t+1}B(r+1,t+1)$ and $B(x,y)$ the Beta
function. The requirement $\langle y\rangle=1$ reduces the number of free
parameters to three.

The spinodal equation (3) can be solved numerically in the condensation (C)
regime. More interesting here is the demixing regime. According to Eq. (5),
for temperatures greater than a threshold temperature $T^{*}$, the demixing
instability changes from being determined by the interaction complexity (RD)
to being governed by the composition (CD), i.e. the distribution of $y$.
Inserting (7) into (5) one finds

$T^{*}=\tilde{s}\rho y_{\max}\sqrt{(Z_{r,t-2}/Z_{r,t})\langle
y^{2}\rangle_{r,t-2}}$ (8)

where $\langle\cdot\rangle_{l,m}$ denotes the average over the Beta
distribution (7) with exponent parameters $l$ and $m$. For $T>T^{*}$ the
spinodal equation (3) can be solved analytically, yielding

$T=\tilde{s}\rho y_{\max}\sqrt{(Z_{r,t-1}/Z_{r,t})\langle y\rangle_{r,t-1}}.$
(9)

Comparing (8) and (9), we see that whenever $Z_{r,t-1}\langle
y\rangle_{r,t-1}>Z_{r,t-2}\langle y^{2}\rangle_{r,t-2}$ the demixing spinodal
is dictated by the mixture composition. This condition is independent of the
total density $\rho$ and controlled only by the shape of $p(y)$. If e.g. we
fix $y_{\min}=0$, this condition for CD reduces to $t>r+3$, meaning that the
upper edge of the distribution has a much longer tail than its lower edge. In
other words, the transition occurs whenever a small fraction of species has
significantly larger density than the average, mirroring our results from
Example 2.

Regarding the nature of the instability, we can exploit the results of Ref.
Lee and Schnelli (2016) to show that the instability direction is delocalized
across species in the RD regime ($T<T^{*}$). In the CD regime, the instability
direction is concentrated on a few dominant species that have entries of
$O(1)$ in $\bm{v}$ and will, therefore, dictate the nature of the spinodal
instability. Translating the results of Lee and Schnelli (2016) further to our
context, the contribution of the highest-density species (denoted by
$\bm{e}_{1}$) to the instability direction is
$|\bm{v}^{\mathsf{T}}\bm{e}_{1}|^{2}=1-(T^{*}/T)^{2}$, which is independent of
$\rho$ along the CD spinodal. For the subsequent dominant species ($j>1$) one
has $|\bm{v}^{\mathsf{T}}\bm{e}_{j}|^{2}\sim
1/(M^{\gamma}T^{2}|\lambda_{j}-\lambda_{1}|^{2})$ for some exponent $\gamma>0$
that depends on the shape of $p(y)$.

The above results are the continuum analog of the ones obtained in the single
dominant species case: by changing the upper edge of the distribution $p(y)$
one can control the nature of instabilities, from delocalized to partially
concentrated onto a few dominant species. Fig. 3 summarizes this example by
showing the complete spinodal line in each regime. To demonstrate the
localization of the instability direction onto the few dominant species, we
also show the inverse participation ratio (IPR) along the spinodal. The IPR is
defined as $\textrm{IPR}=\sum_{\alpha}v_{\alpha}^{4}$ and so expected to be
$O(M^{-1})$ for delocalized instabilities, while it reaches $O(1)$ whenever
the instability is concentrated on a few species.

Here again one observes composition imbalance amplification in the CD regime:
if purely entropic effects were at play, the components of the instability
direction $\bm{v}$ should be distributed according to $p(y)$. This would yield
a much lower IPR ($O(1/M)$ rather than $O(1)$) than we find, cf. Fig. 3. The
effect again requires $s>0$, i.e. complexity in the underlying interactions.

Finally, we illustrate in Fig. 3.c how the presence of systematic interactions
creates further competition between different forms of demixing. To illustrate
this, we consider second virial coefficients of the form
$\epsilon_{\alpha\gamma}=-b+s\eta_{\alpha\gamma}-\sigma_{\alpha}\sigma_{\gamma}$,
with an additional systematic term parameterized by an interaction strength
$\sigma_{\alpha}$ associated with each species. We write
$\langle\sigma\rangle$ and $\textrm{Var}_{\sigma}$ for the mean and variance
of the $\sigma_{\alpha}$ and introduce the parameter
$\delta=\tilde{s}/\sqrt{\textrm{Var}_{\sigma}+\tilde{s}^{2}}$ measuring the
relative strength of the random and systematic interactions. As density is
increased, for sufficiently small $\delta$, the systematic interaction now
dominates the instability, cf. Fig. 3.c. The inclusion of any finite number of
similarly structured interactions into our framework is straightforward and we
discuss the extension of existing models Graf and Machta (2022); Carugno  _et
al._ (2022) in the Suppplemental Material .

To conclude, we have provided a framework for understanding instabilities in
mixtures of arbitrary composition, allowing also for the straightforward
inclusion of volume exclusion and systematic interactions. In particular we
have obtained an exact equation for the spinodal line in the limit of many
components, $M\gg 1$, for mixtures with complex interactions as proposed by
Sear and Cuesta Sear and Cuesta (2003). In simple yet paradigmatic examples we
showed that a small number of higher-density mixture components can strongly
control the nature of instabilities through a surprising interplay between
entropic effects and the complexity of interactions, resulting in a strong
composition imbalance amplification. This new form of instability is the main
physical insight of the present letter. Since in many biological mixtures
different components are present in different amounts, we expect instabilities
to phase separation in such systems to be strongly dictated by the components
with the highest concentration, strongly influencing the structures and phases
found in the steady state. Our results thus also point to a new route for
biological systems to control patterns of phase separation by fine-tuning
mixture composition imbalances.

##  I Acknowledgments

This work was supported by the German Research Foundation (DFG) under grant
numbers SO 1790/1-1 and KR 3844/5-1.

## References

  * Alberti (2017) S. Alberti, Current Biology 27, R1097 (2017).
  * Lafontaine  _et al._ (2020) D. L. Lafontaine, J. A. Riback, R. Bascetin, and C. P. Brangwynne, Nature Reviews Molecular Cell Biology (2020), 10.1038/s41580-020-0272-6.
  * Sear (2007) R. P. Sear, Soft Matter 3, 680 (2007).
  * Berry  _et al._ (2018) J. Berry, C. P. Brangwynne, and M. Haataja, Reports on Progress in Physics 81 (2018), 10.1088/1361-6633/aaa61e.
  * Choi  _et al._ (2020) J. M. Choi, A. S. Holehouse, and R. V. Pappu, Annual Review of Biophysics 49, 107 (2020).
  * Keating (2012) C. D. Keating, Accounts of Chemical Research 45, 2114 (2012).
  * Jacobs (2021) W. M. Jacobs, Physical Review Letters 126, 258101 (2021).
  * Zwicker and Laan (2022) D. Zwicker and L. Laan,  , 1 (2022), arXiv:2201.10898 .
  * Jacobs and Frenkel (2013) W. M. Jacobs and D. Frenkel, Journal of Chemical Physics 139 (2013), 10.1063/1.4812461.
  * Jacobs and Frenkel (2017) W. M. Jacobs and D. Frenkel, Biophysical Journal 112, 683 (2017), arXiv:1703.01223 .
  * Sear and Cuesta (2003) R. P. Sear and J. A. Cuesta, Physical Review Letters 91, 1 (2003), arXiv:0307326 [cond-mat] .
  * Shrinivas and Brenner (2021) K. Shrinivas and M. P. Brenner, Proceedings of the National Academy of Sciences of the United States of America 118, 1 (2021).
  * Sear (2005) R. P. Sear, Journal of Physics Condensed Matter 17 (2005), 10.1088/0953-8984/17/45/052, arXiv:0507517 [cond-mat] .
  * Gasic  _et al._ (2021) A. G. Gasic, A. Sarkar, and M. S. Cheung, Physical Review Research 3 (2021), 10.1103/physrevresearch.3.033220.
  * Neidhardt and Umbarger (1987) F. Neidhardt and H. Umbarger, E. coli and S. typhimurium: Cellular and Molecular Biology (1987).
  * Weber and Brangwynne (2015) S. C. Weber and C. P. Brangwynne, Current Biology 25, 641 (2015).
  * Goehring and Hyman (2012) N. W. Goehring and A. A. Hyman, Current Biology 22, R330 (2012).
  * Riback  _et al._ (2020) J. A. Riback, L. Zhu, M. C. Ferrolino, M. Tolbert, D. M. Mitrea, D. W. Sanders, M. T. Wei, R. W. Kriwacki, and C. P. Brangwynne, Nature 581, 209 (2020).
  * Zhang  _et al._ (2015) H. Zhang, S. Elbaum-Garfinkle, E. M. Langdon, N. Taylor, P. Occhipinti, A. A. Bridges, C. P. Brangwynne, and A. S. Gladfelter, Molecular Cell 60, 220 (2015).
  * Fritsch  _et al._ (2021) A. W. Fritsch, A. F. Diaz-Delgadillo, O. Adame-Arana, C. Hoege, M. Mittasch, M. Kreysing, M. Leaver, A. A. Hyman, F. Jülicher, and C. A. Weber, bioRxiv ., 2021.02.11.430794 (2021).
  * Weber  _et al._ (2019) C. A. Weber, D. Zwicker, F. Jülicher, and C. F. Lee, Reports on Progress in Physics 82 (2019), 10.1088/1361-6633/AB052B.
  * Sollich (2002) P. Sollich, Journal of Physics Condensed Matter 14 (2002), 10.1088/0953-8984/14/3/201, arXiv:0109292 [cond-mat] .
  * Akemann  _et al._ (2011) G. Akemann, J. Baik, and P. Di Francesco,  _The Oxford handbook of random matrix theory_ (Oxford University Press, 2011).
  * Note (1) Only if the entries of $\bm{D}$ are relatively close to each other. Otherwise, the spectrum may split intro multiple bulk pieces.
  * Hwang (2004) S.-G. Hwang, The American mathematical monthly 111, 157 (2004).
  * Benaych-Georges and Nadakuditi (2011) F. Benaych-Georges and R. R. Nadakuditi, Advances in Mathematics 227, 494 (2011), arXiv:0910.2120 .
  * Voiculescu  _et al._ (1992) D. V. Voiculescu, K. J. Dykema, and A. Nica,  _Free random variables_ , 1 (American Mathematical Soc., 1992).
  * Bouchbinder  _et al._ (2021) E. Bouchbinder, E. Lerner, C. Rainone, P. Urbani, and F. Zamponi, Phys. Rev. B 103, 174202 (2021).
  * Moran and Bouchaud (2019) J. Moran and J.-P. Bouchaud, Phys. Rev. E 100, 032307 (2019).
  * Note (2) Sear and Cuesta Sear and Cuesta (2003) omitted the explicit entropic solvent contribution, resulting in the spinodal curve always being linear in $\rho$. Including the solvent entropy as done in other recent studies Jacobs and Frenkel (2017); Jacobs (2021); Shrinivas and Brenner (2021) makes $\theta$ dependent on $T$ and $\rho$ and gives a nonlinear spinodal curve in the condensation regime.
  * Lee and Schnelli (2016) J. O. Lee and K. Schnelli, Probability Theory and Related Fields 164, 165 (2016), arXiv:1310.7057 .
  * Graf and Machta (2022) I. R. Graf and B. B. Machta, Physical Review Research 4, 033144 (2022).
  * Carugno  _et al._ (2022) G. Carugno, I. Neri, and P. Vivo, Physical Biology (2022).