# Parametric simulation studies on the wave propagation of solar radio
emission:
the source size, duration, and position

PeiJin Zhang CAS Key Laboratory of Geospace Environment, School of Earth and
Space Sciences,
University of Science and Technology of China (USTC), Hefei, Anhui 230026,
China  CAS Center for the Excellence in Comparative Planetology, USTC, Hefei,
Anhui 230026, China  ChuanBing Wang CAS Key Laboratory of Geospace
Environment, School of Earth and Space Sciences,
University of Science and Technology of China (USTC), Hefei, Anhui 230026,
China  CAS Center for the Excellence in Comparative Planetology, USTC, Hefei,
Anhui 230026, China  Eduard P. Kontar School of Physics and Astronomy,
University of Glasgow, Glasgow G12 8QQ, UK

(Received October 19, 2020; Accepted xxx)

###### Abstract

The observed features of the radio sources indicate complex propagation
effects embedded in the waves of solar radio bursts. In this work, we perform
ray-tracing simulations on radio wave transport in the corona and
interplanetary region with anisotropic electron density fluctuations. For the
first time, the variation of the apparent source size, burst duration, and
source position of both fundamental emission and harmonic emission at
frequency 35 MHz are simulated as the function of the anisotropic parameter
$\alpha$ and the angular scattering rate coefficient
$\eta=\epsilon^{2}/h_{0}$, where $\epsilon^{2}={\langle\delta
n^{2}\rangle}/{n^{2}}$ is the density fluctuation level and $h_{0}$ is its
correlation length near the wave exciting site. It is found that isotropic
fluctuations produce a much larger decay time than a highly anisotropic
fluctuation for fundamental emission. By comparing the observed duration and
source size with the simulation results in the parameter space, we can
estimate the scattering coefficient and the anisotropy parameter
$\eta=8.9\times 10^{-5}\,\mathrm{km^{-1}}$ and $\alpha=0.719$ with point pulse
source assumption. Position offsets due to wave scattering and refraction can
produce the co-spatial of fundamental and harmonic waves in observations of
some type III radio bursts. The visual speed due to the wave propagation
effect can reach 1.5 $c$ for $\eta=2.4\times 10^{-4}\,\mathrm{km^{-1}}$ and
$\alpha=0.2$ for fundamental emission in the sky plane, accompanying with
large expansion rate of the source size. The visual speed direction is mostly
identical to the offset direction, thus, for the observation aiming at
obtaining the source position, the source centroid at the starting point is
closer to the wave excitation point.

solar radio burst — source size and position — wave propagation effects

††journal: xxx

##  1 Introduction

The imaging and spectroscopy observations of the solar radio bursts can
provide information on the non-thermal electrons associated with the transient
energy release in the solar active region and parameters of the background
plasma. For example, McCauley et al. (2018) inspected the background density
of the solar corona with the interferometric imaging of Type III radio bursts
in the frequency range of 80-240 MHz. High-cadence radio imaging spectroscopy
shows evidence of the existence of particle acceleration by solar flare
termination shock (Chen et al., 2015; Yu et al., 2020). The combined
observations of radio imaging and extreme ultraviolet /white-light imaging
indicate that the particle acceleration occurs at the flank of coronal mass
ejection shock (Chrysaphi et al., 2018; Morosan et al., 2019; Chen et al.,
2014). The imaging spectroscopy study with LOw Frequency Array (LOFAR) reveals
that the velocity dispersion of the electron beams is a key factor for the
duration of type III radio burst (Zhang et al., 2019).

The corona plasma is an in-homogeneous refractive media for solar radio waves,
so the refraction and scattering can cause the deformation of the observed
radio source including the expansion of the source size, the offset of the
visual source position from the wave generation position (Wild et al., 1959;
Kontar et al., 2017; Bisoi et al., 2018), and the duration broadening in the
dynamic spectrum (Zhang et al., 2019). Generally, the refraction can cause the
inward offset of the visual source from the position of wave excitation (Mann
et al., 2018). The visual source could be shifted outward if the scattering is
considered (Stewart, 1972, 1976; Arzner & Magun, 1999; Kontar et al., 2019).
For the wave of solar radio burst with fundamental-harmonic (F-H) pair
structure, the H-emission is generated at a higher height than the F-emission
with the same frequency due to the plasma emission mechanism, while imaging
result indicates that the source of F and H emission can have a similar
position for some bursts (Stewart, 1972; Dulk & Suzuki, 1980). This seems to
indicate that the waves of F-emission and H-emission have experienced
different amounts of refraction and scattering.

For a given frequency, the Langmuir wave excitation responsible for type III
emission is limited to the cross-section of the guiding magnetic field line
for the electrons and the layer with local plasma frequency close to $1$ or
$1/2$ times of the radio wave frequency. The structure of the magnetic field
and the distribution of the background density is stable within the time scale
of electron beam transit time. Thus the excitation site at a fixed frequency
is stable in both size and position. While the visual velocity and expansion
rate of the observed sources can be very large. Kontar et al. (2017) observed
1/4 $c$ radial speed with LOFAR beamformed observation of the fundamental part
of a type IIIb striae, where $c$ is the speed of light. Kuznetsov et al.
(2020) observed 1/3 $c$ at 32 MHz in a drift pair radio burst. The
interferometric imaging with the remote baseline of LOFAR finds about 4 $c$
visual speed for the fundamental part of a type III-IIIb pair event source at
26 MHz (Zhang et al., 2020). The observed source size of the radio burst is
also considerably larger than the size estimated from the bandwidth of the
striae (Kontar et al., 2017; Zhang et al., 2020). The fast variations of the
source is believed to be due to the wave propagation effects.

The ray-tracing simulation is an effective method to investigate the wave
propagation effects, which can help us interpreting the imaging and
spectroscopy observation and diagnosing the solar corona properties from the
observation results (Steinberg et al., 1971; Riddle, 1974). The ray-tracing
method is introduced to the solar radio study by Fokker (1965) to investigate
the source size of type I radio burst due to the scattering effect. Bougeret &
Steinberg (1977) indicated that the over-dense fiber structures in the coronal
should be considered to understand the position of the radio source, the
moving bursts, and the great variety of space-time shapes observed within the
same storm center, and Robinson (1983) introduced fiber in-homogeneity into
the ray tracing simulation and find that the in-homogeneity of the medium can
account for the source displacement. It has been pointed out that the
anisotropic scattering of the wave due to the statistical inhomogeneity is
essential to interpret the observed radio source properties (Arzner & Magun,
1999; Kontar et al., 2019; Bian et al., 2019). Recently, Kontar et al. (2019)
developed a model to include the anisotropic scattering of the wave in the ray
tracing process. It was found that the duration of type III radio burst
decrease with the anisotropy level of the background, and the anisotropic
density fluctuations are required to account for the source sizes and decay
times simulaneously. The theory of anisotropic scattering has been
successfully used to interpret the source properties of the drift pair
bursts(Kuznetsov & Kontar, 2019; Kuznetsov et al., 2020).

These case studies of observation with the corresponding simulation have shown
that the ray-tracing is an effective method to analyze the wave propagation
effect, while the detailed dependency relation between the observed source
characteristics and the density fluctuation properties of the medium is still
unclear. A further understanding of the physical process behind the visual
source variation requires more detailed simulation-observation comparison
study.

In this work, for the first time, we performed a large set of Monte Carlo
simulations of the radio wave propagation by ray tracing to explore a
parameter-space study on the radio source properties with various background
plasma parameters. The paper is arranged as flows: in Section 2, the model
used in the simulation is introduced. Section 3 shows the simulation results
of the radio source size, duration, and position for fundamental waves and
harmonic waves at 35 MHz. The results are discussed and compared in Section 4,
and a brief summary is given in Section 5.

##  2 Simulation model

We implemented the three-dimensional (3D) radio wave ray-tracing simulation
for a point pulse source on basis of the theory and algorithm proposed by
Kontar et al. (2019). The simulation solves the Hamilton ray equation and
nonlinear Langevin equation corresponding to the Fokker-Planck equation, which
can be expressed as

$\displaystyle\frac{\mathrm{d}r_{i}}{\mathrm{d}t}$
$\displaystyle=\frac{\partial\omega}{\partial
k_{i}}=\frac{c^{2}}{\omega}k_{i},$ (1)
$\displaystyle\frac{\mathrm{d}k_{i}}{\mathrm{d}t}$
$\displaystyle=-\frac{\partial\omega}{\partial r_{i}}+\frac{\partial
D_{ij}}{\partial k_{i}}+B_{ij}\xi_{j},$ (2)

where $r_{i}\,(i=x,y,z)$ is the position vector of photons, $k_{i}$ is the
wave vector of the radio wave in Cartesian coordinate. $\omega$ is the angular
frequency of the wave satisfying the dispersion relation of the unmagnetized
plasma:

$\omega^{2}=\omega_{pe}^{2}+c^{2}k^{2}.$ (3)

Here $\omega_{pe}$ is the local plasma frequency expressed as
$\omega_{pe}(\bm{r})=\sqrt{{4\pi e^{2}n(\bm{r})}/{m_{e}}}$ (where $e$ and
$m_{e}$ are respectively the electron charge and mass, and $n$ is the plasma
density). $D_{ij}$ is the diffusion tensor appropriate to scattering, which is
given by:

$D_{ij}=\frac{\pi\omega_{pe}^{4}}{4\omega c^{2}}\int
q_{i}q_{j}S(\bm{q})\delta(\bm{q}\cdot\bm{k})\frac{\mathrm{d}^{3}q}{(2\pi)^{3}},$
(4)

where $\bm{q}$ is the wave-vector of the density fluctuation. $S(\bm{q})$ is
the spectrum of the density fluctuation normalized to the relative density
fluctuation variance:

$\epsilon^{2}=\frac{\langle\delta n^{2}\rangle}{n^{2}}=\int
S(\bm{q})\frac{\mathrm{d}^{3}q}{(2\pi)^{3}},$

where $n$ is the local average plasma density, taken to be a slowly varying
function of position. $B_{ij}$ is a positive-semi-definitive matrix given by
$D_{ij}=(1/2)B_{im}B_{jm}^{T}$. $\bm{\xi}$ is a random vector satisfying the
Gaussian distribution with zero mean and unit variance.

To consider the scattering in a medium with anisotropic density fluctuation,
an anisotropy tensor

$\mathsf{A}=\left({\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&\alpha^{-1}\end{array}}\right)$ (5)

is introduced. Here $\alpha$ is the anisotropic parameter representing the
ratio of density fluctuation wavelength in the direction of anisotropy and its
perpendicular direction. In the model of this work, the density fluctuations
are mostly in the direction perpendicular to the solar radius when $\alpha<1$,
and vice versa.

From Equation 4, one can get the angular scattering rate per unit distance for
a classic Gaussian spectrum of density fluctuations (Steinberg et al., 1971;
Chrysaphi et al., 2018; Kontar et al., 2019),

$\frac{\mathrm{d}\langle\theta^{2}\rangle}{\mathrm{d}x}=\frac{\sqrt{\pi}}{2}\frac{\epsilon^{2}}{h}\frac{\omega^{4}_{pe}}{(\omega^{2}-\omega_{pe}^{2})^{2}}$
(6)

where $h$ is the correlation length of density fluctuations. In this work, the
density fluctuation is characterized as power-law distribution with a cutoff
at the inner scale of $l_{i}=(r/R_{s})\,[\mathrm{km}]$ and outer scale of
$l_{o}=0.25R_{s}(r/R_{s})^{0.82}$ (Coles & Harmon, 1989; Wohlmuth et al.,
2001; Krupar et al., 2018), where $R_{s}$ is the solar radius, $r$ is the
heliocentric distance, the inner scale and outer scale represent the smallest
and largest wave length of the density fluctuations, respectively. The angular
scattering rate per unit length derived from the power-law fluctuation
(Equation (67) in Kontar et al. (2019)) can be expressed in a similar form of
Equation 6 by replacing $h$ with an equivalent scale length $h_{eq}$ given as

$h\equiv h_{eq}=l_{o}^{2/3}l_{i}^{1/3}/\pi^{3/2},$ (7)

and $h$ is a slowly varying function of heliocentric distance $r$.

For a given wave frequency, the angular scattering rate depends on $\epsilon$,
$h$, and the frequency ratio $\omega_{pe}/\omega$. As $h^{-1}$ and
$\omega_{pe}/\omega$ in Equation 6 decreasing with $r$, the scattering rate
decreases quickly when the radio waves propagate outward from their original
site. As a result, the scattering strength experienced by the waves is mainly
determined by the coefficient $\eta$ of the angular scattering rate per unit
distance in Equation 6, which is defined by

$\eta=\epsilon^{2}/h_{0},$ (8)

where $h_{0}=h(r_{0})$ is the scale length of density fluctuations near the
wave exciting site at $r_{0}$. In this work, the function $h(r)$ is given,
therefore the tuning parameter for scattering rate is the density fluctuation
level $\epsilon$ for different simulation cases, but we reiterate that it is
$\eta$ the key parameter.

The intensity decay $I/I_{0}$ due to the Coulomb collision absorption is
described by the optical depth integral along the ray path for each photon:

$\displaystyle I/I_{0}$ $\displaystyle=e^{-\tau_{a}},$ (9)
$\displaystyle\tau_{a}$ $\displaystyle=\int\gamma(\bm{r}(t))\,dt,$ (10)
$\displaystyle\gamma$
$\displaystyle=\frac{4}{3}\sqrt{\frac{2}{\pi}}\frac{e^{4}n_{e}\ln\Lambda}{mv_{Te}^{3}}\frac{\omega_{pe}^{2}}{\omega^{2}},$
(11)

where $v_{Te}$ is the thermal speed and a constant Coulomb logarithm
$\ln\Lambda\simeq 20$ is assumed.

##  3 Simulation results

With the anisotropic scattering model, we can perform ray-tracing simulation
for the photons launched at a given position. After sufficient steps forward,
when the majority of the photons reached 0.95 AU of heliocentric distance, the
final positions ($\bm{r}$) and wave vectors ($\bm{k}$) of the photons near the
direction of observation are collected to reconstruct the radio image. The
observed position of each photon in the sky plane is estimated by back-
tracking the wave vector. The apparent source intensity map $I(x,y)$ is the
weighted photon number density distribution in the sky plane, and the spread
of photon arrival times determines the temporal variation of the source.

In the simulation, the background density function $n(r)$ is assumed to be
spherically symmetric, the same as the model used in Kontar et al. (2019). The
source of radio emission is considered to be a point pulse source, the photons
are launched at the starting point at the same time. The initiate wave vectors
are randomly distributed in an outward direction ($\bm{r}\cdot\bm{k}>0$),
namely, a burst source with isotropic emission is considered.

Figure 1: The simulated source intensity image in the sky plane for
fundamental emission, which is reconstructed from all collected photons. The
frequency of the wave is 35 MHz with $f_{0}/f_{pe}=1.1$, $\alpha=0.72$,
$\epsilon=0.28$, and $\eta\simeq 9\times 10^{-5}\,\textup{km}^{-1}$. A
simulation of $3\times 10^{5}$ photons is conducted. The black circle shows
the outline of the solar disk, the blue line shows the FWHM of the
reconstructed source. The green ‘+’ marks the starting point of the photons.
The upper and right panel shows the histogram of the flux intensity in $x$ and
$y$ direction.

As a case study, Figure 1 shows the flux intensity map reconstructed from all
collected photons for fundamental emission, representing the time integral
imaging of the observed source. Figure 2 shows the time-intensity profile and
the temporal variations of the source position and size, which is obtained
from photons with arrival times within the corresponding time segment. In this
simulation case, the radio waves are fundamental waves initiated at the
position angle $(\theta_{0},\phi_{0})=(20^{\circ},0^{\circ})$, where $\theta$
is the longitude and $\phi$ is the latitude. The density fluctuation level is
$\epsilon=0.28$ and the anisotropy parameter is $\alpha=0.72$. The equivalent
correlation length ($h_{0}$) at the starting point $r_{0}=1.750R_{s}$ is about
860 km, and the corresponding value of the angular scattering rate coefficient
is $\eta\simeq 9\times 10^{-5}\,\textup{km}^{-1}$ near the wave generation
site.

In Figure 1, the source centroid offsets $0.11R_{s}$ in $x$ direction from the
starting point (marked as green ‘+’ in the figure). The offset in $y$
direction is negligible. This is because the centroid offset is in the radial
direction due to the spherical symmetry of the density model used, which is in
the $x$ direction in the sky plane for the case shown in this figure. The blue
line in Figure 1 shows the full width half maximum (FWHM) estimation of the
reconstructed source. The statistical results of time duration, FWHM size, the
visual speed, and size expansion rate can be extracted from the reconstructed
imaging and its variation.

Figure 2: Upper panel: the simulated time-intensity profile. The flux
variation is fitted with a double Gaussian function, shown as a cyan line.
FWHM range of the fitted flux are marked as gray shadow and the peak of the
fitted flux is marked as a black thin line in all panels. Middle panel: the
position offset of the simulated source centroid from its original site, where
red and blue represents X and Y in heliocentric coordinate. Lower panel: the
FWHM size of the simulated source. The linear fitting results of position
offset and size for the rising and decay phase are shown in the middle and
lower panel as orange and purple lines. The simulation parameters are the same
as Figure 1.

The FWHM of the time-intensity profile can be divided into two phases, namely
the rising and decay phase. We fitted the flux variation with a double
Gaussian function (Zhang et al., 2019):

$G(t)=A\exp{\left(-\frac{(t-t_{0})^{2}}{\tau^{2}}\right)},\qquad\tau=\left\\{\begin{aligned}
\tau_{R},&&{t\leq t_{0}}\\\ \tau_{D},&&{t>t_{0}}\end{aligned}\right..$ (12)

The fitted result is shown as cyan line in Figure 2. The peak and the FWHM of
the flux is obtained from the fitted curve. In Figure 2, the FWHM range marked
by gray shadow is divided by a thin vertical line at the peak of flux
($t=1.83$ s). The regime before the peak is the rising phase, and after is the
decay phase. The duration of the rising phase and decay phase are
$\tau_{R}=0.4$ and $\tau_{D}=2.0$ seconds for this case.

The source offset measures the vector distance between the observed centroid
and the starting point of the photons in the sky plane. In this case, the
average offset is $0.105R_{s}$ during the rising phase and $0.115R_{s}$ during
the decay phase, the offset in y direction is negligible (less than
$0.002R_{s}$).

The source size is measured by the FWHM in $x$ and $y$ direction:
$\textup{FWHM}_{x,y}=2\sqrt{2\ln{2}}\sigma_{x,y}$, where $\sigma_{x,y}$ is the
variance of the distribution. In this case, the average FWHM in $x$ direction
is $0.513R_{s}$ $(0.137^{\circ})$ for the rising phase and $0.662R_{s}$
$(0.176^{\circ})$ for the decay phase. The average FWHM in $y$ direction is
$0.538R_{s}$ $(0.143^{\circ})$ for the rising phase and $0.680R_{s}$
$(0.181^{\circ})$ for the decay phase. The source FWHM size in y direction is
slightly larger than that in $x$ direction.

The speed and expansion rate of the source in the sky plane can be obtained
from the linearly fit of the variation of the source position and size with
time. From Figure 2, we can see that, the source behaves differently in
movement and expansion during the rising and decay phase, for this case, the
moving speed and the expansion rate are both larger in the rising phase. The
source size has a positive expansion rate in both $x$ and $y$ direction.

In order to study the properties and variations of the source in different
parameter sets, namely the density fluctuation level and the anisotropic
scale, we run the above simulation process in the parameter space of
$\alpha\in[0.05,0.99]$ and $\epsilon\in[0.03,0.45]$. We uniformly select
$36\times 36$ points in the parameter space. For each case of parameter set
($\epsilon,\alpha$), a simulation of $3\times 10^{5}$ photons is conducted. In
each case, the frequency of the wave is assumed to be 35 MHz, the frequency
ratio between the wave and local plasma frequency at the starting point
($f_{0}/f_{pe}$) is 1.1 for fundamental emission and 2.0 for harmonic
emission. The photons are launched at the center of the solar disk
($\theta_{0}=\phi_{0}=0$), for fundamental emission $r_{0}=1.750R_{s}$ and
$h_{0}=860$ km, for harmonic emission $r_{0}=2.104R_{s}$ and $h_{0}=1010$ km.
The angular scattering rate coefficient $\eta$ changes in the range between
$8.9\times 10^{-7}\,\textup{km}^{-1}$ to $2.4\times 10^{-4}\,\textup{km}^{-1}$
for $\epsilon\in[0.03,0.45]$. The statistical result of the source properties
is presented in the following subsections.

###  3.1 Source size and duration

The source size and duration are the most common information that can be
extracted from the imaging and spectroscopy observations of solar radio
bursts. In this subsection, the source size is measured as the intensity
weighted FWHM size for all collected photons in the simulation as done in
Figure 1. The decay time is measured as the time of intensity decreasing by
$1/e$ from the peak, namely, $\tau_{D}$.

The source size and decay time of the fundamental emission
($f_{0}/f_{pe}=1.1$) is shown in Figure 3. Within the parameter space, the
source size varies from 0.2 to 0.95 solar radii and the decay time varies from
about 0 to 7 seconds. The source size is largely determined by the fluctuation
variance value (or the scattering coefficient $\eta$), the tendency of contour
lines in Figure 3a is nearly parallel to y-axes in the range of $\alpha>0.1$.
The decay time is sensitive to both $\alpha$ and $\epsilon$ (or $\eta$). The
decay time increases with the density fluctuation level and the anisotropy
parameter, thus decrease with the anisotropy degree in the density
fluctuation.

Figure 3: The source size and decay time of the simulation results for the
fundamental emission. Purple lines in the left panel show the contour of
source size, and black lines in the right panel show the contour of the decay
time for the simulation data. The bold lines mark the value of
$D_{\rm{FWHM}}=0.675R_{s}$, $\tau_{D}=2.30\,\rm{s}$ for $f_{0}=35$ MHz in
Equation 13 and 14. The green square marks the cross point of the bold contour
lines at $\epsilon=0.277$ and $\alpha=0.719$. The equivalent scale length of
density fluctuation $h_{0}$ is about 860 km near the wave generation site.

For the harmonic emission with the starting frequency of $f_{0}/f_{pe}=2.0$,
as shown in Figure 4 the source size varies from 0.1 to 0.9 $R_{s}$ and the
decay time varies from 0 to 0.7 seconds within the parameter space. Compared
to the fundamental emission with $f_{0}/f_{pe}=1.1$, the source size of
harmonic emission is slightly smaller than the size of fundamental emission
with a value of about 0.1 $R_{s}$, while the decay time of harmonic emission
is significantly smaller than the fundamental emission. Both the source size
and duration are not sensitive to the variation of fluctuation anisotropy for
$\alpha>0.2$.

Figure 4: The source size and decay time of the simulation results for the
harmonic emission, the specs are the same as Figure 3, and $h_{0}=1010$ km.

The simulation results in the parameter space can be used to constrain the
background plasma properties, compared with the corresponding observation.
According to the previous fitted result of multiple observations for the
source size and decay time of type III radio bursts (Kontar et al., 2019), the
relationship between the observed source size [Degree] and frequency [MHz] can
be expressed as:

$D_{\rm{FWHM}}=(11.78\pm 0.06)\times f^{-0.98\pm 0.05}.$ (13)

The decay time [s] and frequency [MHz] can be expressed as:

$\tau=(72.23\pm 0.05)\times f^{-0.97\pm 0.03}.$ (14)

The observed source size and decay time of the 35 MHz source are
$0.18^{\circ}$ (0.675 $R_{s}$) and 2.30 seconds accordingly. Figure 3 shows
the simulation results of the source size and decay time with different
density fluctuation variance and the anisotropy parameter. The parameter set
$\epsilon=0.277$ and $\alpha=0.719$ (marked as a green square in Figure 3) can
satisfy both the observation estimation of size and decay time at 35 MHz. The
corresponding scattering rate coefficient for $\epsilon=0.277$ in the
simulation is $\eta=8.9\times 10^{-5}\,\textup{km}^{-1}$. The source
properties variations of $\epsilon=0.28$ and $\alpha=0.72$ is shown in Figure
1 and 2, note that the starting point of the photons is chosen to be located
at $\theta_{0}=20^{\circ}$ in Figure 1 and 2 for displaying the source
position offset.

###  3.2 Position offset of the source

In Section 3.1, the starting point is set as in the center of the solar disk,
and the centroid of the reconstructed source is located near $(0,0)$. Thus,
the offset of the source from the starting point is close to zero when the
position angle of the starting point $\theta_{0}$ is 0. When the position
angle of the starting point is away from zero, the offset is not negligible.
Considering the sphere-symmetric of the background parameters, the simulated
observation for a given $\theta_{0}$ can be obtained by rotating the result
wave vector ($\bm{k}$) and ($\bm{r}$) of all photons to the reverse angle
direction of $\theta_{0}$ and redo the photon collecting and image
reconstruction process. For simplification without losing generality, we use
$y_{0}=0$ for the starting point so that the simulated offset is mainly in the
$x$ direction. We measure the offset in the $x$ direction as $\Delta
x=x_{c}^{obs}-x_{0}$. The source centroid $x_{c}^{obs}$ is measured as the
weighted average position in $x$ axis for all collected photons as in Figure
1.

Figure 5: The offset of the reconstructed source centroid from the starting
point of the fundamental emission. The left and right panels show the result
of the source locations at $\theta_{0}=30^{\circ}$ and $60^{\circ}$,
respectively. The positive value represents the outward offset from the solar
disk center, the negative value represents the inward offset in this figure.
The equivalent scale length of density fluctuation $h_{0}$ is about 860 km
near the wave generation site.

Figure 6: The offset from the reconstructed source centroid and the starting
point of the photons for harmonic emission. The specs are the same as Figure
5, and $h_{0}=1010$ km.

Different density fluctuation variance and anisotropic scale can result in
different offset amounts. Figure 5 and 6 shows the source offset ($\Delta x$)
of the fundamental ($f_{0}/f_{pe}=1.1$) and harmonic ($f_{0}/f_{pe}=2.0$)
emission in the parameter space of $\epsilon\in[0.03,0.45]$ and
$\alpha\in[0.05,0.99]$. For the fundamental emission, as shown in Figure 5,
the result of the offset value is positive for a major part of the parameter
space, meaning the direction of the offset is outward. In a small portion of
the parameter space at $\epsilon<0.03$ with weak scattering, the offset is
inward. The source offset $\Delta x$ increases with the density fluctuation
variance $\epsilon$ for the two starting position angles
$\theta_{0}=30^{\circ}$ and $\theta_{0}=60^{\circ}$. The offset of the source
with a starting position angle of $60^{\circ}$ is larger than that of
$30^{\circ}$. The maximum offset of $\theta_{0}=30^{\circ}$ is about
$0.28\,R_{s}$ and the maximum offset of $\theta_{0}=60^{\circ}$ is about
$0.43\,R_{s}$.

The source offset of harmonic emission is largely different from the
fundamental emission, as shown in Figure 6. Considerable portions of the
results in parameter space have inward offset with $\Delta x<0$. The outward
offset can exist when the density fluctuation level is high and the anisotropy
degree is low, for example, $\epsilon>0.3$ and $\alpha>0.3$. The offset value
of $\theta_{0}=30^{\circ}$ varies in $(-0.07,0.11)$ within the parameter
space. The offset value of $\theta_{0}=60^{\circ}$ varies in $(-0.10,0.17)$
within the parameter space. Comparing the offset of the fundamental and
harmonic. The offset value of the harmonic emission is much smaller than that
of the fundamental emission.

Figure 7: The offset ratio ($\Delta x/x_{0}$) of the source centroid from the
starting point of the photons with different starting position angle. The
simulation parameter sets are labeled on the top of each panel.

For further study of the source offset, we selected four cases with the
parameter of $\epsilon=0.102,0.354$, and $\alpha=0.157,0.802$. These four
cases represent the relatively high and low value of angular scattering rate
and a relatively large and small degree of anisotropy. Figure 7 shows the
offset ratio ($\Delta x/x_{0}$) of the source centroid to the starting point.
The four panels of Figure 7 show that for all these four parameter sets, the
offset of fundamental emission is more outward than the harmonic emission. For
relatively small anisotropy parameter $\alpha$ (high degree of anisotropy),
the relative offset for both fundamental and harmonic emission decrease with
the position angle $\theta_{0}$. For relatively larger anisotropy parameter
$\alpha$ (low degree of anisotropy), the relative offset is stable when
$\sin{\theta_{0}}<0.8$, and increase with position angle when
$\sin{\theta_{0}}>0.8$. Comparing the upper and lower column of Figure 7, one
can find that the value of relative offset is larger when the density
fluctuation (or angular scattering rate) is larger.

Assuming the wave of fundamental and harmonic emission with the same given
frequency is generated at the same angular direction from the solar center.
The visual distance between the starting points of the fundamental and the
harmonic is $\Delta R\sin\theta_{0}$ in the sky plane, where $\Delta R$ is the
height difference between the starting points of the fundamental and harmonic
wave. The source offset due to the refraction and scattering will cause the
observed visual distance to deviate from $\Delta R\sin\theta_{0}$. In this
work, for the wave of 35 MHz, the fundamental wave is generated at
$R_{F}=1.750R_{s}$ by assuming $f_{0}/f_{pe}=1.1$, the harmonic wave is
generated at $R_{H}=2.104R_{s}$ by assuming $f_{0}/f_{pe}=2.0$, $\Delta R$ for
this case is $0.354R_{s}$. Figure 8 shows the distance between the
reconstructed source centroids of fundamental and harmonic for 35 MHz. From
Figure 8 we can see that the value of the distance between the reconstructed
sources is significantly smaller than the visual distance of the starting
points of the two waves. In other words, the apparent sources of the
fundamental and harmonic emission will be much closer with each other than
their true source in the sky plane.

Figure 8: The distance of the reconstructed source centroid of fundamental and
harmonic emission (red lines), compared with the geometry distance of the
starting point in the sky plane (black lines).

The influence of the background medium to the visual source position can be
divided into two parts, namely the refraction and scattering, marked as ‘refr’
and ‘scat’ in the following discussion. With the dispersion relation of
unmagnetized cold plasma described in Equation (3), the motion direction of
the photon is aligned with the wave vector as shown in Equation 1. Thus, the
variation of the wave vector leads to the bending of the ray-path of the
photon, and eventually causes the offset of the visual position of the source
centroid. The variation of the wave vector in Equation 2 can be split into
these two parts accordingly:

$\displaystyle\frac{\mathrm{d}k_{i,\mathrm{refr}}}{\mathrm{d}t}$
$\displaystyle=-\frac{\partial\omega}{\partial r_{i}},$ (15)
$\displaystyle\frac{\mathrm{d}k_{i,\mathrm{scat}}}{\mathrm{d}t}$
$\displaystyle=\frac{\partial D_{ij}}{\partial k_{i}}+B_{ij}\xi_{j}.$ (16)

In the ray tracing process, we can measure the cumulative change of the wave
vector $k_{i}$ due to the refraction and scattering by individually
integrating the variation of wave vector due to these two factors for each
photon:

$\displaystyle\Delta K_{i,\mathrm{refr}}$
$\displaystyle=\sum\mathrm{\delta}k_{i,\mathrm{refr}},$ (17)
$\displaystyle\Delta K_{i,\mathrm{scat}}$
$\displaystyle=\sum\mathrm{\delta}k_{i,\mathrm{scat}}.$ (18)

The average change of the wave vector for all photons collected near the
observation site can provide a qualitative estimation of the relative
influence of scattering and refraction effects on the apparent source
position. In the simulation, the reconstructed image shows that the source
offset is mostly in the $x$ direction in the sky plane. Thus, we use $\Delta
K_{x}^{l}$ to represent the total bending influence of the observed source for
the $l\textup{th}$ photon. For the collected photons in each viewpoint, we
calculate the average value $\overline{\Delta K_{x}}=\Sigma_{l=1}^{l=N}\Delta
K_{x}^{l}/N$ as a measure of the influence of the propagation effects on the
shift of observed source position. Figure 9 shows the variation of
$\overline{\Delta K_{x}}$ with the position angle of the starting point, where
the solid lines and dash lines represent the effects of refraction and
scattering, respectively.

Figure 9: The average value $\overline{\Delta K_{x}}$ of all collected photons
due to influences of scattering (dash lines) and refraction (solid lines) for
the waves of fundamental and harmonic emission, respectively.

From Figure 9 one can see that, the contribution of the refraction is in the
positive-$x$ direction for the wave vector, while the contribution of
scattering is in the negative-$x$ direction. The influences of both refraction
and scattering effect are stronger for the fundamental wave than for the
harmonic wave in all four cases. When the scattering is more isotropic
($\alpha=0.802$ shown as cyan lines in the right panel), the influence of
scattering is overwhelmingly larger than the refraction effect for the
fundamental wave, while the influences on the harmonic wave are insignificant.
For high anisotropic scattering cases with small $\alpha$ (left panels in
Figure 9), the value of $\Delta K_{x}$ for refraction is comparable to the
scattering.

It needs to note that projecting the wave vector back to the sky plane, a
positive change of the wave vector in the $x$ direction indicates a negative
offset of the source position, and vice versa. Comparing Figure 8 and 9, one
can get that a larger angular scattering rate and more isotropic scattering
would cause the apparent positions of the fundamental and harmonic emission
being closer to each other.

###  3.3 Source expansion rate and the source speed in the sky plane

As the evolution of radio observation technology, the radio imaging can be
carried out with high time resolution, enabling the study of radio source
variation within a solar radio burst, namely the visual speed and expansion
rate of the source.

For the simulation, the temporal variation of the source properties can be
obtained from the series of the reconstructed source frame sectioned according
to the arrival time of the photons, as shown in Figure 2. For each parameter
set ($\epsilon$, $\alpha$), we linearly fit the time-distance profile of the
source centroid for the frames within the flux FWHM in the $x$ direction, the
resulting slope of the linear-fit yields the visual speed ($V_{x}$). The
$V_{x}$ is calculated in the parameter space of $\epsilon\in[0.03,0.45]$ and
$\alpha\in[0.05,0.99]$ with two given starting position angle of
$\theta_{0}=30^{\circ}$ and $\theta_{0}=60^{\circ}$. Accordingly, the angular
scattering rate coefficient varies in the space of $\eta\in[8.9\times
10^{-7},2.4\times 10^{-4}]\,\textup{km}^{-1}$ for $\epsilon\in[0.03,0.45]$.

As shown in Figure 4 and 6, the duration broadening and position offset due to
the wave propagation effects are small for harmonic emission. Therefore, there
is large uncertainty in linearly-fitting the source variation trend of
harmonic emission. Observation also shows that at a given frequency, the
source of the harmonic wave is relatively stable compared with the fundamental
wave (Zhang et al., 2020). In the following only variation of the fundamental
emission source is presented.

Figure 10: The visual speed of the source for fundamental emission. The left
and right panels show the result of the source locations at
$\theta_{0}=30^{\circ}$ and $\theta_{0}=60^{\circ}$. The positive value
represents the outward motion from the solar disc center, while the negative
value represents the inward motion in this figure. The scale length of density
fluctuation $h_{0}$ is about 860 km near the wave generation site.

Figure 10 shows the visual speed of the source of the fundamental emission.
The scattering makes the apparent source moving outward from the Sun. The
patterns of the visual speed are similar for $30^{\circ}$ and $60^{\circ}$,
and the major part of the parameter space is divided into two parts. The
source tends to have a small visual speed for the more isotropic background
($\alpha>0.5$) and large visual speed requires the background to be highly
anisotropic ($\alpha<0.4$). The visual speed reach the maximum at about
$\alpha=0.2$ and $\epsilon=0.45$ (corresponds to $\eta=2.4\times
10^{-4}\,\textup{km}^{-1}$) within the parameter space for both $30^{\circ}$
and $60^{\circ}$ of starting position angle. The maximum visual speed is about
$0.65\,R_{s}/\textup{s}$ or 1.5 $c$ for $\theta_{0}=30^{\circ}$ and about
$0.5\,R_{s}/\textup{s}$ or 1.2 $c$ for $\theta_{0}=60^{\circ}$. The visual
speed of the source with $30^{\circ}$ starting position angle is larger that
of $60^{\circ}$. Comparing Figure 5 and Figure 10, one can find that direction
of offset is mostly identical to the direction of visual speed within the
parameter space.

Figure 11: The expansion rate of the source size for fundamental emission
located at the solar disk center, measured with the increasing rate of the
FWHM width of the visual source.

For calculating the expansion rate of source size, we use the cases that the
original source is located at the solar disk center. The expansion rate is
measured as the slope of linear-fitting the FWHM width in $x$ axis for the
flux intensity frames, within the time period defined by the FWHM range in the
time-intensity profile. Figure 11 shows the expansion rate. The large
expansion rate mainly gathers in the high anisotropic part ($\alpha<0.4$)
within the parameter space. The expansion rate can reach $1.4R_{s}/\textup{s}$
or $22\,\textup{Arcmin/s}$. Comparing Figure 10 and 11, one can find that the
high expansion rate of the source size is often associated with the large
speed of the source apparent motion.

##  4 Discussion

In this study, the source in the calculation is assumed as point pulse source,
while the real source could have finite size and duration. Assuming there are
no interaction between the electromagnetic waves generated at different time
and position, the observed source can be expressed as the convolution of the
point pulse response and the time spatial distribution of the real source. The
beam electron generation process at energy release site is considered to be
fragmented (Benz, 1994; Bastian et al., 1998), the intrinsic radio sources can
have complex spatial and temporal structures. One should be cautious when the
simulation results presented in Section 3 are directly compared with the
observation results. For a burst element of the fragments like a type III
burst in a given frequency, the decay time is determined by the duration of
the local excitation decay phase, the velocity dispersion of beam electrons,
and the propagation effects (Li et al., 2008; Ratcliffe et al., 2014; Zhang et
al., 2019). The observed size of the source is determined by the real size of
the wave excitation regime and the broadening due to the propagation effect
(Kontar et al., 2017). Thus, the anisotropic degree and angular scattering
rate estimated from the simulation results shown in Figure 3 using Equation
(13) and (14) are the upper limits. The radio burst with fine structures and
the short term narrowband radio bursts can provide more constraint on the
parameters of the background. For example, the wave generation site of the
type IIIb radio burst is believed to be in a compact regime. The apparent
source size may be determined mainly by the propagation effects. However, the
observed source size varies largely from case to case, indicating that the
variation of the corona density fluctuations.

The density fluctuation variance and its length scale determine the angular
scattering rate of the radio wave. In the simulation of this work, the tuning
parameter for the angular scattering rate is $\epsilon$. The length scale is
described by the inner scale ($l_{i}$) and outer scale ($l_{o}$) of the
fluctuation spectrum. Both $l_{i}$ and $l_{o}$ are obtained from empirical
models (Coles & Harmon, 1989; Wohlmuth et al., 2001). The inner scale is
considered to be the inertial length of ion $d=v_{A}/\Omega_{i}$ (Coles &
Harmon, 1989; Spangler & Gwinn, 1990), where $v_{A}$ is the local Alfven
speed, $\Omega_{i}$ is the gyro-frequency of ion. The outer scale represents
the scale of energy containing in fluctuation. The density fluctuation
spectrum and cutoff scales in the corona may be largely different in active
region and quiet region. Observations from multiple methods like the
interplanetary scintillation observations (Chang et al., 2016), remote and in
situ observations can help constrain $l_{i},l_{o}$ and $\epsilon$, as well as
the angular scattering rate, and benefit the simulation model for the radio
wave ray-tracing. Also the radio burst observation combined with simulation
results can help to derive the density fluctuation property of the corona
(Krupar et al., 2020).

For a given frequency, the wave generation position of fundamental and
harmonic emission is at different height according to the plasma emission.
While observational results indicate that the apparent source of fundamental
and harmonic emission have co-spatial relationship for some type III radio
bursts (Dulk & Suzuki, 1980). One possible interpretation on this is that the
waves are generated and confined in an density depleted tube (Duncan, 1979; Wu
et al., 2002). Alternatively, this may be caused by the refraction and
scattering of the waves in their way propagating from the exciting site to the
observer (Steinberg et al., 1971). In this work, we inspect the source offset
and the relative position of the fundamental and harmonic emission of 35 MHz.
The results show that the scattering effect can cause the centroid of the
fundamental emission close to that of the harmonic emission. For the four
parameter sets shown in Figure 8, the visual distance of the reconstructed
sources is closer in the case of larger density fluctuation scale. For the
case of $\epsilon=0.354$ and $\alpha=0.157$, the position difference
$x_{H}-x_{F}$ is negative, meaning the source of fundamental emission can
appear even slightly higher than the harmonic emission.

The position offset due to the scattering and refraction is not negligible for
metric radio bursts generated far away from the solar disk center. The
relative error caused by neglecting the wave propagation effects could be up
to 50% for the fundamental and 10% for the harmonic as shown in Figure 7. The
correction of wave propagation effects is essential for the study concerning
the position of the source centroid, especially for the fundamental emission.

From the LOFAR high-cadence imaging spectroscopy, it is found that the
apparent sources of a solar type IIIb striae move outward quickly from the Sun
with a speed varying from 0.25 to 4 times the speed of light (Kontar et al.,
2017; Zhang et al., 2020). Meanwhile the observed expansion rate of the source
size can be as high as $382\pm 30\,\textup{arcmin}^{2}\textup{s}^{-1}$. A
commonly held belief in plasma emission is that type IIIb bursts are generated
in a source region with high density inhomogeneity. The simulation results in
Figure 10 and 11 show that both the visual speed and the expansion rate can be
very large, when the fluctuation variance $\epsilon$ is large and the
anisotropy parameter $\alpha$ is small. This is in consistent with the
observations and indicates that the density fluctuation may be highly
anisotropic in the source region of type IIIb bursts.

##  5 Summary

In this paper, we performed ray-tracing simulation of wave transport for point
pulse source. For the first time, we explored the parameter space of the
scattering rate coefficient $\eta$ (represented by density fluctuation level
$\epsilon$) and the anisotropy parameter $\alpha$ with massive number of ray
tracing simulation. We analyzed the simulation results to study the influences
of wave propagation on the observed source size, position, expansion rate,
visual speed, and duration of solar radio bursts for both fundamental emission
and harmonic emission. The following are the major conclusions:

  * •

For fundamental emission, both the source size and decay time increase with
the scattering rate coefficient or the density fluctuation variance. The
isotropic fluctuation can produce much larger decay time than a highly
anisotropic fluctuation, while the source size is not sensitive to the level
of anisotropy. For harmonic emission, both the source size and decay time are
largely determined by the density variance. The decay time of harmonic
emission is significantly smaller than that of fundamental emission for the
same background parameters.

  * •

By comparing the source size and decay time derived from simulation results to
the observational statistics, we obtained the estimation of $\eta=8.9\times
10^{-5}\,\textup{km}^{-1}$ and $\alpha=0.719$ near the source region of
fundamental emission at frequency 35 MHz.

  * •

We derived the position offset of the source and analyzed its dominant factor
by decoupling the scattering and refraction in the wave propagation. The
statistical results of source offset show that, the source of the fundamental
emission is more outward shifted than that of harmonic emission by the
propagation effect, which could account for the co-spatial of fundamental and
harmonic emission in some observations.

  * •

The observed source position and size can have significant visual motion and
expansion due to the wave propagation effects. Both the visual speed and the
expansion rate tend to be large for highly anisotropic medium. For fundamental
emission, the visual speed and the expansion rate for FWHM of the source size
can reach respectively about 1.5 $c$ and 22 Arcmin/s at $\eta=2.4\times
10^{-4}\,\mathrm{km^{-1}}$ and $\alpha=0.2$.

A comprehensive comparison of the observed source characteristics and their
corresponding values simulated in the parameter space can help us precisely
diagnose plasma properties near the wave exciting site of the radio bursts.

In this work, the simulation considers the radio emission of a single
frequency (35 MHz) with the frequency ratio of $f_{0}/f_{pe}=1.1$ for
fundamental emission and $f_{0}/f_{pe}=2.0$ for harmonic emission. The future
work of simulation concerning the parameter space of frequency and frequency-
ratio of wave excitation can help to understand the dynamic spectrum of a
complete radio burst event and its beam electron property. The background
density model is spherically symmetric in the present work, while the solar
corona is a non-uniform medium with a number of large structures. More
specific density models with helmet streamers, coronal holes, under-dense or
over-dense flux tubes may be also studied in future simulations.

##  6 Acknowledgements

The research was supported by the National Nature Science Foundation of China
(41974199 and 41574167) and the B-type Strategic Priority Program of the
Chinese Academy of Sciences (XDB41000000). The numerical calculations in this
paper have been done on the supercomputing system in the Supercomputing Center
of University of Science and Technology of China.

## References

  * Arzner & Magun (1999) Arzner, K., & Magun, A. 1999, Astronomy and Astrophysics, 351, 1165
  * Bastian et al. (1998) Bastian, T., Benz, A., & Gary, D. 1998, Annual Review of Astronomy and Astrophysics, 36, 131
  * Benz (1994) Benz, A. 1994, Space science reviews, 68, 135
  * Bian et al. (2019) Bian, N., Emslie, A., & Kontar, E. 2019, The Astrophysical Journal, 873, 33
  * Bisoi et al. (2018) Bisoi, S. K., Sawant, H., Janardhan, P., et al. 2018, The Astrophysical Journal, 862, 65
  * Bougeret & Steinberg (1977) Bougeret, J., & Steinberg, J. 1977, Astronomy and Astrophysics, 61, 777
  * Chang et al. (2016) Chang, O., Gonzalez-Esparza, J., & Mejia-Ambriz, J. 2016, Advances in Space Research, 57, 1307
  * Chen et al. (2015) Chen, B., Bastian, T. S., Shen, C., et al. 2015, Science, 350, 1238, doi: 10.1126/science.aac8467
  * Chen et al. (2014) Chen, Y., Du, G., Feng, L., et al. 2014, The Astrophysical Journal, 787, 59
  * Chrysaphi et al. (2018) Chrysaphi, N., Kontar, E. P., Holman, G. D., & Temmer, M. 2018, The Astrophysical Journal, 868, 79
  * Coles & Harmon (1989) Coles, W., & Harmon, J. 1989, The Astrophysical Journal, 337, 1023
  * Dulk & Suzuki (1980) Dulk, G., & Suzuki, S. 1980, Astronomy and Astrophysics, 88, 203
  * Duncan (1979) Duncan, R. 1979, Solar Physics, 63, 389
  * Fokker (1965) Fokker, A. 1965, Bulletin of the Astronomical Institutes of the Netherlands, 18, 111
  * Kontar et al. (2017) Kontar, E., Yu, S., Kuznetsov, A., et al. 2017, Nature communications, 8, 1515
  * Kontar et al. (2019) Kontar, E. P., Chen, X., Chrysaphi, N., et al. 2019, The Astrophysical Journal, 884, 122
  * Krupar et al. (2018) Krupar, V., Maksimovic, M., Kontar, E. P., et al. 2018, ApJ, 857, 82, doi: 10.3847/1538-4357/aab60f
  * Krupar et al. (2020) Krupar, V., Szabo, A., Maksimovic, M., et al. 2020, The Astrophysical Journal Supplement Series, 246, 57
  * Kuznetsov et al. (2020) Kuznetsov, A. A., Chrysaphi, N., Kontar, E. P., & Motorina, G. 2020, The Astrophysical Journal, 898, 94
  * Kuznetsov & Kontar (2019) Kuznetsov, A. A., & Kontar, E. 2019, Astronomy & Astrophysics, 631, L7
  * Li et al. (2008) Li, B., Cairns, I. H., & Robinson, P. A. 2008, Journal of Geophysical Research: Space Physics, 113
  * Mann et al. (2018) Mann, G., Breitling, F., Vocks, C., et al. 2018, Astronomy & Astrophysics, 611, A57
  * McCauley et al. (2018) McCauley, P. I., Cairns, I. H., & Morgan, J. 2018, Solar Physics, 293, 132
  * Morosan et al. (2019) Morosan, D. E., Carley, E. P., Hayes, L. A., et al. 2019, Nature Astronomy, 3, 452
  * Ratcliffe et al. (2014) Ratcliffe, H., Kontar, E., & Reid, H. 2014, Astronomy & Astrophysics, 572, A111
  * Riddle (1974) Riddle, A. 1974, Solar Physics, 35, 153
  * Robinson (1983) Robinson, R. 1983, in Proceedings of the Astronomical Society of Australia, Vol. 5, 208–211
  * Spangler & Gwinn (1990) Spangler, S. R., & Gwinn, C. R. 1990, The Astrophysical Journal, 353, L29
  * Steinberg et al. (1971) Steinberg, J., Aubier-Giraud, M., Leblanc, Y., & Boischot, A. 1971, Astronomy and Astrophysics, 10, 362
  * Stewart (1972) Stewart, R. 1972, Publications of the Astronomical Society of Australia, 2, 100
  * Stewart (1976) —. 1976, Solar Physics, 50, 437
  * Wild et al. (1959) Wild, J., Sheridan, K., & Trent, G. 1959, in Symposium-International Astronomical Union, Vol. 9, Cambridge University Press, 176–185
  * Wohlmuth et al. (2001) Wohlmuth, R., Plettemeier, D., Edenhofer, P., et al. 2001, Space Science Reviews, 97, 9
  * Wu et al. (2002) Wu, C. S., Wang, C. B., Yoon, P. H., Zheng, H. N., & Wang, S. 2002, The Astrophysical Journal, 575, 1094, doi: 10.1086/341468
  * Yu et al. (2020) Yu, S., Chen, B., Reeves, K. K., et al. 2020, ApJ, 900, 17, doi: 10.3847/1538-4357/aba8a6
  * Zhang et al. (2019) Zhang, P., Yu, S., Kontar, E. P., & Wang, C. B. 2019, The Astrophysical Journal, 885, 140
  * Zhang et al. (2020) Zhang, P., Zucca, P., Sridhar, S. S., et al. 2020, Astronomy & Astrophysics, 639, A115


# Fusion-Based Multi-User Semantic Communications for Wireless Image
Transmission over Degraded Broadcast Channels

Tong Wu, Zhiyong Chen, Meixia Tao, Bin Xia, Wenjun Zhang
Cooperative Medianet Innovation Center, Shanghai Jiao Tong University,
Shanghai, China
Email: {wu_tong, zhiyongchen, mxtao, bxia<EMAIL_ADDRESS>
###### Abstract

Degraded broadcast channels (DBC) are a typical multi-user communications
scenario. There exist classic transmission methods, such as superposition
coding with successive interference cancellation, to achieve the DBC capacity
region. However, semantic communications method over DBC remains lack of in-
depth research. To address this, we design a fusion-based multi-user semantic
communications system for wireless image transmission over DBC in this paper.
The proposed architecture supports a transmitter extracting semantic features
for two users separately, and learns to dynamically fuse these semantic
features into a joint latent representation for broadcasting. The key here is
to design a flexible image semantic fusion (FISF) module to fuse the semantic
features of two users, and to use a multi-layer perceptron (MLP) based neural
network to adjust the weights of different user semantic features for flexible
adaptability to different users channels. Experiments present the semantic
performance region based on the peak signal-to-noise ratio (PSNR) of both
users, and show that the proposed system dominates the traditional methods.

##  I Introduction

In recent years, semantic communications have received significant attention
from both industry and academia. With the help of artificial intelligence
(AI), semantic communications can extract the semantic information from the
original data, and further transmit it, thereby significantly improving
communication efficiency [1]. Therefore, semantic communications have been
considered a promising solution for the sixth-generation (6G) wireless
networks [2].

Several studies have been conducted on semantic communications for different
types of original information, such as text [3, 4], image [5, 6], and video
[7, 8, 9]. For text transmission, a deep learning-based semantic communication
system is proposed in [3], named DeepSC, which has an advantage in the low
signal-to-noise ratio (SNR) regime. For image transmission, a deep learning-
based semantic image coding method is designed in [5] to encode images beyond
pixel level. For video transmission, the end-to-end joint source-channel
coding (JSCC) video transmission scheme is first proposed in [7]. Then, [8]
designs a novel deep joint source-channel coding approach to achieve wireless
video transmission, which can outperform traditional wireless video coded
transmission schemes.

It is worth noting that previous works mainly focus on point-to-point semantic
communications, while research on multiuser semantic communications is
relatively limited. In [10], a heterogeneous semantic and bit communication
framework is designed for multiple access channels that utilizes a method
called semi-nonorthogonal multiple access (NOMA) and achieves better
performance than the classic NOMA system. Meanwhile, a novel joint image
compression and transmission scheme for the multi-user uplink scenario is
presented in [11], which utilizes NOMA and incorporates deep neural networks
(DNNs) into the transmitters. For broadcasting channels, a one-to-many scheme
is proposed in [12] for text transmission, where the transmitter concatenates
these texts together and extracts their semantic features for transmission.
For relay channels, a semantic-and-forward scheme is first designed in [13] to
address the heterogeneous background knowledge problem. Then, a novel deep
joint source-channel coding scheme for image transmission over a half-duplex
cooperative relay channel is presented in [14].

Actually, multi-user semantic communications are not simply point-to-point
semantic communications but require corresponding design for multi-user
channels. Motivated by this, we consider a degraded broadcast channel (DBC) in
this paper, which is a typical multiuser communication scenario. There is a
transmitter and multiple users located in different geographical locations.
The capacity region of DBC is well-known, and there are many traditional
transmission methods, such as superposition coding with successive
interference cancellation, time division (TD) and frequency division (FD), to
achieve the DBC capacity region. However, in semantic communications, can the
semantic information of two users be deeply integrated through AI networks
rather than simply transmitted together using superposition coding?

To address this issue, we propose a fusion-based multi-user semantic
communications system for wireless image transmission over two-user DBC. In
the proposed architecture, a transmitter can extract and fuse the semantic
features as a joint latent representation of both users. The worse user can
only decode its own image from the joint representation, while the better user
reconstructs the other image first and then obtains its own image based on the
reconstructed image. To deal with the fuse of semantic features for two users,
we design a flexible image semantic fusion (FISF) scheme to dynamically
control the weight of two users’ semantic features in the joint latent
representation by using a neural network based on multi-layer perceptron
(MLP). Meanwhile, to adapt the respective semantic features to the respective
channels with different SNRs, the proposed FISF scheme uses the attention
mechanism with channel state information (CSI) to adapt the different channel
condition [15, 16]. Numerical results based on real-world datasets show that
the proposed system can significantly improve the peak signal-to-noise ratio
(PSNR) of the images for both users.

##  II System Model and Strategy Design

In this section, we propose a semantic communication strategy for wireless
image transmission over the degraded broadcast channel. The system consists of
a transmitter and two users, where two image messages $\mathbf{s}_{1}$ and
$\mathbf{s}_{2}$ are expected to be delivered to the two users through
semantic communication.

Figure 1: The structure of the proposed degraded broadcast semantic
communication system

###  II-A System Overview

As shown in Fig. 1, two semantic encoders (SE) $f_{e}(\cdot;{\theta_{1}})$ and
$f_{e}(\cdot;{\theta_{2}})$ can extract the image features
$\mathbf{x}_{1}^{e}$ and $\mathbf{x}_{2}^{e}\in{\mathbb{R}^{l}}$ from the two
source images $\mathbf{s}_{1}$ and $\mathbf{s}_{2}\in{\mathbb{R}^{h\times
n\times 3}}$ respectively, where $l$ is the output dimension of the semantic
encoders, $h$ and $n$ denote the height and width of the image, and 3 is the
color channels $R$, $G$ and $B$. $f_{e}$ is the module strcuture and
$\theta_{i}$, $i=1,2$, is learning parameter. $\mathbf{x}_{1}^{e}$ and
$\mathbf{x}_{2}^{e}$ are then fed into a mutual attention (MA) module
$f_{ma}(\cdot;\zeta)$, which first computes a element-wise weight
$\mathbf{w}\in\mathbb{R}^{l}$. The relationship between $\mathbf{x}_{1}^{e}$
and $\mathbf{x}_{2}^{e}$ is associated with $\mathbf{w}$. Then, the input
features can be fused with $\mathbf{w}$ in the element-wise product as

$\displaystyle\mathbf{x}_{1}^{ma}=\mathbf{x}_{1}^{e}\odot\mathbf{w},$ (1)
$\displaystyle\mathbf{x}_{2}^{ma}=\mathbf{x}_{2}^{e}\odot\mathbf{(1-w)},$ (2)

where $\mathbf{x}_{1}^{ma}$ and $\mathbf{x}_{2}^{ma}$ are the outputs of the
$f_{ma}(\cdot;\zeta)$. We then design FISF module
$f_{if}(\cdot;{\varphi_{1}})$ to fuse
$\mathbf{x}_{1}^{ma},\mathbf{x}_{2}^{ma}\in\mathbb{R}^{l}$ into a joint latent
representation $\mathbf{y}\in\mathbb{R}^{k}$ by using a fusion ratio $\alpha$
for controlling the reconstruction quality of two users, where $k$ is the
number of channel uses. The transmitter also inform FISF with the channel
state information and thus the output $\mathbf{y}$ can fit the degrade
broadcast channel.

We consider two distant users, one with Gaussian noise power $\sigma_{1}^{2}$
and the other with Gaussian noise power $\sigma_{2}^{2}$. Without loss of
generality, we assume that $\sigma_{1}^{2}>\sigma_{2}^{2}$. The received
signals of the two users are $\mathbf{y}_{1}=\mathbf{y}+\mathbf{n}_{1}$ and
$\mathbf{y}_{2}=\mathbf{y}+\mathbf{n}_{2}$ respectively, where
$\mathbf{n}_{1}$ and $\mathbf{n}_{2}$ are Gaussian noise with noise power
$\sigma_{1}^{2}$ and $\sigma_{2}^{2}$, respectively. We call the user with
noise power $\sigma_{1}^{2}$ as the worse user and the other as the better
user. Similar to the traditional DBC, the worse user can only decode its own
message. Therefore, at the worse user, $\mathbf{y}_{1}$ is fed into the de-
fusion(DF) module $g_{df}(\cdot;{\varphi_{2}})$ with $\alpha$ and CSI,
yielding the output $\mathbf{y}_{1}^{df}\in\mathbb{R}^{k}$. The user then
performs a semantic decoder (SD) $g_{d}(\cdot;{\theta_{2}})$ to reconstruct
the image $\mathbf{s}_{1}$ as $\mathbf{z}_{1}\in\mathbb{R}^{h\times n\times
3}$.

The better user first reconstructs the image $\mathbf{s}_{1}$ and then
reconstructs its own image $\mathbf{s}_{2}$ based on $\mathbf{s}_{1}$.
Specifically, upon receiving $\mathbf{y}_{2}$, the user can reconstruct the
image $\mathbf{s}_{1}$ as $\tilde{\mathbf{z}}_{1}$ by performing the DF module
$g_{df}(\cdot;{\varphi_{3}})$, the SD module $g_{d}(\cdot;\theta_{4})$ and the
SE module $g_{e}(\cdot;\theta_{5})$. Therefore, the user can obtain the
features of $\tilde{\mathbf{z}}_{1}$ as
$\tilde{\mathbf{x}}_{1}\in\mathbb{R}^{l}$. Meanwhile, with
$g_{df}(\cdot;{\varphi_{4}})$ and $\tilde{\mathbf{x}}_{1}$, the user can
reconstruct the image $\mathbf{s}_{2}$ as
$\mathbf{z}_{2}\in\mathbb{R}^{h\times n\times 3}$ by performing the SD module
$g_{d}(\cdot;\theta_{6})$.

Figure 2: The structure of flexible image semantic fusion

We note here that neural networks are utilized for the SE module, the MA
module and the SD in this paper. In the following, we detail the design of the
FISF module.

###  II-B Flexible Image Semantic Fusion Module

For DBC, the channel input $\mathbf{y}$ is a joint latent represenatation of
$\mathbf{s}_{1}$ and $\mathbf{s}_{2}$. The component of $\mathbf{y}$ from
$\mathbf{s}_{1}$ is required to fit the worse channel because the worse user
only needs to reconstruct $\mathbf{s}_{1}$ from the received signal
$\mathbf{y}_{1}$. Meanwhile, all the componets of $\mathbf{y}$ is required to
fit the better channel because the better user requires to reconstruct
$\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ from $\mathbf{y}_{2}$. Therefore, how to
flexibly fuse $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ to $\mathbf{y}$ for
adapting to both channel states is crucial to the system design.

In this paper, we develop a flexible image semantic fusion strategy to fuse
$\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ to flexibly adapt the different channel
conditions and dynamic control the weight of $\mathbf{s}_{1}$ and
$\mathbf{s}_{2}$ in $\mathbf{y}$. The architecture of the proposed FISF module
is shown in Fig. 2. First, the output $\mathbf{x}_{1}^{ma}$ of MA is combined
with the SNR of the worse channel and fed into an MLP to generate a vector
$\mathbf{c}_{1}\in\mathbb{R}^{l}$. Likewise, we can obtain a vector
$\mathbf{c}_{2}\in\mathbb{R}^{l}$ from the output $\mathbf{x}_{2}^{ma}$ of MA
for the better channel. As a result, $\mathbf{c}_{1}$ and $\mathbf{c}_{2}$ are
the attention masks that contain the image feature and the channel station
information. We can adjust the image feature by scaling the attention mask to
transmit the image feature in a more robust form in the channel as following

$\displaystyle\mathbf{x}_{1}^{ca}=\mathbf{x}_{1}^{ma}\odot\mathbf{c}_{1},\
\mathbf{x}_{2}^{ca}=\mathbf{x}_{2}^{ma}\odot\mathbf{c}_{2}.$ (3)

Next, to dynamic control the weight of $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$
in $\mathbf{y}$, $\mathbf{x}_{1}^{ca}$ is combined with a fusion ratio
$\alpha$ and then fed into full-connected (FC) layers. The FC can sort the
semantic vector and put the important parts in the front of the vector. The
output of the FC is multiplied by a non-square identity matrix
$\mathbf{D}_{1}=[I_{\lfloor 2\alpha k\rfloor},O]\in\mathbb{R}^{\lfloor\alpha
l\rfloor\times l}$ to produce $\mathbf{x}_{1}^{rc}\in\mathbb{R}^{\alpha l}$,
where $\lfloor\cdot\rfloor$ is round-down function. Here, $I_{j}$ is
j-dimension identity matrix and $O$ is zero matrix. Similarly, we have
$\mathbf{x}_{2}^{rc}\in\mathbb{R}^{(1-\alpha)l}$ based on
$\mathbf{D}_{2}=[I_{\lceil
2(1-\alpha)k\rceil},O]\in\mathbb{R}^{\lceil(1-\alpha)l\rceil\times l}$, where
$\lceil\cdot\rceil$ is round-up function. This process is described as

$\displaystyle\mathbf{x}_{1}^{rc}=\mathbf{D}_{1}*(\mathbf{A}_{1}[(\mathbf{x}_{1}^{ca})^{T},\alpha]+\mathbf{b}_{1}),$
(4)
$\displaystyle\mathbf{x}_{2}^{rc}=\mathbf{D}_{2}*(\mathbf{A}_{2}[(\mathbf{x}_{2}^{ca})^{T},1-\alpha]+\mathbf{b}_{2}),$
(5)

where $\mathbf{A}_{i}$ and $\mathbf{b}_{i}$ are affine function parameter and
their bias of FCs, respectively. Then, $\mathbf{x}_{1}^{rc}$ and
$\mathbf{x}_{2}^{rc}$ are passed through a MLP to generate
$\mathbf{y}\in\mathbb{R}^{k}$. According to the definition of MLP, we have

$\mathbf{y}=\tanh(\mathbf{Q}*[(\mathbf{x}_{1}^{rc})^{T},(\mathbf{x}_{2}^{rc})^{T}]^{T}+\mathbf{d}),$
(6)

where $\mathbf{Q}\in\mathbb{R}^{k\times l}$ and $\mathbf{d}\in\mathbb{R}^{k}$
are learning parameters. Here, we use tanh activation as the activation
function for this MLP. Based on (6), the $i$-th transmitted symbol can be
expressed as

$\displaystyle y_{i}=\tanh(\sum_{j=1}^{\lfloor\alpha
l\rfloor}q_{ij}{x_{1,j}^{rc}}+\sum_{j=\lfloor\alpha
l\rfloor+1}^{l}q_{ij}{x_{2,(j-\lfloor\alpha l\rfloor)}^{rc}}+d_{i}).$ (7)

Finally, we can perform the power normalization for $\mathbf{y}$ and deliver
it over the channel.

######  Remark 1.

Equ. (7) reveals that we can dynamic control the weight of $\mathbf{s}_{1}$
(or $\mathbf{s}_{2}$) in $\mathbf{y}$ by adjusting $\alpha$ to obtain
different decoding performance of two users. We can also see from (7) that it
is different with the superposition coding scheme in DBC. For the
superposition coding, $y_{i}$ should be $q_{i}x_{1,i}^{rc}+w_{i}x_{2,i}^{rc}$,
which means $x_{1,i}^{rc}$ and $x_{2,i}^{rc}$ only transmit one time over the
channel. However, it is shown in (7) that $x_{1,i}^{rc}$ and $x_{2,i}^{rc}$
can be transmitted over the channel multiple times.

##  III LOSS FUNCTION AND TRAINING METHOD

We can observe from Fig.1 that the performance of $\mathbf{z}_{1}$ and
$\mathbf{z}_{2}$ is dependent on the semantic encoder/decoder and the FISF. SE
$f_{e}(\cdot;\theta)$ and SD $g_{d}(\cdot;\theta)$ have been well-researched
in the point-to-point semantic communication, thus we mainly focus on the loss
function design of FISF module in this paper. For the point-to-point system,
$\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ can be encoded as $\mathbf{x}^{e2e}_{1}$
and $\mathbf{x}^{e2e}_{2}$, and can be decoded as $\mathbf{z}_{1}^{e2e}$ and
$\mathbf{z}_{2}^{e2e}$, respectively.

The worse user only needs to reconstruct $\mathbf{s}_{1}$. Thus, in this
paper, we can design the training object to maximize the conditional mutual
information between $\mathbf{z}_{1}$ and $\mathbf{z}_{1}^{e2e}$ given
$\mathbf{x}^{e2e}_{1}$ and $\mathbf{x}^{e}_{1}$, as given by

$\displaystyle\max\
I(\mathbf{z}_{1};\mathbf{z}_{1}^{e2e}|\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1}).$
(8)

It indicates the proposed system tries to output a similar image at the worse
user as the excellent image in the point-to-point system. Based on the work of
[17], this optimization object is hard to achieve. We can achieve a relaxation
object by predicting
$\hat{\mathbf{z}}_{1}^{e2e}=h_{\kappa}(\mathbf{z}_{1},\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1})$
first and then estimating the posterior distribution
$\mathcal{P}_{\kappa}(\mathbf{z}_{1}^{e2e}|\hat{\mathbf{z}}_{1}^{e2e})$. The
relaxation form can be written as

$\displaystyle
I(\mathbf{z}_{1};\mathbf{z}_{1}^{e2e}|\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1})=\sup_{h_{\kappa}}\
\mathbb{E}_{p(\mathbf{x}_{1}^{e2e})}\
[H({p(\mathbf{z}_{1}^{e2e}|\mathbf{x}_{1}^{e2e})})]$
$\displaystyle+\mathbb{E}_{p(\mathbf{s}_{1},\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1})}\
[\log\mathcal{P}_{\kappa}(\mathbf{z}_{1}^{e2e}|\hat{\mathbf{z}}_{1}^{e2e})],$
(9)

where $H(x)$ denotes the entropy of random variable $x$. The first term is
regularization to avoid collapse in the point-to-point system and the second
term is log-likelihood prediction term for target representation. Because the
end-to-end system is well-train, the first term is irrelevant to the training
process.

Therefore, (8) can be solved by deriving the training loss as

$\displaystyle\min\
L_{1}(\theta,\zeta,\varphi)=-\log\mathcal{P}_{\kappa}(\mathbf{z}_{1}^{e2e}|\hat{\mathbf{z}}_{1}^{e2e})$
(10)

where $\theta$, $\zeta$, $\varphi$ are all the parameters of the whole
network. If the estimating posterior distribution $\mathcal{P}$ is Gaussian
distribution, it becomes mean squared error (MSE) loss. If it is Boltzmann
distribution, it becomes softmax cross-entropy loss.

The better user reconstructs $\mathbf{s}_{1}$ first and then reconstructs
$\mathbf{s}_{2}$. Likewise, the training object can writen as

$\displaystyle\max\
I(\tilde{\mathbf{z}}_{1},\mathbf{z}_{2};\tilde{\mathbf{z}}_{1}^{e2e},\mathbf{z}_{2}^{e2e}|\mathbf{x}^{e}_{1},\mathbf{x}^{e}_{2},\mathbf{x}^{e2e}_{1},\mathbf{x}^{e2e}_{2})$
$\displaystyle=\sup_{h_{\kappa}}\
\mathbb{E}_{p(\mathbf{x}_{1}^{e2e},\mathbf{x}_{2}^{e2e})}\
[H({p(\tilde{\mathbf{z}}_{1}^{e2e},\mathbf{z}_{2}^{e2e}|\mathbf{x}_{1}^{e2e},\mathbf{x}_{2}^{e2e})})]$
$\displaystyle+\mathbb{E}_{p(\mathbf{s}_{1},\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1},\mathbf{s}_{2},\mathbf{x}^{e}_{2},\mathbf{x}^{e2e}_{2})}\
[\log\mathcal{P}_{\kappa}(\tilde{\mathbf{z}}_{1}^{e2e},\mathbf{z}_{2}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e},\hat{\mathbf{z}}_{2}^{e2e})].$
(11)

where
$\hat{\mathbf{z}}_{2}^{e2e}=h_{\kappa}(\mathbf{z}_{2},\mathbf{x}^{e}_{2},\mathbf{x}^{e2e}_{2})$.
Similar to (III), the secord term is critial. Furthermore, we can prove it as

$\displaystyle\mathbb{E}_{p(\mathbf{s}_{1},\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1},\mathbf{s}_{2},\mathbf{x}^{e}_{2},\mathbf{x}^{e2e}_{2})}\
[\log\mathcal{P}_{\kappa}(\tilde{\mathbf{z}}_{1}^{e2e},\mathbf{z}_{2}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e},\hat{\mathbf{z}}_{2}^{e2e})]$
$\displaystyle=\mathbb{E}_{p(\mathbf{s}_{1},\mathbf{x}^{e}_{1},\mathbf{x}^{e2e}_{1},\mathbf{s}_{2},\mathbf{x}^{e}_{2},\mathbf{x}^{e2e}_{2})}\
[\log\mathcal{P}_{\kappa}(\tilde{\mathbf{z}}_{1}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e},\hat{\mathbf{z}}_{2}^{e2e})$
$\displaystyle+\log\mathcal{P}_{\kappa}(\mathbf{z}_{2}^{e2e}|\tilde{\mathbf{z}}_{1}^{e2e},\hat{\tilde{\mathbf{z}}}_{1}^{e2e},\hat{\mathbf{z}}_{2}^{e2e})]$
$\displaystyle=\mathbb{E}_{p(\mathbf{s}_{1},\mathbf{x}^{e}_{1},\mathbf{x}^{e2r}_{1},\mathbf{s}_{2},\mathbf{x}^{e}_{2},\mathbf{x}^{e2e}_{2})}\
[\log\mathcal{P}_{\kappa}(\tilde{\mathbf{z}}_{1}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e})$
(12)
$\displaystyle+\log\frac{\mathcal{P}_{\kappa}(\hat{\mathbf{z}}_{2}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e},\tilde{\mathbf{z}}_{1}^{e2e})}{\mathcal{P}_{\kappa}(\hat{\mathbf{z}}_{2}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e})}+\log\mathcal{P}_{\kappa}(\mathbf{z}_{2}^{e2e}|\tilde{\mathbf{z}}_{1}^{e2e},\hat{\tilde{\mathbf{z}}}_{1}^{e2e},\hat{\mathbf{z}}_{2}^{e2e})].$

The first term indicates that $\mathbf{s}_{1}$ is required to be reconstructed
first without any information about $\mathbf{s}_{2}$ and the third term shows
$\mathbf{s}_{2}$ should be reconstructed under the condition that
$\mathbf{s}_{1}$ has been reconstructed. The structure of the designed user
corresponds to formulation that we first reconstruct $\mathbf{s}_{1}$ as
$\tilde{\mathbf{z}}_{1}$ and then based on $\tilde{\mathbf{z}}_{1}$,
$\mathbf{z}_{2}$ is reconstructed. The second term indicates that
$\mathbf{s}_{1}$ has effects on $\mathbf{s}_{2}$. We design the mutual
attention module to address the effects. Therefore, when designing the loss
function, the second term is omitted and the third term only contains
$\hat{\tilde{z}}_{2}^{e2e}$. $\beta$ is given to balance the importance of the
two terms. The training loss is designed as

$\displaystyle\min\ L_{2}(\theta,\zeta,\varphi)$
$\displaystyle=-\log\mathcal{P}_{\kappa}(\mathbf{z}_{2}^{e2e}|\hat{\mathbf{z}}_{2}^{e2e})$
$\displaystyle-\beta\log\mathcal{P}_{\kappa}(\tilde{\mathbf{z}}_{1}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e}).$
(13)

Finally, the whole goal of the proposed system is to minimize $L_{1}$ and
$L_{2}$ at the same time. It also becomes a multi-criterion problem that aims
to find the Pareto optimal points. Scalarization is a standard technique for
finding Pareto optimal points. The final problem can be solved by deriving the
training loss as

$\displaystyle\min\ L(\theta,\zeta,\varphi)=L_{1}+\lambda
L_{2}=-\log\mathcal{P}_{\kappa}(\mathbf{z}_{1}^{e2e}|\hat{\mathbf{z}}_{1}^{e2e})$
$\displaystyle-\lambda\log\mathcal{P}_{\kappa}(\mathbf{z}_{2}^{e2e}|\hat{\mathbf{z}}_{2}^{e2e})-\lambda\beta\log\mathcal{P}_{\kappa}(\tilde{\mathbf{z}}_{1}^{e2e}|\hat{\tilde{\mathbf{z}}}_{1}^{e2e}),$
(14)

where $\lambda$ is the scalarization parameter. In this paper, we consider the
posterior distribution $\mathcal{P}_{\kappa}$ is Gaussian distribution and
therefore the loss $L$ can be computed as

$\displaystyle
L=-||\mathbf{z}_{1}-\mathbf{s}_{1}||_{2}-\lambda||\mathbf{z}_{2}-\mathbf{s}_{2}||_{2}-\lambda\beta||\tilde{\mathbf{z}}_{1}-\mathbf{s}_{1}||_{2}.$
(15)

When training the model, $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ should be
generated from the same dataset $S$ individually. A copy of the dataset
$\tilde{S}$ is loaded with shuffling. $\mathbf{s}_{1}$ comes from the batch of
$S$ and $\mathbf{s}_{2}$ comes from the batch $\tilde{S}$. For the FISF
module, the SNR of the worse channel ($SNR_{1}$) is randomly set in a given
range and we then randomly set $\gamma>0$ so that the SNR of better channel
($SNR_{2}$) is $SNR_{1}+\gamma$. Fusion rate $\alpha$ is also randomly
selected in the range between 0 and 1 with step 0.1. When it is 0 or 1, which
means only one source is expected to be delivered, the model degrades to an
point-to-point model. The whole system takes $L$ as the loss function and all
the parameters are updated jointly according to the loss $L$. The whole
training procedures are described in Algorithm 1.

Input: Training set $S$, hyper-parameter $\lambda$ and $\beta$.
Output: The trained model with one transmitter and two users.

Algorithm 1 Training algorithm

1: Copy a dataset as $\tilde{S}$ and shuffle it.

2: while the training stop condition is not met do

3: Take a batch $\mathbf{s}_{1}$ from the set $S$.

4: Take another batch $\mathbf{s}_{2}$ from the dataset $\tilde{S}$.

5: Randomly sample $SNR_{1}$, $\gamma$, $\alpha$ individually.

6: Set $SNR_{2}=SNR_{1}+\gamma$.

7: Semantic encode and fuse $\mathbf{s}_{1}$ and $\mathbf{s}_{2}$ as
$\mathbf{y}$ (tramsmitter).

8: Transmit $\mathbf{y}$ over the two channels and users get $\mathbf{y}_{1}$,
$\mathbf{y}_{2}$.

9: Semantic de-fuse and decode $\mathbf{y}_{1}$ as $\mathbf{z}_{1}$ (worse
user).

10: Compute loss $L_{1}$ (worse user).

11: Semantic de-fuse and decode $\mathbf{y}_{2}$ as $\tilde{\mathbf{z}}_{1}$
and $\mathbf{z}_{2}$ (better user).

12: Compute $L_{2}$ (better user).

13: if $\alpha==0$ then

14: Set $\beta=0$ and then compute loss $L=L_{2}$.

15: else if $\alpha==1$ then

16: Compute loss $L=L_{1}$.

17: else

18: Compute loss $L=L_{1}+\lambda L_{2}$.

19: end if

20: Update all the parameters to minimize $L$.

21: end while

##  IV EXPERIMENTS

In this section, we evaluate the performance of the proposed semantic
communication scheme for DBC to transmit the image by using CIFAR-10 dataset.
We use Adam optimizer to train the system for 100 epochs with a learning rate
of $1\times 10^{-4}$ and then train for another 50 epochs with a learning rate
$1\times 10^{-5}$. The batch size is 128. For the scalarization parameters, we
set $\lambda=6$ and $\lambda\beta=0.1$. Without loss of generality, we use
$l=2k$, $h=n=32$ and the bandwidth ratio $\frac{k}{h\times n\times 3}=0.25$ in
the experiment.

(a) Raw image $\mathbf{s}_{1}$

(b) Proposed

(c) TD

(d) PA (e) Raw image $\mathbf{s}_{2}$ (f) Proposed (g) TD (h) PA

]

Figure 3: Raw iamges and images reconstructed by different methods with
CIFAR10 dataset.

In the experiment, we consider the power allocation (PA) scheme and the TD
scheme as benchmarks. For both benchmarks, $\mathbf{s}_{i}$ can pass the
semantic encoder and then be fed into its own MLP, yielding the transmitted
symbols for one user. For the PA scheme, the superposition coding with
successive interference cancellation is then exploited to produce the
transmission symbols $\mathbf{y}$. For the TD scheme, the transmission symbols
for different users are transmitted over their assigned time slots.

Fig. 3 depicts the visible results of the reconstructed images based on
different methods. We use $\alpha=0.5$, $SNR_{1}=5\leavevmode\nobreak\ dB$ and
$SNR_{2}=10\leavevmode\nobreak\ dB$. The upper image is for the worse user,
while the lower image is for the better user. It can be observed that the
proposed scheme produces the most clear recovered images among those based on
the TD and PA schemes.

Figure 4: Semantic performance region for DBC with different methods. Here,
the SNR of the worse/better user is $-5\leavevmode\nobreak\
dB$/$0\leavevmode\nobreak\ dB$. Figure 5: PSNR vs. SNR of the better user with
the proposed schme. Here, $\alpha$ is 0.5. Figure 6: PSNR vs. SNR of the worse
user with the proposed schme. Here, $\alpha$ is 0.5.

Next, we evaluate the performance of DBC by using PSNR. The PSNR of a single
user cannot reflect the comprehensive performance of DBC. Therefore, we can
describe the achievable PSNR groups of both users, which form a region called
the _semantic performance region_. Fig. 4 shows the semantic performance
region with different schemes. For the proposed scheme, we can obtain
different PSNR groups by adjusting fusion ratio $\alpha$, as shown in Fig. 4.
It is seen that with the proposed scheme, the PSNR of the worse user increases
with $\alpha$, but leads to the decrease of PSNR of the better user. This
result matches (7). We also clearly see that the semantic performance region
of the proposed scheme strictly contains the region of other bechmarks. This
fact shows that the proposed scheme can achieve the best performance for both
users in DBC compared to traditional methods. The PA scheme without channel
adaptive (CA) has the smallest region, which indicates that the power
allocation scheme is not suitable in the semantic communications system. The
gap between the power allocation scheme without CA and that of with CA reveals
the CA gain. The gap between the power allocation scheme with CA and the
proposed scheme reveals the fusion gain.

Fig. 5 and Fig. 6 show PSNR vs. SNR for the better user and the worse user
based on the proposed scheme, respectively.We can observe that the proposed
FISF module incorporating CSI into semantic features to adapt the channel can
provide significant performance gain. For example, for the better user at
$SNR=-5$ dB the system training at -5 dB achieve the best performance than
that of other SNR-fixed training, e.g. 0 dB. However, it is still 1.29 dB
lower in the PSNR performance than proposed scheme.

##  V CONCLUSION

In this paper, we have proposed a novel semantic communications system for
wireless image transmission over two-user degraded broadcast channels. The
transmitter can extract the semantic features of two users’ images and fuse
these semantic features into a joint latent representation for broadcasting.
We have designed a flexible image semantic fusion scheme that dynamically
controls the weight of semantic features in the joint latent representation
and adapts the respective semantic features to the respective channels with
different SNRs. Experimental results have shown that the proposed system
significantky dominates the traditional methods, such as TD and PA, for
wireless image transmission over two-user degraded broadcast channels.

## References

  * [1] E. Bourtsoulatze, D. Burth Kurka, and D. Gündüz, “Deep joint source-channel coding for wireless image transmission,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 5, no. 3, pp. 567–579, 2019.
  * [2] W. Yang, H. Du, Z. Q. Liew, W. Y. B. Lim, Z. Xiong, D. Niyato, X. Chi, X. S. Shen, and C. Miao, “Semantic communications for future internet: Fundamentals, applications, and challenges,” _IEEE Communications Surveys & Tutorials_, pp. 1–1, 2022.
  * [3] H. Xie, Z. Qin, G. Y. Li, and B.-H. Juang, “Deep learning enabled semantic communication systems,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 2663–2675, 2021.
  * [4] Q. Zhou, R. Li, Z. Zhao, C. Peng, and H. Zhang, “Semantic communication with adaptive universal transformer,” _IEEE Wireless Communications Letters_ , vol. 11, no. 3, pp. 453–457, 2022.
  * [5] D. Huang, F. Gao, X. Tao, Q. Du, and J. Lu, “Toward semantic communications: Deep learning-based image semantic coding,” _IEEE Journal on Selected Areas in Communications_ , vol. 41, no. 1, pp. 55–71, 2023.
  * [6] H. Zhang, S. Shao, M. Tao, X. Bi, and K. B. Letaief, “Deep learning-enabled semantic communication systems with task-unaware transmitter and dynamic data,” _IEEE Journal on Selected Areas in Communications_ , vol. 41, no. 1, pp. 170–185, 2023.
  * [7] T.-Y. Tung and D. Gündüz, “Deepwive: Deep-learning-aided wireless video transmission,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 9, pp. 2570–2583, 2022.
  * [8] S. Wang, J. Dai, Z. Liang, K. Niu, Z. Si, C. Dong, X. Qin, and P. Zhang, “Wireless deep video semantic transmission,” _IEEE Journal on Selected Areas in Communications_ , vol. 41, no. 1, pp. 214–229, 2023.
  * [9] P. Jiang, C.-K. Wen, S. Jin, and G. Y. Li, “Wireless semantic communications for video conferencing,” _IEEE Journal on Selected Areas in Communications_ , vol. 41, no. 1, pp. 230–244, 2023.
  * [10] X. Mu, Y. Liu, L. Guo, and N. Al-Dhahir, “Heterogeneous semantic and bit communications: A semi-noma scheme,” _IEEE Journal on Selected Areas in Communications_ , vol. 41, no. 1, pp. 155–169, 2023.
  * [11] S. F. Yilmaz, C. Karamanli, and D. Gündüz, “Distributed deep joint source-channel coding over a multiple access channel,” _https://arxiv.org/abs/2211.09920_ , 2022.
  * [12] H. Hu, X. Zhu, F. Zhou, W. Wu, R. Q. Hu, and H. Zhu, “One-to-many semantic communication systems: Design, implementation, performance evaluation,” _IEEE Communications Letters_ , vol. 26, no. 12, pp. 2959–2963, 2022.
  * [13] X. Luo, B. Yin, Z. Chen, B. Xia, and J. Wang, “Autoencoder-based semantic communication systems with relay channels,” in _2022 IEEE International Conference on Communications Workshops (ICC Workshops)_ , 2022, pp. 711–716.
  * [14] C. Bian, Y. Shao, and D. Gündüz, “Deep joint source-channel coding over cooperative relay networks,” _https://arxiv.org/abs/2211.06705_ , 2022.
  * [15] J. Xu, B. Ai, W. Chen, A. Yang, P. Sun, and M. Rodrigues, “Wireless image transmission using deep source channel coding with attention modules,” _IEEE Transactions on Circuits and Systems for Video Technology_ , vol. 32, no. 4, pp. 2315–2328, 2022.
  * [16] H. Wu, Y. Shao, K. Mikolajczyk, and D. Gündüz, “Channel-adaptive wireless image transmission with ofdm,” _IEEE Wireless Communications Letters_ , vol. 11, no. 11, pp. 2400–2404, 2022.
  * [17] W. Su, X. Zhu, C. Tao, L. Lu, B. Li, G. Huang, Y. Qiao, X. Wang, J. Zhou, and J. Dai, “Towards All-in-one Pre-training via Maximizing Multi-modal Mutual Information,” _https://arxiv.org/abs/2211.09807_ , 2022.