# Signature of topology via heat transfer analysis in the Su–Schrieffer–Heeger
(SSH) model

Vipul Upadhyay  Department of Physics, Indian Institute of Technology Delhi,
Hauz Khas 110 016, New Delhi, INDIA.  M. Tahir Naseem  Faculty of Engineering
Science, Ghulam Ishaq Khan Institute of Engineering Sciences and Technology,
Topi 23640, Khyber Pakhtunkhwa, Pakistan,  Özgür E. Müstecaplıoğlu
<EMAIL_ADDRESS>Department of Physics, Koç University, 34450 Sariyer,
Istanbul, Türkiye,  TÜBİTAK Research Instute for Fundamental Sciences, 41470
Gebze, Türkiye.  Rahul Marathe<EMAIL_ADDRESS>Department of
Physics, Indian Institute of Technology Delhi, Hauz Khas 110 016, New Delhi,
INDIA.

###### Abstract

In this work, we explore the potential of thermodynamics as a tool for
identifying the topological phase transition. Specifically, we focus on a one-
dimensional Su–Schrieffer–Heeger (SSH) chain sandwiched between two fermionic
baths. To investigate distinctive thermodynamic signatures associated with the
topological phase, we employ heat flow analysis. Our results, derived using a
global master equation, unveil a significant suppression of heat flow as we
transition from the trivial to the topological phase. This decline in heat
flow can be attributed to the reduction in transmission coefficients of non-
zero energy modes within the topological phase. It may serve as an indicator
of a phase transition. Furthermore, we investigate the heat flow asymmetry to
search for phase transition indicators. Interestingly, no asymmetry is
observed when employing fermionic baths. However, upon substituting fermionic
baths with bosonic ones, we report a non-zero heat flow asymmetry. For the SSH
model with a few fermionic sites, this asymmetry is more pronounced in the
topological phase compared to the trivial phase. Therefore, the observed
behavior of the heat diode provides an additional means of distinguishing
between the topological and trivial phases. Finally, we delve into the
contributions from both bulk and edge effects in heat flow and rectification
to explore the impact of small system sizes on our findings.

##  I Introduction

In classical physics, phase transitions are identified through the examination
of thermodynamic signatures, manifesting as abrupt changes in properties like
specific heat or magnetization. However, quantum phase transitions, distinct
from classical ones, occur solely due to quantum mechanical properties,
particularly at absolute zero temperatures. These transitions give rise to
exotic quantum phases, emphasizing the importance of identifying and
understanding them. This unique realm of condensed matter physics has become a
focal point of research Azzouz (2006); Zhang and Song (2015); Miao  _et al._
(2017); Roy  _et al._ (2019); Liu  _et al._ (2020); Luo (2022), demanding
innovative approaches to detect and comprehend these transitions, as they hold
the key to unlocking novel quantum materials and technological advancements.

Transport studies in topological materials have been a subject of intense
research Molignini  _et al._ (2017); Benenti  _et al._ (2009); Bandyopadhyay
_et al._ (2021); Mondal  _et al._ (2022); van Caspel  _et al._ (2019). This is
because the separation of edge and bulk in the topological phase leads to
exotic transport properties for a system in topological phase Rivas and
Martin-Delgado (2017); Molignini  _et al._ (2017). In mathematics, topology is
a geometrical property for any system Kane and Moore (2011) which is robust
against continuous deformation. Depending upon the context of the specific
system, a topological invariant may be mathematically defined. In quantum
mechanics, the topology appears due to the independence in defining the phase
of the wavefunction. As such we can define the topology of energy band
structures in a quantum system. It is a property that is used to distinguish
between two phases of a Hamiltonian which can’t be smoothly transformed into
each other without crossing the energy gap between ground and bulk. These are
also called non-adiabatically connected regions Kane and Moore (2011); Asbóth
_et al._ (2016). For isolated systems at absolute zero, winding number or
Berry phase Asbóth  _et al._ (2016) plays the role of topological invariant
helping us to distinguish between the phases of the system. Coupling the
system to heat baths will usually not preserve the topological invariants and
new topological invariants like the Ulham phase are needed Viyuela  _et al._
(2014); Kempkes  _et al._ (2016). In 1-D quantum systems topological phase
transition results in localization of wavefunction at edges. This localization
can then be utilized for altering transport properties and indicators of the
phase of the system. A lot of studies have been dedicated to establishing heat
transport as an indicator of the topological phase. Experiments for such
systems are possible in cold atoms, and nanowires Meier  _et al._ (2016a);
Drost  _et al._ (2017); Kanungo  _et al._ (2022); Meier  _et al._ (2016b);
Gröning  _et al._ (2018).

In this study, we set out to address two main inquiries: firstly, the
viability of utilizing heat transport for identifying the topological phase
within the SSH model; and secondly, the prospect of harnessing the asymmetry
in heat transport across different phases to engineer a specialized device,
notably a heat diode Tanaka  _et al._ (2022); Fluckey  _et al._ (2022);
Karabassov  _et al._ (2023); Li  _et al._ (2017), whose operational efficiency
hinges upon the specific topological phase encountered. Our investigation
yields affirmative responses to both questions, subject to defined parameter
regimes. Additionally, we analyze thermodynamic intricacies associated with
topological phase transitions. This analysis entails the identification of an
effective temperature corresponding to a specific energy mode of the system,
following the analytical framework proposed by previous studies Kempkes  _et
al._ (2016). Notably, our analysis differs from prior works by adopting a
global master equation, thus ensuring thermodynamic consistency and obviating
the reliance on the local master equation He and Chien (2023); Nava  _et al._
(2023) which may violate the second law of thermodynamics Levy and Kosloff
(2014); Naseem  _et al._ (2018). Furthermore, our method does not require
external chemical potentials at each site, distinguishing it from previous
studies. For the rectification analysis Landi  _et al._ (2014); Jiang  _et
al._ (2010); Ordonez-Miranda  _et al._ (2017); Karg ı  _et al._ (2019a); Segal
and Nitzan (2005), we replace fermionic baths with bosonic counterparts and
exploit inherent bath-system interaction asymmetry to induce rectification.
While this technique echoes prior studies in inducing rectification Segal and
Nitzan (2005), its application as a potent harbinger of phase transitions
within a one-dimensional SSH chain represents uncharted territory within our
study.

The structure of this paper is as follows: In Sec. II, we present the model
and introduce the global master equation. Sec. III is dedicated to the
analysis of heat flow results. In Sec. IV, we define the effective temperature
and explore the thermodynamic characteristics of our system. We replace the
fermionic baths with bosonic baths and investigate heat current rectification
in Sec. V. To conclude, in Sec. VI, we summarize our findings and offer
concluding remarks. In the appendices A and B we give detailed derivations of
some of the results used in the main text.

##  II The SSH Model

Figure 1: Schematic of the model system. Double and single bonds signify the
relative strengths of interaction between sites with double bonds having more
interaction strength than the single bond.

In our model, we consider a one-dimensional system consisting of fermions
occupying $N$ sites. The interaction strengths of these fermions vary
depending on whether the bonds between sites are odd or even. It is important
to note that $N$ is always chosen to be an even number. The Hamiltonian
governing this model is expressed as follows Su  _et al._ (1979); Sirker  _et
al._ (2014):

$\displaystyle\hat{H}_{S}$
$\displaystyle=-t\sum_{n}^{N-1}(1+(-1)^{n}\delta)\hat{c}_{n+1}^{\dagger}\hat{c_{n}}+H.C$
(1)

where, $v_{1}=-t(1-\delta)$, $v_{2}=-t(1+\delta)$ are the odd and even bond
strengths, respectively. This represents the simplest 1-D model permitting a
topological phase. It is well known Su  _et al._ (1979); Sirker  _et al._
(2014) that when $v_{1}<v_{2}$, the system enters the topological phase,
whereas the phase characterized by $v_{2}<v_{1}$ is the trivial phase.
Consequently, $\delta$ can be interpreted as a parameter indicative of
topological order, wherein $\delta<0$ signifies the trivial phase, while
$\delta>0$ corresponds to the topological phase. The primary distinction
between the trivial and the topological phase is the existence of two zero-
energy modes within the latter Asbóth  _et al._ (2016). This difference leads
to the localization of energy eigenvectors corresponding to these zero energy
modes at the two ends of the chain. This system is coupled at its two ends
with two distinct fermionic baths that do not overlap. Each of these baths is
governed by its respective Hamiltonian

$\displaystyle\hat{H}_{B}=\sum_{n}\omega^{L}_{n}\hat{b}_{n}^{\dagger
L}\hat{b}^{L}_{n}+\sum_{n}\omega^{R}_{n}\hat{b}^{\dagger
R}_{n}\hat{b}^{R}_{n},$ (2)

here $\hat{b}_{n}^{\dagger L}(\hat{b}^{L}_{n})$ are the creation (destruction)
operator for the $n$-th mode of the left bath. The system bath coupling is
given by

$\displaystyle\hat{H}_{SB}$
$\displaystyle=(\hat{c}_{1}^{\dagger}+\hat{c}_{1})\otimes\sqrt{\kappa_{L}}\sum_{k}(g^{L}_{k}\hat{b}^{L}_{k}+g^{*L}_{k}\hat{b}_{k}^{\dagger
L})$
$\displaystyle+(\hat{c}_{N}^{\dagger}+\hat{c}_{N})\otimes\sqrt{\kappa_{R}}\sum_{k}(g^{R}_{k}\hat{b}^{R}_{k}+g^{*R}_{k}\hat{b}_{k}^{\dagger
R}),$ (3)

where $\sqrt{\kappa_{L}}(\sqrt{\kappa_{R}})$ is the system-bath interaction
strength with the left (right) bath.

###  II.1 System eigenvalues and eigenvector

The eigenvalues of the SSH model are presented as Sirker  _et al._ (2014);
Ikramov (1999)

$\displaystyle\lambda_{k}=\pm\sqrt{2}t\sqrt{1+\delta^{2}+(1-\delta^{2})\cos{\theta_{k}}}$
(4)

where $\theta_{k}$ can be found by solving the equation

$\displaystyle(1-\delta)T(\theta_{k},N/2)+(1+\delta)T(\theta_{k},N/2-1)=0,$
(5)

where,

$\displaystyle
T(\theta_{k},n)=\frac{\sin{[(n+1)\theta_{k}]}}{\sin{\theta_{k}}}.$ (6)

The equation above has $N$ real solutions when $\delta<0$. However, with
$\delta>0$, its behavior changes: it includes ($N-2$) real solutions along
with 2 complex roots. These complex roots correspond to the zero energy modes
of the SSH model, which are linked to the complex solutions of $\theta_{k}$.
The elements of the eigenvector $\bm{\Phi}^{k}$ associated with the eigenvalue
$\lambda_{k}$ are given as Sirker  _et al._ (2014):

$\displaystyle\bm{\Phi}_{2n-1}^{k}$
$\displaystyle=\frac{1+\delta}{1-\delta}T(\theta_{k},n-2)+T(\theta_{k},n-1),$
$\displaystyle\bm{\Phi}_{2n}^{k}$
$\displaystyle=\pm\frac{\lambda_{k}}{(1-\delta)}T(\theta_{k},n-1).$ (7)

Considering that $t$ serves as a mere multiplication factor for our matrix
$\bm{\Phi}^{k}$, the eigenvectors should remain independent of the value of
$t$. Moreover, owing to the symmetries inherent in the Hamiltonian, we can
expect additional simplifications in our eigenvectors Nava  _et al._ (2023).
As the Hamiltonian in Eq. (1) exhibits symmetry under the transformation $n\to
N+1-n$, we find that $|\bm{\Phi}^{k}_{1}|^{2}=|\bm{\Phi}^{k}_{N}|^{2}$ for non
near-zero energy modes. For zero-energy modes, wavefunction localization
occurs at either end, rendering this symmetry inapplicable. Finally, as a
result of particle-hole symmetry, we have
$|\bm{\Phi}^{k}_{n}|^{2}=|\bm{\Phi}^{-k}_{n}|^{2}$.

###  II.2 Master equation and heat flow

Utilizing the eigenvectors mentioned above and constructing a unitary matrix
$U$, we can diagonalize the system’s Hamiltonian and express it in the energy
basis as (see Appendix A for details)

$\displaystyle\hat{H}_{S}=\sum_{k}\omega_{k}\hat{a}^{\dagger}_{k}\hat{a}_{k},$
(8)

where, $\omega_{k}$ represents the energy associated with the $k$-th mode,
while $\hat{a}^{\dagger}_{k}$ and $\hat{a}_{k}$ denote the corresponding
creation and destruction operators, respectively. The derivation of the master
equation for our model adheres to the conventional Born-Markov and secular
approximations. In the Schrödinger picture, it is given by Breuer and
Petruccione (2002)

$\frac{d}{dt}\hat{\rho}(t)=-i[H_{S},\rho]+\mathcal{L}_{L}\hat{\rho}(t)+\mathcal{L}_{R}\hat{\rho}(t),$
(9)

where, $\hat{\rho}$ is the density matrix of the N site SSH model, with (see
Appendix A),

$\displaystyle\mathcal{L}_{L}\hat{\rho}(t)$
$\displaystyle=\kappa^{L}_{j}\sum_{j}G(\omega_{j},T_{L})\mathcal{D}(\hat{a}_{j})[\hat{\rho}]+G(-\omega_{j},T_{L})\mathcal{D}(\hat{a}_{j}^{\dagger})[\hat{\rho}],$
$\displaystyle\mathcal{L}_{R}\hat{\rho}(t)$
$\displaystyle=\kappa^{R}_{j}\sum_{j}G(\omega_{j},T_{R})\mathcal{D}(\hat{a}_{j})[\hat{\rho}]+G(-\omega_{j},T_{R})\mathcal{D}(\hat{a}_{j}^{\dagger})[\hat{\rho}],$
(10)

Liouvillian super-operators of the left and right baths, respectively. The
index $j$ in the sum runs over all the eigenvalues,
$\kappa^{L}_{j}=\kappa_{L}|U_{1,j}|^{2}$, and
$\kappa^{R}_{j}=\kappa_{R}|U_{N,j}|^{2}$. Considering that $|U_{i,j}|^{2}$
represents the square of the $i$-th component of the $j$-th eigenvector, we
can express $\kappa^{L}_{j}$ $(\kappa^{R}_{j})$ in terms of the eigenvectors
defined in Eq. (II.1). Specifically,
$\kappa^{L}_{j}=\kappa_{L}|\bm{\Phi}^{j}_{1(N)}|^{2}$ and
$\kappa^{R}_{j}=\kappa_{R}|\bm{\Phi}^{j}_{1(N)}|^{2}$. Here, we have leveraged
the symmetry property as outlined in Sec. II.1. The Lindblad dissipators given
in Eq. (II.2) are defined as

$\displaystyle\mathcal{D}(\hat{A})[\hat{\rho}]=2\hat{A}\hat{\rho}\hat{A}^{\dagger}-\hat{A}^{\dagger}\hat{A}\hat{\rho}-\hat{\rho}\hat{A}^{\dagger}\hat{A},$
(11)

and bath spectral density is given by

$\displaystyle
G(\omega,T_{i})=\begin{cases}J(\omega)(1-F_{i}(\omega))&\omega>0,\\\
J(|\omega|)F_{i}(|\omega|)&\omega<0.\end{cases}$ (12)

In this work, unless stated otherwise, we assume flat spectra for the baths,
implying that $J(\omega)=J$ holds true for all $\omega$ values. The Fermi-
Dirac distribution is given as

$\displaystyle F_{i}(\omega)=\frac{1}{e^{\omega/T_{i}}+1}.$ (13)

Finally, the heat current from the left and right baths can be calculated
using Rivas and Martin-Delgado (2017)

$\displaystyle I_{L(R)}=\text{Tr}(\mathcal{L}_{L(R)}H_{S}).$ (14)

##  III Results

The expression for the heat current, as computed using Eq. (14) for the SSH
model, is given by (see Appendix B)

$\displaystyle
I_{L}=4\sum_{\omega_{j}>0}\mathcal{T}_{j}\omega_{j}(F_{L}(\omega_{j})-F_{R}(\omega_{j})).$
(15)

We observe that the expression for heat current presented above adheres to the
Landauer form. This representation unveils that the total heat current arising
from the summation of heat currents flowing through independent energy
channels with frequencies $\pm\omega_{j}$. Additionally, each energy channel
$\omega_{j}$ is characterized by its respective transmission coefficient,
denoted as $\mathcal{T}_{j}$ Datta (2005)

$\displaystyle\mathcal{T}_{j}=\frac{k^{L}_{j}k^{R}_{j}}{k^{L}_{j}+k^{R}_{j}}=|\bm{\Phi}^{j}_{1(N)}|^{2}\frac{\kappa_{L}\kappa_{R}}{\kappa_{L}+\kappa_{R}}.$
(16)

Figure 2: (a) Energy spectrum with $\delta$, (b) Variation of heat current
with $\delta$ for various $t$’s, (c) Transmission coefficient for all the bulk
energy modes (the near zero energy modes are not shown) with $\delta$ and (d)
First and last component of eigenvector of a typical non-zero energy mode. For
all the case, if not specified otherwise, the values of the system parameters
are $N=200,t=1,T_{L}=1,T_{R}=0.1,\kappa_{L}=0.1,\kappa_{R}=0.1$.

In Fig. 2(a), we present the energy spectrum of the SSH model as a function of
the order parameter, $\delta$. Notably, a substantial portion of the energy
eigenvalues lies above the $`2\delta^{\prime}$ threshold in both phases. The
primary spectral distinction between the trivial and topological phases is the
presence of two nearly zero-energy modes in the topological phase. Looking at
the variation of heat current with $\delta$ in Fig. 2(b), It becomes apparent
that the heat current in the topological phase experiences significant
suppression compared to the trivial phase. Notably, the current does not
approach zero as the bond strength significantly weakens in the trivial phase;
instead, it saturates at $2t(F_{L}(2t)-F_{R}(2t))$. This is attributed to the
breakdown of the secular approximation in the derivation of the global master
equation, particularly when dealing with a weak intrasystem coupling strength
Hofer  _et al._ (2017). Therefore, we limit our analysis to scenarios in which
$\kappa<<t(1\pm\delta)$ to ensure the continued validity of the secular
approximation. For convenience, we will generally work in the range
$-0.5<\delta<0.5$.

To gain deeper insights into the substantial disparity in heat transport
between the trivial and topological phases, we perform a thorough transmission
coefficient analysis below.

###  III.1 Fall in heat current via transmission coefficient analysis:

Examining the heat current expression in Eq. (15), we observe two distinct
factors influencing its variation at fixed bath temperatures. The first factor
pertains to alterations in energy levels, while the second factor relates to
modifications in transmission coefficients. Upon closer examination of Fig.
2(a), it becomes evident that all heat conduction channels on either side are
the same, with the exception of two near-zero energy channels within the
topological phase. However, it’s worth noting that these near-zero energy
channels make only negligible contributions to heat current in large systems.
Consequently, this difference in energy spectrum diminishes as we scale up the
system size. Therefore, changes in energy levels alone cannot account for the
substantial decrease in heat current observed in the topological region. Now
we analyze the transmission coefficients associated with energy levels greater
than $2\delta$, as defined in Eq. $\eqref{transmission_coefficients}$. As we
go deeper into the topological phase, the wavefunctions of the two zero-energy
modes gradually localize toward the two ends of the chain. Consequently, the
first and $N$-th components of non-zero eigenvectors begin approaching zero.
The further we advance into the topological region, the more pronounced this
effect becomes. Consequently, the transmission coefficients for all bulk
energy channels gradually tend toward zero, resulting in a significant
reduction in heat current within the topological region. This observation
aligns with previous findings reported in Nava  _et al._ (2023).

In Fig. 2(c), we present the bulk transmission coefficients for all available
channels. In the trivial phase, the average of all these transmission
coefficients remains almost the same, but as we transition into the
topological phase, it starts to decrease noticeably. In Fig. 2(d), we plot the
first and last components of a non-zero energy eigenvector. Intriguingly, both
of these components are equal and tend to zero once we are deep in the
topological phase. This tendency of near-zero energy modes to localize at the
ends is an inherent feature of the topological property characterizing the SSH
model. Therefore, we can confidently assert that the reduction in heat current
observed in the topological region is a direct consequence of the system’s
topological characteristics.

###  III.2 Heat current with system size

Figure 3: Variation of heat current with $N$ for (a) Topological Phase, (b)
Trivial phase. (c) Comparison of fitting quality of $F_{1}$ and $F_{2}$ for
variation of heat current with $N$ for topological and trivial (inset) phase
and (d) Variation of edge and bulk current with $\delta$. The edge current in
fit $F_{2}$ has been scaled by a factor of ‘$\frac{1}{50}$’ for visual
convenience. All the fits are done for data in the range $N=\\{6,50\\}$.

Examining the behavior of heat current concerning the system size N in Figs.
3(a) and 3(b). When $N$ exceeds approximately 20, heat transport appears to
become independent of system size. This observation implies that, for large
system sizes, heat transport follows a ballistic pattern and deviates from
Fourier’s law. In contrast, for considerably smaller system sizes, the current
exhibits a non-monotonic decline with $N$, attributed to finite size effects.
Notably, this drop in heat current with $N$ is more pronounced in the
topological phase compared to the trivial phase. To gain a quantitative
understanding of this phenomenon, we employ two models to fit the behavior of
heat current with system size. Both models are selected to ensure that heat
current saturates in the thermodynamic limit, maintaining a ballistic
transport. The first model, $F_{1}$, involves dividing the total current into
ballistic and diffusive components Landi  _et al._ (2022)

$\displaystyle F_{1}:=I_{L}(N)=J_{B}+J_{E}/N,$ (17)

where $J_{B}$ is the bulk current and $J_{E}$ is the edge current. The second
model $F_{2}$ we use is given by

$\displaystyle F_{2}:=I_{L}(N)=J_{B}+J_{E}e^{-N}.$ (18)

Examining Fig. 3(c), we observe that the fit is significantly better for the
exponentially decaying model. However, in qualitative terms, both models
reveal the same underlying property: as we enter the topological phases, the
contribution from the edges increases rapidly, reaches a maximum, and then
gradually declines.

A practical approach to distinguish between edge and bulk current can be to
isolate the current contributions originating from the two near-zero energy
modes and compare them to the currents associated with energy modes at
energies $\geq 2\delta$. By defining the contribution to current from energy
modes below $2\delta$ as $J_{Ep}$, we can draw insights from Fig. 3(d), where
we observe that this contribution qualitatively aligns with the fitting
analysis we have presented.

##  IV Effective temperature and steady state density matrix

The master equation (9) can be solved by defining an effective temperature
Rivas and Martin-Delgado (2017)

$\displaystyle
T^{j}_{eff}=\frac{|\omega_{j}|}{\log\Big{(}\frac{\kappa^{L}_{j}e^{\frac{|\omega_{j}|}{T_{L}}}\big{(}e^{\frac{|\omega_{j}|}{T_{R}}}+1\big{)}+\kappa^{R}_{j}e^{\frac{|\omega_{j}|}{T_{R}}}\big{(}e^{\frac{|\omega_{j}|}{T_{L}}}+1\big{)}}{\kappa^{L}_{j}\big{(}e^{\frac{|\omega_{j}|}{T_{R}}}+1\big{)}+\kappa^{R}_{j}\big{(}e^{\frac{|\omega_{j}|}{T_{L}}}+1\big{)}}\Big{)}}.$
(19)

By defining this effective temperature, Eq. (9) reduces to

$\displaystyle\frac{d\hat{\rho}}{dt}$
$\displaystyle=-i[H_{S},\rho]+\sum_{j}\frac{1}{\kappa^{L}_{j}+\kappa^{R}_{j}}\Big{(}\kappa^{L}_{j}G(\omega_{j},T^{j}_{eff})+\kappa^{R}_{j}G(\omega_{j},T^{j}_{eff})\Big{)}\mathcal{D}(\hat{a}_{j})[\hat{\rho}]$
$\displaystyle+\frac{1}{\kappa^{L}_{j}+\kappa^{R}_{j}}\Big{(}\kappa^{L}_{j}G(-\omega_{j},T^{j}_{eff})+\kappa^{R}_{j}G(-\omega_{j},T^{j}_{eff})\Big{)}\mathcal{D}(\hat{a}_{j}^{\dagger})[\hat{\rho}].$
(20)

We can see that if $T_{L}=T_{R}$, $T_{eff}^{j}=T_{L}(T_{R})$. Given that the
global master equation ensures the thermalization of the density matrix at
equilibrium, our system’s steady state is consequently established as

$\displaystyle\hat{\rho}_{SS}=\frac{1}{Z}e^{-\sum_{j}{\omega_{j}}/{T_{eff}^{j}}\hat{a}_{j}^{\dagger}\hat{a}_{j}},$
(21)

where the partition function is given by

$\displaystyle
Z=Tr[e^{-\sum_{j}{\omega_{j}}/{T_{eff}^{j}}\hat{a}_{j}^{\dagger}\hat{a}_{j}}].$
(22)

This indicates that we can conceptualize our system as an assembly of $N$
independent two-state systems, each in equilibrium with a bath at temperature
$T_{eff}^{j}$. The energy gap between the states in each of these systems is
determined by $\omega_{j}$.

Examining the behavior of $\kappa^{L}_{j}$ and $\kappa^{R}_{j}$ for the two
near-zero energy modes, it becomes apparent that, owing to the localization of
these modes at opposite ends, one tends toward zero while the other converges
to one,

$\displaystyle\lim\limits_{\kappa^{L}_{j}(\kappa^{R}_{j})\to
1,\kappa^{R}_{j}(\kappa^{L}_{j})\to 0}T_{eff}^{j}=T_{L}(T_{R}).$ (23)

When we apply this limit to the effective temperature expression in Eq. (19),
we observe that one of the near-zero energy modes attains an effective
temperature $T_{L}$, while the other reaches an effective temperature $T_{R}$.
Additionally, given that $\kappa^{L}_{j}=\kappa^{R}_{j}$ for bulk energy
modes, the effective temperature $T_{eff}^{j}$ for bulk modes remains
unaffected by the wave function and relies solely on the energy of the modes
as well as the temperatures of the left and right baths. Plotting the
effective temperature against the mode energy in Figs. 4(a) and 4(b) and 4(b),
we observe that the effective temperature remains the same for bulk modes in
both the trivial and topological phases. The primary distinction lies in the
two edge modes within the topological phase. Notably, we observe a gradual
shift in the temperatures of the near-zero energy modes toward $T_{L}$ and
$T_{R}$.

###  IV.1 Thermodynamics of bulk and edge contributions

Figure 4: Variation of effective temperature with energy $E$ for (a)
Topological Phase, (b) Trivial phase. Variation of edge and bulk contribution
with $\delta$ for (c) entropy $S$ and (d) specific heat $C_{T}$.

The effective partition function for our system can be written as

$\displaystyle Z=\Pi_{j=1}^{j=N}(1+e^{-\beta_{j}\omega_{j}}),$ (24)

where, $\beta_{j}=\frac{1}{T_{eff}^{j}}$. With this partition function, we can
define the system’s effective free energy

$\displaystyle
F=-\sum_{j}\frac{1}{\beta_{j}}\log(1+e^{-\beta_{j}\omega_{j}}).$ (25)

Similarly, we can define the entropy of our system as

$\displaystyle
S=-\frac{dF}{dT}=\sum_{j}\log(1+e^{-\beta_{j}\omega_{j}})+\frac{\omega_{j}\beta_{j}e^{-{\beta_{j}\omega_{j}}}}{1+e^{-\beta_{j}\omega_{j}}}.$
(26)

Utilizing the approach detailed in Kempkes  _et al._ (2016), we formulate the
total entropy as a summation of edge entropy and bulk entropy

$\displaystyle S=S_{E}+S_{B}*N,$ (27)

where $S_{E}$ and $S_{B}$ are edge and bulk entropies, respectively.
Similarly, we can analyze the specific heat,

$\displaystyle C_{T}=-T\frac{d^{2}F}{dT^{2}},$ (28)

and using Eq. (25), the specific heat evaluates to

$\displaystyle
C_{T}=\sum_{j}\frac{\beta_{j}^{2}\omega_{j}^{2}e^{-\beta_{j}\omega_{j}}}{1+e^{-\beta_{j}\omega_{j}}}-\frac{\beta_{j}^{2}\omega_{j}^{2}e^{-2\beta_{j}\omega_{j}}}{(1+e^{-\beta_{j}\omega_{j}})^{2}}.$
(29)

Again distributing total specific heat into edge and bulk part using Hill’s
thermodynamics Hill (1994), we have

$\displaystyle C_{T}=C_{E}+N*C_{B}$ (30)

Upon performing data fitting, we observe an excellent fit. As illustrated in
Figs. 4(c) and 4(d), when transitioning from the trivial phase to the
topological phase, there is a notable increase in the contribution from the
edge to both entropy and specific heat. This observation underscores the
heightened significance of the edge in the topological state, aligning with
findings presented in Kempkes  _et al._ (2016) for an equilibrium system.

##  V Heat diode with bosonic baths

It is reported in previous studies that in a spin-boson system, the asymmetry
in the system bath coupling for a spin system sandwiched between two bosonic
baths results in heat flow rectification Segal and Nitzan (2005); Karg ı  _et
al._ (2019b); Upadhyay  _et al._ (2021). In our case, we observe from the heat
current expression in Eq. (2) that, regardless of any type of asymmetry in
system parameters, the heat current remains the same upon exchanging the bath
temperatures. Therefore, imposing interaction asymmetry between the system and
the left and right baths does not result in heat rectification within this
system. To explore further, we opt for bosonic baths at the two ends instead
of fermionic ones and reiterate our analysis. This yields the subsequent
expression for heat current

$\displaystyle
I_{L}=4\sum_{\omega_{j}>0}\omega_{j}\frac{|\bm{\Phi}^{j}_{1(N)}|^{2}\kappa_{L}\kappa_{R}\big{(}N_{L}(\omega_{j})-N_{R}(\omega_{j})\big{)}}{\kappa_{L}(1+2N_{L}(\omega_{j}))+\kappa_{R}(1+2N_{R}(\omega_{j}))},$
(31)

where,

$\displaystyle N_{i}(\omega)=\frac{1}{e^{\omega/T_{i}}-1}.$ (32)

Once again, we note that the heat current is directly proportional to
$|\bm{\Phi}^{j}_{1(N)}|^{2}$. Consequently, the heat current exhibits a
significant decrease in the topological phase.

###  V.1 Heat rectification

Figure 5: Variation of left and right current with (a) $\delta$, (b) $\chi$
and variation of rectification with (c) $\delta$, (d) $\chi$. Variation of
rectification with (e) system size $N$. (f) Variation of edge and bulk
contribution to rectification with $\delta$. If not specified otherwise
$\chi=-0.5$.

When we exchange the temperatures of the bath in equation Eq. (31), the
denominator undergoes a change Segal and Nitzan (2005), while the numerator
remains constant. This implies the potential for heat rectification when
$\kappa_{L}\neq\kappa_{R}$. Considering this, we assume

$\displaystyle\kappa_{L}=\kappa(1+\chi),$
$\displaystyle\kappa_{R}=\kappa(1-\chi),$ $\displaystyle\kappa=0.1.$ (33)

The difference in heat current on the exchange of temperature is given bySegal
and Nitzan (2005)

$\displaystyle\Delta I$ $\displaystyle=I_{L}(T_{L},T_{R})+I_{L}(T_{R},T_{L})$
$\displaystyle=\sum_{\omega_{j}>0}\frac{4\omega_{j}|\bm{\Phi}^{j}_{1(N)}|^{2}\kappa\chi(1-\chi^{2})(N_{L}(\omega_{j})-N_{R}(\omega_{j}))^{2}}{(1+N_{L}(\omega_{j})+N_{R}(\omega_{j}))^{2}-\chi^{2}(N_{L}(\omega_{j})-N_{R}(\omega_{j}))^{2}}.$
(34)

Sign Convention: We adopt a positive sign convention for the left-to-right
bath current for $T_{L}>T_{R}$, and vice versa. This convention also dictates
the sign of rectification, with ”positive” indicating that the current flows
more from left to right, and vice versa. The rectification factor can be
calculated by Karg ı  _et al._ (2019b); Upadhyay  _et al._ (2021)

$\mathcal{R}=\frac{I_{L}(T_{L},T_{R})+I_{L}(T_{R},T_{L})}{\text{Max}[|I_{L}(T_{L},T_{R})|,|I_{L}(T_{R},T_{L})|]}.$
(35)

We observe that the current is suppressed from left to right when the left
bath is more strongly connected to the system, and vice versa.

In Fig. 5 (a), we observe that the behavior of heat current with respect to
the topological order parameter $\delta$ is qualitatively similar to the
fermionic case. There is a sharp decline in the heat flow as we transition
from the trivial phase into the topological phase. This decline is attributed
to the heat flow dependency on the value of $|\bm{\Phi}^{j}_{1(N)}|^{2}$.
Analyzing Fig. 5(b), we observe that the introduction of coupling strength
asymmetry results in the suppression of current from the strongly coupled
site, consistent with the findings of previous studies Segal and Nitzan
(2005). Examining the rectification in relation to $\delta$ in Fig. 5(c), we
find that based on the asymmetry of heat flow it is not possible to
differentiate between the trivial and topological phases. Lastly, Fig. 5(d)
reveals that reversing the sign of the asymmetry parameter $\chi$ results in a
reversal of the current rectification direction.

###  V.2 Diode with a few fermionic sites

The preceding analysis indicates that, for large systems, our diode’s
rectification ability is not a robust indicator for distinguishing between the
topological and trivial phases. However, when we plot rectification against
system size in Fig. 5(e), a consistent pattern emerges for small $N$:
rectification is consistently greater in the topological phase. As we
incrementally increase the system size, the rectification becomes equal in
both phases. This can be directly attributed to the presence of near-zero
energy modes in the topological phase. Accordingly, for very small system
size, our diode proves to be a reliable indicator of the system’s topological
phase. When examining the rectification curve with respect to system size $N$,
we again make a fit using the following function

$\displaystyle\mathcal{R}=\mathcal{R}_{B}+\frac{\mathcal{R}_{E}}{N}.$ (36)

We observe a sharp increase in the contribution from the edge as we enter the
topological phase. Similar to the heat flow analysis concerning bulk and edge
contributions (as shown in Fig. 3), we differentiate between the heat flow
linked to near-zero energy modes and that associated with modes at energies
$>2\delta$. We plot the contribution of edge and bulk heat flow asymmetry in
the rectification factor $\mathcal{R}$ as a function of the order parameter
$\delta$ in Fig. 5(f). It becomes evident from the results that in the trivial
phase, only bulk heat flow contributes to heat rectification. In contrast, in
the topological phase, the contribution of edge heat flow asymmetry dominates
over bulk heat flow. This observation could serve as an indicator for
distinguishing between the trivial and topological phases in the SSH model.

##  VI Conclusions

In this work, we have explored the potential of thermodynamics as a tool for
identifying topological phase transitions in the one-dimensional
Su–Schrieffer–Heeger (SSH) model. Our investigation revealed intriguing
thermodynamic signatures associated with the transition from the trivial to
the topological phase, emphasizing the unique transport properties that emerge
in the SSH model. The thermodynamically consistent global master equation
served as a powerful framework for our analysis, allowing us to precisely
quantify the suppression of heat flow as we enter the topological phase. This
decline in heat flow was linked to the reduction in transmission coefficients
of non-zero energy modes within the topological phase. We have demonstrated
that this reduction can be a valuable indicator of a phase transition in the
SSH model. Moreover, our investigation of heat flow asymmetry provided
additional insights. We observed that when employing fermionic baths, no
asymmetry was detected. However, upon substituting fermionic baths with
bosonic ones, a non-zero heat flow asymmetry emerged. Remarkably for systems
with a few fermionic sites, this asymmetry was more pronounced in the
topological phase compared to the trivial phase. Consequently, for such
systems the behavior of the heat diode serves as an additional means of
distinguishing between these two phases. While such asymmetry has been used
for rectification in earlier systems, its role as an indicator of phase
transition in SSH chains has not been thoroughly studied before.

Our findings contribute to the broader understanding of the intriguing
behavior exhibited by the SSH model, providing insights that may have
implications for various fields, including condensed matter physics and
quantum information science.

##  VII Acknowledgements

RM gratefully acknowledges financial support from Sci- ence and Engineering
Research Board (SERB), India, under the Core Research Grant (Project No.
CRG/2020/000620).

##  Appendix A Master Equation

The master equation derivation is similar to the one given in Rivas and
Martin-Delgado (2017). We begin by writing the total Hamiltonian,

$\displaystyle\hat{H}_{T}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{SB}$ (37)

where the system Hamiltonian is,

$\displaystyle\hat{H}_{S}=\hat{V}+\hat{H}_{0}=\langle C^{\dagger}|M|C\rangle$
(38)

where $M$ is a matrix,

$\displaystyle|C\rangle=(|\hat{c}_{1},\hat{c}_{2},\hat{c}_{3},\hat{c}_{4}\ldots\hat{c}_{N}\rangle)^{T}$
(39)

We can diagonalise the system Hamiltonian using the unitary transformation:

$\displaystyle M_{D}$ $\displaystyle=U^{\dagger}MU$ $\displaystyle M$
$\displaystyle=UM_{D}U^{\dagger}$ $\displaystyle M_{i,j}$
$\displaystyle=\sum_{k,l}U_{i,l}(M_{D})_{l,k}U^{\dagger}_{k,j}=\sum_{k}U_{i,k}(M_{D})_{k,k}U^{\dagger}_{k,j}$
(40)

So,

$\displaystyle H_{S}$ $\displaystyle=\langle
C^{\dagger}|M|C\rangle=\sum_{i,j}\hat{c}^{\dagger}_{i}M_{i,j}\hat{c}_{j}$
$\displaystyle=\sum_{i,j,k}\hat{c}^{\dagger}_{i}U_{i,k}(M_{D})_{k,k}U^{\dagger}_{k,j}\hat{c}_{j}$
(41)

If,

$\displaystyle\hat{a}_{k}$
$\displaystyle=\sum_{j}U^{\dagger}_{k,j}\hat{c}_{j},$
$\displaystyle\hat{a}^{\dagger}_{k}=\sum_{j}U_{j,k}\hat{c}^{\dagger}_{j}$
$\displaystyle|A\rangle$ $\displaystyle=U^{\dagger}|C\rangle$
$\displaystyle|C\rangle=U|A\rangle$ (42)

Studying the commutation relation between $\hat{a}^{\prime}s$. We know,

$\displaystyle\\{\hat{c}_{i},\hat{c}_{j}\\}=0,$
$\displaystyle\\{\hat{c}^{\dagger}_{i},\hat{c}^{\dagger}_{j}\\}=0,$
$\displaystyle\\{\hat{c}_{i},\hat{c}^{\dagger}_{j}\\}=\delta_{i,j}.$ (43)

Then,

$\displaystyle\\{\hat{a}_{i},\hat{a}_{j}\\}=0,$
$\displaystyle\\{\hat{a}^{\dagger}_{i},\hat{a}^{\dagger}_{j}\\}=0,$
$\displaystyle\\{\hat{a}_{i},\hat{a}^{\dagger}_{j}\\}=\delta_{i,j}.$ (44)

So, we can write the system Hamiltonian in the following basis:

$\displaystyle
H_{S}=\sum_{k}(M_{D})_{k,k}\hat{a}^{\dagger}_{k}\hat{a}_{k}=\sum_{k}\omega_{k}\hat{a}^{\dagger}_{k}\hat{a}_{k}$
(45)

Writing the system bath interaction Hamiltonian in the diagonal basis of the
system Hamiltonian:

$\displaystyle\hat{c}_{1}=\sum_{k}U_{1,k}\hat{a}_{k},$
$\displaystyle\hat{c}_{1}^{\dagger}=\sum_{k}U^{\dagger}_{k,1}\hat{a}^{\dagger}_{k}$
(46)

So,

$\displaystyle H^{L}_{SB}$
$\displaystyle=(\hat{c}_{1}^{\dagger}+\hat{c}_{1})\otimes\sum_{k}(g^{L}_{k}\hat{b}_{k}+g^{L}_{k}\hat{b}_{k}^{\dagger})$
$\displaystyle=\sum_{j}(U_{1,j}\hat{a}_{j}+U^{\dagger}_{j,1}\hat{a}^{\dagger}_{j})\otimes\sum_{k}(g^{L}_{k}\hat{b}^{\dagger
L}_{k}+g^{*L}_{k}\hat{b}^{L}_{k})$ (47)

where,

$\displaystyle B_{j,k}=(g^{L}_{k}\hat{b}^{\dagger
L}_{k}+g^{*L}_{k}\hat{b}^{L}_{k})$ (48)

and,

$\displaystyle\hat{Y}_{j}=\hat{Y}_{L}=\sum_{k}B_{j,k},$
$\displaystyle\hat{X}^{L}_{j}=U_{1,j}\hat{a}_{j}+U^{\dagger}_{j,1}\hat{a}^{\dagger}_{j}$
(49)

The Redfield master equation in the interaction picture is given as Breuer and
Petruccione (2002),

$\frac{d\hat{\rho}^{I}(t)}{dt}=-\int_{0}^{\infty}ds\text{
}Tr_{B}\\{[\hat{H}_{SB}^{I}(t),[\hat{H}_{SB}^{I}(t-s),\hat{\rho}_{S}^{I}(t)\otimes\hat{\rho}_{B}^{I}(0)]]\\}$
(50)

Simplifying,

$\displaystyle\frac{d\hat{\rho}^{I}(t)}{dt}$
$\displaystyle=\mathcal{D}_{L}[\hat{\rho}]+\mathcal{D}_{R}[\hat{\rho}]$
$\displaystyle=\sum_{k,l}\int_{0}^{\infty}d\tau\Gamma_{l,k}(\tau,T)[\hat{X}_{k}(t-\tau)\rho,\hat{X}_{l}(t)]+H.C.$
(51)

and,

$\displaystyle\Gamma_{k,l}(\tau,T)=Tr_{B}(\rho_{B}\hat{Y}_{k}(\tau)\hat{Y}_{l})$
(52)

Since, $\hat{Y}_{j}=\hat{Y}_{L}=\sum_{k}B_{j,k}$,

$\displaystyle\Gamma_{k,l}(\tau,T)=\Gamma(\tau,T)=Tr_{B}(\rho_{B}\hat{Y}(\tau)\hat{Y})$
(53)

$\displaystyle\mathcal{D}_{L}[\hat{\rho}]$
$\displaystyle=\sum_{k}\int_{0}^{\infty}d\tau\Gamma(\tau,T)[\hat{X}^{L}_{k}(t-\tau)\rho,\hat{X}^{k}_{l}(t)]+H.C.$
$\displaystyle=\sum_{k}|U_{1,k}|^{2}\int_{0}^{\infty}d\tau\Gamma(\tau,T)e^{i\tau\omega_{k}}[\hat{a}_{k}\rho,\hat{a}^{\dagger}_{k}]+H.C.$
$\displaystyle+\sum_{k}|U_{1,k}|^{2}\int_{0}^{\infty}d\tau\Gamma(\tau,T)e^{-i\tau\omega_{k}}[\hat{a}^{\dagger}_{k}\rho,\hat{a}_{k}]+H.C.$
(54)

Putting in the value of $\Gamma(\tau,T)$ and using the secular approximation,
we get the master equation given in Eq. (A) in the Schrodinger picture.

$\displaystyle\frac{d\rho}{dt}$
$\displaystyle=-i[H_{S},\rho]+\sum_{j}\Big{(}\kappa^{L}_{j}G(\omega_{j},T_{L})+\kappa^{R}_{j}G(\omega_{j},T_{R})\Big{)}\mathcal{D}(\hat{a}_{j})[\rho]$
$\displaystyle+\Big{(}\kappa^{L}_{j}G(-\omega_{j},T_{L})+\kappa^{R}_{j}G(-\omega_{j},T_{R})\Big{)}\mathcal{D}(\hat{a}_{j}^{\dagger})[\rho]$
(55)

##  Appendix B Heat Current Expression

The heat current,

$\displaystyle I_{L}$ $\displaystyle=Tr(\mathcal{L}_{L}H_{S})$
$\displaystyle=$ $\displaystyle
2\sum_{j}-\omega_{j}\kappa^{L}_{j}G(\omega_{j},T_{L})\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle+\omega_{j}\kappa^{L}_{j}G(-\omega_{j},T_{L})(1-\langle\hat{a}^{\dagger}_{j}\hat{a}_{j}\rangle)$
$\displaystyle=$ $\displaystyle
2\sum_{j}\frac{|\omega_{j}|(k^{L}_{j}k^{R}_{j})}{k^{L}_{j}+k^{R}_{j}}(F_{L}(|\omega_{j}|)-F_{R}(|\omega_{j}|))$
(56)

The value of $\langle\hat{a}_{k}^{\dagger}\hat{a}_{k}\rangle$ in the above
equation has been substituted from the following master equation (A), we get,

$\displaystyle\frac{d\langle\hat{a}^{\dagger}_{k}\hat{a}_{k}\rho\rangle}{dt}$
$\displaystyle=-i\langle\hat{a}^{\dagger}_{k}\hat{a}_{k}[H_{S},\rho]\rangle$
$\displaystyle+\Big{\langle}\sum_{j}\Big{(}\kappa^{L}_{j}G(\omega_{j},T_{L})+\kappa^{R}_{j}G(\omega_{j},T_{R})\Big{)}\hat{a}^{\dagger}_{k}\hat{a}_{k}\mathcal{D}(\hat{a}_{j})[\rho]$
$\displaystyle+\Big{(}\kappa^{L}_{j}G(-\omega_{j},T_{L})+\kappa^{R}_{j}G(-\omega_{j},T_{R})\Big{)}\hat{a}^{\dagger}_{k}\hat{a}_{k}\mathcal{D}(\hat{a}_{j}^{\dagger})[\rho]\Big{\rangle}$
(57)

Solving for the steady state,

$\displaystyle\langle\hat{a}^{\dagger}_{k}\hat{a}_{k}\rangle={\frac{\Big{(}\kappa^{L}_{j}G(-\omega_{j},T_{L})+\kappa^{R}_{j}G(-\omega_{j},T_{R})\Big{)}}{\Big{(}\kappa^{L}_{j}G(\omega_{j},T_{L})+\kappa^{R}_{j}G(\omega_{j},T_{R})\Big{)}+\Big{(}\kappa^{L}_{j}G(-\omega_{j},T_{L})+\kappa^{R}_{j}G(-\omega_{j},T_{R})\Big{)}}}$
$\displaystyle=\frac{\Big{(}\kappa^{L}_{j}G(-\omega_{j},T_{L})+\kappa^{R}_{j}G(-\omega_{j},T_{R})\Big{)}}{\kappa^{L}_{j}+\kappa^{R}_{j}}$
(58)

## References

  * Azzouz (2006) Mohamed Azzouz, “Field-induced quantum criticality in low-dimensional heisenberg spin systems,” Phys. Rev. B 74, 174422 (2006).
  * Zhang and Song (2015) G. Zhang and Z. Song, “Topological characterization of extended quantum ising models,” Phys. Rev. Lett. 115, 177204 (2015).
  * Miao  _et al._ (2017) Jian-Jian Miao, Hui-Ke Jin, Fu-Chun Zhang, and Yi Zhou, “Exact solution for the interacting kitaev chain at the symmetric point,” Phys. Rev. Lett. 118, 267701 (2017).
  * Roy  _et al._ (2019) Saptarshi Roy, Titas Chanda, Tamoghna Das, Debasis Sadhukhan, Aditi Sen(De), and Ujjwal Sen, “Phase boundaries in an alternating-field quantum xy model with dzyaloshinskii-moriya interaction: Sustainable entanglement in dynamics,” Phys. Rev. B 99, 064422 (2019).
  * Liu  _et al._ (2020) Zi-An Liu, Tian-Cheng Yi, Jin-Hua Sun, Yu-Li Dong, and Wen-Long You, “Lifshitz phase transitions in a one-dimensional gamma model,” Phys. Rev. E 102, 032127 (2020).
  * Luo (2022) Qiang Luo, “Analytical results for the unusual grüneisen ratio in the quantum ising model with dzyaloshinskii-moriya interaction,” Phys. Rev. B 105, L060401 (2022).
  * Molignini  _et al._ (2017) Paolo Molignini, Evert van Nieuwenburg, and R. Chitra, “Sensing floquet-majorana fermions via heat transfer,” Phys. Rev. B 96, 125144 (2017).
  * Benenti  _et al._ (2009) Giuliano Benenti, Giulio Casati, Toma ž Prosen, Davide Rossini, and Marko Žnidarič, “Charge and spin transport in strongly correlated one-dimensional quantum systems driven far from equilibrium,” Phys. Rev. B 80, 035110 (2009).
  * Bandyopadhyay  _et al._ (2021) Souvik Bandyopadhyay, Sourav Bhattacharjee, and Diptiman Sen, “Driven quantum many-body systems and out-of-equilibrium topology,” Journal of Physics: Condensed Matter 33, 393001 (2021).
  * Mondal  _et al._ (2022) Saikat Mondal, Diptiman Sen, and Amit Dutta, “Disconnected entanglement entropy as a marker of edge modes in a periodically driven kitaev chain,” Journal of Physics: Condensed Matter 35, 085601 (2022).
  * van Caspel  _et al._ (2019) Moos van Caspel, Sergio Enrique Tapias Arze, and Isaac Pérez Castillo, “Dynamical signatures of topological order in the driven-dissipative Kitaev chain,” SciPost Phys. 6, 026 (2019).
  * Rivas and Martin-Delgado (2017) Ángel Rivas and Miguel A. Martin-Delgado, “Topological heat transport and symmetry-protected boson currents,” Scientific Reports 7, 6350 (2017).
  * Kane and Moore (2011) Charles Kane and Joel Moore, “Topological insulators,” Physics World 24, 32 (2011).
  * Asbóth  _et al._ (2016) Já nos K. Asbóth, László Oroszlány, and András Pályi,  _A Short Course on Topological Insulators_ (Springer International Publishing, 2016).
  * Viyuela  _et al._ (2014) O. Viyuela, A. Rivas, and M. A. Martin-Delgado, “Uhlmann phase as a topological measure for one-dimensional fermion systems,” Physical Review Letters 112, 130401 (2014).
  * Kempkes  _et al._ (2016) S. N. Kempkes, A. Quelle, and C. Morais Smith, “Universalities of thermodynamic signatures in topological phases,” Scientific Reports 6, 38530 (2016).
  * Meier  _et al._ (2016a) Eric J. Meier, Fangzhao Alex An, and Bryce Gadway, “Observation of the topological soliton state in the su–schrieffer–heeger model,” Nature Communications 7, 13986 (2016a).
  * Drost  _et al._ (2017) Robert Drost, Teemu Ojanen, Ari Harju, and Peter Liljeroth, “Topological states in engineered atomic lattices,” Nature Physics 13, 668–671 (2017).
  * Kanungo  _et al._ (2022) S. K. Kanungo, J. D. Whalen, Y. Lu, M. Yuan, S. Dasgupta, F. B. Dunning, K. R. A. Hazzard, and T. C. Killian, “Realizing topological edge states with rydberg-atom synthetic dimensions,” Nature Communications 13, 972 (2022).
  * Meier  _et al._ (2016b) Eric J. Meier, Fangzhao Alex An, and Bryce Gadway, “Observation of the topological soliton state in the su–schrieffer–heeger model,” Nature Communications 7, 13986 (2016b).
  * Gröning  _et al._ (2018) Oliver Gröning, Shiyong Wang, Xuelin Yao, Carlo A Pignedoli, Gabriela Borin Barin, Colin Daniels, Andrew Cupo, Vincent Meunier, Xinliang Feng, Akimitsu Narita, Klaus Müllen, Pascal Ruffieux, and Roman” Fasel, “Engineering of robust topological quantum phases in graphene nanoribbons,” Nature 560, 209–213 (2018).
  * Tanaka  _et al._ (2022) Yukio Tanaka, Bo Lu, and Naoto Nagaosa, “Theory of giant diode effect in $d$-wave superconductor junctions on the surface of a topological insulator,” Phys. Rev. B 106, 214524 (2022).
  * Fluckey  _et al._ (2022) Stephen P. Fluckey, Sabyasachi Tiwari, Christopher L. Hinkle, and William G. Vandenberghe, “Three-dimensional-topological-insulator tunnel diodes,” Phys. Rev. Appl. 18, 064037 (2022).
  * Karabassov  _et al._ (2023) Tairzhan Karabassov, Emir S. Amirov, Irina V. Bobkova, Alexander A. Golubov, Elena A. Kazakova, and Andrey S. Vasenko, “Superconducting diode effect in topological hybrid structures,” Condensed Matter 8 (2023).
  * Li  _et al._ (2017) Jiaqian Li, Xiaofeng Zhou, Jing Li, Lufeng Che, Jun Yao, Glen McHale, Manoj K. Chaudhury, and Zuankai Wang, “Topological liquid diode,” Science Advances 3, eaao3530 (2017).
  * He and Chien (2023) Yan He and Chih-Chun Chien, “Particle and thermal transport through one dimensional topological systems via lindblad formalism,” Physics Letters A 473, 128826 (2023).
  * Nava  _et al._ (2023) Andrea Nava, Gabriele Campagnano, Pasquale Sodano, and Domenico Giuliano, “Lindblad master equation approach to the topological phase transition in the disordered su-schrieffer-heeger model,” Phys. Rev. B 107, 035113 (2023).
  * Levy and Kosloff (2014) Amikam Levy and Ronnie Kosloff, “The local approach to quantum transport may violate the second law of thermodynamics,” Europhys. Lett. 107, 20004 (2014).
  * Naseem  _et al._ (2018) M. Tahir Naseem, André Xuereb, and Özgür E. Müstecaplıoğlu, “Thermodynamic consistency of the optomechanical master equation,” Phys. Rev. A 98, 052123 (2018).
  * Landi  _et al._ (2014) Gabriel T. Landi, E. Novais, Mário J. de Oliveira, and Dragi Karevski, “Flux rectification in the quantum $xxz$ chain,” Phys. Rev. E 90, 042142 (2014).
  * Jiang  _et al._ (2010) J. W. Jiang, J. S. Wang, and B. Li, “Topology-induced thermal rectification in carbon nanodevice,” Europhys. Lett. 89, 46005 (2010).
  * Ordonez-Miranda  _et al._ (2017) Jose Ordonez-Miranda, Younès Ezzahri, and Karl Joulain, “Quantum thermal diode based on two interacting spinlike systems under different excitations,” Phys. Rev. E 95, 022128 (2017).
  * Karg ı  _et al._ (2019a) Cahit Karg ı, M. Tahir Naseem, Tomá š Opatrný, Özgür E. Müstecaplıoğlu, and Gershon Kurizki, “Quantum optical two-atom thermal diode,” Phys. Rev. E 99, 042121 (2019a).
  * Segal and Nitzan (2005) Dvira Segal and Abraham Nitzan, “Spin-boson thermal rectifier,” Phys. Rev. Lett. 94, 034301 (2005).
  * Su  _et al._ (1979) W. P. Su, J. R. Schrieffer, and A. J. Heeger, “Solitons in polyacetylene,” Phys. Rev. Lett. 42, 1698–1701 (1979).
  * Sirker  _et al._ (2014) J Sirker, M Maiti, N P Konstantinidis, and N Sedlmayr, “Boundary fidelity and entanglement in the symmetry protected topological phase of the SSH model,” Journal of Statistical Mechanics: Theory and Experiment 2014, P10032 (2014).
  * Ikramov (1999) Khakim D. Ikramov, “Shin’s formulas for eigenpairs of symmetric tridiagonal 2-toeplitz matrices,” Bulletin of the Australian Mathematical Society 59, 119–120 (1999).
  * Breuer and Petruccione (2002) H. P. Breuer and F. Petruccione,  _The theory of open quantum systems_ (Oxford university press, 2002).
  * Datta (2005) Supriyo Datta,  _Quantum Transport: Atom to Transistor_ (Cambridge University Press, 2005).
  * Hofer  _et al._ (2017) Patrick P Hofer, Martí Perarnau-Llobet, L David M Miranda, Géraldine Haack, Ralph Silva, Jonatan Bohr Brask, and Nicolas Brunner, “Markovian master equations for quantum thermal machines: local versus global approach,” New Journal of Physics 19, 123037 (2017).
  * Landi  _et al._ (2022) Gabriel T. Landi, Dario Poletti, and Gernot Schaller, “Nonequilibrium boundary-driven quantum systems: Models, methods, and properties,” Rev. Mod. Phys. 94, 045006 (2022).
  * Hill (1994) T.L. Hill,  _Thermodynamics of Small Systems_, Dover Books on Chemistry (Dover Publications, 1994).
  * Karg ı  _et al._ (2019b) Cahit Karg ı, M. Tahir Naseem, Tomá š Opatrný, Özgür E. Müstecaplıoğlu, and Gershon Kurizki, “Quantum optical two-atom thermal diode,” Phys. Rev. E 99, 042121 (2019b).
  * Upadhyay  _et al._ (2021) Vipul Upadhyay, M. Tahir Naseem, Rahul Marathe, and Özgür E. Müstecaplıoğlu, “Heat rectification by two qubits coupled with dzyaloshinskii-moriya interaction,” Phys. Rev. E 104, 054137 (2021).


Typical X-ray Outburst Light Curves of Aql X-1–Typical X-ray Outburst Light
Curves of Aql X-1

# Typical X-ray Outburst Light Curves of Aql X-1

Ömer Faruk Çobanand Ünal Ertan
Sabancı University E-mail<EMAIL_ADDRESS>Orhanlı  Tuzla 34956
Istanbul  Turkey

(Accepted XXX. Received YYY; in original form ZZZ; 2023)

###### Abstract

We show that a typical X-ray outburst light curve of Aql X-1 can be reproduced
by accretion onto the neutron star in the frame of the disc instability model
without invoking partial accretion or propeller effect. The knee and the
subsequent sharp decay in the X-ray light curve can be generated naturally by
taking into account the weak dependence of the disc aspect ratio, $h/r$, on
the disc mass-flow rate, $\dot{M}_{\mathrm{in}}$, in the X-ray irradiation
flux calculation. This $\dot{M}_{\mathrm{in}}$ dependence of $h/r$ only
slightly modifies the irradiation temperature profile along the hot disc in
comparison to that obtained with constant $h/r$. Nevertheless, this small
difference has a significant cumulative effect on the hot disc radius leading
to a much faster decrease in the size of the hot disc, and thereby to a
sharper decay in the X-ray outburst light curve. The same model also produces
the long-term evolution of the source consistently with its observed outburst
recurrence times and typical light curves of Aql X-1. Our results imply that
the source accretes matter from the disc in the quiescent state as well. We
also estimate that the dipole moment of the source $\mu\lesssim 2\times
10^{26}$ G cm3.

###### keywords:

accretion – accretion discs – instabilities – X-rays: binaries – stars:
individual: Aql X-1

##  1 Introduction

Low-mass X-ray binaries (LMXBs) are systems containing either a neutron star
(NS) or a black hole (BH) and a low-mass companion star with mass
$M<1M_{\odot}$. In LMXBs, mass is transferred from the companion to the
compact object through Roche-lobe overflow. Due to its angular momentum, the
matter cannot fall directly onto the star. Instead, the material forms a
geometrically thin accretion disc around the NS (Frank et al., 2002). The
matter in the disc moves in Keplerian orbits with speed
$\upsilon_{\mathrm{K}}=(GM/r)^{1/2}$ where $G$ is the gravitational constant,
$M$ is the mass of the NS, and $r$ is the radial distance from the centre.
Through viscous processes, angular momentum is carried outwards while the
matter flows inwards toward the compact object.

A small fraction of these systems are persistent X-ray sources, while the
remaining larger fraction consists of transient systems, with both neutron
stars (NSXTs) and black holes (BHXTs) that show X-ray outbursts with large
ranges of recurrence and duration time-scales (see e.g. Psaltis, 2006). The
inner disc in NSXTs is cut at a radius, $r_{\mathrm{in}}$, depending on the
mass-flow rate, $\dot{M}_{\mathrm{in}}$, of the disc and the dipole moment
$\mu=BR_{*}^{3}$ of the NS, where $B$ and $R_{*}$ are the dipole field
strength at the equator and the radius of the star respectively. At the co-
rotation radius $r_{\mathrm{co}}=(GM/\Omega_{*}^{2})^{1/3}$, the angular speed
of the NS, $\Omega_{*}$, equals the Keplerian angular speed,
$\Omega_{\mathrm{K}}$, of the disc matter. If $r_{\mathrm{in}}\leq
r_{\mathrm{co}}$, the matter can flow from the inner disc along the closed
field lines onto the NS. In the case that $r_{\mathrm{in}}>r_{\mathrm{co}}$
the matter is estimated to be thrown out from the inner disc (propeller
effect; Illarionov & Sunyaev, 1975). For a given $\mu$, the critical
$\dot{M}_{\mathrm{in}}$ level for the onset of the propeller mechanism is not
well known. For spherical accretion onto a non-rotating magnetised star,
magnetic pressure balances the ram pressure of the inflowing matter at
Alfv$\acute{\mathrm{e}}$n radius,
$r_{\mathrm{A}}\simeq(GM)^{-1/7}\mu^{4/7}\dot{M}_{\mathrm{in}}^{-2/7}$
(Davidson & Ostriker, 1973; Lamb et al., 1973). In the case of disc accretion,
$r_{\mathrm{in}}$ is conventionally estimated by equating the viscous and
magnetic stresses and found to be $r_{\xi}=\xi r_{\mathrm{A}}$ with $\sim
0.5<\xi<1$ (Ghosh & Lamb, 1979; Kluźniak & Rappaport, 2007).

It was proposed that there is a minimum critical X-ray luminosity,
$L_{\mathrm{X}}$, for a persistent LMXB such that the entire disc is kept hot
above the hydrogen ionisation temperature by the X-ray irradiation flux with
an $L_{\mathrm{X}}$ above this critical level (van Paradijs, 1996). When
$L_{\mathrm{X}}$ decrease below the critical rate, the outer cold regions of
the disc enter a low viscosity state with an inefficient mass transfer from
the outer to the inner disc with a rate lower than the mass-flow rate,
$\dot{M}_{\mathrm{c}}$, from the companion to the outer disc. This leads to a
gradual increase in the surface densities and temperatures of the outer disc.
When the outer disc temperature exceeds the critical level at a radius, a
heating front propagates to inner and outer radii, which could take most of
the disc into the "hot state" with higher viscosities enhancing the mass
transfer to the inner disc. The resultant abrupt increase in
$\dot{M}_{\mathrm{in}}$ at the inner disc is observed as an X-ray outburst. In
the outburst phase, efficient mass transfer from the outer radii decreases the
surface densities and temperatures of the outer disc. Eventually, starting
from the outer disc, propagation of a cooling front inward could take the
entire disc into the "cold state" with much lower $\dot{M}_{\mathrm{in}}$ and
$L_{\mathrm{X}}$ (quiescent state). To sum up, the outer disc behaves like a
mass reservoir in the quiescent state, and this mass is released as a result
of thermal-viscous disc instability producing an X-ray outburst (see e.g.
Lasota 2001; Frank et al. 2002 for details).

The $L_{\mathrm{X}}$ peak in outbursts and the recurrence time-scale,
$\tau_{\mathrm{rec}}$, depend on the size of the disc together with the mass
transfer rate from the companion, as well as the X-ray irradiation strength
during the outburst. The X-ray irradiation slows down the inward propagation
of the cooling front during the outburst state. Therefore, the outer disc of a
system with relatively high $L_{\mathrm{X}}$ is evacuated more efficiently
during outbursts, which needs a longer $\tau_{\mathrm{rec}}$ to refill the
disc. For instance, $\tau_{\mathrm{rec}}\sim$ weeks for dwarf novae (DNs),
while for NSXTs and BHXTs $\tau_{\mathrm{rec}}$ varies from months to decades
(Frank et al., 2002).

In the disc instability models (DIMs), the disc diffusion equation is usually
solved using the $\alpha$-prescription of the kinematic viscosity $\nu=\alpha
c_{\mathrm{s}}h$ (Shakura & Sunyaev, 1973), where $c_{\mathrm{s}}$ is the
sound speed and $h$ is the pressure scale-height of the disc. The basic
characteristics of DN light curves can be reproduced with
$\alpha_{\mathrm{h}}\simeq 0.1$ and $\alpha_{\mathrm{c}}\simeq 0.01-0.05$ for
the hot and cold viscosity states, respectively. DIMs can also account for the
outburst light curves and recurrence times of transient LMXBs with
$\alpha_{\mathrm{h}}$ and $\alpha_{\mathrm{c}}$ values similar to those used
for DNs, provided that the effect of X-ray irradiation on the disc dynamics
during the outburst states is included in the models (King, 1998; King &
Ritter, 1998; Dubus et al., 1999; Dubus et al., 2001).

A characteristic "knee" followed by a steepening in the X-ray outburst light
curves observed from a variety of magnetised stars in binaries including some
NSXTs was proposed to be the signature of the onset of the propeller effect at
a critical $L_{\mathrm{X}}$ corresponding to $r_{\xi}=\xi
r_{\mathrm{A}}=r_{\mathrm{co}}$ with $\xi\simeq 0.5$ (Zhang et al., 1998;
Campana et al., 1998; Campana et al., 2018). Among these NSXTs, Aql X-1 is a
prototypical source, accreting matter from a K-type companion. The spin period
of the source $P=1.8~{}\mathrm{ms}$ (Casella et al., 2008). In quiescence,
$L_{\mathrm{X}}\simeq 5\times 10^{33}~{}$erg s-1 for a distance $d=4.5$ kpc
(Campana et al., 2014; Galloway et al., 2008). In the outburst states, the
$L_{\mathrm{X}}$ peak can exceed $10^{37}~{}$erg s-1 within a few days
(Campana et al., 2013). After the peak of a typical outburst of Aql X-1, the
X-ray light curve shows a smooth and slow exponential decay lasting a few
weeks followed by a sharp turn downward, forming a knee-like morphology, and
subsequently $L_{\mathrm{X}}$ decreases more steeply down to the quiescent
level within several days. The outbursts are repeated with a recurrence time
$\tau_{\mathrm{rec}}\sim 1$ yr.

A steepening of the X-ray light curve during the decay phase is a behaviour
estimated in the DIM resulting from the inward propagation of the cooling
front (see e.g. King & Ritter 1998). Recently, Lipunova et al. (2022) showed
through detailed numerical analysis that the knee in the decay curve of Aql
X-1 could be produced as a result of the shrinking size of the hot inner disc
during the decay phase of the outburst. Initially, the X-ray light curve is
similar to that of a purely viscous relaxation (without disc instability). The
presence of a cold-hot border is communicated to the innermost disc after a
viscous time across the hot inner disc. While $r_{\mathrm{h}}$ is decreasing,
the viscous time-scale of the hot disc also decreases. Below a critical level,
this leads to a steeper $\dot{M}_{\mathrm{in}}$ and $L_{\mathrm{X}}$ decay
compared to those in purely viscous decay. This sequence of events governed by
the dynamics and viscous time-scale of the inner hot disc cause a transition
from a slower to a faster decay in $L_{\mathrm{X}}$, producing a knee-like
morphology in the X-ray outburst light curve. In their model, to account for
the sharp decay, Lipunova et al. (2022) included a partial accretion regime
when $r_{\mathrm{in}}>r_{\mathrm{co}}$ with a varying accretion fraction
depending on the fastness parameter,
$(r_{\mathrm{in}}/r_{\mathrm{co}})^{3/2}$, adopting $r_{\mathrm{in}}=r_{\xi}$
(see the discussion in section 2).

In this work, we also investigate the X-ray outburst light curves of Aql X-1
together with its recurrence characteristics in the frame of DIMs. In our
model, the physical reason for the formation of the knee is similar to that
found by Lipunova et al. (2022), while the functional form of the sharp decay
after the knee is produced in a way rather different from the mechanism
adopted in Lipunova et al. (2022). Without including partial accretion or a
strong propeller in the model, we show that the sharp $L_{\mathrm{X}}$ decay
can be accounted for by taking the weak $\dot{M}_{\mathrm{in}}$ dependence of
the disc aspect ratio, $h/r$, into account in the X-ray irradiation flux
calculations. In the model, repeating outbursts with similar X-ray outburst
light curves and with a recurrence time in agreement with the observations can
also be reproduced. In section 2, we describe our model briefly. We discuss
the results of the model calculations in section 3 and summarise our
conclusions in section 4.

##  2 Model

We use the numerical code employed earlier by Ertan & Alpar (2002) to simulate
the fast-rise-exponential-decay (FRED) type outburst light curves of BHXTs.
Here, we briefly describe the model calculations.

We represent the mass distribution of the disc before the X-ray outburst with
a Gaussian surface density distribution
$\Sigma(r,t=0)=\Sigma_{0}\exp\bigl{[}-\bigl{\\{}(r-r_{\mathrm{circ}}/\Delta
r\bigl{\\}}^{2}\bigr{]}$ centred at the circularisation radius,
$r_{\mathrm{circ}}$. Through numerical calculations, we solve the diffusion
equation

$\frac{\partial\Sigma}{\partial t}=\frac{3}{r}\frac{\partial}{\partial
r}\left[r^{1/2}\frac{\partial}{\partial r}(\nu\Sigma r^{1/2})\right]$ (1)

(Frank et al., 2002) to calculate $\Sigma(r,t)$ together with corresponding
$\dot{M}_{\mathrm{in}}$ evolution. We use the $\alpha$ prescription of the
kinematic viscosity, $\nu=\alpha c_{\mathrm{s}}h$ (Shakura & Sunyaev, 1973)
where $h=c_{\mathrm{s}}\Omega_{\mathrm{K}}$ is the pressure scale height of
the disc, $c_{\mathrm{s}}=\sqrt{k_{\mathrm{B}}T_{\mathrm{c}}/\mu
m_{\mathrm{p}}}$ is the sound speed, $m_{\mathrm{p}}$ is the proton mass,
$\mu$ is the mean molecular weight in units of $m_{\mathrm{p}}$, and
$T_{\mathrm{c}}$ is the mid-plane temperature of the disc. For a better
resolution of the inner disc, we solve the diffusion equation substituting
$x=2r^{1/2}$ and $S=x\Sigma$ (Frank et al., 2002). We divide the disc into 800
radial grids equally spaced in $x$-space.

The viscous dissipation rate per unit disc area (from the mid-plane to the
surface of the disc) can be written as

$D(r)=\frac{9}{8}\nu\Sigma\frac{GM}{r^{3}}=\frac{4\sigma}{3\tau}T_{\mathrm{c}}^{4}$
(2)

where $\sigma$ Stephan-Boltzmann constant, and $\tau=\kappa\Sigma$ is the
vertically integrated optical depth of the disc. For an optically thick,
geometrically thin disc, $\tau\gg 1$, and the dissipated energy is emitted
locally in the form of blackbody radiation. We use the Rosseland mean
opacities from the opacity tables for population I stars (Alexander & Ferguson
1994, for $\log T\leq 3.7$ and Iglesias & Rogers 1996, for $\log T>3.7$).

In addition to viscous dissipation, the disc is also heated by the X-rays
produced by the accretion onto the NS. The X-ray irradiation flux for a point
source is given by

$F_{\mathrm{irr}}=\sigma
T_{\mathrm{irr}}^{4}=\frac{\eta\dot{M}_{\ast}c^{2}(1-\epsilon)}{4\pi
r^{2}}\frac{h}{r}\Biggl{(}\frac{\mathrm{d\>ln}h}{\mathrm{d\>ln}r}-1\Biggr{)}=C\frac{\dot{M}_{\ast}c^{2}}{4\pi
r^{2}}$ (3)

(Shakura & Sunyaev, 1973; Dubus et al., 1999) where $\eta$ is the efficiency
of the conversion of rest mass energy into radiation and $\epsilon$ is the
X-ray albedo of the disc surface, and

$C\equiv\eta(1-\epsilon)\frac{h}{r}\Biggl{(}\frac{\mathrm{d\>ln}h}{\mathrm{d\>ln}r}-1\Biggr{)}$
(4)

is the irradiation parameter which is estimated to be in the range of
$10^{-4}-10^{-3}$ from the model fits to the optical and X-ray spectra of
LMXBs (de Jong et al., 1996; Dubus et al., 1999; Ertan & Alpar, 2002). Note
that if $F_{\mathrm{irr}}$ is defined in terms of
$L_{\mathrm{X}}=\eta\dot{M}_{\ast}c^{2}$, then the irradiation parameter
becomes greater than $C$ given by equation (4) by a factor of $1/\eta$.

The effective temperature of the disc
$T_{\mathrm{eff}}\simeq\left[(D+F_{\mathrm{irr}})/\sigma\right]^{1/4}$. At the
inner disc regions, $D$ dominates $F_{\mathrm{irr}}$, while outside a radius
($\sim 10^{9}$ cm) the X-ray irradiation is the dominant heating mechanism.
Since the $h/r$ ratio in the equation (4) is not very sensitive to
$\dot{M}_{\mathrm{in}}$ and
$r\>\bigl{(}h/r\propto\dot{M}_{\mathrm{in}}^{3/20}r^{1/8}\bigr{)}$, it seems
reasonable to take $C$ constant in the calculations (Frank et al., 2002).
Nevertheless, in our simulations, we noticed that the variation of $h/r$ with
$\dot{M}_{\mathrm{in}}$ and the resultant change in $F_{\mathrm{irr}}$
significantly modify the slope of the X-ray light curve in the sharp decay
phase after the knee. Including the weak $\dot{M}_{\mathrm{in}}$ and $r$
dependence of $h/r$, equation (4) can be written as

$C=C_{0}\>\dot{M}_{\mathrm{in},17}^{3/20}\>r_{10}^{1/8}$ (5)

where $\dot{M}_{\mathrm{in},17}=(\dot{M}_{\mathrm{in}}/10^{17}~{}$g s-1),
$~{}r_{10}=(r/10^{10}~{}\mathrm{cm})$, and $C_{0}$ will be a free parameter of
our model.

It was shown earlier through detailed analysis that the disc structure is not
affected by the X-ray irradiation (Dubus et al., 1999; Lasota, 2001). This
means that the irradiation flux does not change the $h/r$ ratio and its
$\dot{M}_{\mathrm{in}}$ dependence along the hot portion of the disc.
Nevertheless, the irradiation significantly affects the stability criteria of
the disc by reducing the minimum critical surface densities for the transition
to the cold state (King, 1998; King & Ritter, 1998; Dubus et al., 1999; Dubus
et al., 2001; Lasota, 2001).

At the peak of the X-ray outburst, a large portion of the disc is irradiated.
With decreasing $L_{\mathrm{X}}$ in the decay phase, the radius of the hot
disc, $r_{h}$, decreases as well. The irradiation flux, $F_{\mathrm{irr}}$,
determines the current position of the hot disc radius, $r_{h}$, and governs
the speed of its inward propagation. In this phase, the X-ray light curve
characteristics are produced mainly by the dynamics of the hot inner disc.
Meanwhile, the outer cold regions of the disc that remains outside $r_{h}$
evolve with purely viscous decay, without any significant effect on the decay
light curve. In other words, the $\dot{M}_{\mathrm{in}}$ dependence of the
disc aspect ratio enters the model calculations across the hot disc only which
could be approximated by a steady disc profile corresponding to the current
mass-flow rate of the disc.

Following the conventional approach, we calculate the kinematic viscosity with
two different $\alpha$ parameters, $\alpha_{\mathrm{h}}$ and
$\alpha_{\mathrm{c}}$, for the hot and cold viscosity states, respectively.
For a given $r$ in the cold (hot) state, there is a maximum (minimum) critical
surface density $\Sigma_{\mathrm{max}}$ ($\Sigma_{\mathrm{min}}$) for the
transition to the hot (cold) state. We adopt $\Sigma_{\mathrm{max}}$ and
$\Sigma_{\mathrm{min}}$ equations obtained by Dubus et al. (2001) through
detailed vertical disc analyses. Interpolating the critical $\Sigma$ values
found for different radii of the disc, Dubus et al. (2001) obtained

$\Sigma_{\mathrm{max}}=(10.8-10.3\>\xi)\>\alpha_{\mathrm{c}}^{-0.84}~{}M_{1}^{-0.37+0.1\xi}~{}r_{10}^{1.11-0.27\xi}~{}\mathrm{g\,cm^{-2}},$
(6)
$\Sigma_{\mathrm{min}}=(8.3-7.1\>\xi)\>\alpha_{\mathrm{h}}^{-0.77}M_{1}^{-0.37}~{}r_{10}^{1.12-0.23\xi}~{}\mathrm{g\,cm^{-2}}$
(7)

where $\xi=(T_{\mathrm{irr}}/10^{4}~{}\mathrm{K})^{2}$, and $M_{1}$ is the
mass of the NS in solar masses. For our initial mass distribution,
$L_{\mathrm{X}}\simeq 0$ until the surface densities increase at the inner
disc, which is not realistic, since $L_{\mathrm{X}}$ does not decrease below
the quiescence level. Considering these situations, in addition to the
critical surface densities, we also impose the $\alpha=\alpha_{\mathrm{h}}$
condition for radii at which $T_{\mathrm{eff}}>6000$ K, independent of
$T_{\mathrm{irr}}$.

The X-ray luminosity produced by accretion onto the star is related to the
total disc luminosity through

$L_{\mathrm{acc}}\simeq\frac{GM\dot{M}_{\ast}}{R_{*}}\simeq\frac{2r_{\mathrm{in}}}{R_{*}}~{}L_{\mathrm{disc}}$
(8)

when $r_{\mathrm{in}}$ is not close to $R_{*}$. A large fraction of the disc
luminosity is emitted from the innermost disc regions. In our numerical
calculations, we take $r_{\mathrm{in}}=r_{\mathrm{co}}$, and
$\dot{M}_{\ast}=\dot{M}_{\mathrm{in}}$ throughout the outburst state. The
exact position of $r_{\mathrm{in}}$ does not significantly affect the
$\dot{M}_{\mathrm{in}}$ evolution. During the transition to the propeller
phase, $L_{\mathrm{X}}$ decreases rapidly by a factor greater than $\sim
2r_{\mathrm{co}}/R_{*}\simeq 5$, since $r_{\mathrm{in}}$ should be greater
than $r_{\mathrm{co}}$ in the propeller phase. Our results show that such a
sharp modification of the model curve corresponding to accretion/propeller
transition is not required to account for the observed outburst behaviour of
the source.

The maximum inner disc radius at which a steady propeller mechanism can be
established could be written as

$R_{\mathrm{in,max}}^{25/8}~{}\Big{|}1-R_{\mathrm{in,max}}^{-3/2}\Big{|}~{}\simeq~{}1.26~{}\alpha_{-1}^{2/5}~{}M_{1.4}^{-7/6}~{}\dot{M}_{\mathrm{in,16}}^{-7/20}~{}\umu_{30}~{}P^{-13/12}$
(9)

(Ertan, 2017, 2018), where
$R_{\mathrm{in,max}}=r_{\mathrm{in,max}}/r_{\mathrm{co}}$,
$M_{1.4}=(M/1.4M_{\odot})$,
$\dot{M}_{\mathrm{in,16}}=(\dot{M}_{\mathrm{in}}/10^{16}~{}$g s-1),
$\umu_{30}=(\umu/10^{30}$ G cm3), $\alpha_{-1}=(\alpha/0.1)$, and $P$ is the
spin period of the NS. Ertan (2017) derived this formula through analytical
calculations adopting the basic results of the disc-magnetosphere interaction
model proposed by Lovelace et al. (1995), which was developed later to
simulate the propeller phase as well (Lovelace et al., 1999; Ustyugova et al.,
2006). In this model, the inner disc settles down at a radius where the field
lines can force the innermost disc matter into co-rotation. The inner disc-
field interaction takes place in a narrow boundary with radial width $\Delta
r<r$. The field lines interacting with matter inflate and open up within an
interaction time-scale
$\tau_{\mathrm{int}}=\left|\Omega_{\mathrm{K}}-\Omega_{\ast}\right|^{-1}$. In
the propeller phase, the matter can be thrown out along the open field lines.
Subsequently, the lines reconnect on a time-scale similar to
$\tau_{\mathrm{int}}$. The field lines outside the boundary are open and
decoupled from the disc. With these conditions, Ertan (2017) found that the
inflowing disc matter can be stopped at $r_{\mathrm{A}}$, but cannot be thrown
out from the system efficiently. The resultant pile-up pushes the inner disc
inwards to the radius at which the field is strong enough to expel all the
inflowing mass from the inner disc, continuously evacuating the inner
boundary. Equation (9) gives the maximum radius at which this strong propeller
condition is satisfied. In the propeller phase, $r_{\mathrm{in}}=r_{\eta}$
does not scale with and in some cases significantly smaller than
$r_{\mathrm{A}}$ (see Ertan 2017, 2018, 2021 for details).

The actual inner disc radius is estimated to be close to
$r_{\mathrm{in,max}}$, and can be written as $r_{\eta}=\eta
r_{\mathrm{in,max}}$ where $\eta$ is a parameter close to unity. The
$L_{\mathrm{X}}$ level corresponding to strong propeller/accretion transition
can be estimated from equation (9). The predictions of this model for the
critical $L_{\mathrm{X}}$ and the torque variations are in agreement with the
properties of transitional millisecond pulsars (tMSPs; see Papitto & de
Martino 2022 for a recent review) during their transitions between the LMXB
and the radio millisecond pulsar (RMSP) states (Ertan, 2017, 2018). The model
was recently developed to include all the rotational phases of the NSs in
LMXBs (Ertan, 2018), and applied to the torque-reversal properties of LMXBs as
well (Gençali et al., 2022). We use equation 9 to estimate the dipole field
strength of Aql X-1 that is consistent with the results of our X-ray light
curve analysis.

We also perform simulations to analyse and compare the repeating outburst
light curves produced in the model with the observed recurrent outbursts of
Aql X-1. To represent the mass flow from the companion, we add mass into the
radial grid at $r_{\mathrm{circ}}$ at each time step. For Aql X-1,
$r_{\mathrm{co}}\simeq 2.5\times 10^{6}~{}\mathrm{cm}$,
$r_{\mathrm{circ}}=5.0\times 10^{10}~{}\mathrm{cm}$, and we take
$r_{\mathrm{out}}=1.8\times 10^{11}~{}\mathrm{cm}$ in our calculations. We
have repeated the numerical calculations with different mass-flow rates,
$\dot{M}_{\mathrm{c}}$, from the companion until we obtain typical outburst
morphology and recurrence time of Aql X-1. In section 3, we discuss the
results of the model calculations summarised above.

##  3 Results and Discussion

###  3.1 The 2010 Outburst of Aql X-1

The model curve in Figure 1 is obtained with $\alpha_{\mathrm{h}}=0.1$,
$\mu=0.6$ in the hot state, and $\alpha_{\mathrm{c}}=5\times 10^{-3}$,
$\mu=0.9$ in the cold state. The model can reproduce the basic features of the
observed X-ray light curve with $C_{0}\simeq 1.3\times 10^{-4}$. The
fluctuations seen in the decay phase are not addressed in our model.
Therefore, a seemingly good fit to data is sufficient for our purpose in this
work.

The sharp decay in $L_{\mathrm{X}}$ starts when $r_{\mathrm{h}}$ decreases to
about $10^{10}~{}\mathrm{cm}$. After the knee, the model cannot account for
the functional form of the sharp decay curve without including the local
$\dot{M}_{\mathrm{in}}$ dependence of the $h/r$ ratio and its effect on
$F_{\mathrm{irr}}$. We smooth out the realistic opacities to prevent numerical
fluctuations. The $L_{\mathrm{X}}$ fluctuations during the quiescent state are
due to continuous mass flow toward the inner region of the cold outer disc and
the resultant local instabilities continually taking part of this region into
the hot state. We calculate the entire $L_{\mathrm{X}}$ curve with
$\dot{M}_{\ast}=\dot{M}_{\mathrm{in}}$. That is, $L_{\mathrm{X}}$ is produced
by accretion onto the star without invoking any partial accretion or strong
propeller. In this case, we can estimate an upper limit to the dipole field
strength $B$ on the surface of the star. Using equation (9) and the strong
propeller condition $r_{\eta}>1.26r_{\mathrm{co}}$ (see Ertan 2021 for
details), the minimum $\dot{M}_{\mathrm{in}}$ level of the source ($\sim
3\times 10^{13}~{}$g s-1) gives $B<2.3\times 10^{8}$ G.

A model including the $\dot{M}_{\mathrm{in}}$ dependence of the disc aspect
ratio in the $F_{\mathrm{irr}}$ calculation produces a much sharper decay
curve after the knee in comparison with a model neglecting this dependence.
Figure 1 shows the difference in the model curves obtained for these two
cases.

How does a very weak $\dot{M}_{\mathrm{in}}$ and $r$ dependence of $h/r$ yield
the difference in $L_{\mathrm{X}}$ curves seen in Figure 1? This can be
understood by comparing hot disc radii, $r_{\mathrm{h}}^{\prime}$ and
$r_{\mathrm{h}}$, obtained with constant and varying $h/r$, respectively. From
equation 3, it can be shown that
$r_{\mathrm{h}}^{\prime}\propto\dot{M}_{\mathrm{in}}^{1/2}$. Using C given in
equation 5, we obtain $r_{\mathrm{h}}\propto\dot{M}_{\mathrm{in}}^{46/75}$.
During the decay phase, both $r_{\mathrm{h}}$ and $r_{\mathrm{h}}^{\prime}$
decrease, while the $r_{\mathrm{h}}^{\prime}/r_{\mathrm{h}}$ ratio increases.
For comparison, it is seen in Figure 1 that $r_{\mathrm{h}}\simeq
r_{\mathrm{h}}^{\prime}$ when $L_{\mathrm{X}}\simeq 10^{37}~{}$erg s-1. With
two orders decrease in $L_{\mathrm{X}}$ to $\sim 10^{35}~{}$erg s-1
$r_{\mathrm{h}}^{\prime}/r_{\mathrm{h}}$ increases to $1.7$.

It is the dynamics of the hot inner disc that determines the X-ray outburst
light curve morphology. The matter at the outermost regions of the hot disc
contributes to accretion onto the star after one viscous time across the disc.
Relatively fast inward propagation of $r_{\mathrm{h}}$ also reduces the amount
of hot matter more rapidly, which produces a cumulative effect leading to a
sharper decline in $L_{\mathrm{X}}$, and thus to an even faster-shrinking size
of the hot disc. For the case with constant $C$, the disc evolves with a
larger hot disc size, that is, with more matter that will accrete to the star
during the outburst state. This produces an outburst with relatively high
$L_{\mathrm{X}}$ and long duration.

The model agrees with the observed $L_{\mathrm{X}}$ level in quiescence as
well. The observed $L_{\mathrm{X}}$ is not smooth in quiescence. There are
flares with ten-fold variations (Coti Zelati et al., 2014; Bernardini et al.,
2013). These $L_{\mathrm{X}}$ fluctuations are also seen in our model curve.
Nevertheless, correct modelling of this behaviour is not easy due to the
resolution problem for the innermost disc for a disc size greater than
$10^{11}$ cm, which is beyond the scope of this work. We should also note that
there are differences in the observed $L_{\mathrm{X}}$ levels after different
outbursts (Campana et al., 2014).

Figure 1: X-ray outburst light curve of Aql X-1 taken from Campana et al.
(2014). The solid model curve is obtained with $\alpha_{\mathrm{c}}=5\times
10^{-3},\alpha_{\mathrm{h}}=0.1,\mu_{\mathrm{c}}=0.9,\mu_{\mathrm{h}}=0.6$,
and $C_{0}=1.3\times 10^{-4}$. For the $L_{\mathrm{X}}$ calculation, we take
$\dot{M}_{\ast}=\dot{M}_{\mathrm{in}}$ for the entire model curve. For
comparison, the dotted curve is obtained with the same parameters, but with a
constant irradiation parameter $C=1.3\times 10^{-4}$.

###  3.2 The Quiescent State and Recurrent Outbursts

Using the same model parameters producing the model curve in Figure 1, we have
continuously added mass into the radial grid at $r_{\mathrm{circ}}=5\times
10^{10}$ cm to represent the mass flow into the outer disc from the companion.
The outburst light curves seen in Fig 2 are obtained with
$\dot{M}_{\mathrm{c}}=4\times 10^{16}~{}$g s-1 which produces outbursts with
$\tau_{\mathrm{rec}}\sim$ 1 year. It is interesting that the model light
curves of the second and fourth outbursts also seemingly fit well to the X-ray
outburst data of Aql X-1. For instance, Figure 3 zooms in on the second
outburst seen in Figure 2. Simulations with higher values of
$\alpha_{\mathrm{c}}$ produce outburst light curves with different
morphologies, including those similar to the light curve seen in Figure 1.

Figure 2: The same model parameters are used for the solid curve in Figure 1.
This model curve is obtained by injecting mass into the radial grid at
$r=r_{\mathrm{circ}}$ with a rate of $4\times 10^{16}~{}$g s-1. The recurrence
time of the outbursts is approximately one year. The second outburst seen in
this figure is enlarged in Figure 3

In the quiescent state of Aql X-1, $L_{\mathrm{X}}\simeq 5\times
10^{33}~{}$erg s-1. For comparison, tMSPs are estimated to accrete matter from
the disc down to $L_{\mathrm{X}}$ levels of $\sim 10^{33}~{}$erg s-1
(Archibald et al., 2009; Papitto & de Martino, 2022; Bassa et al., 2014).
Below this level, X-ray pulses are switched off, and they show RMSP behaviour.
Dipole field strengths of tMSPs seem to be similar to that indicated by the
upper limit $B<2.3\times 10^{8}$ G estimated for Aql X-1 in this work.
Considering also the similarities in the rotational properties, Aql X-1 could
have a critical $\dot{M}_{\mathrm{in}}$ rate for termination of accretion and
transition to the propeller phase that is similar to that of tMSPs. This might
be the reason for the lack of radio pulsations from Aql X-1 during its
quiescent state. It is important to note here that the minimum pulsed
$L_{\mathrm{X}}$ level of tMSPs is about two orders of magnitude smaller than
the rate corresponding to $r_{\xi}=r_{\mathrm{co}}$.

Figure 3: The second outburst seen in Figure 2.

##  4 Conclusions

We have shown that the characteristic X-ray outburst light curve features of
Aql X-1 can be reproduced by the mass accretion onto the neutron star along
the entire X-ray outburst phase without requiring a transition into the
propeller phase. The knee at the end of the initial slow exponential decay
comes about as the cold/hot border of the disc approaches the inner disc,
which is in line with the results obtained by Lipunova et al. (2022). The
sharp decay phase after the knee can be accounted for by taking the weak
$\dot{M}_{\mathrm{in}}$ and $r$ dependence of the $h/r$ ratio into account in
the $F_{\mathrm{irr}}$ calculation. Throughout the outburst phase, we take
$\dot{M}_{\mathrm{in}}=\dot{M}_{\mathrm{acc}}$ without including any partial
accretion regime or propeller effect in the $L_{\mathrm{X}}$ calculations.
Injecting mass into the grid at $r_{\mathrm{circ}}$ with a rate of $\sim
4\times 10^{16}~{}$g s-1, the model can also generate recurrent X-ray
outbursts that mimic typical outburst light curves of Aql X-1 consistently
with the observed $\tau_{\mathrm{rec}}\sim 1$ yr. We have also estimated that
$B<2.3\times 10^{8}$ G corresponding to $r_{\mathrm{in}}\leq r_{\mathrm{co}}$
for the lowest $L_{\mathrm{X}}$ level of the source.

## Acknowledgements

We acknowledge research support from TÜBİTAK (The Scientific and Technological
Research Council of Turkey) through grant 120F329 and from Sabancı University.
We thank Ali Alpar for the useful comments that improved our manuscript.

## Data Availability

No new data were analysed in support of this letter.

## References

  * Alexander & Ferguson (1994) Alexander D. R., Ferguson J. W., 1994, ApJ, 437, 879
  * Archibald et al. (2009) Archibald A. M., et al., 2009, Science, 324, 1411
  * Bassa et al. (2014) Bassa C. G., et al., 2014, MNRAS, 441, 1825
  * Bernardini et al. (2013) Bernardini F., Cackett E. M., Brown E. F., D’Angelo C., Degenaar N., Miller J. M., Reynolds M., Wijnands R., 2013, MNRAS, 436, 2465
  * Campana et al. (1998) Campana S., Stella L., Mereghetti S., Colpi M., Tavani M., Ricci D., Dal Fiume D., Belloni T., 1998, ApJ, 499, L65
  * Campana et al. (2013) Campana S., Coti Zelati F., D’Avanzo P., 2013, MNRAS, 432, 1695
  * Campana et al. (2014) Campana S., Brivio F., Degenaar N., Mereghetti S., Wijnands R., D’Avanzo P., Israel G. L., Stella L., 2014, MNRAS, 441, 1984
  * Campana et al. (2018) Campana S., Stella L., Mereghetti S., de Martino D., 2018, A&A, 610, A46
  * Coti Zelati et al. (2014) Coti Zelati F., Campana S., D’Avanzo P., Melandri A., 2014, MNRAS, 438, 2634
  * Davidson & Ostriker (1973) Davidson K., Ostriker J. P., 1973, ApJ, 179, 585
  * Dubus et al. (1999) Dubus G., Lasota J.-P., Hameury J.-M., Charles P., 1999, MNRAS, 303, 139
  * Dubus et al. (2001) Dubus G., Hameury J. M., Lasota J. P., 2001, A&A, 373, 251
  * Ertan (2017) Ertan Ü., 2017, MNRAS, 466, 175
  * Ertan (2018) Ertan Ü., 2018, MNRAS, 479, L12
  * Ertan (2021) Ertan Ü., 2021, MNRAS, 500, 2928
  * Ertan & Alpar (2002) Ertan Ü., Alpar M. A., 2002, A&A, 393, 205
  * Frank et al. (2002) Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics: Third Edition
  * Galloway et al. (2008) Galloway D. K., Muno M. P., Hartman J. M., Psaltis D., Chakrabarty D., 2008, ApJS, 179, 360
  * Gençali et al. (2022) Gençali A. A., et al., 2022, A&A, 658, A13
  * Ghosh & Lamb (1979) Ghosh P., Lamb F. K., 1979, ApJ, 234, 296
  * Iglesias & Rogers (1996) Iglesias C. A., Rogers F. J., 1996, ApJ, 464, 943
  * Illarionov & Sunyaev (1975) Illarionov A. F., Sunyaev R. A., 1975, A&A, 39, 185
  * King (1998) King A. R., 1998, MNRAS, 296, L45
  * King & Ritter (1998) King A. R., Ritter H., 1998, MNRAS, 293, L42
  * Kluźniak & Rappaport (2007) Kluźniak W., Rappaport S., 2007, ApJ, 671, 1990
  * Lamb et al. (1973) Lamb F. K., Pethick C. J., Pines D., 1973, ApJ, 184, 271
  * Lasota (2001) Lasota J.-P., 2001, New A Rev., 45, 449
  * Lipunova et al. (2022) Lipunova G., Malanchev K., Tsygankov S., Shakura N., Tavleev A., Kolesnikov D., 2022, MNRAS, 510, 1837
  * Lovelace et al. (1995) Lovelace R. V. E., Romanova M. M., Bisnovatyi-Kogan G. S., 1995, MNRAS, 275, 244
  * Lovelace et al. (1999) Lovelace R. V. E., Romanova M. M., Bisnovatyi-Kogan G. S., 1999, ApJ, 514, 368
  * Papitto & de Martino (2022) Papitto A., de Martino D., 2022, in Bhattacharyya S., Papitto A., Bhattacharya D., eds, Astrophysics and Space Science Library Vol. 465, Astrophysics and Space Science Library. pp 157–200 (arXiv:2010.09060), doi:10.1007/978-3-030-85198-9_6
  * Psaltis (2006) Psaltis D., 2006, in , Vol. 39, Compact stellar X-ray sources. pp 1–38
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Ustyugova et al. (2006) Ustyugova G. V., Koldoba A. V., Romanova M. M., Lovelace R. V. E., 2006, ApJ, 646, 304
  * Zhang et al. (1998) Zhang S. N., Yu W., Zhang W., 1998, ApJ, 494, L71
  * de Jong et al. (1996) de Jong J. A., van Paradijs J., Augusteijn T., 1996, A&A, 314, 484
  * van Paradijs (1996) van Paradijs J., 1996, ApJ, 464, L139