# Non-Markovian stochastic Gross-Pitaevskii equation for the exciton-polariton
Bose-Einstein condensate

A. D. Alliluev  V. I. Il’ichev Pacific Oceanological Institute, Far East
Branch of the Russian Academy of Sciences, Vladivostok 690041, Russia  D.V.
Makarov  V. I. Il’ichev Pacific Oceanological Institute, Far East Branch of
the Russian Academy of Sciences, Vladivostok 690041, Russia  N.A. Asriyan
N.L. Dukhov Research Institute of Automatics (VNIIA), Moscow 127030, Russia
A.A. Elistratov  N.L. Dukhov Research Institute of Automatics (VNIIA), Moscow
127030, Russia  Yu. E. Lozovik  Institute for Spectroscopy RAS, Troitsk
108840, Moscow, Russia  MIEM, National Research University Higher School of
Economics, Moscow 101000, Russia

###### Abstract

In this paper, a non-Markovian Gross-Pitaevskii equation is proposed to
describe the formation of a condensate in an exciton-polariton system under
incoherent pumping. By introducing spatially delta-correlated noise terms, we
observe a transition from a spatially ordered phase to a disordered one as the
temperature increases. In course of this process, the population of the
condensate is significantly reduced. Irregularly located separate dense spots
of condensate above the transition temperature are revealed. Using the Gabor
transform, it is shown that, with increasing temperature, the condensate
decoheres, that is accompanied by the transition from narrowband to broadband
spectral density.

††preprint: APS/123-QED

##  I Introduction

Open quantum systems attract much interest in recent years, which is largely
due to the rapid development of quantum computing technologies. A common
approach for studying such systems is the Markov approximation i.e. using the
Lindblad equation for the density matrix or the stochastic Schrödinger
equation [1, 2].

When dealing with a condensate of quasiparticles, we can consider it as an
open quantum system linked to a reservoir of noncondensed particles. For a
theoretical description, the stochastic Gross-Pitaevskii equation can be used
[3, 4]. However, due to the significant development of experimental methods,
the reservoir may have a relatively low temperature as well as a narrow energy
spectrum. This makes Markovian approximation invalid for such systems. Non-
Markovian evolution infers existence of memory effects in the system. The
corresponding evolution equation is an integro-differential one with time non-
local terms [5]. For instance, Nakajima-Zwanzig [6, 7] equation or non-
Markovian versions of the Lindblad equation [8, 9] are used for the density
matrix. Though the non-Markovian Schrödinger equation for non-interacting
particles was first introduced in [10], its generalisation for systems with
interparticle interaction is still a challenge. Significant progress has been
achieved in [11], where, using the Keldysh diagram technique, a system of non-
Markovian equations for an exciton-polariton Bose-Einstein condensate in a
Fabry-Perot optical microresonator has been derived.

An exciton-polariton is a bound state of an exciton and a photon, the lifetime
of which is significantly affected by the quality factor of the
microresonator. Exciton-polariton condensation arises as a result of their
thermalization stimulated by laser pumping. Exciton polaritons are attractive
for technological applications largely because of the prospects for their use
in quantum computing applications [12], as well as in the context of the
polariton laser [13].

In [11] it was assumed that the condensate corresponds to a macroscopically
populated state with zero momentum. Therefore, the resulting system of
equations can be treated a discrete non-Markovian Gross-Pitaevskii equation.
Unfortunately, a rigorous derivation of a similar Gross-Pitaevskii equation
with spatial dynamics taken into account remains challenging. However, as a
reasonable approximation, we can use the evolution equation from [11] a small
neighborhood of zero momentum. It corresponds to a smoothly varying spatial
profile of the macroscopic wave function of the condensate. Under this
assumption one can propose a non-Markovian version of the Gross-Pitaevskii
equation for the exciton-polariton condensate.

The paper is organized as follows. In the next section we introduce the non-
Markovian stochastic Gross-Pitaevskii equation. The process of condensate
formation is studied in Section III by means of numerical simulation. In
Conclusion we summarize and discuss the results obtained.

##  II Evolution equations

When dealing with exciton-polariton condensation one may usually neglect the
population of the upper spectral branch (upper polariton states) [14]. It
leads to the following form of the Gross-Pitaevskii equation:

$\displaystyle i\hbar\frac{\partial\psi(\mathbf{r},t)}{\partial
t}=\hat{H}_{0}\psi(\mathbf{r},t)+P_{\mathrm{coh}}(\mathbf{r},t)+$ (1)
$\displaystyle\hat{D}_{\text{cav}}\psi(\mathbf{r},t)+\hat{D}_{\text{ex}}\psi(\mathbf{r},t),$

where $\hat{H}_{0}$ is an unitary operator including kinetic and potential
energies, and energy corrections due to interparticle interactions,

$\hat{H}_{0}=-\frac{\hbar^{2}}{2m^{*}}\nabla^{2}+V(\mathbf{r})+\alpha_{\mathrm{c}}|\psi(\mathbf{r},t)|^{2}+\alpha_{\mathrm{r}}\rho_{\text{r}}(\mathbf{r},t).$
(2)

Here $m^{*}$ stands for the effective polariton mass, which is approximately
$10^{-5}$–$10^{-4}$ $m_{\text{e}}$, $V(\mathbf{r})$ is the external potential,
$\alpha_{\text{c}}$ is the interaction constant for condensate polaritons,
$\alpha_{\text{r}}$ describes condensate-reservoir interaction,
$\rho_{\text{r}}(\mathbf{r},t)$ is the exciton reservoir density. The function
$P_{\mathrm{coh}}(\mathbf{r},t)$ in (1) describes coherent laser pumping of
microresonator photons.

Non-Hermitian operators $\hat{D}_{\text{cav}}$ and $\hat{D}_{\text{ex}}$
correspond to interactions with the photon and exciton reservoirs
respectively. Hereafter we assume that the interaction with the photon
reservoir is Markovian. This assumption leads to

$\hat{D}_{\mathrm{cav}}\psi=-i\gamma_{\mathrm{cav}}\psi+\eta_{\mathrm{cav}}(\mathbf{r},t).$
(3)

Using the truncated Wigner approximation, one can derive the following
expression for the autocorrelation function:

$\braket{\eta^{*}_{\mathrm{cav}}(\mathbf{r},t)\eta_{\mathrm{cav}}(\mathbf{r^{\prime}},t^{\prime})}=\frac{\gamma_{\mathrm{cav}}}{\Delta
x\Delta y}\delta(\mathbf{r}-\mathbf{r^{\prime}})\delta(t-t^{\prime}),$ (4)

with $\Delta x$ and $\Delta y$ being the grid cell sizes.

The non-Markovian behaviour of the system is mainly associated with the
interaction with the exciton reservoir. Considering only the low-momentum
states as a condensate, we can use the following approximation

$\hat{D}_{\text{ex}}\psi(\mathbf{r})\simeq\hbar\int\limits_{0}^{t}dt^{\prime}\Sigma^{\text{R}}(t,t^{\prime})\psi(\mathbf{r},t^{\prime})+\eta_{\text{ex}}(\mathbf{r},t),$
(5)

i.e. treat the retarded self-energy term $\Sigma^{\text{R}}(t,t^{\prime})$ as
a spatially homogeneous one. We use expressions derived in [11]

$\Sigma^{\text{R}}(t,t^{\prime})=i\frac{\rho_{\text{r}}^{2}\alpha_{\text{c}}^{2}}{\hbar^{2}}\frac{e^{-\gamma_{\text{ex}}(t-t^{\prime})}}{1+\left[\frac{k_{\text{B}}T}{\hbar}(t-t^{\prime})\right]^{2}}\theta(t-t^{\prime}),$
(6)

with $\theta(t)$ being the Heaviside function, $\gamma_{\text{ex}}$ –
condensate exciton decay rate. The corresponding expression for the Keldysh
component of the self-energy term is as follows:

$\Sigma^{\text{K}}(t,t^{\prime})=-i\frac{\rho_{\text{r}}^{2}\alpha_{\text{c}}^{2}}{\hbar^{2}}\frac{e^{-\gamma_{\text{ex}}|t-t^{\prime}|}}{1+\left[\frac{k_{\text{B}}T}{\hbar}(t-t^{\prime})\right]^{2}}.$
(7)

This expression defines the temporal autocorrelation function of exciton
fluctuations:

$\braket{\eta^{*}(\mathbf{r},t)\eta(\mathbf{r^{\prime}},t^{\prime})}=i\hbar^{2}\delta(\mathbf{r},\mathbf{r^{\prime}})\Sigma^{\text{K}}(t,t^{\prime}).$
(8)

When modeling the spatial structure of exciton fluctuations, we utilize above
a simple approximation, considering them as delta-correlated. Condensate
formation with correlated noise was studied in [15]. Thus, equation (1) takes
the form

$\displaystyle i\hbar\frac{\partial\psi(\mathbf{r},t)}{\partial
t}=\hat{H}_{0}\psi(\mathbf{r},t)-\frac{i\hbar\gamma_{\text{cav}}}{2}\psi(\mathbf{r},t)+P_{\mathrm{coh}}(\mathbf{r},t)$
$\displaystyle+\eta_{\text{cav}}(\mathbf{r},t)+\eta_{\text{ex}}(\mathbf{r},t)+\hbar\int\limits_{0}^{t}dt^{\prime}\Sigma^{\text{R}}(t,t^{\prime})\psi(\mathbf{r},t^{\prime})$
(9)

In our model, the evolution of the density of the exciton reservoir is
described by the equation

$\displaystyle\frac{\partial\rho(\mathbf{r},t)}{\partial t}$
$\displaystyle{=}\frac{1}{\hbar}P_{\text{incoh}}(\mathbf{r},t){-}\gamma_{\text{exR}}\rho(\mathbf{r},t){-}\frac{2}{\hbar}\text{Im}\left[\psi^{*}(\mathbf{r},t)\eta(\mathbf{r},t)\right]-$
$\displaystyle-2\text{Im}\left[\psi^{*}(\mathbf{r},t^{\prime})\int\,dt^{\prime}\Sigma^{\text{R}}(t,t^{\prime})\psi(\mathbf{r},t^{\prime})\right],$
(10)

where the term $P_{\text{incoh}}(\mathbf{r},t)$ describes the incoherent
pumping of the reservoir, $\gamma_{\text{exR}}$ is the decay rate of the
reservoir excitons. The third term on the right-hand side describes polariton
exchange between the condensate and the reservoir via the fluctuations.

In [16] we have demonstrated that at low temperatures an exponential
approximation can be used both for (6) and (7):

$\displaystyle\Sigma^{\text{R}}(t,t^{\prime})$ $\displaystyle\simeq
i\frac{\rho_{\text{r}}^{2}\alpha_{\text{c}}^{2}}{\hbar^{2}}e^{-\gamma_{\text{eff}}(t-t^{\prime})}\theta(t-t^{\prime}),$
(11) $\displaystyle\Sigma^{\text{K}}(t,t^{\prime})$
$\displaystyle\simeq-i\frac{\rho_{\text{r}}^{2}\alpha_{\text{c}}^{2}}{\hbar^{2}}e^{-\gamma_{\text{eff}}|t-t^{\prime}|}.$
(12)

Moreover, for low temperatures, $\gamma_{\text{eff}}$ is linearly dependent on
temperature [15]. In this case, the autocorrelation function (8) corresponds
to the complex-valued stochastic Ornstein-Uhlenbeck process, when
$\tilde{\eta}(\mathbf{r},t)$ is the solution of the following Langevin
equation:

$\frac{d\tilde{\eta}(\mathbf{r},t)}{dt}=-\gamma_{\text{eff}}\tilde{\eta}(\mathbf{r},t)+\sqrt{2\gamma_{\text{eff}}}\xi(\mathbf{r},t).$
(13)

Here $\xi(\mathbf{r},t)$ is a complex white noise with unit variance,

$\braket{\xi^{*}(\mathbf{r},t)\xi(\mathbf{r^{\prime}},t^{\prime})}=\delta(t-t^{\prime})\delta(\mathbf{r}-\mathbf{r^{\prime}}).$
(14)

The exponential form of the memory kernel of the evolution equation (9) makes
if possible to transform it into an equivalent Markov equation by introducing
an auxiliary wave function

$\phi(\mathbf{r},t)=\psi_{0}(\mathbf{r})e^{-\gamma_{\text{eff}}t}+\int\limits_{0}^{t}\,dt^{\prime}e^{-\gamma_{\text{eff}}(t-t^{\prime})}\psi(\mathbf{r},t^{\prime}),$
(15)

where $\psi_{0}(\mathbf{r})=\psi(\mathbf{r},t=0)$. This technique is known as
Markov embedding [5, 17]. The equations (9) and (10) are now expressed as:

$\displaystyle i\hbar\frac{\partial\psi(\mathbf{r},t)}{\partial t}$
$\displaystyle=\hat{H}_{0}\psi(\mathbf{r},t)-\frac{i\hbar\gamma_{\text{cav}}}{2}\psi(\mathbf{r},t)+P_{\mathrm{coh}}(\mathbf{r},t)+\eta(\mathbf{r},t)+i\frac{\alpha_{\text{c}}^{2}\rho_{\text{r}}^{2}(\mathbf{r},t)}{\hbar\gamma_{\text{eff}}}\left[\phi(\mathbf{r},t)-\psi_{0}(\mathbf{r})e^{-\gamma_{\text{eff}}t}\right],$
(16) $\displaystyle\frac{\partial\rho_{text{r}}(\mathbf{r},t)}{\partial t}$
$\displaystyle=\frac{1}{\hbar}P_{{\text{incoh}}}(\mathbf{r})-\gamma_{\text{exR}}\rho_{\text{r}}(\mathbf{r},t)-2\frac{\alpha_{\text{c}}^{2}\rho_{\text{r}}^{2}(\mathbf{r},t)}{\hbar^{2}}\text{Re}\left\\{\psi^{*}(\mathbf{r},t)\left[\phi(\mathbf{r},t)-\psi_{0}(\mathbf{r})e^{-\gamma_{\text{eff}}t}\right]\right\\}-\frac{2}{\hbar}\text{Im}\left\\{\psi^{*}(\mathbf{r},t)\eta(\mathbf{r},t)\right\\}.$
(17)

These equations have to be supplemented with the evolution equation for the
auxiliary wave function $\phi$:

$\frac{\partial\phi(\mathbf{r},t)}{\partial
t}=\gamma_{\text{eff}}\left[\psi(\mathbf{r},t)-\phi(\mathbf{r},t)\right].$
(18)

##  III Numerical simulation

Figure 1: Spatial distribution of condensate density at time $t=130$ ps.
Temperature values: (a) 5 K, (b) 20 K, (c) 35 K, (d) 50 K. All data presented
correspond to a single realization of the fluctuation field.

In the present paper, the system of equations (16)–(18) is used to study
condensate formation in a polariton gas exposed to incoherent pumping, i.e.,
in the absence of direct coherent photon pumping into the microcavity,
$P_{\text{coh}}=0$. We consider the case of constant pumping:

$P_{\text{incoh}}(\mathbf{r})=\gamma_{\text{exR}}\rho_{0}w(\mathbf{r}),$ (19)

with the function $w(\mathbf{r})$ defining the spot profile

$w(\mathbf{r})=\exp\left[-\left(\frac{\mathbf{r}-\mathbf{r_{\text{c}}}}{\sigma_{\text{{r}}}}\right)^{2}\right].$
(20)

We place the pump spot center at the origin, $\mathbf{r_{\text{c}}}=0$. The
value $\rho_{0}$ represents the maximum reservoir density in the equilibrium
state. We use $\rho_{0}=0.5\times 10^{12}$ $cm^{-2}$ for numerical simulation.
The parameter $\sigma_{\text{r}}=20$ $\mu m$ sets the half-width of the pump
spot. The values of $\gamma_{\text{cav}}$, $\gamma_{\text{ex}}$ and
$\gamma_{\text{exR}}$ and are determined by the corresponding lifetimes:
$\tau_{\text{cav}}=1/\gamma_{\text{cav}}=3.8$ ps,
$\tau_{\text{ex}}=1/\gamma_{\text{ex}}=1$ ps,
$\tau_{\text{exR}}=1/\gamma_{\text{exR}}=10$ ps. The interexciton interaction
constant $\alpha_{\text{c}}$ is set to $6\cdot 10^{-14}$ eV$\cdot$cm2, the
size of the grid cell $\Delta x=\Delta y=0.5$ $\mu m$. We consider the case of
zero initial conditions

$\psi(t=0)=\phi(t=0)=\rho_{\text{r}}(t=0)=0,$ (21)

i.e., photon and exciton fluctuations play the role of the ”seeds” for
condensate nucleation. In this case the condensate formation develops a
process of self-organization accompanied by gradual increasing of condensate
healing length. It is reasonable to assume that the dynamic memory in the
Gross-Pitaevskii equation should contribute to this process. Since memory time
is determined by the reservoir temperature, we expect significant influence of
the latter on the dynamics of the system.

Figure 2: Spatial distribution of the condensate phase absolute value at time
$t=130$ ps. The temperature values are: (a) 5 K, (b) 20 K, (c) 35 K, (d) 50 K.
All data presented correspond to a single realization of the fluctuation
field.

Results of numerical simulation confirm the latter assumption. Figure 1 shows
condensate density distributions for different temperatures and individual
realisations of $\eta_{\text{cav}}(\mathbf{r},t)$ and
$\eta_{\text{ex}}(\mathbf{r},t)$. All cases correspond to $t=130$ ps. By this
time, all transients have finished, and the system reaches a quasi-equilibrium
state.

As follows from the data presented in Fig. 1, the structure of this state is
qualitatively temperature dependent. At $5$ and $20$ K, the condensate
occupies a vast region, with the spatial density distribution roughly
repeating the intensity distribution for the incoherent pump beam. At $35$ and
$50$ K, however, we observe the condensate as individual spots of irregular
shape. This indicates that dynamic memory, the duration of which increases
with decreasing temperature, plays a significant role in establishing the
spatial correlations of the condensate. The phenomenon of condensate
fragmentation with increasing temperature was previously discussed in [15].

Destruction of spatial correlations of the condensate with increasing
temperature is more clearly demonstrated in Fig. 2, which depicts spatial
condensate phase distributions for the same realizations as in Fig. 1. At
$T=5$ K, the lines of constant phase take the form of concentric rings,
somewhat distorted due to the presence at the periphery of several singular
points corresponding to the vortex cores. Such a phase configuration indicates
the presence of matter waves being emitted from the pump region with a fairly
regular wave front shape. At $T=20$ K the wave front undergoes significant
spatial distortions due to the influence of spatial condensate fluctuations.
With a further increase in temperature, the contribution of fluctuations
increases, which leads eventually to the destruction of the ring phase
configuration. At $T=35$ K traces of phase coherence are preserved only for
individual condensate fragments, and at $T=50$ K they almost completely
disappear.

The behavior of individual statistical realizations sheds light on the
qualitative dependence of the condensate dynamics on the temperature, though
statistical modeling is yet necessary.

As one of the indicators of the condensate state, let us consider condensate
density defined as

$\bar{\rho}=\frac{1}{\pi r_{\text{c}}^{2}}\int
d\mathbf{r}f(\mathbf{r})|\psi(\mathbf{r})|^{2}\,$ (22)

where the function $f(\mathbf{r})$ extracts the regions of the pump spot with
the highest density,

$f(\mathbf{r})=\Biggl{\\{}\Biggr{.}\begin{aligned}
1,\quad|\mathbf{r}-\mathbf{r}_{\text{c}}|\leq\sigma,\\\
0,\quad|\mathbf{r}-\mathbf{r}_{\text{c}}|>\sigma.\end{aligned}$ (23)

Figure 3: Time dependencies of ensemble-averaged density $\bar{\rho}$ (panel
(a)) and the scintillation index (panel (b)) for different temperature values.

We use also the scintillation index as a characteristics of the spatial
condensate distribution

$\text{SI}=\frac{\braket{\bar{\rho}^{2}}}{\braket{\bar{\rho}}^{2}}-1.$ (24)

The time dependence of $\bar{\rho}$ and scintillation index are presented in
Fig. 3. First of all, we note that the low-temperature phase-ordered states (5
K and 20 K) are described by significantly high densities and very weak
spatial fluctuations. The sharp decrease in density that occurs with
increasing temperature is accompanied by enhancement of fluctuations: we
observe a pattern corresponding to the Berezinskii–Kosterlitz–Thouless(BKT)
phase transition. At temperature $T=35$ K the condensate exists in the form of
separate spots of high density on a turbulent background. This leads to high
values of scintillation index, $\text{SI}\simeq 2$. At $T=50$ K the
contribution of such spots significantly decreases. As a consequence, the
scintillation index approaches unity, which corresponds to the statistical
saturation regime [20]. Also note that even after the condensate density
stabilizes, it continues to exhibit slow oscillations. This may be due to the
so-called relaxation oscillations [16, 21, 22, 23].

The transition of the condensate to a quasi-equilibrium regime with constant
or slowly varying density corresponds to the establishment of a balance
between pumping from the reservoir and decay. The presence of such a balance
may indicate the emergence of a state with PT (parity-time) symmetry [24]. In
this case the energy spectrum becomes real despite the non-Hermiticity of
exciton-polaritons. The appearance of PT symmetry in a system with memory has
been considered earlier in [25].

Figure 4: Optical radiation spectrum of the condensate obtained using the
Gabor transform. Temperature values: (a) 5 K, (b) 20 K, (c) 35 K, (d) 50 K.

In order to trace changes in the spectrum during condensate formation, we
studied the optical radiation of the condensate. To simplify the analysis, we
get rid of the spatial dependence of this signal by using a weighted
integration, with the weight function being given by the pumping profile:

$\Psi_{\text{eff}}(t)=\frac{\int
d\mathbf{r}w(\mathbf{r})\psi(\mathbf{r},t)}{\int d\mathbf{r}w(\mathbf{r})}.$
(25)

To obtain the spectrum of the signal, we use the Gabor transform, which is a
variation of the Fourier window transform with the Gaussian window:

$g_{W}(\nu,t)=\int_{-\infty}^{\infty}dt^{\prime}\Psi_{\text{eff}}(t^{\prime})\chi_{\text{W}}(\nu,t^{\prime}-t),$

where

$\chi_{\text{W}}(\nu,t)=R_{\text{W}}(t)\exp{(2\pi i\nu t)}$ (26)

and

$R_{\text{W}}(t)=\frac{1}{\sqrt{W\sqrt{2\pi}}}\exp{\left(-\frac{t^{2}}{4W^{2}}\right)}.$
(27)

The window is centered at $t=t^{\prime}$, and $W$ sets its size. By moving the
center of the window along the time axis, one may trace the time dependence of
the frequency spectrum of $\Psi_{\text{eff}}(t)$. Fig. 4 shows the spectral
density of the signal averaged over 100 realizations. One may observe that at
temperatures 5 K and 20 K the condensate is concentrated in low-energy states
within a narrow spectral band. The mean frequency of this band corresponds to
the equilibrium chemical potential. Its increase with increasing temperature
is associated with enhancement of the condensate inhomogeneity and, as a
consequence, higher contribution of kinetic energy. At 35 K, the condensate
spectrum is still given by a fairly narrow band, but the shape and intensity
of this band changes with time. Comparing the evolution of the spectrum with
the data presented in Fig. 3(a), one may conclude that the decay is not
completely compensated by pumping. This may be due to the metastability of
condensate spots arising at 35 K. The spectrum at 50 K is qualitatively
different from the previous cases: instead of excitation of a narrow spectral
band, we observe a broadband spectrum with irregular variability with time.
Such behavior is consistent with the pattern of disordered turbulent
condensate dynamics.

##  IV Conclusion

Let’s summarize the main results of the paper. We have presented a non-
Markovian stochastic Gross-Pitaevskii model for the exciton-polariton Bose-
Einstein condensation. We used it to study the influence of the temperature on
the condensate formation process in the absence of coherent pumping. The
results obtained reveal a the complex structure of the transition from the
coherent to the disordered phase. In particular, one may note the appearance
of dense condensate spots on the turbulent background, detected at 35 K. Due
to the presence of these spots, the spectrum of the optical signal emitted by
the condensate is a narrow-band one. Further work in this field will be
devoted to the improvement of the model used to describe the condensate and
the reservoir.

###### Acknowledgements.

This work was supported by the Basic Research Program of the Institute of
Experimental Problems of the Far East Branch of the Russian Academy of
Sciences, project no. 121021700341-2 and by the foundation for the advancement
of theoretical physics and mathematics “Basis”. The work of Yu. E. Lozovik and
N.A. Asriyan is supported by RFBR projects 20-02-00410 and 21-52-12038.

## References

  * [1] Heinz-Peter Breuer, F. Petruccione, The Theory of Open Quantum Systems ,(2002).
  * [2] A. R. Kolovsky, D. L. Shepelyansky, Ann. Phys., 531, 1900231 (2019).
  * [3] H. T. C. Stoof, M. J. Bijlsma, J. Low. Temp. Phys., 124, 431 (2001).
  * [4] S. P. Cockburn, N. P. Proukakis, Laser Physics, 19, 558 (2009).
  * [5] I. De Vega, D. Alonso, Rev. Mod. Phys., 89, 015001 (2017).
  * [6] S. Nakajima, Progr. Theor. Phys., 20, 948 (1958).
  * [7] R. Zwanzig, J. Chem. Phys., 83, 1338 (1960).
  * [8] H.-P. Breuer, Phys. Rev. A, 70, 012106 (2004).
  * [9] H.-P. Breuer, E.-M. Laine, J. Piilo, B. Vacchini, Rev. Mod. Phys., 88, 021002 (2016).
  * [10] L. Diósi, W. T. Strunz, Phys. Lett. A, 235, 569 (1997).
  * [11] A. A. Elistratov, Yu. E. Lozovik, Phys. Rev. B, 93, 014525 (2018).
  * [12] Y. Xue, I. Chestnov, E. Sedov, E. Kiktenko, A. K. Fedorov, S. Schumacher, X. Ma, A. Kavokin, Phys. Rev. Res., 3, 013099 (2021).
  * [13] S. Kim, B. Zhang, Z. Wang, J. Fischer, S. Brodbeck, M. Kamp, C. Schneider, S. Höfling, H. Deng, Phys. Rev. X, 6, 011026 (2016).
  * [14] S. S. Gavrilov, Physics-Uspekhi, 63 (2): 123 (2020).
  * [15] A. D. Alliluev, D. V. Makarov, N. A. Asriyan, A. A. Elistratov, Yu. E. Lozovik, Phys. Lett. A, 453, 128492 (2022).
  * [16] D. V. Makarov, A. A. Elistratov, Yu. E. Lozovik, Phys. Lett. A, 384, 126942 (2020).
  * [17] X. Li, Phys. Lett. A, 387, 127036 (2021).
  * [18] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. Keeling, F. Marchetti, M. Szymańska, R. André, J. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, Le Si Dang, Nature, 443, 409 (2006).
  * [19] V.B. Timofeev, FTP, 46, 865 (2012).
  * [20] R. Dashen, J. Math. Phys., 20, 894 (1979).
  * [21] M. De Giorgi, D. Ballarini, P. Cazzato, G. Deligeorgis, S.I. Tsintzos, Z. Hatzopou- los, P.G. Savvidis, G. Gigli, F.P. Laussy, D. Sanvitto, Phys. Rev. Lett., 112, 113602 (2014).
  * [22] A. Opala, M. Pieczarka, M. Matuszewski, Phys. Rev. B, 98, 195312 (2018).
  * [23] C. Tian, L. Chen, Y. Zhang, L. Zhu, W. Hu, Y. Pan, Z. Wang, F. Zhang, L. Zhang, Nano Lett., 22, 3026 (2022).
  * [24] I.Yu. Chestnov, S.S. Demirchyan, A.P. Alodjants, Yu.G. Rubo, A.V. Kavokin, Sci. Rep., 22, 3026 (2022).
  * [25] Z.A. Cochran, A. Saxena, Y.N. Joglekar, Phys. Rev. Res., 3, 013135 (2021).


# Inferred Linear Stability of Parker Solar Probe Observations using One- and
Two-Component Proton Distributions

K.G. Klein University of Arizona, Tucson, AZ, USA  J.L. Verniero  Space
Sciences Laboratory, University of California, Berkeley, CA 94720-7450, USA
B. Alterman  Space Science and Engineering, Southwest Research Institute, San
Antonio, TX, USA  S. Bale Physics Department, University of California,
Berkeley, CA 94720-7300, USA  Space Sciences Laboratory, University of
California, Berkeley, CA 94720-7450, USA  The Blackett Laboratory, Imperial
College London, London, SW7 2AZ, UK  School of Physics and Astronomy, Queen
Mary University of London, London E1 4NS, UK  A. Case  Smithsonian
Astrophysical Observatory, Cambridge, MA, USA  J.C. Kasper  Smithsonian
Astrophysical Observatory, Cambridge, MA, USA  K. Korreck  Smithsonian
Astrophysical Observatory, Cambridge, MA, USA  D. Larson  Space Sciences
Laboratory, University of California, Berkeley, CA 94720-7450, USA  E. Lichko
University of Arizona, Tucson, AZ, USA  R. Livi  Space Sciences Laboratory,
University of California, Berkeley, CA 94720-7450, USA  M. McManus  Space
Sciences Laboratory, University of California, Berkeley, CA 94720-7450, USA
M. Martinović  University of Arizona, Tucson, AZ, USA  A. Rahmati  Space
Sciences Laboratory, University of California, Berkeley, CA 94720-7450, USA
M. Stevens  Smithsonian Astrophysical Observatory, Cambridge, MA, USA  P.
Whittlesey  Space Sciences Laboratory, University of California, Berkeley, CA
94720-7450, USA

###### Abstract

The hot and diffuse nature of the Sun’s extended atmosphere allows it to
persist in non-equilibrium states for long enough that wave-particle
instabilities can arise and modify the evolution of the expanding solar wind.
Determining which instabilities arise, and how significant a role they play in
governing the dynamics of the solar wind, has been a decades-long process
involving in situ observations at a variety of radial distances. With new
measurements from Parker Solar Probe (PSP), we can study what wave modes are
driven near the Sun, and calculate what instabilities are predicted for
different models of the underlying particle populations. We model two hours-
long intervals of PSP/SPAN-i measurements of the proton phase-space density
during PSP’s fourth perihelion with the Sun using two commonly used
descriptions for the underlying velocity distribution. The linear stability
and growth rates associated with the two models are calculated and compared.
We find that both selected intervals are susceptible to resonant
instabilities, though the growth rates and kind of modes driven unstable vary
depending on if the protons are modeled using one or two components. In some
cases, the predicted growth rates are large enough to compete with other
dynamic processes, such as the nonlinear turbulent transfer of energy, in
contrast with relatively slower instabilities at larger radial distances from
the Sun.

solar wind — plasmas — instabilities — Sun: corona

##  1 Introduction

Wave-particle interactions are suspected of affecting the evolution of the
solar wind as it is accelerated from the Sun’s surface and expands into the
heliosphere; c.f. reviews in Matteini et al. (2012); Yoon (2017); Verscharen
et al. (2019). Such instabilities are driven by departures from local
thermodynamic equilibrium (LTE) that are frequently modeled using velocity
distributions with anisotropic temperatures $T_{\perp,j}$ and
$T_{\parallel,j}$ with respect to local magnetic field $\mathbf{B}$, relative
field-aligned drifts between constituent plasma populations $\Delta
v_{i,j}=(\mathbf{V}_{i}-\mathbf{V}_{j})\cdot\mathbf{B}/|\mathbf{B}|$, and
temperature disequilibrium between species $T_{i}\neq T_{j}$. The simultaneous
effects of multiple sources of free energy can complicate a simple linear
analysis; for instance, it has been found that the free energy contributions
to unstable behavior from different ion and electron species can be non-
negligible (Chen et al., 2016). To address this difficulty, previous works
have applied a numerical implementation of the Nyquist instability criterion
(Nyquist, 1932; Klein et al., 2017) to selected solar wind observations from
the Wind (Klein et al., 2018) and Helios (Klein et al., 2019) missions,
finding that a majority of intervals were unstable, including many intervals
that simple parametric models accounting for a single source of free energy
would have predicted to be stable.

Given the complexity of phase-space distributions typically found in weakly
collisional plasmas, a number of different schemes for modeling the underlying
velocity-space structure are frequently used; for instance, it is common to
treat the protons as a single, anisotropic bi-Maxwellian or kappa
distribution, or as a linear combination of core and relatively drifting beam
distributions, each with distinct parallel and perpendicular temperatures; see
the introduction of Alterman et al. (2018) for a review of solar wind
observations of secondary ion populations.

In this work, we select two hours-long time intervals observed by the SPAN-i
instrument from the SWEAP instrument suite (Kasper et al., 2015) on Parker
Solar Probe (PSP) (Fox et al., 2015) during its fourth encounter with the Sun,
where significant ion-scale wave activity is observed, similar to activity
previously reported in Bowen et al. (2020) and Verniero et al. (2020). We
generate both a one-component and two-component model for each measurement of
the proton velocity distribution, calculating and comparing the associated
linear stability. Using the different models produces significantly different
instabilities, either in the robustness of the associated growth rates or the
kinds of waves driven unstable. The two-component model generally predicts
ion-scale waves with characteristics more in line with the observed wave
activity than models using a single proton component. This suggests that using
overly simplistic models for ion distributions may neglect essential kinetic-
scale processes responsible for the generation of these waves, even if these
models capture macroscopic departures from LTE.

##  2 Data and Methodology

###  2.1 Parker Solar Probe Data

We select two hours-long sections from the outbound pass of PSP’s fourth
encounter with the Sun, when SPAN-i had sufficient coverage of the proton
velocity distribution to model $f_{p}(\mathbf{v})$, specifically Selection A:
2020/01/30 11:00-13:30 (SA, Fig. 1) and Selection B: 2020/02/01 00:10-02:00
(SB, Fig. 2). During both selections, ion-scale electromagnetic waves are
observed by the FIELDS instrument suite (Bale et al., 2016). Figs. 1 and 2
show the vector magnetic field components, as well as the trace power spectral
density normalized to an ansatz power-law distribution for the background
turbulent spectrum of $f^{-5/3}$, and the polarization of the transverse
components of the magnetic fields, where red (blue) indicates right-handed
(left-handed) circular polarization in the spacecraft frame. In SA, we see an
abundance of power above a $f^{-5/3}$ spectrum persist for several hours near
$3$ Hz. At the same frequencies, we see a clear signature (red) of right-hand
polarization persist for nearly the entire duration of the more than two-hour
selection.

Unlike in SA, in SB there is not a persistent signature at a nearly constant
frequency of ion-scale waves of a single handedness; both left-handed (blue)
and right-handed (red) polarized waves are observed. There are also times
during SB where no enhanced wave activity near ion frequencies is observed.

Figure 1: Magnetic field characteristics observed by FIELDS/PSP during
Selection A, 2020/01/30 11:00-13:30. Top row: vector components of
$\mathbf{B}$. Second row: Trace power spectral density normalized by
$k^{-5/3}$ power law. Third row: Polarization of transverse magnetic field
components, where red indicates right-handed circular polarization in the
spacecraft frame.

Figure 2: Magnetic field characteristics observed by FIELDS/PSP during
Selection B, 2020/02/01 00:10-02:00, organized in the same fashion as Fig. 1.

###  2.2 One- and Two-Component Proton Distributions

For each $\approx 7$ second measurement where a significant fraction of the
thermal proton distribution is in the SPAN-i field of view, a two-component
fit of the observed proton energy and angle spectra is attempted, modeling the
protons as a combination of two relatively drifting bi-Maxwellian
distributions,

$f_{p}^{\textrm{2-comp.}}(v_{\perp},v_{\parallel})=\sum\limits_{j=c,b}\frac{n_{j}}{\pi^{3/2}w_{\perp,j}^{2}w_{\parallel,j}}\exp\left[-\frac{v_{\perp}^{2}}{w_{\perp,j}^{2}}-\frac{\left(v_{\parallel}-V_{j}\right)^{2}}{w_{\parallel,j}^{2}}\right].$
(1)

Parallel and perpendicular are defined with respect to the local mean-magnetic
field direction, $n_{j}$ is the component density, $V_{j}$ the component bulk
speed, and $w_{\perp,\parallel;j}=\sqrt{2T_{\perp,\parallel;j}/m_{j}}$ the
component thermal velocities. This fit represents our two-component model. To
mitigate the partial FOV coverage of SPAN-i, all fitted densities were
calibrated to QTN densities. All calculations using this model are performed
in the proton center-of-mass frame.

For a model with the same macroscopic thermodynamic quantities, i.e. total
proton density as well as parallel and perpendicular thermal pressures, that
are used in a linear instability calculation that does not represent the beam-
and-core structure of the protons observed in the inner heliosphere, we
construct a one-component model as

$f_{p}^{\textrm{1-comp.}}(v_{\perp},v_{\parallel})=\frac{n_{p}}{\pi^{3/2}w_{\perp,p}^{2}w_{\parallel,p}}\\\
\exp\left[-\frac{v_{\perp}^{2}}{w_{\perp,p}^{2}}-\frac{v_{\parallel}^{2}}{w_{\parallel,p}^{2}}\right].$
(2)

Here, the proton density is $n_{p}=n_{c}+n_{b}$ and the total thermal
velocities are $w_{\perp,\parallel;p}=\sqrt{2T_{\perp,\parallel;p}/m_{p}}$. We
have defined the perpendicular proton temperature as

$T_{\perp,p}=\frac{n_{c}T_{\perp,c}+n_{b}T_{\perp,b}}{n_{c}+n_{b}}$ (3)

and the parallel proton temperature as

$T_{\parallel,p}=\frac{n_{c}T_{\parallel,c}+n_{b}T_{\parallel,b}+\left(\frac{n_{c}n_{b}}{n_{c}+n_{b}}\right)m_{p}\Delta
v_{cb}^{2}}{n_{c}+n_{b}}.$ (4)

We emphasize that this is not equivalent to fitting the measured proton VDF
with a single bi-Maxwellian distribution. Our method is employed so that both
models have the same macroscopic perpendicular and parallel proton pressures,
which would not necessarily be the case for a single bi-Maxwellian fit of
protons with a significant secondary population.

The parameters from both models, along with measurements of the magnetic field
strength averaged to the SPAN-i measurement cadence, are combined into the
dimensionless parameters used as inputs for the Nyquist instability analysis.
We will see that the significant differences in the underlying proton phase-
space densities for the two models lead to significant differences in the
predicted unstable behavior.

###  2.3 Instability Analysis

We employ a numerical implementation of the Nyquist instability criterion
(Nyquist, 1932; Klein et al., 2017) for the hot plasma dispersion relation for
an arbitrary number of relatively drifting bi-Maxwellian components as
determined by the PLUME numerical dispersion solver (Klein & Howes, 2015). The
Nyquist criterion determines the stability of a linear system of equations
through a conformal mapping of the contour integral of a dispersion relation
$\mathcal{D}(\omega,\mathbf{k},\mathcal{P})$ over the upper-half of the
complex frequency plane. This integral counts the number of normal mode
solutions that are unstable, having $\gamma>0$, for a specific wavevector
$\mathbf{k}$ and set of dimensionless parameters $\mathcal{P}$;
$\omega_{\textrm{r}}$ and $\gamma$ are the real and imaginary components of
the complex frequency $\omega$. Iterating this process for multiple contours
with increasing values of $\gamma$ enables the determination of the maximum
growth rate and associated characteristics of the fastest growing mode
supported by a particular $\mathbf{k}$. We have set $\gamma=10^{-4}\Omega_{p}$
as the minimum growth rate for a wavevector to be considered unstable. We
repeat this process over a log-spaced grid in wavevector space
$k_{\perp}\rho_{p}\in[10^{-3},3]$ and $k_{\parallel}\rho_{p}\in[10^{-2},3]$,
enabling the determination of the fastest growing mode for all wavevectors
given a particular parameter set $\mathcal{P}$.

For the one-component model, the set of dimensionless plasma parameters is

$\mathcal{P}_{\textrm{1-comp}}=\left(\beta_{\parallel,p},\frac{w_{\parallel,p}}{c},\frac{T_{\perp,p}}{T_{\parallel,p}}\right)$
(5)

while for the two-component model, the dimensionless plasma parameters are

$\displaystyle\mathcal{P}_{\textrm{2-comp}}=$
$\displaystyle\left(\beta_{\parallel,c},\frac{w_{\parallel,c}}{c},\frac{T_{\perp,c}}{T_{\parallel,c}},\frac{T_{\perp,b}}{T_{\parallel,b}},\right.$
$\displaystyle\left.\frac{n_{b}}{n_{c}},\frac{T_{\parallel,b}}{T_{\parallel,c}},\frac{\Delta
v_{b,c}}{v_{Ac}}\right),$

where we define the thermal-to-magnetic pressure ratio
$\beta_{\parallel,j}=8\pi n_{j}T_{\parallel,j}/B^{2}$, the core-proton Alfvén
velocity as $v_{Ac}=B/\sqrt{4\pi m_{p}n_{c}}$, and the speed of light $c$.
Frequencies are normalized to the proton gyrofrequency
$\Omega_{p}=q_{p}B/m_{p}c$. For this study, we neglect the contribution of
alphas and other minor ions and treat the electrons as a single isotropic
distribution with density and velocity necessary to enforce quasi-neutrality
and zero net current. The impact of the non-proton components on stability
will be the focus of future study.

Figure 3: SPAN-i observation of the proton velocity distribution for the
interval under analysis in Fig. 4 as a function of $v_{z}$ and $v_{r}$ (top)
and $v_{y}$ and $v_{r}$ (bottom), in SPAN-i instrument coördinates where
$v_{r}=\sqrt{v_{x}^{2}+v_{y}^{2}}$. Diamonds represent the central values of
the instrument’s velocity space bins, color the proton distribution phase-
space density, and the arrow magnetic field orientation with the length
representing the Alfvén speeed.  Figure 4: Comparison of linear stability and
resonances for the one- and two-component models, left and right columns,
associated with the SPANi observation shown in Fig. 3. Top row: Fastest
growing mode calculated by the Nyquist method as a function of
$k_{\perp}d_{p}$ and $k_{\parallel}d_{p}$. Second row: linear dispersion
relation $\omega_{\textrm{r}}(k_{\parallel}d_{p})/\Omega_{p}$ for the four
weakly damped, parallel propagating linear modes. Third: Normalized growth
(solid) or damping (dashed) rates $\gamma/\omega_{\textrm{r}}$ for the same
modes. Fourth: Cyclotron resonant velocities normalized to $v_{A}$. Bottom:
Illustration of phase (dashed vertical) and resonant (solid) velocities for
the four modes at the wavevector associated with the maximum growth rate (dot-
dashed line in middle panels), the associated curves of constant energy in the
wave-frame, and the phase-space densities associated with the one- and two-
component models (grey-scale).

Given an example SPAN-i measurement of $f_{p}(\mathbf{v})$, shown in Fig. 3,
both the one-component and two-component models are constructed, producing the
sets of dimensionless parameters $\mathcal{P}_{\textrm{1-comp}}$ and
$\mathcal{P}_{\textrm{2-comp}}$. For the selected example, starting at
11:16:22 on 01/30/2020, these sets are:

$\displaystyle\mathcal{P}_{\textrm{1-comp}}=$
$\displaystyle\left(\beta_{\parallel,p}=1.0646,\frac{w_{\parallel,p}}{c}=2.786\times
10^{-4}\right.,$ (7)
$\displaystyle\left.\frac{T_{\perp,p}}{T_{\parallel,p}}=0.389\right)$

and

$\displaystyle\mathcal{P}_{\textrm{2-comp}}=$
$\displaystyle\left(\beta_{\parallel,c}=0.410,\frac{w_{\parallel,c}}{c}=1.861\times
10^{-4},\right.$
$\displaystyle\left.\frac{T_{\perp,c}}{T_{\parallel,c}}=0.770,\frac{T_{\perp,b}}{T_{\parallel,b}}=0.620,\frac{n_{b}}{n_{c}}=0.157,\right.$
$\displaystyle\left.\frac{T_{\parallel,b}}{T_{\parallel,c}}=2.465,\frac{\Delta
v_{b,c}}{v_{Ac}}=-1.350\right).$

Given these sets, we calculated
$\gamma^{\textrm{max}}(\mathbf{k}d_{p})/\Omega_{p}$ using the Nyquist method,
shown in the top two panels in Fig. 4, which in turn allows the calculation of
$\gamma^{\textrm{max}}/\Omega_{p}$ over the entire wavevector range, as well
as the associated $\omega_{\textrm{r}}^{\textrm{max}}/\Omega_{p}$,
$k^{\textrm{max}}d_{p}$, $\theta_{kB}^{\textrm{max}}$, and other
eigenfunctions of the unstable modes. For this measurement and associated
models, $\gamma^{\textrm{max}}/\Omega_{p}$ is significantly larger for the
two-component model and the wavevector region supporting unstable modes is
broader compared to the one-component model, though both models predict the
same mode, the parallel propagating firehose/fast-magnetosonic wave, to be
linearly unstable.

For validation, we compare these predicted properties to the normal mode
solutions for the forward and backward parallel propagating Alfvén and fast-
magnetosonic waves numerically calculated using the PLUME dispersion
solver(Klein & Howes, 2015). The central rows of Fig. 4 show the real
component of the normal mode frequency
$\omega_{\textrm{r}}(k_{\parallel}d_{p})/\Omega_{p}$ for fixed
$k_{\perp}d_{p}=10^{-3}$, the normalized growth or damping rates
$\gamma(k_{\parallel}d_{p})/|\omega_{\textrm{r}}|$, and the normalized $n=\pm
1$ cyclotron resonant velocities,

$\frac{v_{\textrm{res}}(k_{\parallel})}{v_{A}}=\frac{\omega_{\textrm{r}}(k_{\parallel})-n\Omega_{p}}{k_{\parallel}v_{A}}$
(9)

where the choice of sign of $n$ is determined by the wave’s polarization and
direction of propagation; $n=+1$ for the forward Alfvén and backwards fast
modes and $n=-1$ for the backwards Alfvén and forward fast modes. For these
nearly parallel modes, there is no significant $n=0$ contribution to the wave-
particle interaction. We find good agreement with the kinds of modes and
region of wavevectors predicted to be stable and unstable from both the
Nyquist and traditional dispersion calculation.

Both models are unstable to the parallel firehose instability for this
interval, but there are significant differences—illustrated in the bottom
panels of Fig. 4— in the resonant coupling between the protons and the
electric field. The wave-phase velocity for each of the four parallel
propagating modes at a fixed wavevector $k_{\parallel}d_{p}$, set to be
$|k^{\textrm{max}}|d_{p}$ for the one- or two-component model, is illustrated
as a dashed vertical line compared to the model phase-space density
$f_{p}(v_{\perp}/v_{A},v_{\parallel}/v_{A})$. The $n=\pm 1$ cyclotron resonant
velocity is shown as a solid vertical line, and contours of constant energy in
the wave-frame are illustrated as colored half-circles. The sign of the pitch
angle gradient of $f_{p}$ where the resonant velocity meets the contours of
constant energy determines if energy is transferred from the wave to the
protons, leading to damping of the wave, or from the protons to the wave,
leading to excitation and instability. For this interval, the fitting of a
secondary proton population leads to the suppression of the unstable anti-
beam-aligned fast mode and the enhancement of the beam-aligned fast mode’s
growth rate. The beam component also significantly increases the damping rate
of the anti-beam aligned Alfvén mode, leading it to switch propagation
directions at $k_{\parallel}d_{p}\approx 0.3$.

##  3 Inferred Stability Across Selections

The Nyquist instability analysis described in §2.3 is performed over the
entirety of SA, Fig. 5, and SB, Fig. 6, for both the one- and two-component
models (red and blue).

Figure 5: Dimensionless parameters from one- and two-component (red and blue)
models for SA, (a-f) and calculated instability characteristics, (g-j). (a):
thermal-to-magnetic pressure ratio $\beta_{\parallel,p}$ or
$\beta_{\parallel,c}$, (b): thermal-speed ratio $w_{\parallel,p}/c$ or
$w_{\parallel,c}/c$, (c): temperature anisotropy $T_{\perp,p}/T_{\parallel,p}$
or $T_{\perp,c}/T_{\parallel,c}$ ($T_{\perp,b}/T_{\parallel,b}$ in teal), (d):
temperature disequilibrium $T_{\parallel,b}/T_{\parallel,c}$, (e): density
ratio $n_{b}/n_{c}$, (f): relative drift velocity $\Delta v_{bc}/v_{A}$. (g):
maximum growth rate $\gamma^{\textrm{max}}/\Omega_{p}$ (h): normal mode real
frequency $\omega_{\textrm{r}}^{\textrm{max}}/\Omega_{p}$ (i & j): Amplitude
and angle, $|k|^{\textrm{max}}d_{p}$ and $\theta_{kB}^{\textrm{max}}$ of the
wavevector associated with fastest growing mode.  Figure 6: Dimensionless
parameters and calculated instability characteristics from one- and two-
component (red and blue) models for SB organized in the same format as Fig. 5.

For both selections, we see different predicted unstable behavior for the two
models. Using the one-component model for SA, only $40.5\%$ of the intervals
are found to be unstable, and of those most have relatively weak growth rates,
with a median value of
$\bar{\gamma}^{\textrm{max}}_{1-\textrm{comp}}=2.33^{3.41}_{1.57}\times
10^{-4}\Omega_{p}$. The sub- and super-scripts represent the
$25^{\textrm{th}}$ and $75^{\textrm{th}}$ percentiles of the unstable mode
growth rate distribution. These are parallel firehose instabilities, where
sufficiently extreme parallel-to-perpendicular thermal pressure ratios,
manifest in a one-component proton distribution, change the sign of the
velocity gradient at the cyclotron resonant velocity such that energy is
extracted from the protons to drive an unstable fast-magnetosonic mode. Due to
the symmetry of the one-component model, both forward and backward propagating
modes are driven. No other kinds of unstable modes are supported by the one-
component model during SA.

For the two-component model, $99.9\%$ of the intervals in SA are found to be
unstable, with a median growth rate of
$\bar{\gamma}^{\textrm{max}}_{2-\textrm{comp}}=2.54^{3.24}_{1.92}\times
10^{-2}\Omega_{p}$, two orders of magnitude larger than for the one component
model. All of the unstable intervals are associated with parallel propagating
fast-magnetosonic modes with $|k|^{\textrm{max}}d_{p}\approx 0.5$. Unlike the
symmetrically emitted unstable waves from the one-component model, the
unstable modes from the two-component model only propagate in the same
direction as the secondary proton population.111We define the radial component
of our coördinate system to align with the mean magnetic field. In both SA and
SB, PSP was in a region of Sunward magnetic polarity, meaning that the anti-
Sunward propagating secondary proton populations have a negative velocity with
respect to the primary proton population. The maximum growth rate of the
unstable fast mode is enhanced due to an increased phase-space density
associated with the secondary proton population, while the anti-beam aligned
fast-mode resonance is effectively starved of protons with which to interact,
leading to damping rather than instability for this mode.

We find differences in the kinds of instabilities predicted for the two models
in SB. Ninety-nine percent of the intervals are predicted to be linearly
unstable to the parallel propagating firehose instability for the one-
component model, with a median growth rate of
$\bar{\gamma}^{\textrm{max}}_{1-\textrm{comp}}=9.26^{13.3}_{6.12}\times
10^{-4}\Omega_{p}$. This is not the case for the two-component model. The
median growth rate for the two-component model is similar,
$\bar{\gamma}^{\textrm{max}}_{2-\textrm{comp}}=6.43^{15.6}_{3.30}\times
10^{-4}\Omega_{p}$, however only $55.7\%$ of the intervals are found to be
unstable and the associated fastest growing mode oscillates between a beam-
aligned, parallel propagating firehose mode and an oblique instability. This
demonstrates that fitting a secondary component does not universally enhance
the predicted growth rate and that more sophisticated treatments of velocity-
space structure can lead to the generation of different kinds of unstable
modes.

As seen in Fig. 7, $\gamma^{\textrm{max}}/\Omega_{p}$ is generally larger for
the two-component model than for the one-component model for SA. This is not
the case for SB, where more of the one-component intervals are unstable, while
the variance in the growth rate for the two-component model is larger. When
re-normalized to the normal mode frequency
$\omega_{\textrm{r}}^{\textrm{max}}$, Fig 7b, we see an enhancement in the
growth rates for the two-component model in SB, while the the other growth
rates remain relatively unaffected.

Other time scales of potential interest include an estimate for the non-linear
cascade rate at the wavevector of fastest growth,

$\displaystyle\gamma^{\textrm{max}}\tau_{nl}=$
$\displaystyle\left(\frac{\gamma^{\textrm{max}}}{v_{A}}\right)\left(k_{\textrm{break}}\right)^{-1/3}(|\mathbf{k}^{\textrm{max}}|)^{-2/3}$
(10) $\displaystyle=$
$\displaystyle\left(\frac{\gamma^{\textrm{max}}}{\Omega_{p}}\right)(\frac{2\pi
f_{\textrm{break}}}{\Omega_{p}}\frac{v_{A}}{v_{sw}})^{-1/3}(|\mathbf{k}^{\textrm{max}}d_{p}|)^{-2/3}$

where we approximate the transition from the injection to the inertial ranges
of turbulence as $k_{\textrm{break}}=2\pi f_{\textrm{break}}/v_{sw}$ with
$f_{\textrm{break}}$ found to be approximately $10^{-3}$ Hz when constructing
trace power-spectral density curves for either SA or SB, not shown. These
values are in rough agreement with the results reported in Chen et al. (2020).
The cascade time is estimated as the critically balanced nonlinear cascade
rate, $\tau_{nl}\sim\omega_{\textrm{Alfv\'{e}n}}^{-1}$(Goldreich & Sridhar,
1995; Mallet et al., 2015). Previous analysis between 0.3 and 0.7 au (Klein et
al., 2019) found that $\gamma^{\textrm{max}}$ never exceeded the estimated
nonlinear cascade rate, though the two rates were found to be within an order
of magnitude, with $50\%$ of the intervals having
$\gamma^{\textrm{max}}\tau_{nl}\gtrsim 0.2$. For the two-component model in
SA, the maximum growth rate is of the same order as $\tau_{nl}^{-1}$, with a
median value of
$\bar{\gamma}^{\textrm{max}}_{2-\textrm{comp}}\tau_{nl}=0.618^{0.813}_{0.463}$,
indicating that these predicted instabilities operate on similar timescales as
the nonlinear transport of energy through these spatial scales. Importantly,
while $\bar{\gamma}^{\textrm{max}}_{2-\textrm{comp}}\sim\tau_{nl}^{-1}$, the
median value of $\bar{\gamma}^{\textrm{max}}_{1-\textrm{comp}}\tau_{nl}$ is
$4.70^{6.90}_{3.14}\times 10^{-3}$ for the same interval. This emphasizes that
our choice of different models for the proton phase-space density will lead to
drastically different interpretations of the importance of different physical
processes. The impact of these instabilities, especially when the ions are
modeled as multiple components, on the turbulent transport of energy must be
considered in future modeling efforts. The median values of
$\gamma^{\textrm{max}}\tau_{nl}$ are comparable for SB, with
$\bar{\gamma}^{\textrm{max}}_{1-\textrm{comp}}\tau_{nl}=2.08^{3.04}_{1.45}\times
10^{-2}$ and
$\bar{\gamma}^{\textrm{max}}_{2-\textrm{comp}}\tau_{nl}=3.55^{6.81}_{1.74}\times
10^{-2}$, again showing that the two-component model does not universally
enhance growth rates compared to the one-component model. To remove variations
associated with the normalization by $\Omega_{p}$ due to changes in
$|\mathbf{B}|$ as a function of time, we also plot the growth rate in Hertz,
Fig 7d, and see a distribution of growth rates similar to that seen in panel
a.

Figure 7: Comparison of normalized growth rates for the two models for SA
(teal) and SB (gold), with the abscissa and ordinate mapping the one- and two-
component rates. In panels a,b,c, and d, $\gamma^{\textrm{max}}$ is normalized
to $\Omega_{p}$, $\omega_{\textrm{r}}^{\textrm{max}}$, $\tau_{nl}^{-1}$ and
$1$ Hz respectively. Black dots and bars correspond to medians and
$25^{\textrm{th}}$ and $75^{\textrm{th}}$ percentiles associated with the
unstable intervals.

By design, the one- and two-component models have the same parallel and
perpendicular thermal pressures for a given interval, which can be
characterized by the firehose (Kunz et al., 2015)

$\Lambda_{F}=\frac{\beta_{\parallel}-\beta_{\perp}}{2}+\frac{\sum_{j}n_{j}m_{j}|\Delta\tilde{v}_{j}|^{2}}{\sum_{j}(n_{j}m_{j})v_{A}^{2}}$
(11)

or mirror (Hellinger, 2007)

$\Lambda_{M}=\sum_{j}\beta_{\perp,j}\left(\frac{\beta_{\perp,j}}{\beta_{\parallel,j}}-1\right)-\frac{\left(\sum_{j}q_{j}n_{j}\frac{\beta_{\perp,j}}{\beta_{\parallel,j}}\right)^{2}}{2\sum_{j}\frac{(q_{j}n_{j})^{2}}{\beta_{\parallel,j}}}$
(12)

criterion, where $\Delta\tilde{v}_{j}$ is the difference between the bulk
speed of component $j$ and the center of mass velocity. When these criterion
exceed unity, large-scale firehose or mirror instabilities are generated. For
both SA and SB, the amplitude of neither criteria exceeds $\sim 0.5$ for
either model; therefore, it is the resonances between the proton distribution
and the associated electromagnetic fields and not the excess macroscopic
parallel or perpendicular pressures that drives the predicted unstable wave
modes.

Slight changes in the relative drift speed between the two proton populations
and their densities can have a significant impact on the kind of unstable mode
predicted to be generated. This is illustrated in Fig. 8, where nine
sequential illustrations of contours of constant
$\gamma^{\textrm{max}}(\mathbf{k})$ are shown for the one- and two-component
models for SPAN-i observations from near the beginning of SB. Throughout these
two minutes both the maximum growth rates and regions of unstable wavevectors
are largely unchanged for the one-component model. This is expected given that
$T_{\perp,p}/T_{\parallel,p}$ and $\beta_{\parallel,p}$ are relatively
constant over this time, remaining consistent with a parallel propagating
firehose instability. For the two-component model, oblique modes are initially
driven. A minute into the sequence, the maximum growth rate transitions to a
parallel propagating wavevector, and then transitions back to an oblique
instability. These transitions correspond to a temporary dip in the relative
density of the beam component and an increase in the relative drift speed.
Given that many kinds of waves are observed in this section of data, it
appears plausible that these transitions between parallel and oblique
instabilities may be real, but are not properly accounted for in overly
simplistic models of the protons as a single anisotropic distribution, which
only drive one kind of unstable mode.

Figure 8: Left: Contours of constant $\gamma^{\textrm{max}}/\Omega_{p}$ as a
function of $k_{\perp}\rho_{c}$ and $k_{\parallel}\rho_{c}$ for the one- and
two-component models (red and blue) for nine intervals at the start of SB.
Right: Temporal variation of $T_{\perp,p}/T_{\parallel,p}$ and
$\beta_{\parallel,p}$ from the one-component model (top) and of the relative
drifts and densities of the two-component model (bottom).

We note that there is not a simple parametric function dependence only on
$n_{b}/n_{c}$ and $\Delta v_{b,c}/v_{A}$ that divides the parallel unstable
modes from the oblique modes. In Fig. 9, we plot the angle of the fastest
growing mode $\theta_{kB}^{\textrm{max}}$ for the two-component model for SB
as a function of these two parameters. Generally, the larger the relative
drift, the more likely the model is predicted to generate an oblique unstable
mode, with the transition between parallel and oblique modes arising at lower
drifts for larger relative beam densities. However, we find many stable
intervals with very similar drifts and densities to the intervals unstable to
the generation of both parallel and oblique unstable modes. This can be
understood by recalling that the variation of the temperatures and
anisotropies of the individual proton components will have a significant
impact on the predicted stability of the system that is not captured in this
reduced parameter space. Due to this complexity, we do not attempt to offer a
simple parametric prescription for this transition between parallel and
oblique instabilities in this work, but do note again that if this
distribution is treated as a single proton population, the only instability
supported is the parallel-propagating, fast/magnetosonic firehose instability.

Figure 9: Wavevector angle $\theta_{kB}^{\textrm{max}}$ of fastest growing
mode during SB calculated using the two-component model, indicated by color,
as a function of relative densities $n_{b}/n_{c}$ and drift velocities $\Delta
v_{b,c}/v_{A}$. Grey squares indicate intervals predicted to be linearly
stable.

##  4 Conclusions

In this work, we have selected two hours-long intervals where in situ
measurements of the local plasma conditions have been made during PSP’s fourth
perihelion orbit. These measurements coincide with significant ion-scale wave
activity as observed by the FIELDS magnetometers. The proton phase-space
densities have been modeled as either a single anisotropic population, or two
relatively drifting anisotropic populations. The linear stability of both
models was calculated, with strikingly different predictions for the supported
linear modes. In the first selection, both models produce the same kind of
unstable mode, but the two-component model drives instabilities that grow
nearly two orders of magnitude faster, fast enough to potentially act on the
same timescales as the local nonlinear turbulent transfer of energy.
Additionally, the two-component model for SA only drives instabilities
propagating in a single direction, as opposed to the one-component model where
waves are driven both Sunward and anti-Sunward due to the enforced symmetry of
the simplified description of the protons. For the second selection, modeling
the protons using two components does not make the plasma more unstable, but
does change the kind of unstable modes driven, leading to an oscillation
between the production of parallel and oblique propagating waves.

As future lines of inquiry, we intend on extending this work to investigate
the predicted growth rates and waves concurrently observed with other plasma
parameters and solar wind conditions, such as intervals where the total
parallel proton pressure is exceeded by the total perpendicular pressure. We
will also include additional sources of free energy associated with minor ions
and electrons, to determine if they act to enhance or stabilize these growing
modes. This work will help to ascertain under what conditions which models may
suffice to properly describe kinetic processes. Importantly, as the
instabilities under consideration are resonant, we must also consider the
impact of departures from bi-Maxwellian distributions, either using other
analytic prescriptions, e.g. kappa (Livadiotis, 2015) or flattop
distributions(Klein & Chandran, 2016; Wilson et al., 2020), or via a direct
numerical integration of the observed phase-space density (Verscharen et al.,
2018).

The SWEAP Investigation and this publication are supported by the PSP mission
under NASA contract NNN06AA01C. K.G.K. is supported by NASA ECIP Grant
80NSSC19K0912. An allocation of computer time from the UA Research Computing
High Performance Computing at the University of Arizona is gratefully
acknowledged.

## References

  * Alterman et al. (2018) Alterman, B. L., Kasper, J. C., Stevens, M. L., & Koval, A. 2018, Astrophys. J., 864, 112
  * Bale et al. (2016) Bale, S. D., Goetz, K., Harvey, P. R., et al. 2016, Space Sci. Rev., doi:10.1007/s11214-016-0244-5
  * Bowen et al. (2020) Bowen, T. A., Mallet, A., Huang, J., et al. 2020, Astrophys. J. Supp., 246, 66
  * Chen et al. (2016) Chen, C. H. K., Matteini, L., Schekochihin, A. A., et al. 2016, Astrophys. J. Lett., 825, L26
  * Chen et al. (2020) Chen, C. H. K., Bale, S. D., Bonnell, J. W., et al. 2020, Astrophys. J. Supp., 246, 53
  * Fox et al. (2015) Fox, N. J., Velli, M. C., Bale, S. D., et al. 2015, Space Sci. Rev., doi:10.1007/s11214-015-0211-6
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, Astrophys. J., 438, 763
  * Hellinger (2007) Hellinger, P. 2007, Phys. Plasmas, 14, 082105
  * Kasper et al. (2015) Kasper, J. C., Abiad, R., Austin, G., et al. 2015, Space Sci. Rev., 1
  * Klein et al. (2018) Klein, K. G., Alterman, B. L., Stevens, M. L., Vech, D., & Kasper, J. C. 2018, Phys. Rev. Lett., 120, 205102
  * Klein & Chandran (2016) Klein, K. G., & Chandran, B. D. G. 2016, Astrophys. J., 820, 47
  * Klein & Howes (2015) Klein, K. G., & Howes, G. G. 2015, Phys. Plasmas, 22, 032903
  * Klein et al. (2017) Klein, K. G., Kasper, J. C., Korreck, K. E., & Stevens, M. L. 2017, Journal of Geophysical Research (Space Physics), 122, 9815
  * Klein et al. (2019) Klein, K. G., Martinović, M., Stansby, D., & Horbury, T. S. 2019, Astrophys. J., 887, 234
  * Kunz et al. (2015) Kunz, M. W., Schekochihin, A. A., Chen, C. H. K., Abel, I. G., & Cowley, S. C. 2015, Journal of Plasma Physics, 81, 325810501
  * Livadiotis (2015) Livadiotis, G. 2015, Journal of Geophysical Research (Space Physics), 120, 1607
  * Mallet et al. (2015) Mallet, A., Schekochihin, A. A., & Chandran, B. D. G. 2015, Mon. Not. Roy. Astron. Soc., 449, L77
  * Matteini et al. (2012) Matteini, L., Hellinger, P., Landi, S., Trávníček, P. M., & Velli, M. 2012, Space Sci. Rev., 172, 373
  * Nyquist (1932) Nyquist, H. 1932, Bell system technical journal, 11, 126
  * Verniero et al. (2020) Verniero, J. L., Larson, D. E., Livi, R., et al. 2020, Astrophys. J. Supp., 248, 5
  * Verscharen et al. (2018) Verscharen, D., Klein, K. G., Chandran, B. D. G., et al. 2018, Journal of Plasma Physics, 84, 905840403
  * Verscharen et al. (2019) Verscharen, D., Klein, K. G., & Maruca, B. A. 2019, Living Rev. Solar Phys., 16, 5
  * Wilson et al. (2020) Wilson, Lynn B., I., Chen, L.-J., Wang, S., et al. 2020, Astrophys. J., 893, 22
  * Yoon (2017) Yoon, P. H. 2017, Reviews of Modern Plasma Physics, 1, 4