# A note on an open-source toolbox for simulation of acoustic waves: inclusion
of time-varying source

Ashkan Javaherian<EMAIL_ADDRESS>Department of Medical Physics &
Biomedical Engineering,
University College London, London, UK. WC1E 6BT

(Date: January 2023)

## Abstract

Simulating propagation of acoustic waves via solving a system of three-coupled
first-order linear differential equations using a k-space pseudo-spectral
method is popular for biomedical applications, firstly because of availability
of an open-source toolbox for implementation of this numerical approach, and
secondly because of its efficiency. The k-space pseudo-spectral method is
efficient, because it allows coarser computational grids and larger time steps
than finite difference and finite element methods for the same accuracy. The
goal of this study is to compare this numerical wave solver with an analytical
solution to the wave equation using the Green’s function for computing
propagation of acoustic waves in homogeneous media. This comparison is done in
the frequency domain. Using the k-Wave solver, a match to the Green’s function
is obtained after modifying the approach taken for including mass source in
the linearised equation of continuity (conservation of mass) in the associated
system of wave equations.

##  1\. Introduction

A numerical approach based on a k-space pseudo-spectral method was proposed
for solving the coupled first-order linear differential equations for wave
propagation [1, 2]. The proposed approach is unconditionally stable for the
sound speeds below a nominal sound speed used in the k-space scheme, and
accurate for weakly heterogeneous media like soft tissues [1, 2]. Using
pseudo-spectral approaches, the spatial gradient of a field is approximated
using the information from the entire field [1]. This is in contrast with
finite difference or finite element methods, in which the spatial gradient at
each point is computed based on only the information on nearby points [3, 4].
Furthermore, a k-space correction to the spatial gradients allows using larger
time steps than pseudo-spectral methods for the same accuracy [1, 2, 5].

An open-source toolbox for a numerical implementation of this approach is
popular in the fields related to biomedical acoustics [5, 6]. The goal of this
study is to numerically investigate accuracy of the k-space pseudo-spectral
approach for solving the three-coupled first-order linear wave equations in
homogeneous media. For measuring accuracy of this numerical approach, the
acoustic waves approximated using this approach will be compared to an
analytical solution to the wave equation using the Green’s function [7, 8].

In general, the acoustic sources can be instant or time varying. One example
of instant acoustic sources are formation of sound waves following light
absorption of a material sample, and is used in photoacoustic tomography [9,
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21] or designing ultrasound
transducers [22]. Time-varying ultrasound sources are widely used. Some
important applications of time-varying sources include quantitative
reconstruction of acoustic properties of tissues using full-Wave inversion
[23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38], B-mode
reflection imaging [39, 21, 40] or focused ultrasound [41, 42, 43, 44].

The analysis in this mauscript will be done for approximating propagation of
acoustic waves induced by an excitation of a point source using a time-varying
ultrasonic pulse. The k-Wave toolbox (version 1.4) is used for performing the
simulation experiments in 2D and 3D cases and for different excitation pulses
and Courant-Friedrichs-Lewy (CFL) numbers [5, 6]. The results show that the
numerical approach proposed in [1, 2] is accurate for homogeneous media, i.e.,
the acoustic waves approximated using this numerical approach match the
analytic solution to the wave equation using the Green’s function for
homogeneous media. More importantly, for obtaining this match, the approach
taken for including source in the equation of continuity (conservation of
mass) in the system of three-coupled first-order linear wave equations in the
k-Wave toolbox is revisited and modified. (cf. Eq. 2.19 in the k-Wave user
manual [45].)

Section 2 first describes the system of three-coupled first order linear wave
equations for propagation of acoustic waves, and then a numerical
implementation of this system using a k-space pseudo-spectral method is
described, and an inclusion of time-varying sources to this system of
equations is modified and explained. An analytic solution to the wave equation
using the Green’s function, which will be used as the benchmark, is also
explained in this section. The numerical results will be given in section 3. A
discussion of the results follows in section 4.

##  2\. Wave equation for homogeneous media

###  2.1. Wave equation in the continuous domain

This section explains formula describing the acoustic wave propagation in
homogeneous media. Let $\bm{x}=(x^{1},...,x^{d})$ denote a spatial position in
$\mathbb{R}^{d}$ with $d$ the number of dimensions. ($d\in\\{2,3\\}$.) An
emission element is excited by an ultrasonic pulse $s(\bm{x},t)$ in times
$t\in(0,T_{s})$. Here, the emission element is assumed isotropic and a point
at position $\bm{x}_{e}$ for simplicity, but an extension to an emission
element with finite volume is straightforward. The propagation of the induced
pressure can be defined using the second-order wave equation in the form

$\displaystyle\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}p(\bm{x},t)-\nabla^{2}p(\bm{x},t)=s(\bm{x},t),$ (1)

where $p$ denotes the acoustic pressure, and $c$ is the sound speed which is
assumed homogeneous. By ignoring the dependence on $\bm{x}$ and $t$, $p$ is
the scalar acoustic pressure with units of $\text{kg}\ \text{m}^{2-d}\
\text{s}^{-2}$, and the excitation pulse $s$ is a scalar with units of
$\text{kg}\text{m}^{-d}\text{s}^{-2}$. The three-coupled first-order linear
acoustic propagation equations can be defined in the form

$\displaystyle\begin{split}\frac{\partial}{\partial
t}\bm{u}(\bm{x},t)&=-\frac{1}{\rho_{0}}\nabla p(\bm{x},t)\\\
\frac{\partial}{\partial
t}\bm{\rho}(\bm{x},t)&=-\rho_{0}\nabla\bm{u}(\bm{x},t)+s_{M}(\bm{x},t)\\\
p(\bm{x},t)&=c(\bm{x})^{2}\bm{\rho}(\bm{x},t),\end{split}$ (2)

Here, $\rho_{0}$ is the ambient density with units of
$\text{kg}\text{m}^{-d}$, and is assumed homogeneous. By ignoring the
dependence on $\bm{x}$ and $t$, $\bm{\rho}$ and $\bm{u}$ denote the acoustic
density and particle velocity vectors with components in units of
$\text{kg}\text{m}^{-d}$ and $\text{m}\text{s}^{-1}$, respectively. Also,
$s_{M}(\bm{x},t)$ denotes the time-varying mass source which represents the
time rate of the input of mass per unit volume, and is in units of
$\text{kg}\text{m}^{-d}\text{s}^{-1}$. The relation between the mass source
and the excitation pulse satisfies

$\displaystyle s(\bm{x},t)=\frac{\partial s_{M}(\bm{x},t)}{\partial t}.$ (3)

Here, for a point source at position $\bm{x}_{e}$, the excitation is defined
as

$\displaystyle s(\bm{x},t)=\delta(\bm{x}-\bm{x}_{e})\ \mathfrak{s}(t),$ (4)

where $\mathfrak{s}(t)$ has a time duration $T_{s}$. Also, the pressure field
$p(\bm{x},t)$ is recorded for time $T$ with $T\gg T_{s}$.

###  2.2. Wave equation in the discretised domain

This section explains a discretisation of the wave equation (2). Let
$\bm{X}=(X^{1},...,X^{d})$ denote the position of the grid points for a
$d$-dimensional grid with $\bm{X}^{\zeta}$ the position $\bm{X}$ along the
Cartesian coordinate $\zeta\in\\{1,...,d\\}$. We also use
$n\in\\{1,...,N_{s},...N_{T}\\}$ for the sampled times within the times
$t\in\\{0,...,T_{s},...,T\\}$. (The same notations will be used for quantities
in the continuous and discretised domains, and we hope this is not confusing
for the readers.) Accordingly, the excitation pulse in the discretised domain
is defined by $s(\bm{X},n)$.

In this manuscript, compared to the k-Wave [5, 6, 45], two corrections are
applied for including source in the acoustic wave equation. At first, by
ignoring the dependence on time, in an ideal case when a point source at
position $\bm{x}_{0}$ matches the point $\bm{X}_{0}$ on the grid, the relation
between the point source in the continuous and discretised domains satisfies

$\displaystyle s(\bm{x})\delta(\bm{x}-\bm{x}_{0})\equiv\frac{1}{(\Delta
x)^{d}}s(\bm{X}_{0}),$ (5)

where $\Delta x$ is the grid spacing. The scaling factor accounts for spatial
sampling of dirac delta function in a $d$-dimensional space.

In addition, using (3), the relation between the discretised directional
components of the mass source and the discretised excitation pulse satisfies

$\displaystyle s_{M}^{\zeta}(\bm{X},n)=\frac{\Delta
t}{d}\sum_{n^{\prime}=1}^{n}s(\bm{X},n^{\prime}),$ (6)

A discretisation of the system of equations (2) on a grid staggered in space
and time is outlined in Algorithm 1.

Algorithm 1 Temporal iterations for the system of coupled first-order wave
equations

1:input: $c,\rho_{0},\Delta t,\Lambda^{\zeta}\
(\zeta\in\\{1,..,d\\}),s_{M}^{\zeta}(\bm{X},n+\frac{1}{2})\
(n\in\\{1,...,N_{T}-1\\})$

2:initialise: $p(\bm{X},1)=0,\
\bm{\rho}(\bm{X},1)=0,\bm{u}(\bm{X},\frac{1}{2})=0,\ \forall\bm{X}$
$\triangleright$ Set initial parameters

3:Start temporal iterations at $n=1$ $\triangleright$ Iterate

4:Record $p(\bm{X},n)$ on receivers

5:$\bm{u}^{\zeta}(\bm{X},n+\frac{1}{2})\leftarrow\Lambda^{\zeta}\Big{[}\Lambda^{\zeta}\bm{u}^{\zeta}(\bm{X},n-\frac{1}{2})-\frac{\Delta
t}{\rho_{0}}\
\frac{\partial}{\partial\zeta}p(\bm{X},n)\Big{]}$$\triangleright$ Update
$\bm{u}$

6:$\bm{\rho}^{\zeta}(\bm{X},n+1)\
\leftarrow\Lambda^{\zeta}\Big{[}\Lambda^{\zeta}\bm{\rho}^{\zeta}(\bm{X},n)-\Delta
t\
\rho_{0}\frac{\partial}{\partial\zeta}\bm{u}_{\zeta}(\bm{X},n+\frac{1}{2})\Big{]}+\Delta
t\ s_{M}^{\zeta}(\bm{X},n+\frac{1}{2})$$\triangleright$ Update $\bm{\rho}$

7:$p(\bm{X},n+1)\ \ \ \leftarrow c^{2}\
\sum_{\zeta=1}^{d}\bm{\rho}^{\zeta}(\bm{X},n)$$\triangleright$ Update $p$

8:Terminate temporal iterations at $n=N_{T}-1$

9:Record $p(\bm{X},N_{T})$ on receivers

10:output: $p(\bm{X},n),n\in\\{1,...,N_{T}\\}$ on receivers

In Algorithm 1, the spatial derivatives are computed using a k-space pseudo-
spectral method [2, 45], and $s_{M}^{\zeta}(\bm{X},n+1/2)$, which is a
discretisation of the mass source shifted by $\Delta t/2$, can be computed by
enforcing a k-space correction step [45]. Also, $\Lambda^{\zeta}$ is an
operator representing the perfectly matched layers (PMLs) along the Cartesian
coordinate $\zeta$.

######  Remark 1.

In the k-Wave solver, the relation between the discretised directional
components of the mass source and the discretised excitation pulse
$s(\bm{X},n)$ (source.p) is defined by $\displaystyle
s_{M}^{\zeta}(\bm{X},n)=\frac{s(\bm{X},n)}{c^{2}d}\frac{2c}{\Delta\bm{x}}.$
(7) (cf. Eq. (2.19) in the k-Wave user manual [45], or Algorithm 1, line 5,
the second term in the right-hand-side in [23].) However, Eq. (7), in which
$s$ has been assumed pressure in
$\text{kg}\text{m}^{2-d}\text{s}^{-2}(\text{Pa})$, has not been understood by
the author. (cf. Table 3.3 in [45].)

###  2.3. Analytic solution to the wave equation

This section describes an analytic solution to the system of wave equations
(2) using the Green’s function for homogeneous media. This approach will be
used as the benchmark for measuring accuracy of the k-Wave solver for
homogeneous media.

A frequency-domain version of (2) at a single angular frequency $\omega$ can
be described using a Helmholtz equation of the form

$\displaystyle\left((\frac{\omega}{c})^{2}+\nabla^{2}\right)p(\omega,\bm{x})=-s(\omega,\bm{x}).$
(8)

The analytic solution to Eq. (8) in terms of the Green’s function can be
written in the form

$\displaystyle p(\omega,\bm{x};\bm{x}^{\prime})=\int
g(\omega,\bm{x};\bm{x}^{\prime})\ s(\omega;\bm{x}^{\prime})\
d\bm{x}^{\prime},$ (9)

where $g(\omega,\bm{x};\bm{x}^{\prime})$ is the Green’s function at point
$\bm{x}$ with respect to point $\bm{x}^{\prime}$ at the single angular
frequency $\omega$. Eq. (9) is used for an emission element with finite
volume. As $s$ is a point source here, this equation becomes a scaling

$\displaystyle
p(\omega,\bm{x};\bm{x}_{e})=g(\omega,\bm{x};\bm{x}_{e})s(\omega,\bm{x}_{e}).$
(10)

For a 2D homogeneous medium, the Green’s function is in the form [8]

$\displaystyle
g_{0,2D}(\omega,\bm{x};\bm{x}^{\prime})\approx\frac{1}{\left(8\pi\phi_{0}(\bm{x};\bm{x}^{\prime})\right)^{1/2}}e^{{i\mkern
1.0mu}(\phi_{0}(\bm{x};\bm{x}^{\prime})+\pi/4)}.$ (11)

For a 3D homogeneous medium, the Green’s function is in the form [8]

$\displaystyle
g_{0,3D}(\omega,\bm{x};\bm{x}^{\prime})=\frac{1}{4\pi|\bm{x}-\bm{x}^{\prime}|}e^{{i\mkern
1.0mu}\phi_{0}(\bm{x};\bm{x}^{\prime})}.$ (12)

Here, $\phi_{0}(\bm{x};\bm{x}^{\prime})$ denotes the phase of the Green’s
function, and satisfies for a homogeneous medium

$\displaystyle\phi_{0}(\bm{x};\bm{x}^{\prime})=k|\bm{x}-\bm{x}^{\prime}|,$
(13)

where $k=\omega/c$, and is the wavenumber.

##  3\. Numerical results

This section will describe numerical experiments demonstrating the accuracy of
the k-Wave and its corrected version, in which Eq. (7) has been replaced by
Eqs. (5) and (6), as discussed in section 2.2. It must be reminded that the
latter equation enforces a dimension-dependent scaling on the source to
account for spatial sampling. The k-Wave (version 1.4) was used in this study.
For measuring accuracy of the k-Wave, an analytic solution to the wave
equation using Eq. (9) was used as the benchmark.

###  3.1. 2D case

The k-Wave simulations were performed on a grid of size $0.4$ mm. One emitter
and 256 receivers were placed on a ring with a radius of 9.5 cm. The sound
speed inside the ring was set 1500 $\text{ms}^{-1}$ and homogeneous. The
emitter was placed on the position $\bm{x}_{e}=\big{[}0,-9.5\big{]}\
\text{cm}$. Figure 1 shows the emitter (red circle) and receivers (blue
circles). The interpolations between the position of the transducers and the
grid points were performed using the k-Wave toolbox [6, 46].

Figure 1. Transducers: (a) 2D case, (b) 3D case.

The k-Wave simulations were performed using three different excitation pulses
for a time $T=163.26\mu s$. Figures 2, 2 and 2 show the excitation pulses 1, 2
and 3, respectively. Here, a low-pass filter has been applied on each
excitation pulse using the filterTimeSeries.m function in the k-Wave to ensure
that the frequency spectrum of the pulse be smaller than the maximum frequency
supported by the computational grid, and the maximum absolute amplitude of the
excitation pulses has then been normalised. The maximum frequency supported by
the grid for the k-Wave simulation can be determined by the Shannon-Nyquist
limit [2]. For a homogeneous medium with sound speed $c$, the maximum
supported frequency is [45]

$\displaystyle f_{\text{max}}=\frac{c}{2\Delta x}.$ (14)

Here, $\Delta x$ is the grid spacing, and $c$ is the sound speed in water
($1500\text{m}{s}^{-1}$), but here, the low-pass filter was chosen more
conservatively using a lower sound speed ($1470\text{m}{s}^{-1}$) in (14) for
computing $f_{\text{max}}$. The computed $f_{\text{max}}$ was $1.84$ MHz.

The simulation 1 was performed using a CFL number 0.1, which provides a time
spacing $\Delta t=25.25$ $ns$, which is equivalent to a sampling rate of
$39.6$ MHz. The simulation 2 was performed using a CFL number 0.2, which gives
a time spacing $50.51$ ns, which is equivalent to a sampling rate of 19.80
MHz. The simulation 3 was performed using a CFL number 0.3, which is
equivalent to a time spacing $75.76$ ns and a sampling rate of 13.20 MHz.
Figures 2, 2 and 2 show the amplitude and phase of the excitation pulse in the
frequency domain for the simulations 1, 2 and 3, respectively.

For each simulation using the k-Wave and modified k-Wave, the pressure field
was computed in time and was recorded on the receivers at the same discretised
times as used for the k-Wave simulation. The recorded signal on each receiver
was then decomposed into the frequency components (phase and amplitude) using
a Fourier Transform operator. As a benchmark, the pressure field was also
calculated analytically in the frequency domain using Eq. (10), in which the
Green’s function is approximated using (11).

Figures 3, 3 and 3 show the amplitude of the pressure field induced by an
excitation of the emitter (the red circle in figure 1) and recorded on the
receiver 100 (the green circle in figure 1) as a function of frequency, when
the excitation pulses 1, 2 and 3 were used, respectively. As shown in these
figures, the amplitudes of the pressure field approximated using the modified
k-Wave match the amplitudes analytically calculated using the Green’s function
for all frequencies, but the amplitudes approximated using the k-Wave are
different from the analytically calculated amplitudes. In the same way,
figures 3, 3 and 3 show the phases of the pressure field recorded on the
receiver 100, when the excitation pulses 1, 2 and 3 were used, respectively.
The phases were wrapped to $\big{[}-\pi,\pi\big{]}$ rad. As shown in the these
figures, the phases approximated by the modified k-Wave and phases
analytically calculated using the Green’s function match for all frequencies
below 1.6 MHz, but the phases approximated by the k-Wave are different from
the analytically calculated phases. Using the excitation pulses 2 and 3, a
mismatch between the modified k-Wave and green’s function was found for
frequencies above 1.6 MHz, as shown in figures 3 and 3, respectively. This
mismatch is because the amplitude of the excitation pulses 2 and 3 at these
frequencies is very small, and therefore, simulation of propagation of
acoustic waves at these frequencies has been dominated by numerical errors.

Figure 2. 2D case. The excitation pulse in the time domain: a) 1, b) 2 and c)
3. The amplitude and phase of the excitation pulse in the frequency domain: d)
1, e) 2 and f) 3.

Figure 3. 2D case. The pressure field induced by the emitter and recorded on
receiver 100 in the frequency domain. Amplitude for simulations: a) 1, b) 2
and c) 3. Phase for simulations: d) 1, e) 2 and f) 3. The phases were wrapped
to $\big{[}-\pi,\pi\big{]}$ rad. The amplitude of the excitation pulses in (b)
and (c) is very small at frequencies above 1.6 MHz, and therefore, the
approximation of phases by k-Wave has been dominated by numerical errors.

Figures 4, 4 and 4 show the amplitudes of the pressure field recorded on all
the receivers as a function of distance of receivers to the emitter, when the
excitation pulses 1, 2 and 3 were used, respectively. In these figures, the
amplitudes were computed at the single frequency $f_{\text{max}}/2=0.92$ MHz.
As shown in these figures, for all the excitation pulses, the amplitudes
computed using the modified k-Wave simulation (red plot) match the analytic
solution to the wave equation using the Green’s function (green plot), but the
amplitudes computed by the k-Wave (blue plot) are different from the
analytically calculated amplitudes by almost the same factor.

In the same way, figures 4, 4 and 4 show the phase of the pressure field
recorded on all the receivers, when the excitation pulses 1, 2 and 3 were
used, respectively. The phases were wrapped to $\big{[}-\pi,\pi\big{]}$ rad.
As shown in these figures, the phases approximated by the modified k-Wave
match the analytically calculated phases on all the receivers, but the phases
approximated using the k-Ware are different from the analytically calculated
phases .

Figure 4. 2D case. The pressure field at the single frequency
$f_{\text{max}}/2$ on the receivers as a function of distance of the receivers
to the emitter. Amplitude for simulations: a) 1, b) 2 and c) 3. Phase for
simulations: d) 1, e) 2 and f) 3. The phases were wrapped to
$\big{[}-\pi,\pi\big{]}$ rad.

###  3.2. 3D case

The k-Wave simulations were performed on a grid of size $0.5$ mm. One emitter
and 256 receivers were placed on a hemisphere with a radius of 12.35 cm, as
shown in figure 1. The emitter was placed at the bottom of the hemisphere at
position $\bm{x}_{e}=\big{[}1.54,0,-12.25\big{]}\text{cm}$. Figure 1 shows the
emitter by the red colour and the receivers by the blue colour. The
interpolations between the position of the transducers and the grid points
were performed using the k-Wave toolbox [6, 46]. The sound speed inside the
hemisphere was set 1500 $\text{ms}^{-1}$ and homogeneous. The k-Wave
simulations were performed using three excitation pulses for a time
$T=212.2\mu s$. Figure 5, 5 and 5 show the excitation pulses 1, 2 and 3 in the
time domain, respectively. Here, a low-pass filter has been applied on each
excitation pulse using the filterTimeSeries.m function in the k-Wave to ensure
that the frequency spectrum of the pulse be smaller than the maximum frequency
supported by the computational grid, and then the maximum absolute amplitude
of the excitation pulses was normalised. The low pass filter was chosen
conservatively using a sound speed of 1470 $\text{ms}^{-1}$ for computing
$f_{\text{max}}$ in (14). The computed $f_{\text{max}}$ was 1.47 MHz.

For the simulation using the excitation pulse 1, the CFL number was chosen
0.1, which is equivalent to a time spacing $\Delta t=31.57$ ns and a sampling
rate of 31.68 MHz. The simulation 2 was performed using the CFL number 0.2,
which provides a time spacing 63.13 MHz, which is equivalent to a sampling
rate 15.84 MHz. The simulation using the excitation pulse 3 was performed
using the CFL number 0.3, which is equivalent to a time spacing $94.70$ ns and
a sampling rate of 10.56 MHz. Figures 5, 5 and 5 show the amplitude and phase
of the excitation pulse in the frequency domain for the simulations 1, 2 and
3, respectively.

Figure 5. 3D case. The excitation pulse in the time domain: a) 1, b) 2 and c)
3. The amplitude and phase of the excitation pulse in the frequency domain: d)
1, e) 2 and f) 3.

For each simulation using the k-Wave and modified k-Wave, the pressure field
was computed in time and was recorded on the receivers at the same discretised
times as used for the k-Wave simulation. The recorded signal on each receiver
was then decomposed into the frequency components (phase and amplitude) using
a Fourier Transform operator. The pressure field was also calculated
analytically in the frequency domain using Eq. (10), in which the Green’s
function is approximated using (12).

Figures 6, 6 and 6 show the amplitude of the pressure field recorded on
receiver 100 (the green circle in figure 1) after being produced by the
emitter (the red circle in figure 1) in the frequency domain for the
simulation using the excitation pulse 1, 2 and 3, respectively. Because the
amplitudes computed by the k-Wave (the blue plot) are much smaller than the
amplitudes calculated using the Green’s function (green plot), a variant of
the amplitudes which are enlarged by a factor 666.7 was also shown by the
magenta colour for a better visualisation. As shown in these figures, the
amplitudes approximated using the modified k-Wave match the amplitudes
calculated analytically using the Green’s function for all frequencies. In
addition, considering Eq. (12), the amplitude of the Green’s function in the
3D case depends only on distance, and is not changed with frequency. Note that
for a nonabsorbing medium, the amplitude of the Green’s function represents
only geometrical attenuation. (cf. [7] and [8].) Therefore, the fact that the
amplitude of the Green’s function in (12) is independent from frequency
implies that the shape of the ultrasonic pulse should remain unchanged after
propagation in the medium. A comparison between the amplitudes of the
excitation pulse in figures 5, 5 and 5 and the approximated amplitudes of the
pressure field on the receiver 100 in the respective figures 6, 6 and 6 show
that the geometrical attenuation approximated by the k-Wave has changed with
frequency, but the attenuation approximated using the modified k-Wave has
remained unchanged with frequency. For example, for simulation 1, the maximum
amplitude of the excitation pulse in the frequency domain, the amplitude
analytically calculated using the Green’s function, and the amplitude
approximated using the modified k-Wave in the frequency domain are all close
to 0.6 MHz, but the maximum amplitude approximated using the k-Wave has been
shifted to about 0.9 MHz, as shown in figure6 by the magenta plot.

Figures 6, 6 and 6 show the phases of the pressure field recorded on the
receiver 100, when the excitation pulses 1, 2 and 3 were used, respectively.
The phases were wrapped to $\big{[}-\pi,\pi\big{]}$ rad. As shown in these
figures, the phases approximated by the modified k-Wave and the phases
calculated using the Green’s function match for all frequencies, but the
phases approximated using the k-Wave are different.

Figure 6. 3D case. The pressure field induced by the emitter and recorded on
receiver 100 in the frequency domain. Amplitude for simulations: a) 1, b) 2
and c) 3. Phase for simulations: d) 1, e) 2 and f) 3. The phases were wrapped
to $\big{[}-\pi,\pi\big{]}$ rad.

Figures 7, 7 and 7 show the amplitudes of the pressure field recorded on all
the receivers as a function of distances of receivers to the emitter, when the
excitation pulses 1, 2 and 3 were used, respectively. In these figures, the
amplitudes were computed at the single frequency $f_{\text{max}}/2=0.73$ MHz.
Because the amplitudes computed by the k-Wave (blue plot) are much smaller
than the Green’s function (green plot), they are also enlarged by a factor
666.7 for a better visualisation (magenta plot). As shown in these figures,
using all the chosen excitation pluses, the amplitudes computed using the
modified k-Wave (red plot) match the analytically calculated amplitudes (green
plot), but the amplitudes computed by the k-Wave are different from the
Green’s function by almost the same factor.

In the same way, figures 7, 7 and 7 show the phases, which have been wrapped
to $\big{[}-\pi,\pi\big{]}$ rad. As shown in these figures, the phases
approximated by the modified k-Wave match the analytically calculated phases
for all the receivers, but the phases approximated using the k-Wave are
different.

Figure 7. 3D case. The pressure field at the single frequency
$f_{\text{max}}/2$ on the receivers as a function of distance of the receivers
to the emitter. Amplitude for simulations: a) 1, b) 2 and c) 3. Phase for
simulations: d) 1, e) 2 and f) 3. The phases were wrapped to
$\big{[}-\pi,\pi\big{]}$ rad.

##  4\. Discussion

This manuscript showed numerically that the k-space pseudo-spectral method
proposed in [1, 2] for solving the three-coupled first order wave equations is
almost exact for homogeneous media. A Matlab toolbox for a numerical
implementation of this approach is available and was used in this study [5, 6,
45], but the approach taken for including source in the system of differential
equations was revisited and modified. For measuring accuracy of this wave
solver, an analytic solution to the second-order wave equation using the
Green’s function was used as the benchmark. After applying the corrections
discussed in section 2.2 for including the source in the equation of
continuity (conservation of mass) in the system of three-coupled linearised
wave equations, a good match between the wave solver and analytic solution to
the wave equation was obtained.

In this study, the transducers were assumed points which are positioned off
the grid and on the surface of a ring or hemisphere. The match between Eq. (9)
(the Green’s relation) and Algorithm 1 (the wave solver) for a point source
implies that an approximation of the pressure field in 2D (resp. 3D) media
using real transducers and with finite size requires an integration of point
sources over area (resp. volume) of the transducers.

## References

  * [1] T. D. Mast, L. P. Souriau, D. -L. D. Liu, M. Tabei, A. I. Nachman and R. C. Waag, A k-space method for large-scale models of wave propagation in tissue, IEEE Trans. Ultrason. Ferroelectr. Freq., Vol. 48, no. 2, pp. 341-354, March 2001, doi: 10.1109/58.911717.
  * [2] M. Tabei, T. D. Mast, and R. C. Waag, “A k-space method for coupled first-order acoustic propagation equations,” J. Acoust. Soc. Am. vol. 111, pp. 53–63, 2002.
  * [3] S. Finette, Computational methods for simulating ultrasound scattering in soft tissue, IEEE Trans. Ultrason., Ferroelect. Freq. Contr., vol. 34, no. 3, pp. 283–292, 1987.
  * [4] M. J. Bencomo and W. W. Symes, Discretization of multipole sources in a finite difference setting for wave propagation problems, Journal of Computational Physics, vol. 386, pp. 296-322, 2019.
  * [5] B. E. Treeby and B. T. Cox, k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields J. Biomed. Opt. vol. 15, no. 2, 021314, 2010.
  * [6] www.k-Wave.org. (version 1.4)
  * [7] A. Javaherian and B. Cox, Ray-based inversion accounting for scattering for biomedical ultrasound tomography, Inverse Problems, vol. 37, no.11, 115003, 2021.
  * [8] A. Javaherian, Hessian-inversion-free ray-born inversion for high-resolution quantitative ultrasound tomography, arXiv:2211.00316. https://arxiv.org/abs/2211.00316.
  * [9] J. Qian, P. Stefanov, G. Uhlmann and H. Zhao, An Efficient Neumann series–based algorithm for thermoacoustic and photoacoustic tomography with variable sound speed, SIAM Journal on Imaging Sciences, Vol. 4, no. 3, 2011. doi: 10.1137/100817280.
  * [10] R. Kowar and O. Scherzer. Attenuation Models in Photoacoustics. In: Ammari, H. (eds) Mathematical Modeling in Biomedical Imaging II. Lecture Notes in Mathematics, Vol. 2035, 2012. Springer, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-22990-94.
  * [11] X. L. Deán-Ben, A. Buehler, V. Ntziachristos and D. Razansky, Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography, IEEE transactions on medical imaging, Vol. 31, no. 10, pp. 1922-1928, 2012.
  * [12] S. Arridge, M. Betcke, B. Cox, F. Lucka and B E Treeby, On the adjoint operator in photoacoustic tomography, Inverse Problems, vol. 32, 115012, 2016.
  * [13] M. Haltmeier and L. V. Nguyen, Analysis of iterative methods in photoacoustic tomography with variable sound speed, SIAM Journal on Imaging Sciences, Vol. 10, no. 2, 2017. doi: 10.1137/16M1104822.
  * [14] M. Mozaffarzadeh, A. Mahloojifar, M. Orooji, S. Adabi and M. Nasiriavanaki, Double-stage delay multiply and sum beamforming algorithm: Application to linear-array photoacoustic imaging, IEEE Transactions on Biomedical Engineering, Vol. 65, no. 1, pp. 31-42, 2018.
  * [15] A. Javaherian and S. Holman, A continuous adjoint for photo-acoustic tomography of the brain, Inverse Problems, vol. 34, no. 8, 085003, 2018.
  * [16] S. Antholzer, M. Haltmeier and J. Schwab, Deep learning for photoacoustic tomography from sparse data, Inverse Problems in Science and Engineering, Vol. 27, no. 7, pp. 987-1005, 2019. doi: 10.1080/17415977.2018.1518444.
  * [17] A. Javaherian and S. Holman, Direct quantitative photoacoustic tomography for realistic acoustic media, Inverse Problems, vol. 35, no. 8, 084004, 2019.
  * [18] S. Guan, A. A. Khan, S. Sikdar and P. V. Chitnis, Fully dense UNet for 2-D sparse photoacoustic tomography artifact removal, IEEE Journal of Biomedical and Health Informatics, vol. 24, no. 2, pp. 568-576, Feb. 2020, doi: 10.1109/JBHI.2019.2912935.
  * [19] S. Na and L. V. Wang, Photoacoustic computed tomography for functional human brain imaging, Biomed. Opt. Express Vol. 12, pp. 4056-4083, 2021.
  * [20] J. Grohl, M. Schellenberg, K. Dreher and L. Maier-Hein, Deep learning for biomedical photoacoustic imaging: A review, Photoacoustics, Vol. 22, 100241, 2021.
  * [21] E. Park, S. Park, H. Lee, M. Kang, C. Kim and J. Kim, Simultaneous dual-modal multispectral photoacoustic and ultrasound macroscopy for three-dimensional whole-body imaging of small animals, Photonics, Vol. 8, pp. 13, 2021.
  * [22] J. A. Guggenheim, J. Li, T. J. Allen, R. J. Colchester, S. Noimark, O. Ogunlade, I. P. Parkin, I. Papakonstantinou, A. E. Desjardins, E. Z. Zhang and P. C. Beard, Ultrasensitive plano-concave optical microresonators for ultrasound sensing. Nature Photon, Vol. 11, pp. 714–719, 2017. https://doi.org/10.1038/s41566-017-0027-x.
  * [23] F. Lucka, M. Pérez-Liva, B. E. Treeby and B. T. Cox, High resolution 3D ultrasonic breast imaging by time-domain full waveform inversion, Inverse Problems, Vol. 38, 025008, 2022.
  * [24] J. Wiskin, D. T. Borup, S. A. Johnson and M. Berggren, Non-linear inverse scattering: high resolution quantitative breast tissue tomography, J. Acoust. Soc. Am., vol. 131, no. 5, pp. 3802-13, 2012.
  * [25] K. Wang, T. Matthews, F. Anis, C. Li, N. Duric, and M. A. Anastasio, Waveform inversion with source encoding for breast sound speed reconstruction in ultrasound Computed Tomography, IEEE T ULTRASON FERR, vol. 62, no. 3, 2015.
  * [26] J. W. Wiskin, D. T. Borup, E. Iuanow, J. Klock, M. W. Lenox, 3-D Nonlinear acoustic inverse scattering: algorithm and quantitative results, IEEE T ULTRASON FERR, Vol. 64, no. 3, 2017.
  * [27] A. V. Goncharsky and S. Y. Romanov, Iterative methods for solving coefficient inverse problems of wave tomography in models with attenuation, Inverse Problems, Vol. 33, pp. 025003, 2017.
  * [28] T. P. Matthews, K. Wang, C. Li, N. Duric, and M. A. Anastasio, Regularized dual averaging image reconstruction for full-wave ultrasound computed tomography, IEEE T ULTRASON FERR, vol. 64, no. 5, 2017.
  * [29] T. P. Matthews and M. A. Anastasio, Joint reconstruction of the initial pressure and speed of sound distributions from combined photoacoustic and ultrasound tomography measurements, Inverse Problems, vol. 33, pp. 124002, 2017.
  * [30] M. Pérez-Liva, J. L. Herraiz, J. M. Udias, E. Miller, B. T. Cox, and B. E. Treeby, Time domain reconstruction of sound speed and attenuation in ultrasound computed tomography using full wave inversion, J. Acoust. Soc. Am., Vol. 141, no. 3, pp. 1595-1604, 2017. doi: 10.1121/1.4976688.
  * [31] E. Merćep, J. L. Herraiz, X. L. Deán-Ben and D. Razansky, Transmission-reflection optoacoustic ultrasound (TROPUS) computed tomography of small animals. Light Sci Appl. Vol 8: 18. 2019, doi: 10.1038/s41377-019-0130-5.
  * [32] A. Siahkoohi, M. Louboutin, and F. J. Herrmann, The importance of transfer learning in seismic modeling and imaging, Geophysics, Vol. 84, no. 6., pp. A47–A52, 2019. doi:10.1190/geo2019-0056.1.
  * [33] H. Aghamiry, A. Gholami, and S. Operto, 2019, Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers, Geophysics, Vol. 84, no. 1, R139–R162, 2019.
  * [34] L. Guasch, O. Calderón Agudo, M. Tang, P. Nachev, and M. Warner, Full-waveform inversion imaging of the human brain. Nature Digital Medicine, vol. 3, 28, 2020.
  * [35] M. van Driel , C. Boehm, L. Krischer and M. Afanasiev, Accelerating numerical wave propagation using wavefield adapted meshes. Part I: forward and adjoint modelling, Geophys. J. Int., Vol. 221, pp. 1580–1590, 2020.
  * [36] E. Bachmann and J. Tromp, Source encoding for viscoacoustic ultrasound computed tomography J. Acoust. Soc. Am., vol. 147, 3221–35, 2020.
  * [37] M. Feigin , D. Freedman, and B. W. Anthony, A Deep learning framework for single-sided sound speed inversion in medical ultrasound, IEEE Transactions on Biomedical Engineering, Vol. 67, no. 4, pp. 1142-1151, 2020.
  * [38] W. W. Symes, Error bounds for extended source inversion applied to an acoustic transmission inverse problem, Inverse Problems, Vol. 38, no. 11., 2022, 115002 (21pp), 2022.
  * [39] D. Perdios, A. Besson, M. Arditi and J. -P. Thiran, A deep learning approach to ultrasound image recovery, 2017 IEEE International Ultrasonics Symposium (IUS), Washington, DC, USA, 2017, pp. 1-4, doi: 10.1109/ULTSYM.2017.8092746.
  * [40] V. H. J. van Hal, H. De Hoop, J. -W. Muller, M. R. H. M. van Sambeek, H. -M. Schwab and R. G. P. Lopata, Multiperspective Bistatic Ultrasound Imaging and Elastography of the Ex Vivo Abdominal Aorta, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, Vol. 69, no. 2, pp. 604-616, 2022, doi: 10.1109/TUFFC.2021.3128227.
  * [41] J. K. Mueller, L. Ai, P. Bansal and W. Legon, Numerical evaluation of the skull for human neuromodulation with transcranial focused ultrasound, J. Neural Eng., Vol. 14, 066012 (19pp), 2017.
  * [42] S. Jiménez-Gambin, N. Jiménez, J. M. Benlloch and F. Camarena, Holograms to focus arbitrary ultrasonic fields through the skull, Physical Review Applied, Vol. 12, 014016, 2019.
  * [43] A. N. Pouliopoulos, S-Y Wu, M. T. Burgess, M. E. Karakatsani, H. A. S. Kamimura, E. E. Konofagou, A clinical system for non-invasive blood–brain barrier opening using a neuronavigation-guided single-element focused ultrasound transducer, Ultrasound in Medicine & Biology, Vol. 46, no. 1, pp. 73-89, 2020.
  * [44] J. A. Cain, S. Visagan, M.A., Johnson, J. Crone, R. Blades, N. M. Spivak, D. W. Shattuck and M. M. Monti, Real-time and delayed effects of subcortical low intensity focused ultrasound. Sci Rep, Vol. 11, 6100, 2021. https://doi.org/10.1038/s41598-021-85504-y.
  * [45] k-Wave user manual, A Matlab toolbox for the time domain simulation of acoustic wave fields, Version 1.1, 27th August 2016 (the last version).
  * [46] E. S. Wise, B. T. Cox, J. Jaros, B. E. Treeby, Representing arbitrary acoustic source and sensor distributions in Fourier collocation methods, J. Acoust. Soc. of Am., vol. 146, no. 1, pp. 278-288, 2019.


# Non-equilibrium steady states of long-range coupled harmonic chains

Francesco Andreucci  SISSA and INFN, Sezione di Trieste, Via Bonomea 265,
I-34136 Trieste, Italy  Stefano Lepri  Consiglio Nazionale delle Ricerche,
Istituto dei Sistemi Complessi, Via Madonna del Piano 10 I-50019 Sesto
Fiorentino, Italy.  Istituto Nazionale di Fisica Nucleare, Sezione di Firenze,
via G. Sansone 1, I-50019 Sesto Fiorentino, Italy  Stefano Ruffo  SISSA and
INFN, Sezione di Trieste, Via Bonomea 265, I-34136 Trieste, Italy  Consiglio
Nazionale delle Ricerche, Istituto dei Sistemi Complessi, Via Madonna del
Piano 10 I-50019 Sesto Fiorentino, Italy.  Andrea Trombettoni  Department of
Physics, University of Trieste, Strada Costiera 11, I-34151 Trieste, Italy
SISSA and INFN, Sezione di Trieste, Via Bonomea 265, I-34136 Trieste, Italy
CNR-IOM DEMOCRITOS Simulation Center, Via Bonomea 265, I-34136 Trieste, Italy

###### Abstract

We perform a numerical study of transport properties of a one-dimensional
chain with couplings decaying as an inverse power $r^{-(1+\sigma)}$ of the
inter-site distance $r$ and open boundary conditions, interacting with tho
heat reservoirs. Despite its simplicity, the model displays highly nontrivial
features in the strong long-range regime, $-1<\sigma<0$. At weak coupling with
the reservoirs, the energy flux departs from the predictions of perturbative
theory and displays anomalous superdiffusive scaling of the heat current with
the chain size. We trace back this behavior to the transmission spectrum of
the chain, which displays a self-similar structure with a characteristic
sigma-dependent fractal dimension.

##  I Introduction

The main task of statistical mechanics is to relate the microscopic
interactions of a given system to its macroscopic properties. One typical
instance is the context of heat transfer. Suppose we apply a temperature
gradient $\nabla T$ to a system, after a while the system will reach a
stationary state characterized by the presence of a heat flux $\mathcal{J}$.
The thermal conductivity $\kappa$ is defined in terms of these quantities as:

$\mathcal{J}=-\kappa\nabla T,$ (1)

In the case of diffusive transport, Fourier’s law holds and $\kappa$ does not
depend on the size of the system $N$ in the thermodynamic limit. This is
typically the case for three-dimensional systems with short-range
interactions. We remark, however, that there is currently no generic way,
given the microscopic properties of a system, to know whether Fourier’s law
holds or not.

A case in which Fourier’s law is systematically violated is the case of
harmonic interactions. For instance, for the harmonic crystal each phonon
propagates freely and the transport is ballistic. This was showed for the
first time for a chain with nearest-neighboors interactions in the seminal
paper by Rieder, Lebowitz and Lieb Rieder  _et al._ (1967). They found that
the thermal conductivity $\kappa$ diverges as $\kappa\propto N$, $N$ being the
number of particles in the chain. Moreover, the bulk temperature profile is
flat, while Fourier’s law would lead to a linear one. The non-equilibrium
properties of quantum harmonic lattices have also been considered in the last
decades Zurcher and Talkner (1990); Dhar and Roy (2006); Saito  _et al._
(2000); Asadian  _et al._ (2013); Freitas and Paz (2014b).

Generally speaking, in harmonic lattices transport features are dictated by
the spectral properties of both the thermal reservoirs and the system itself.
For instance in the case of disordered lattices displaying Anderson
localization, the conductivity (or energy flux) depends on the localization
lengths, but also on the boundary conditions Visscher (1971), the spectral
density of the baths at low frequencies Dhar (2001) as well as on the
distribution and correlations of the random disorder Herrera-González  _et
al._ (2015); Ash  _et al._ (2020). For more general, non homogeneous harmonic
networks, the spectral properties can be accounted by random matrix theory and
can describe also current fluctuations Schmidt  _et al._ (2013). This is even
more striking for active (non-equilibrium) baths that can lead to non-trivial
transport regimes even for the ordered harmonic chain Santra and Basu (2022).

It became progressively become clear that in one (and two) dimensions there
are violations of Fourier’s law also for nonlinear systems Lepri  _et al._
(2003); Dhar (2008); Lepri  _et al._ (2016); Benenti  _et al._ (2020, 2023),
such as the Fermi-Pasta-Ulam-Tsingou (FPUT) chain. In one dimension, these
violations manifest themselves as a power-law divergence of the thermal
conductivity $\kappa$ with the system’s size $\kappa\propto N^{\alpha}$.
Transport in these cases is called anomalous. It is now clear that
superdiffusive transport is a generic feature of non-linear one-(and
two)-dimensional non-integrable systems conserving momentum, energy and
stretch. There are both numerical and analytical evidences that the exponent
$\alpha$ can be used to identify different universality classes Benenti  _et
al._ (2020). For weakly non-integrable models the scenario may be more
involved since quasi-particles may have very large mean-free paths Di Cintio
_et al._ (2018); Lepri  _et al._ (2020).

A further element of interest is represented by the presence of forces that
are not strictly local. Indeed, much less is known about systems with long-
range interactions, that is, systems in which the inter-particle interaction
scales with the particle distance $r$ as $V(r)\propto r^{-d-\sigma}$. Several
physical systems are characterized by long-range interaction, both classical
(gravity, pure plasmas, $2d$ hydrodynamics) and quantum (dipolar systems and
trapped atoms). As a concrete experimental instance we mention trapped ion
chains, where ions are confined in periodic arrays and interact with external
reservoirs Bermúdez  _et al._ (2013); Ramm  _et al._ (2014). On a macroscale,
effective long-range forces arise for tailored macroscopic systems like chain
of coupled magnets Molerón  _et al._ (2019) and the effects of fluctuations
and nonlinearity may be relevant.

Long-range systems received considerable attention in the last years, for
reviews see for example Campa  _et al._ (2009) and Defenu  _et al._ (2021) for
classical and quantum systems respectively. For what we are going to be
concerned with in this paper, we remind that, at equilibrium, the universality
class of a one-dimensional long-range system depends on the value of $\sigma$.
Indeed, for $-1<\sigma<0$, the critical exponents are the mean field ones,
that is the ones that we obtain by putting $\sigma=-1$. Then, there exists a
non universal value $\sigma^{*}$ such that for $\sigma>\sigma^{*}$ we recover
the critical exponents of the short-range case $\sigma=\infty$. Typically
$\sigma^{*}>0$. Furthermore, excitations in long-range systems can propagate
at diverging velocity Torcini and Lepri (1997); Métivier  _et al._ (2014) and
therefore we can expect some form of superdiffusive transport. There are
already several, mainly numerical, studies of heat transport in long-range
interacting systems that confirm these expectations. On the classical side,
the heat transport was analyzed for the long-range XY model Olivares and
Anteneodo (2016); Iubini  _et al._ (2018) the FPUT chain in Bagchi (2017);
Iubini  _et al._ (2018); Di Cintio  _et al._ (2019); Wang  _et al._ (2020);
Bagchi (2021) and the lattice $\varphi^{4}$ theory Iubini  _et al._ (2022). In
all cases Fourier’s law is violated in different ways according to the value
of $\sigma$. Scaling analysis of equilibrium correlations also suggests that
hydrodynamics is non-standard Di Cintio  _et al._ (2019); Iubini  _et al._
(2022). Thus, one may interpret transport as a fractional diffusion process
with energy carriers performing Lévy flights, with jump statistics controlled
by the exponent $\sigma$.

A classical harmonic long-range model with a stochastic dynamics was studied
analytically in Tamaki and Saito (2020); Suda (2022) and the heat flux and
temperature profile for a mean-field chain were computed in Defaveri  _et al._
(2022); Andreucci  _et al._ (2022). The same system was studied in the quantum
regime in Andreucci  _et al._ (2022) and a hydrodynamic approach to study
transport in quantum magnets was proposed in Schuckert  _et al._ (2020). We
refer again to Defenu  _et al._ (2021) for more references on the study of
dynamics and transport in quantum long-range systems. However, in the
literature there is not yet a detailed study of the plain harmonic chain with
power-law interaction, and this contribution aims at filling this gap. We will
show that the results are far from trivial in the strong long-range case and
deserve careful analysis.

More precisely, in this paper we study numerically heat transport in a
quadratic chain with a power-law interaction by coupling the first and last
site of the system to two heat baths at different temperature. We focus on
computing the heat flux in the stationary state with different approaches. In
section II we introduce the model and the main methods that we will use to
compute the heat flux. In section III-V we report an analysis based on the
spectral properties of the nonequilibrium Green’s function and the
transmission spectra and we discuss them. Finally, we draw our conclusions in
section VI.

##  II Model and methods

###  II.1 The long-range coupled harmonic chain

We consider a one-dimensional chain of particles with a power-law interaction:

$H=\frac{1}{2}\sum_{i}p_{i}^{2}+\frac{1}{2}\sum_{ij}x_{i}\Phi_{ij}x_{j},$ (2)

where the interaction matrix $\Phi$ is given by:

$\Phi_{ij}=\left(2\delta_{ij}-\frac{1}{N_{\sigma}}\frac{1}{|i-j|^{1+\sigma}}\right),\quad
N_{\sigma}=\sum_{l=1}^{N}l^{-\sigma},$ (3)

where $N_{\sigma}$ is the usual Kac factor introduced to guarantee extensivity
of the energy, chosen as site-independent. The matrix correctly reduces to the
discrete Laplacian for large $\sigma$. Note that definition (3) corresponds to
open boundary conditions, which are the ones appropriate for our problem due
to the presence of the baths. For long-ranged systems we expect that the role
of boundary conditions can have very important consequences, even more than
for short-ranged systems, and we focus on this natural choice for simplicity.

In the case of open boundary conditions the spectrum of matrix $\Phi$ is, to
the best of our knowledge, not known analytically. The usual standing waves
are not eigenvectors and the matrix cannot be diagonalized exactly. Even in
the continuum limit, this would correspond to solving the spectral problem for
the fractional Laplacian in a finite domain, which is notoriously not
straightforward Zoia  _et al._ (2007).

For comparison, it is useful to recall the solvable case for periodic boundary
condition where the proper definition of $\Phi$ is:

$\Phi_{ij}=\left(2\delta_{ij}-\frac{1}{N_{\sigma}}\frac{1}{min(N-|i-j|^{1+\sigma},|i-j|^{1+\sigma})}\right).$
(4)

Here the spectrum is known, see for example Defenu (2021). Due to
translational invariance, the eigenvectors are plane waves of wavenumber $k$.
The nature of the eigenfrequency spectra strongly depends on whether $\sigma$
is positive or negative. In the first case, the system has a proper continuum
limit and for low momenta $k$ the squared frequencies $\omega^{2}$ of the
plane waves behave as:

$\omega_{k}^{2}\approx\begin{cases}|k|^{\sigma},\quad&0<\sigma<2,\\\
k^{2},\quad&\sigma>2.\end{cases}$ (5)

Thus, for $\sigma>0$ one has the standard acoustic dispersion and a finite
group velocities while the group velocity diverges as
$|k|^{\frac{\sigma-2}{2}}$ in the first case. This result can also be derived
from the continuum limit, corresponding to a fractional wave equation in the
infinite domain Tarasov (2006). On the other hand, if $\sigma<0$ the spectrum
remains discrete even in the thermodynamic limit and contains a countable
infinite number of frequencies that accumulate at the band edge Defenu (2021).

To simulate the non-equilibrium steady state, we follow the usual procedure
and connect the first and last sites of the system to two Langevin heat baths
at temperatures $T_{L}$ and $T_{R}$, respectively. The coupling with the baths
introduces both noise and dissipation in the dynamics of the system. The
resulting equations of motion are:

$\ddot{x}_{i}=-\sum_{j}\Phi_{ij}x_{j}+\delta_{i1}\left(\xi_{L}-\lambda\dot{x}_{i}\right)+\delta_{iN}\left(\xi_{R}-\lambda\dot{x}_{i}\right),$
(6)

where the $\xi$’s are Gaussian noises that satisfy the fluctuation-dissipation
relation:

$\Braket{\xi_{a}(t)\xi_{a}(t^{\prime})}=2T_{a}\lambda\delta(t-t^{\prime}),\quad
a=L,R.$ (7)

After a transient, the system reaches a stationary state: we are interested in
the heat flux of the chain in this state. To compute this quantity, we will
employ three different methods.

###  II.2 RLL approach

The first method was introduced long time ago in this context in Rieder  _et
al._ (1967). It consists in solving the many-body Fokker-Planck equation
related to (6) (in the following we will refer to this method as the RLL
method). In particular, defining the vector
$y=(x_{1},...x_{N},p_{1},...p_{N})$, and denoting by $P(y,t)$ its probability
at time $t$, the aforementioned equation reads as:

$\frac{\partial P(y,t)}{\partial t}=A_{ij}\frac{\partial}{\partial
y_{i}}(y_{j}P)+\frac{1}{2}D_{ij}\frac{\partial^{2}P}{\partial y_{i}\partial
y_{j}},\quad$ (8)

where the drift and diffusion matrices are

$A=\begin{pmatrix}\mathbb{O}&-\mathbb{I}\\\
-\Phi&\lambda\mathcal{R}\end{pmatrix},\quad
D=\begin{pmatrix}\mathbb{O}&\mathbb{O}\\\ \mathbb{O}&2k_{B}\lambda
T(\mathcal{R}+\eta\mathcal{S})\end{pmatrix},$ (9)

where

$\displaystyle\begin{cases}T=\dfrac{T_{L}+T_{R}}{2},\\\
\eta=\dfrac{T_{L}-T_{R}}{T}\end{cases},\quad\mathcal{R}_{ij}=\delta_{ij}(\delta_{i1}+\delta_{iN}),\quad$
(10) $\displaystyle\mathcal{S}_{ij}=\delta_{ij}(\delta_{i1}-\delta_{iN}).$
(11)

The solution of equation (8) is a multi-variate Gaussian whose covariance
matrix is given by the matrix of correlations among the canonical coordinates:

$\displaystyle
P(y,t)\propto\exp{\left\\{-\frac{1}{2}C^{\,-1}_{ij}y_{i}y_{j}\right\\}},\quad
C=\begin{pmatrix}\Braket{x_{i}x_{j}}&\Braket{x_{i}p_{j}}\\\
\Braket{p_{i}x_{j}}&\Braket{p_{i}p_{j}}\end{pmatrix}.$ (12)

By plugging (12) in the Fokker-Planck equation (8) we get:

$\partial_{t}C=D-AC-CA^{T}.$ (13)

Furthermore, in the stationary state $\partial_{t}C=0$, so we get the so-
called (continuous) Lyapunov equation:

$AC+CA^{T}=D,$ (14)

which has to be solved numerically. Knowing the various correlators, we can
then express the heat flux in the stationary state as the difference between
the temperature of the left bath and the temperature of the first site:

$\mathcal{J}=\lambda\left(T_{L}-T_{1}\right),\quad
T_{i}=\frac{1}{2}\Braket{p_{i}^{2}}.$ (15)

###  II.3 Nonequilibrium Green’s function

The second method consists in writing the exact solution to (6) in terms of
the Green’s function $G(\omega)$, which is possible due to the linearity of
the equations. The details of this method are explained in refs. Dhar and Roy
(2006); Dhar and Saito (2016); Dhar (2008). Since we are interested in the
stationary state, we work directly in frequency space:

$\displaystyle\tilde{x}_{l}(\omega)=\sum_{ln}G_{ln}(\omega)(\tilde{\xi}_{L,n}(\omega)+\tilde{\xi}_{R,n}(\omega)),$
(16) $\displaystyle
G(\omega)=\left(-\omega^{2}\mathbb{I}+\Phi+i\lambda\omega\mathcal{R}\right)^{-1},$
(17)

where the tilde indicates the Fourier transform and $\mathcal{R}$ is the
matrix defined in Eqs.(9, 11). As explained in Dhar (2008), we can express the
heat flux in the stationary state as:

$\mathcal{J}=\frac{2\Delta
T\lambda^{2}}{\pi}\int_{0}^{\infty}d\omega\,\omega^{2}|G_{1N}(\omega)|^{2}.$
(18)

###  II.4 Generalized eigenvalue method

There is in the literature another approach to the Green’s function method,
called generalized eigenvalues method, which we briefly outline below (for a
more detailed explanation see Freitas and Paz (2014b); Freitas  _et al._
(2015); Tisseur and Meerbergen (2001)). Let $G^{L}(s)$ be the Green’s function
defined in Laplace’s space:

$G^{L}(s)=\left(-s^{2}\mathbb{I}+\Phi+\lambda sR\right)^{-1},$ (19)

and introduce the $2N$ complex numbers $\\{s_{a}\\}_{a=1}^{2N}$ and the $2N$
vectors $\\{\bm{r}_{a}\\}_{a=1}^{2N}$ as defined by the following linear
problem:

$G^{L}(s_{a})\bm{r}_{a}=0.$ (20)

Then, the Green’s function (19) can be written as Tisseur and Meerbergen
(2001):

$G^{L}(s)=\sum_{a=1}^{2N}\dfrac{s_{a}}{s-s_{a}}\bm{r}_{a}\bm{r}_{a}^{\dagger}.$
(21)

Note that the $s_{a}$ come in complex conjugate pairs. We now recall that we
can obtain the Green’s function in frequency space $G(\omega)$ via a Wick
rotation $G(\omega)=G^{L}(-is)$. Then we can compute the integral in (18) with
a contour integration one finding Freitas and Paz (2014b):

$\mathcal{J}=2\Delta
T\lambda^{2}\sum_{a,b=1}^{2N}\dfrac{s_{a}^{3}s_{b}}{s_{a}+s_{b}}r_{a,1}r_{a,N}^{\\*}r_{b,N}r_{b,1}^{\\*}.$
(22)

Formula (22) gives yet another way of computing the heat flux and extract the
scaling exponents.

###  II.5 Comments

Before proceeeding, let us comment on the numerical issues connected with the
above approaches. The numerical implementation of the RLL method is rather
straightforward resorting to the numerical routines available to solve the
Lyapunov equation based on the Bartels-Stewart algorithm, as implemented for
instance in the SciPy library (Virtanen et al., 2020). Indeed, one can easily
reach sizes of $N\sim 10^{3}$. Some convergence issues may arise in the case
of strong degeneracies Andreucci  _et al._ (2022). The numerical
implementation of the Green’s function method can be more involved than the
one of the RLL method. Indeed, we need to numerically invert the matrix in the
definition of the Green’s function (17) in the range of $\omega$ where the
transmission is non-vanishing in order to be able to compute the integral in
(18). Furthermore, the sampling over $\omega$ has to be fine enough to ensure
accuracy, especially if the transmission coefficient oscillates rapidly. This
difficulty does occur in our model, as it will be clear in what follows. In
practice, it is difficult to study lattices larger than $N\sim 10^{2}$ using
this method. The generalized eigenvalues method has the advantage of reducing
the problem to the calculation of the eigenvalues and eigenvectors of a
$2N\times 2N$ matrix Freitas  _et al._ (2015), which can be done by standard
linear algebra routines, the main limitation being memory storage and accuracy
of very small eigenvalues and avoiding the sampling problem.

##  III Heat Flux

In the short-range case, $\sigma=\infty$, two of the methods outlined above
have been used to obtain exact analytical results for the heat flux in the
thermodynamic limit Rieder  _et al._ (1967); Dhar (2008). This is possible
because the matrix of the interactions $\Phi$ reduces to the discrete
Laplacian, which is a tridiagonal matrix. In our case the matrix $\Phi$ is
dense, and we are unable to either solve analytically the Lyapunov equation or
to exactly compute the Green function. Nonetheless, it is possible to obtain a
certain amount of informations about the heat flux numerically.

###  III.1 Small coupling

If the coupling with baths $\lambda$ is small, a perturbative calculation of
the steady-state current is possible in terms of the eigenvalues and
eigenvectors of the isolated harmonic chain. This approach yields the so-
called Matsuda-Ishii’s formula, whereby $\mathcal{J}\approx\mathcal{J}_{MI}$
to the leading order in the coupling constant Matsuda and Ishii (1970); Lepri
_et al._ (2003), with $\mathcal{J}_{MI}$ given by

$\mathcal{J}_{MI}=\lambda\Delta
T\sum_{k}\frac{\psi^{2}_{k,1}\psi^{2}_{k,N}}{\psi^{2}_{k,1}+\psi^{2}_{k,N}}$
(23)

where $\Delta T=T_{L}-T_{R}$ and $\psi_{k,,n}$ denotes the $n$ component of
the $k$th eigenvector of the matrix $\Phi$ defined in (3). For the model we
consider here (which is homogeneous and mirror-symmetric, i.e. the first and
last component of each eigenvector are equal $\psi_{k,1}=\psi_{k,N}$ for
$k=1...N$) the above expression simplifies to

$\mathcal{J}_{MI}=\dfrac{\lambda\Delta
T}{2}\sum_{k}\psi_{k,1}^{2}=\dfrac{\lambda\Delta T}{2},$ (24)

where in the last step we used the property of completeness of the set of
eigenvectors. Note that eq. (24) expresses the fact that the chain is a
ballistic conductor.

Typically, in the short-range case $\sigma\rightarrow\infty$, this result
applies for $\lambda\ll\lambda_{0}\approx\mathcal{O}(1)$. In the our long-
range case, however, the situation is more complicated. In Fig. 1, we compare
formula (24) and the numerical solution of the Lyapunov equation. As we can
see, (24) holds for $\lambda$ smaller than a certain threshold
$\lambda_{0}(\sigma,N)$, that depends both on $N$ and on $\sigma$. More
specifically, $\lambda_{0}$ decreases with $\sigma$ and with $N$. On the other
hand, for $\sigma>0$ the perturbative approximation holds well in the
considered range.

Figure 1: Plots of the ratio between the heat flux $\mathcal{J}$, computed
numerically with the RLL method, and the Matsuuda-Ishii heat flux (24) versus
the system size $N$ for several values of $\sigma$ and $\lambda$ in the weak
coupling regime.

To have some insight into these deviations we may perform some further checks.
Usually the perturbative approach is justified assuming that the separation of
the unperturbed normal mode frequencies is smaller than the typical
dissipation caused by the coupling with the baths (Freitas and Paz, 2014b).
This assumption can actually be checked by examining the poles $s_{a}$. In
particular, we compare the spacings between the imaginary parts of consecutive
poles $Im(s_{a+1}-s_{a})$ and the real parts $Re(s_{a})$ . As we can see from
Fig. 2, the former is always much larger than the latter, therefore this
assumption is justified. This suggests that the observed deviations from the
Matsuda-Ishii formula may have a different origin.

Figure 2: Plots of the spacing between the imaginary parts of the poles of the
Green’s function $Im(s_{k+1})-Im(s_{k})$ (circles) and the real parts of the
poles $Re(s_{k})$ (crosses) for $\sigma=-0.5$. Different colors correspond to
different system’s size: $N=256,512,1024$ in blue, orange, green,
respectively.

###  III.2 Strong coupling

We now want to understand how the flux scales with the system size $N$ for not
too weak coupling $\lambda$. In order to so, we computed the heat flux using
the RLL method for several values of $N$ and $\sigma$ for $\lambda=1$ (and we
will set $\lambda=1$ for the rest of the paper) As shown in Fig.3 the data can
be fitted with a power law $\mathcal{J}\propto N^{-\gamma}$.

Although the direct computation of the Green’s function is numerically
cumbersome, we can easily compute its poles, compute the heat flux according
to (22) and fit a power law as we did before. In panel $b)$ of Fig. 4 we
report both the exponents fitted with the generalized eigenvalues method and
with the RLL method. As we can see, they are qualitatively in agreement.

The results of fits using the two methods are reported in Fig. 4. We can
identify three regions. The region close to the mean-field case $\sigma=-1$
and the one close to the short-range case $\sigma>1$, where finite-size
effects are almost absent, and an intermediate region in which finite-size
effects are quite strong. We also note that $\gamma$ seems to be converging to
the short-range value $\gamma=0$ while $\sigma$ goes to $1$. Summarizing, even
if we are not able to extract the exact values of the exponents, it is clear
that the flux scales with some nontrivial power of the system’s size $N$.

Figure 3: Log-log plot of the heat flux $\mathcal{J}$ versus the system’s size
$N$ for $\lambda=1$ and different values of the long-range exponent $\sigma$.
The flux is computed using the RLL method as described in the text. Figure 4:
Plot of the scaling exponent of the flux $\gamma$, defined as
$\mathcal{J}\propto N^{-\gamma}$, $(a)$, we report the exponents obtained by
fitting a power law on the heat flux obtained with the RLL method. To check
the finite-size effects, each data set corresponds to a fit over different
length ranges, $50\leq N\leq 1600$ (circles), $500\leq N\leq 2000$ (squares),
$1500\leq N\leq 7500$ (triangles). Panel $(b)$, comparison between the
exponents obtained by the RLL method (circles) and the generalized eigenvalue
method (triangles).

##  IV Transmission spectra

To understand the origin of the nontrivial dependence of the flux on the size,
let us investigate the transmission spectrum of the chain. We begin by
plotting the transmission coefficient, namely the integrand in (18) as a
function of the frequency $\omega$. In Fig. 5 we report its plot for several
values of $\sigma$. We can see that it is characterized by a rather
complicated peak structure which consists of $N-2$ peaks (as can be checked
numerically).

A manin point we want to make and explore is that the structure of such
resonances determines the scaling of the current. Notice that a change of sign
in $\omega$ in (17) is equivalent to the complex conjugation of $G(\omega)$.
Since the transmission coefficent depends on the square modulus of $G(\omega)$
it is an even function of $\omega$ and we can therefore restrict ourselves to
study positive frequencies. Let us denote by $\omega_{k}$, $k=1,2\ldots$ the
location of the peak frequencies for positive $\omega$. The peaks accumulate
at a band-edge frequency ${\omega_{B}}<2$, i.e $\omega_{k}\to{\omega_{B}}$ for
$k$ large. Furthermore, upon approaching ${\omega_{B}}$, the width of the
peaks decreases. Notice that this is the reason why it is important to finely
sample the Green’s function in $\omega$, especially in the proximity of the
band edge. Indeed, we used a logarithmic sampling in order to increase the
sampling points near ${\omega_{B}}$. The integrand is thus a much more
complicated function of $\omega$ with respect to the mean-field case
$\sigma=-1$ Defaveri  _et al._ (2022); Andreucci  _et al._ (2022), where only
the first peak is present. It can be checked numerically that the first few
peaks are Lorentzian with amplitude $\Delta_{k}\approx N^{-1}$, exactly like
the peak in mean-field case. The subsequent peaks are too narrow to be
resolved. For positive values of $\sigma$ the situation becomes even more
complicated, as a curve emerges below the peaks, as we can see in Fig. 5 for
$\sigma=0.5$.

Figure 5: Panels $a)$, $b)$, $c)$, $d)$ : transmission spectra (the integrand
of the heat flux expression (18)) for different values of the range exponent
$\sigma=-0.7,0.5,0.1,0.5$ and for a chain with $N=100$. Only positive
frequencies are reported. Panels $e)$, $f)$, $g)$, $h)$: rescaled cumulative
function $N^{\gamma}F(\omega)$, for $N=80,100,120,140$ and
$\sigma=-0.7,-0.5,-0.1,0.5$ in panels $a)$, $b)$, $c)$, $d)$ respectively. The
values of $\gamma$ are taken from the blue points in Figure 4. The abrubt
increase of the cumulative function in panels $e)$ and $f)$ at $\omega\approx
1.3$ is due to the dominant contribution of the first peak in panels $a)$ and
$b$. The subsequent, smaller, jumps are due to the contributions of the other
peaks.

For the reasons outlined above, it seems more convenient to consider the
cumulative function $F(\omega)$, that is, the integral (18) performed up to
frequency $\omega$. In the rightmost panels of Fig. 5 we report the function
$F(\omega)$ for several values of $N$ of order $10^{2}$ and $\sigma$, rescaled
by $N^{\gamma}$, where $\gamma$ is the exponent obtained with the RLL method
for values of $N$ of order $10^{2}:10^{3}$. As we can see, the curves nicely
collapse for $\sigma=-0.7,-0.5$, but for higher values of $\sigma$, such as
$\sigma=-0.3$, the collapse is not as good due to the finite-size effects, as
expected. Regardless of the lack of further quantitative progress in the
computation of the exponents, the qualitative information about the peak
structure will be crucial in our understanding of the model, as we will see
later.

##  V Poles of the Green’s function

In view of the numerical difficulties encountered above and for comparison, we
also performed a study of poles of the Green’s function. These are computed
through the generalized eigenvalue method described above.

The main advantage of the analysis is that we gain a new perspective on the
peak structure discussed before. Indeed, the positions $\omega_{k}$ of the
peaks in Fig. 5 are given by the absolute value of the imaginary part of
$s_{a}$, while the absolute value of the imaginary part should be proportional
to their widths $\Delta_{k}$.

Figure 6: Real parts of the poles of the Green’s functions $s_{a}$ versus the
distance of their imaginary parts from the band-edge. Leftmost panel:
$\sigma<0$, vertical axis rescaled by $N^{\delta}$ with $\delta\approx
1-\sigma$. The inset in panel (a) demonstrates the different scaling the
collapse for the widths of the first peaks (the first $80,160,320,640,1280$
for $N=256,512,1024,2048,4096$, respectively) rescaled by $N$. For the other
values of $\sigma$, we get the same scaling for the first peaks. Rightmost
panel: same for $\sigma>0$, with vertical axis rescaled by $N$. Note that such
scaling works for the whole spectrum in this case.

In particular, we consider all the peaks as Lorentzian – for simplicity, but
also because all the peaks that we were able to resolve are actually very well
approximated by a Lorentzian – with width given by $\Delta_{k}(N)=Re(s_{k})$.
In this approximation, as far as scaling with the size is concerned, the heat
flux can be estimated as the sum of the widths of the peaks $\Delta_{k}(N)$.
Furthermore, the height of each peak can be shown to be equal to
$\lambda^{2}/4$ (indeed, note that in Fig. 5, in which $\lambda=1$, the
heights of the peaks are all the same and equal to $1/4$). Thus, we replace
the integrand in eq. (18) with a sum of normalized Lorentzians, and we get:

$\dfrac{\mathcal{J}(N)}{\Delta
T}\approx\int_{-\infty}^{\infty}\dfrac{d\omega}{\pi}\sum_{k=1}^{N-2}\dfrac{\lambda^{2}\Delta_{k}(N)^{2}/4}{(\omega-\omega_{k})^{2}+\Delta_{k}(N)^{2}}=\frac{\lambda^{2}}{4}\sum_{k=1}^{N-2}\Delta_{k}(N).$
(25)

The relevant information should thus be contained in the dependence of the
$\Delta_{k}$ on $k$ and $N$. Physically, this is the effective damping of
plane waves due to the coupling with the thermal reservoirs.

The dependence of $\Delta_{k}$ on $N$ is reported in Fig. 6, where we plot
(parametrically) the real parts of the poles as a function of the imaginary
ones, for negative and positive values of $\sigma$, respectively. Since the
resonances accumulates at the band-edges, it is convenient to report the
frequencies as a function of their relative distance from $\omega_{B}$. Let us
focus on the case of negative $\sigma$, to begin with. From the leftmost
panels of Fig. 6, it is seen that the poles can be grouped in two sets, each
having different dependencies on $\omega_{k}$ and $N$. Empirically, this is
accounted for by the following scaling:

$\Delta_{k}(N)\approx\begin{cases}d_{k}/N,\quad k<k_{o}\\\
d_{k}/N^{\delta},\quad k>k_{o},\end{cases}$ (26)

where $k_{o}<<N$ and $d_{k}$ do not depend $N$. We do not have an a-priori
theoretical estimate of $\delta$, but we find that there is a good collapse
upon choosing $\delta\approx 1+|\sigma|$. It is interesting to point out that
the exponent $\delta$ can be interpreted as the fractal dimension of area
below the graphs in Fig. 5. Indeed, if we increase the system’s size $N$ new
peaks emerge with progressively shrinking area and, in a putative
$N\rightarrow\infty$ limit we would have an infinite number of peaks with
vanishing area.

In addition, there are a few poles whose widths do not follow this scaling and
fall consistently well outside the collapsed curve. It actually turns out that
there are two degenerate eigenvalues between the $s_{a}$s that do not follow
the scaling law. However this is inconsequential, as one can check that the
contribution of the these eigenvalues to (22) vanishes. Heuristically, this is
because, as one can check, the eigenvectors related to these eigenvalues are
localized at the endpoints of the chain and therefore do not contribute to
transport. This also explains why the peaks in Fig. 5 are $N-2$ instead of
$N$. We can therefore infer the following scaling law for the heat flux (22)
plugging (26) into (25):

$\mathcal{J}\approx\dfrac{\sum_{k=1}^{k_{o}}d_{k}}{N}+\dfrac{\sum_{k=k_{o}}^{N}d_{k}}{N^{\delta}}\propto
N^{1-\delta}.$ (27)

The first term scales as $N^{-1}$, since $k_{o}$ does not scale with $N$ (as
can be inferred from Fig.6). On the other hand, the second term scale as
$N^{1-\delta}$ since each $d_{k}$ is of order $1$ and thus their sum scales as
$N$. Finally, since $\delta>0$, we get the reported scaling for the heat flux.
For positive $\sigma$, the scaling of $\Delta_{k}$ is reported in the right-
most panels of Fig. 6: as we can see in this case $\Delta_{k}\approx N^{-1}$,
over the entire spectrum. Therefore, the estimate the heat flux yields

$\mathcal{J}\approx\dfrac{\sum_{k=1}^{N}d_{k}}{N}\approx\mathcal{O}(1).$ (28)

So the heat flux for positive $\sigma$ behaves as the heat flux for
$\sigma=\infty$ (the nearest-neighboors case), that is, it does not scale with
$N$.

To summarize, according to approximation (25) and the numerical estimate of
$\delta$ extracted from the data, we find that the heat flux scale as:

$\mathcal{J}\propto
N^{-\tilde{\gamma}},\quad\quad\tilde{\gamma}\approx\begin{cases}1-\delta,\quad\sigma<0,\\\
0,\quad\sigma>0.\end{cases}$ (29)

As we already mentioned, see Fig. 6, we found a good collapse of the imaginary
part of the poles of the Green’s functions for $\delta\approx 1-|\sigma|$. So
this yields

$\tilde{\gamma}\approx-\sigma$ (30)

for negative $\sigma$. Admittedly, this estimate accounts only qualitatively
for the behavior of the exponents as given in Fig. 4. The deviations are
sizeable and, in addition the dependence of $\gamma$ on $\sigma$ appears to be
non-linear. While this could be due to the aforementioned finite-size effects,
the discrepancy is present even for values of $\sigma$ for which the exponent
$\gamma$ has basically converged (for example $\sigma=-0.7,-0.5$). Another
possibility, which seems more likely, is that, while the widths of the peaks
of Fig. 5 are indeed related to the real parts of $s_{a}$ on general grounds,
they are not exactly equal. On the other hand, we point out that, since the
$s_{a}$ are related to the widths of the peaks, the transition in the scaling
of the $\Delta_{k}$s at $\sigma=0$ suggests that the scaling of the heat-flux
between the short-range and the long-range behaviour has to occur at
$\sigma=0$.

##  VI Conclusions

Heat transport in short-range linear systems has been widely studied Lepri
_et al._ (2003). On the contrary, the behaviour of linear oscillators with
long-range power-law couplings is not yet well understood beyond the mean-
field (fully-coupled) case Defaveri  _et al._ (2022); Andreucci  _et al._
(2022). In this paper, we have made a step forward along this direction by
applying three different methods Rieder  _et al._ (1967); Dhar (2008);
freitas2014analytic that allow to compute numerically both the heat flux and
its scaling with the system’s size. All the methods give a clear scaling of
the current with a power-law in the system’s size. This scaling interpolates
between the short-range behaviour, where the current is constant in the
system’s size, and the mean-field behaviour, where the current is inversely
proporational to the system’s size. However, the fitted scaling exponents show
significant finite-size effects for all the three methods. The method of ref.
Rieder  _et al._ (1967) which consists in solving a matricial equation is
straightforwardly applicable to the long-range case. The Green’s function
approach allows to express the current as an integral over frequencies, which
cannot be solved analytically. However, the integrand has the interesting
property of showing a sequence of peaks that accumulate near the band edges of
the spectrum. Further properties of these peaks can be inferred using the
third method, which allows to compute the poles of the Green’s function.
Indeed, the real and the imaginary part of these poles are related to the
position and the width of the peaks, respectively. We find a sharp transition
in the scaling of the real parts of the poles at the value of the long-range
coupling exponent $\sigma=0$ corresponding to the transition between the long-
range and the short-range behaviour of the system. The crucial problem is now
the dependence on $\sigma$ of the scaling exponent of the current. Assuming
that all of the peaks of the integrand are well-separated Lorentzians and that
their widths are exactly given by the real parts of the poles, we might
conclude that the heat current scales as $\mathcal{J}\propto N^{-|\sigma|}$
for $-1<\sigma<0$. cin agreement with the one derived directly from the fit of
the current, which is anyway affected – at least for small values of
$|\sigma|$ – by significant finite-size effects. The disagreement between
these two scaling exponents remains to be explored, even though our analysis
of the scaling of the real part of the poles of the Green’s function clearly
supports the presence of a transition at $\sigma=0$ from the long-range to the
short-range behaviour.

###### Acknowledgements.

We gratefully thank Celia Anteneodo and Lucianno Defaveri for useful
discussions. SL and SR acknowledge partial support from project MIUR-PRIN2017
Coarse-grained description for non-equilibrium systems and transport phenomena
(CO-NEST) n. 201798CZL.

## References

  * Rieder  _et al._ (1967) Z. Rieder, J. L. Lebowitz, and E. Lieb, Journal of Mathematical Physics 8, 1073 (1967).
  * Zurcher and Talkner (1990) U. Zurcher and P. Talkner, Physical Review A 42, 3278 (1990).
  * Dhar and Roy (2006) A. Dhar and D. Roy, Journal of Statistical Physics 125, 801 (2006)).
  * Saito  _et al._ (2000) K. Saito, S. Takesue, and S. Miyashita, Physical Review E 61, 2397 (2000).
  * Asadian  _et al._ (2013) A. Asadian, D. Manzano, M. Tiersch, and H. Briegel, Physical Review E 87, 012109 (2013).
  * Visscher (1971) W. M. Visscher, Progress of Theoretical Physics 46, 729 (1971).
  * Dhar (2001) A. Dhar, Physical review letters 86, 5882 (2001).
  * Herrera-González  _et al._ (2015) I. Herrera-González, F. Izrailev, and L. Tessieri, Europhysics Letters 110, 64001 (2015).
  * Ash  _et al._ (2020) B. Ash, A. Amir, Y. Bar-Sinai, Y. Oreg, and Y. Imry, Phys. Rev. B 101, 121403(R) (2020).
  * Schmidt  _et al._ (2013) M. Schmidt, T. Kottos, and B. Shapiro, Physical Review E 88, 022126 (2013).
  * Santra and Basu (2022) I. Santra and U. Basu, SciPost Phys. 13, 041 (2022).
  * Lepri  _et al._ (2003) S. Lepri, R. Livi, and A. Politi, Physics Reports 377, 1 (2003).
  * Dhar (2008) A. Dhar, Adv. Phys. 57, 457 (2008).
  * Lepri  _et al._ (2016) S. Lepri, R. Livi, and A. Politi, in  _Thermal transport in low dimensions: From statistical physics to nanoscale heat transfer_ (Springer, 2016) pp. 1–37.
  * Benenti  _et al._ (2020) G. Benenti, S. Lepri, and R. Livi, Frontiers in Physics 8, 292 (2020).
  * Benenti  _et al._ (2023) G. Benenti, D. Donadio, S. Lepri, and R. Livi, La Rivista del Nuovo Cimento 46, 105 (2023).
  * Di Cintio  _et al._ (2018) P. Di Cintio, S. Iubini, S. Lepri, and R. Livi, Chaos, Solitons & Fractals 117, 249 (2018).
  * Lepri  _et al._ (2020) S. Lepri, R. Livi, and A. Politi, Phys. Rev. Lett. 125, 040604 (2020).
  * Bermúdez  _et al._ (2013) A. Bermúdez, M. Bruderer, and M. B. Plenio, Physical review letters 111, 040601 (2013).
  * Ramm  _et al._ (2014) M. Ramm, T. Pruttivarasin, and H. Haeffner, New Journal of Physics 16, 063062 (2014).
  * Molerón  _et al._ (2019) M. Molerón, C. Chong, A. J. Martínez, M. A. Porter, P. G. Kevrekidis, and C. Daraio, New Journal of Physics 21, 063032 (2019).
  * Campa  _et al._ (2009) A. Campa, T. Dauxois, and S. Ruffo, Phys. Rep. 480, 57 (2009).
  * Defenu  _et al._ (2021) N. Defenu, T. Donner, T. Macrì, G. Pagano, S. Ruffo, and A. Trombettoni, (2021), arXiv:2109.01063 [cond-mat.quant-gas] .
  * Torcini and Lepri (1997) A. Torcini and S. Lepri, Phys. Rev. E 55, R3805 (1997).
  * Métivier  _et al._ (2014) D. Métivier, R. Bachelard, and M. Kastner, Phys. Rev. Lett. 112, 210601 (2014).
  * Olivares and Anteneodo (2016) C. Olivares and C. Anteneodo, Phys. Rev. E 94, 042117 (2016).
  * Iubini  _et al._ (2018) S. Iubini, P. Di Cintio, S. Lepri, R. Livi, and L. Casetti, Phys. Rev. E 97, 032102 (2018).
  * Bagchi (2017) D. Bagchi, Phys. Rev. E 95, 032102 (2017).
  * Di Cintio  _et al._ (2019) P. Di Cintio, S. Iubini, S. Lepri, and R. Livi, Journal of Physics A: Mathematical and Theoretical 52, 274001 (2019).
  * Wang  _et al._ (2020) J. Wang, S. V. Dmitriev, and D. Xiong, Physical Review Research 2, 013179 (2020).
  * Bagchi (2021) D. Bagchi, Physical Review E 104, 054108 (2021).
  * Iubini  _et al._ (2022) S. Iubini, S. Lepri, and S. Ruffo, Journal of Statistical Mechanics: Theory and Experiment 2022, 033209 (2022).
  * Tamaki and Saito (2020) S. Tamaki and K. Saito, Phys. Rev. E 101, 042118 (2020).
  * Suda (2022) H. Suda, Nonlinearity 35, 2288 (2022).
  * Defaveri  _et al._ (2022) L. Defaveri, C. Olivares, and C. Anteneodo, Phys. Rev. E 105, 054149 (2022).
  * Andreucci  _et al._ (2022) F. Andreucci, S. Lepri, S. Ruffo, and A. Trombettoni, SciPost Phys. Core 5, 036 (2022).
  * Schuckert  _et al._ (2020) A. Schuckert, I. Lovas, and M. Knap, Phys. Rev. B 101, 020416(R) (2020).
  * Zoia  _et al._ (2007) A. Zoia, A. Rosso, and M. Kardar, Physical Review E 76, 021116 (2007).
  * Defenu (2021) N. Defenu, Proceedings of the National Academy of Sciences 118, e2101785118 (2021).
  * Tarasov (2006) V. E. Tarasov, Journal of Physics A: Mathematical and General 39, 14895 (2006).
  * Dhar and Saito (2016) A. Dhar and K. Saito, in  _Thermal Transport in Low Dimensions_ (Springer, 2016) pp. 39–105.
  * Freitas and Paz (2014b) N. Freitas and J. P. Paz, Phys. Rev. E 90, 042128 (2014b).
  * Freitas  _et al._ (2015) N. Freitas, E. A. Martinez, and J. P. Paz, Physica Scripta 91, 013007 (2015).
  * Tisseur and Meerbergen (2001) F. Tisseur and K. Meerbergen, SIAM Review 43, 235 (2001), https://doi.org/10.1137/S0036144500381988 .
  * Virtanen et al. (2020) P. S. Virtanen et al., Nature Methods 17, 261 (2020).
  * Matsuda and Ishii (1970) H. Matsuda and K. Ishii, Progress of Theoretical Physics Supplement 45, 56 (1970).