# Offloading and Quality Control for AI Generated Content Services in Edge
Computing Networks

Yitong Wang, Chang Liu, Jun Zhao  School of Computer Science and Engineering
Nanyang Technological University, Singapore
{yitong002<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

AI-Generated Content (AIGC), as a novel manner of providing Metaverse services
in the forthcoming Internet paradigm, can resolve the obstacles of immersion
requirements. Concurrently, edge computing, as an evolutionary paradigm of
computing in communication systems, effectively augments real-time interactive
services. In pursuit of enhancing the accessibility of AIGC services, the
deployment of AIGC models (e.g., diffusion models) to edge servers and local
devices has become a prevailing trend. Nevertheless, this approach faces
constraints imposed by battery life and computational resources when tasks are
offloaded to local devices, limiting the capacity to deliver high-quality
content to users while adhering to stringent latency requirements. So there
will be a tradeoff between the utility of AIGC models and offloading decisions
in the edge computing paradigm. This paper proposes a joint optimization
algorithm for offloading decisions, computation time, and diffusion steps of
the diffusion models in the reverse diffusion stage. Moreover, we take the
average error into consideration as the metric for evaluating the quality of
the generated results. Experimental results conclusively demonstrate that the
proposed algorithm achieves superior joint optimization performance compared
to the baselines.

###### Index Terms:

Edge computing, diffusion model, AIGC-as-a-service, resource allocation.

##  I Introduction

With GPT models capturing the spotlight, Generative Artificial Intelligence
(GAI), as a transformative field within the broader landscape of machine
learning and artificial intelligence, has changed the way people interact with
and understand the digital world [1]. The demonstration of the capabilities
inherent in Generative Artificial Intelligence (GAI) models is referred to as
AI-generated content (AIGC). With the development of AIGC techniques, multiple
AIGC models (e.g., diffusion models) can be employed to generate outputs with
diverse forms, including text-to-speech, text-to-image, and image-to-image
[2]. Therefore, AIGC-as-a-Service (AaaS) architecture is proposed to offer the
generated content, repair corrupted images, or alter inputted images,
resulting in providing Metaverse users with immersive AIGC services.

Edge computing, as a novel computing paradigm, not only solves the latency
concerns of cloud computing but also improves the security and quality of the
communication networks involved by integrating with other key technologies
including AI, Blockchain, and digital twin[3]. The convergence of the AIGC
models and edge computing paradigm becomes the focus of the future research
direction, especially in the Metaverse field[4] which emphasizes the immersion
of users.

Nonetheless, generic diffusion models (e.g., Open AI’s DALL-E 2 [5] and
Google’s Imagen [6]) require substantial memory storage, so deploying large
AIGC models is challenging due to the large volume of parameters [7], which
brings obstacles for real-time applications and devices with constrained
computational resources. Though lightweight diffusion models (like the text-
to-image model with 860M UNet and 123M text encoder proposed in [8]) are
employed on consumer mobile devices, the quality of generated content will
also inevitably be affected. Thus, the convergence of diffusion models and
edge computing systems has become the important direction of the future
research area.

Challenges: In light of the resource constraints inherent in mobile devices,
only lightweight diffusion models can be deployed on the local devices,
resulting in the locally generated content typically exhibiting a lower
quality level. While larger models can be deployed on the edge server, it
remains incapable of simultaneously managing all computational tasks while
generating high-quality content. Consequently, the first challenge centers on
the offloading decisions and quality level of the generated content. Secondly,
different locations to process the computational task means the computation
time is always different. As mobile users possess diverse and stringent
latency requirements, how to balance the tradeoff between the offloading
decisions and computation time is the second challenge. Furthermore, the
subjective nature of evaluating the quality level of generated content
necessitates the establishment of a mathematical correlation between quality
levels and diffusion models, thus constituting the third challenge.

Related Work: Recently, multiple studies have reviewed the state-of-the-art
research and development in diffusion models [9, 10, 11, 12, 13, 14]. Notably,
within the research of diffusion steps and image quality, [13] focuses on the
forward diffusion dynamics to bridge the gap between the ideal and the
simulated by adopting smaller diffusion times. Furthermore, [15] treats the
number of reverse diffusion steps as a variable to find the optimal reverse
diffusion step to balance the tradeoff between the image quality and the
diffusion time. For the convergence of the AIGC models and edge computing
paradigm, [16] introduces the diffusion model-based AI-generated optimal
decision (AGOD) algorithm to provide the optimal strategy for the selection of
AIGC service providers (AGPs). Unlike the previous works, we propose a
resource allocation scheme for the AGOD by correlating the image quality to
the diffusion steps and taking the computation time and utility of AIGC models
into consideration.

Contributions: The main contributions in this work are:

  * •

We are pioneering the implementation of the resource allocation for diffusion
models in edge computing systems while guaranteeing the quality of the AIGC,
which could enhance the performance and experience of the system and users
respectively.

  * •

We quantify the quality level of AICG by taking the reverse diffusion steps
within the diffusion process into consideration while considering the average
error of computation results.

  * •

The proposed optimization algorithm jointly optimizes the tradeoff between the
offloading decisions of the computational tasks and the utility of AIGC
models.

##  II System Model

We consider a mobile edge computing communication system with $N$ mobile user
equipments (UEs) and one edge server. Assume in this system that each UE
requests to access AIGC services and content of corresponding quality level is
generated at the chosen point based on the offloading decision and allocated
reverse diffusion steps.

AI-generated content model. To generalize the types of AIGC services (e.g.
text-to-image, text-to-video, and image-to-image generation), we employ score-
based diffusion models within our proposed system as the exemplar of AIGC
models. Assume that there has been a standard sampling and diffusion process
$S$ in the forward stage. AIGC services are then mainly provided by gradually
reversing the diffusion process, step by step. Consequently, the quality of
AIGC computation results is intricately related to the high quantity of
reverse diffusion steps during the reverse diffusion stage.

When the generation requests of AIGC services are evaluated by the monitor,
different offloading decisions and diffusion strategies of computational tasks
are allocated to the edge server or local UEs based on the current energy
conditions of the system. We denote the allocated offloading decisions and
reverse diffusion steps of computation tasks requested by UE $n$ in the
reverse stage as the binary variable $a_{n}\in\\{0,1\\}$ and discretization
variable $s_{n}$ respectively. Specifically, $a_{n}=0$ indicates that UE $n$’s
task will be processed locally on the mobile device, while $a_{n}=1$ indicates
that the computation task will be offloaded to the edge server.

For the reverse diffusion step $s_{n}$, as elucidated in the work by Du et al.
[16], it exhibits a positive correlation with the associated energy
expenditure. Given the finite energy resources inherent in the practical
system, it is thus noteworthy that $s_{n}$ is restricted as each step of the
reverse diffusion process necessitates energy consumption, primarily
associated with the execution of a neural network for Gaussian noise removal.
Therefore, we transform the total energy limitation into the constraint of
total supported reverse diffusion steps at different servers. To enhance the
generality of the proposed algorithm, we denote the total supported reverse
diffusion steps at the edge server as $S_{e}^{\text{max}}$ with the constraint
illustrated as:

$\displaystyle\sum\nolimits_{n\in\mathcal{N}}(a_{n}s_{n})\leq
S_{e}^{\text{max}},$ (1)

where $\mathcal{N}=\\{1,2,...,N\\}$ is the set of UEs. Furthermore, we denote
$S_{n,0}^{\text{max}}$ as the total supported reversed diffusion steps at each
local UE $n$. In order to mitigate the potential for specific computation
tasks to greedily consume computational resources on the edge server, a
constraint is imposed wherein the maximum number of reverse diffusion steps
allocated to each computation task is restricted. This constraint is denoted
as the maximum reverse diffusion steps limit for each UE denoted by
$S_{n,1}^{\text{max}}$, where $n$ signifies the $n$-th UE with $a_{n}=1$.
Therefore, another constraint considering the limitation of energy including
UE $n$ and the edge server is introduced as
$(1-a_{n})(s_{n}-S_{n,0}^{\text{max}})+a_{n}(s_{n}-S_{n,1}^{\text{max}})\leq
0$, and after mathematical transformations, this constraint is expressed as:

$\displaystyle
s_{n}-(1-a_{n})S_{n,0}^{\text{max}}-a_{n}S_{n,1}^{\text{max}}\leq 0,\forall
n\in\mathcal{N}.$ (2)

###  II-A Cost Functions

Computation time: As per the findings reported in Page 6 of the paper [15], it
is elucidated that the overall computation time can be expressed as an affine
function with respect to the reverse diffusion steps. Consequently, the
cumulative computational delays for tasks with different offloading decisions
are determined as:

$\displaystyle T_{n,0}(s_{n})=s_{n}\Delta
t_{n,0};~{}T_{n,1}(s_{n})=s_{n}\Delta t_{n,1},$ (3)

where $\Delta t_{n,0}$ and $\Delta t_{n,1}$ are designated as a constant
temporal interval per individual step within the process at UE $n$ and the
edge server. Given that the system is modeled for edge computing and the
durations associated with transmitting the computed results through the
downlink channels are negligible compared with the processing time of the edge
server, it is pertinent to note that we shall exclude further consideration of
downlink transmission time.

Average error of computation results: In order to mitigate the influence of
subjective factors on image quality and increase the alignment of the AIGC
content with the users’ request, we assess content quality by means of average
error metrics. Building upon the Eq.(16) provided in the work [15], we derive
that the reverse conditional diffusion pathway exhibits exponential error
reduction and proposed a modified version which could be defined as:

$\displaystyle\bar{\epsilon}_{n}(s_{n})=\bar{\epsilon}^{\text{fwd}}(S)e^{-s_{n}C_{1,n}},$
(4)

where function $\bar{\epsilon}^{\text{fwd}}(S)$ is to determine the quality of
processed content by adding Gaussian noise within the context of forward
diffusion process, which can be modeled mathematically by a convex function
related to the forward process $S$, and $C_{1,n}>0$, as the attenuation
factor, represents the recovering ability of the AIGC model. When
$S\rightarrow\infty$, it signifies that the forward diffusion has been
sufficiently close to the unknown and simple noise distribution. Furthermore,
when $S\rightarrow 0$, there is no diffusion forward process and the average
error converges to 0. Hence, function $\bar{\epsilon}^{\text{fwd}}(S)$
increases as the forward diffusion process $S$ increases.

Total energy consumption: Based on the relationship between CPU frequency
(cycles/s) and data size (bits), and drawing from relevant work [17], the
energy consumption associated with different offloading modes can be
formulated as: When $a_{n}=0$, the energy consumption for UE $n$, denoted as
$E_{n,0}(s_{n})$, can be expressed as:

$\displaystyle E_{n,0}(s_{n})=k_{n}T_{n,0}(s_{n})f_{n}^{3},$ (5)

where $f_{n}$ represents the CPU frequency of UE $n$, and $k_{n}$ is the
coefficient reflecting the power efficiency of UE $n$. Conversely, when
$a_{n}=1$, the energy consumption for UE $n$, referred to as $E_{n,1}(s_{n})$,
can be represented as

$\displaystyle E_{n,1}(s_{n})=k_{e}T_{n,1}(s_{n})g_{n}^{3},$ (6)

with $g_{n}$ representing the allocated computing capacity at the edge server
for the computational task requested by UE $n$. Note that $k_{e}$ represents
the analogous coefficient related to the power efficiency of the edge server.

Cost functions: Based on the previously delineated cost considerations, it can
be deduced that, in the scenario where $a_{n}=0$, the cost function denoted as
$\bar{R}_{n,0}(s_{n})$ is amenable to expression as
$\bar{R}_{n,0}(s_{n})=c_{1}T_{n,0}(s_{n})+c_{2}\bar{\epsilon}_{n}(s_{n})+c_{3}E_{n,0}(s_{n})$.
Conversely, in the case where $a_{n}=1$, the cost function assumes the form of
$\bar{R}_{n,1}(s_{n})=c_{1}T_{n,1}(s_{n})+c_{2}\bar{\epsilon}_{n}(s_{n})+c_{3}E_{n,1}(s_{n})$.
Herein, $c_{1}$, $c_{2}$, and $c_{3}$ denote the weighting coefficients for
each cost component.

###  II-B Utility of Reverse Diffusion Steps

We additionally contemplate the utility of reverse diffusion steps $s_{n}$ in
the reverse diffusion process (i.e., the alignment of the generated content
with the request). It is intuitive that the larger the reverse diffusion
steps, the higher the diversity of content generated and the stronger the
alignment with user requests. Consequently, the utility function should
exhibit a non-decreasing relationship with respect to $s_{n}$. Besides, there
is a marginal effect on the generation of the content (i.e., the content
generation approaches saturation as $s_{n}$ increases), so the function should
be concave. Then the utility function is defined as follows which is provided
in [18] for edge system:

$\displaystyle U_{n}(s_{n})=1-e^{-s_{n}C_{2,n}},$ (7)

where $C_{2,n}\geq 0$ is the constant parameter. Variations in the parameter
$s_{n}$ exert influence on the quality of the generated content, whereas
diminished values of $s_{n}$ entail trade-offs in terms of energy consumption
and utility. Therefore, it is imperative to establish a state of equilibrium.

##  III Joint Optimization of Cost and Utility

In this section, we build the original optimization problem and make it convex
by introducing auxiliary variables and adopting the penalized joint policy.

###  III-A Problem Formulation

We define $\bm{a}:=[a_{1},a_{2},\ldots,a_{N}]$ and
$\bm{s}:=[s_{1},s_{2},\ldots,s_{N}]$. In general, throughout this paper, for a
vector $\bm{v}$, we denote its $i$-th dimension by $v_{i}$. A joint
optimization problem, incorporating cost expenditure, offloading
determinations, and utility is formulated as problem $\mathbb{P}_{1}$:

$\displaystyle(\mathbb{P}_{1})\min_{\bm{a},\bm{s}}~{}\sum\nolimits_{n\in\mathcal{N}}\bigr{[}(1-a_{n})R_{n,0}(s_{n})+a_{n}R_{n,1}(s_{n})\bigr{]}$
(8) $\displaystyle\text{subject to}:a_{n}\in\\{0,1\\},\forall
n\\!\in\\!\mathcal{N},$ (9) Inequalities (1) and (2).

where $R_{n,r}(s_{n})=\omega_{1}\bar{R}_{n,r}(s_{n})-\omega_{2}U_{n}(s_{n})$,
$r\in\\{0,1\\}$, signifies the amalgamation of cost functions and utility
functions, and $\omega_{1}$, $\omega_{2}$ serve as weight parameters
specifically designated to modulate the magnitudes of the cost and utility
components. $a_{n}$ and $s_{n}$ are optimization variables. Constraint (9)
ensures the enforcement of binary offloading decisions. Constraint (1) imposes
an upper bound on the total reverse diffusion steps $S_{e}^{\text{max}}$ for
the edge server, thus addressing the constraint on available energy resources.
Additionally, $S_{n,0}^{\text{max}}$ in constraint (2) serves as a mechanism
for regulating the energy consumption of UE $n$, while $S_{n,1}^{\text{max}}$
acts as a constraint to prevent the excessive energy consumption associated
with multiple computation tasks on the edge server.

Then, we derived that the objective function of $\mathbb{P}_{1}$ is not
jointly convex, as the binary value of the optimization variable $a_{n}$.
Therefore, the optimization problem is stuck to be solved. More importantly,
the discrete values of the variable $a_{n}$ and the coupling of the decision
variables in the constraint (1) make this problem become an intractable NP-
hard problem.

###  III-B Transformation of the Non-Convex Problem

Derived from the formulated objective function, it becomes evident that the
terms associated with $a_{n}$ are amenable to extraction through standard
factorization procedures. Henceforth, we employ a strategic approach involving
the decoupling of the optimization variables, subsequently leading to the
transformation of the initial problem $\mathbb{P}_{1}$ into a set of sub-
problems [19] that can be addressed iteratively.

Handling objective function: To tackle the coupling of optimization variables,
surrogate functions of objective function can be constructed by introducing
auxiliary variables [20]. In the $(k+1)$-th iteration, the surrogate function
is expressed as:

$\displaystyle
G(\bm{a},\bm{s}|\bm{u}^{(k)},\bm{v}^{(k)})=\sum\nolimits_{n\in\mathcal{N}}g_{n}(a_{n},s_{n}|u_{n}^{(k)},v_{n}^{(k)}),$
(10)

where $u_{n}^{(k)}=\frac{(1-a_{n}^{(k)})}{2R_{n,0}(s_{n}^{(k)})}$ and
$v_{n}^{(k)}=\frac{a_{n}^{(k)}}{2R_{n,1}(s_{n}^{(k)})}\in\mathbb{R}^{+}$ are
the introduced auxiliary variables related to the convergence of our algorithm
which will be analyzed in Section III-F. Function
$g_{n}(a_{n},s_{n}|u_{n}^{(k)},v_{n}^{(k)})$ for UE $n$ is defined as:

$\displaystyle
g_{n}(\\!a_{n},\\!s_{n}\\!|u_{n}^{(k)}\\!\\!,\\!v_{n}^{(k)})\\!\\!=\\!\\!R_{n,0}^{2}(s_{n})u_{n}^{(k)}\\!\\!\\!+\\!\frac{(1-a_{n})^{2}}{4u_{n}^{(k)}}\\!\\!+\\!\\!R_{n,1}^{2}(s_{n})v_{n}^{(k)}\\!\\!\\!+\\!\\!\frac{a_{n}^{2}}{4v_{n}^{(k)}}.$

Handling constraint (1): In order to decouple the constraint (1), we adopt a
similar approach by introducing an auxiliary variable
$\bm{z}=[z_{1},...,z_{N}]$ to facilitate the formulation of the surrogate
function for each individual sub-problem. Consequently, in the $(k+1)$-th
iteration, the constraint undergoes a transformation as follows:

$\displaystyle-
S_{e}^{\text{max}}+\sum\nolimits_{n\in\mathcal{N}}(s_{n}^{2}z_{n}^{(k)}+\frac{a_{n}^{2}}{4z_{n}^{(k)}})\leq
0,$ (11)

where $z_{n}^{(k)}=\frac{a_{n}^{(k)}}{2s_{n}^{(k)}}$ constrains the variables
in each loop and converges with iterations.

Noting that the coupling among the optimization variables has been addressed,
it is pertinent to acknowledge that the optimization sub-problems remain non-
convex due to the discrete nature of the variables $a_{n}$.

###  III-C Successive Convex Approximation

In this section, we focus on using the penalized joint policy and successive
convex approximation (SCA) technique to handle the discrete variables $a_{n}$.

Handling constraint (9): To solve the discrete variable $a_{n}$ and without
loss of equivalence, it can be rewritten as:

$\displaystyle~{}~{}~{}~{}a_{n}\in[0,1],n\in\mathcal{N},$ (12)
$\displaystyle\sum\nolimits_{n\in\mathcal{N}}a_{n}(1-a_{n})\leq 0.$ (13)

Note that the optimization problem has been transitioned into a continuous
optimization problem, resulting in a notable reduction in computational
complexity when contrasted with the direct resolution of the original discrete
variable $a_{n}$. Nonetheless, the function
$\sum_{n\in\mathcal{N}}a_{n}(1-a_{n})$ in constraint (13) is a concave
function. To further facilitate the solution, we adopt a method that
introduces a penalty term for this concave constraint into the function
$G(\bm{a},\bm{s}|\bm{u}^{(k)},\bm{v}^{(k)})$, represented as:

$\displaystyle
G(\bm{a},\bm{s}|\bm{u}^{(k)},\bm{v}^{(k)})-\tau\cdot\sum\nolimits_{n\in\mathcal{N}}a_{n}(a_{n}-1),$
(14)

where $\tau$ is the penalty parameter with $\tau>0$. Then the objective
function becomes concave due to the concavity of the second term.
Simultaneously, given the second term is differentiable, we utilize the first-
order Taylor series to linearize it at each iteration. Specifically, at the
$(i+1)$-th iteration,we approximate $\sum_{n\in\mathcal{N}}a_{n}(a_{n}-1)$
with
$\sum_{n\in\mathcal{N}}a_{n}^{(i)}(a_{n}^{(i)}-1)+(2a_{n}^{(i)}-1)(a_{n}-a_{n}^{(i)})$
which is denoted as $H(\bm{a}|\bm{a}^{(i)})$, where $\bm{a}^{(i)}$ is defined
as the optimal solution of the $i$-th sub-problem. Consequently, the objective
function is converted to:

$\displaystyle G(\bm{a},\bm{s}|\bm{u}^{(k)},\bm{v}^{(k)})-\tau\cdot
H(\bm{a}|\bm{a}^{(i)}).$ (15)

Then, the $(i+1)$-th sub-problem in the $(k+1)$-th iteration is transformed
equivalently to $\mathbb{P}_{2}$:

$\displaystyle(\mathbb{P}_{2})\min_{\bm{a},\bm{s}}~{}G(\bm{a},\bm{s}|\bm{u}^{(k)},\bm{v}^{(k)})-\tau\cdot
H(\bm{a}|\bm{a}^{(i)})$ (16) $\displaystyle\text{subject
to}:a_{n}\in[0,1],\forall n\in\mathcal{N},$ (16a)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\sum\nolimits_{n\in\mathcal{N}}(s_{n}^{2}z_{n}^{(k)}+\frac{a_{n}^{2}}{4z_{n}^{(k)}})-S_{e}^{\text{max}}\leq
0,$ (16b)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}s_{n}-(1-a_{n})S_{n,0}^{\text{max}}-a_{n}S_{n,1}^{\text{max}}\leq
0,\forall n\in\mathcal{N}.$ (16c)

Until now, the original optimization problem $\mathbb{P}_{1}$ can be solved by
using $\mathbb{P}_{2}$ iteratively, and the process of solving the intra-sub-
problem is listed in Algorithm 1. Given introduced conditions and the
objective function are convex, $\mathbb{P}_{2}$ is convex. Thus, we can adopt
Karush-Kuhn-Tucker (KKT) conditions of $\mathbb{P}_{2}$ to obtain the optimal
solutions.

Algorithm 1 Intra-Sub-Problem Programming

1: Initialize: $k=0$; optimization variable
$\bm{\mathcal{X}}^{(0)}\\!=\\![\bm{a}^{(0)}\\!,\bm{s}^{(0)}]$;

2: Calculate and derive auxiliary variable space:
$\bm{\mathcal{A}}^{(0)}=[\bm{u}^{(0)}\\!\\!,\bm{v}^{(0)}\\!\\!,\bm{z}^{(0)}]$,
where $u_{n}^{(0)}\\!=\frac{(1-a_{n}^{(0)})}{2R_{n,0}(s_{n}^{(0)})}$,
$v_{n}^{(k)}\\!=\frac{a_{n}^{(0)}}{2R_{n,1}(s_{n}^{(0)})}$, and
$z_{n}^{(0)}=\frac{a_{n}^{(0)}}{2s_{n}^{(0)}}$, $\forall n\in\mathcal{N}$;

3: repeat

4: Obtain the optimal variable $\bm{\mathcal{X}}^{(k+1)}$ of $((k+1))$-th
iteration by adopting the Algorithm 2 when given auxiliary variable space
$\bm{\mathcal{A}}^{(k)}$;

5: Update
$\bm{\mathcal{A}}^{(k+1)}=[\bm{u}^{(k+1)}\\!\\!,\bm{v}^{(k+1)}\\!\\!,\bm{z}^{(k+1)}]$
with given $\bm{\mathcal{X}}^{(k+1)}$;

6: $k\leftarrow k+1$;

7: until Convergence or reach the max iteration number $K$.

###  III-D KKT Conditions for Problem $\mathbb{P}_{2}$

We first write down the Lagrange function of $\mathbb{P}_{2}$ by introducing
multipliers for the constraints:

$\displaystyle
L\\!=\\!G(\bm{a},\bm{s}|\bm{u}^{(k)}\\!\\!,\bm{v}^{(k)})\\!-\\!\tau\\!\\!\cdot\\!\\!H(\bm{a}|\bm{a}^{(i)})\\!+\\!\\!\sum\nolimits_{n\in\mathcal{N}}\\![\beta_{n}(\\!-a_{n})\\!+\\!\gamma_{n}(a_{n}\\!\\!-\\!\\!1)]\\!+\\!\delta\cdot$
$\displaystyle\bigr{[}\\!\\!\sum_{n\in\mathcal{N}}\\!\\!(s_{n}^{2}z_{n}^{(k)}\\!+\\!\frac{a_{n}^{2}}{4z_{n}^{(k)}}\\!)\\!-\\!S_{e}^{\text{max}}\bigr{]}\\!\\!+\\!\\!\\!\sum_{n\in\mathcal{N}}\\!\\!\zeta_{n}\bigr{[}s_{n}\\!\\!-\\!(\\!1\\!-\\!a_{n}\\!)S_{n,0}^{\text{max}}\\!-\\!a_{n}S_{n,1}^{\text{max}}\bigr{]},$
(17)

After applying KKT conditions, we get:

Stationarity:

$\displaystyle\frac{\partial L}{\partial
a_{n}}=D_{n}(a_{n},\zeta_{n},\delta)-\beta_{n}+\gamma_{n}=0,$ (18a)
$\displaystyle\frac{\partial L}{\partial s_{n}}=\frac{\partial
G(\bm{a},\bm{s}|\bm{u}^{(k)},\bm{v}^{(k)})}{\partial s_{n}}+2\delta
s_{n}z_{n}^{(k)}+\zeta_{n}=0.$ (18b)

where
$D_{n}(a_{n},\zeta_{n},\delta)=(\frac{1}{2u_{n}^{(k)}}+\frac{1}{2v_{n}^{(k)}}+\frac{\delta}{2z_{n}^{(k)}})a_{n}-\frac{1}{2u_{n}^{(k)}}-\tau(2a_{n}^{(i)}-1)+\zeta_{n}(S_{n,0}^{\text{max}}-S_{n,1}^{\text{max}})$
relating to the variable $a_{n}$.

Complementary slackness:

$\begin{split}&\text{(\ref{kktC}a)}\\!:\\!\beta_{n}(\\!-a_{n}\\!)\\!=\\!0;\hskip
7.0pt\text{(\ref{kktC}c)}\\!:\\!\delta\bigr{[}\\!\\!\sum_{n\in\mathcal{N}}\\!\\!(s_{n}^{2}z_{n}^{(k)}\\!+\\!\frac{a_{n}^{2}}{4z_{n}^{(k)}}\\!)\\!-\\!S_{e}^{\text{max}}\bigr{]}\\!=\\!0;\\\
&\text{(\ref{kktC}b)}\\!:\\!\gamma_{n}(\\!a_{n}\\!-\\!1\\!)\\!=\\!0;\hskip
2.0pt\text{(\ref{kktC}d)}\\!:\\!\zeta_{n}\bigr{[}s_{n}\\!\\!-\\!(\\!1\\!-\\!a_{n}\\!)S_{n,0}^{\text{max}}\\!-\\!a_{n}S_{n,1}^{\text{max}}\bigr{]}\\!=\\!0.\end{split}$
(19)

Primal Feasibility: (16a), (16b), (16c).

Dual Feasibility:

$\begin{split}\text{(\ref{kktD}a)-(\ref{kktD}d)}:\beta_{n},\gamma_{n},\delta,\zeta_{n}\geq
0,\forall n\in\mathcal{N}.\end{split}$ (20)

Under the aforementioned conditions, we then proceed to seek the optimal
solutions by analyzing the KKT conditions and employing the proposed algorithm
which is listed in Algorithm 2.

Algorithm 2 Solve KKT Conditions

1: Given the auxiliary variable space
$\bm{\mathcal{A}}=[\bm{u},\bm{v},\bm{z}]$;

2: for $n\leftarrow 1\text{ to }N$ do

3: Obtain $\widehat{a}_{n}(\zeta_{n},\delta)$ by assuming
$D_{n}(a_{n},\zeta_{n},\delta)$ on condition (18a) equal to 0;

4: Obtain $\widetilde{s}_{n}(\zeta_{n},\delta)$ based on condition (18b);

5: Obtain $\widehat{\zeta}_{n}(\delta)$ when assuming
$I_{n}(\zeta_{n},\delta)=0$ in the condition (LABEL:kktCd);

6: end for

7:
$\delta^{*}\\!\\!\leftarrow\\!\\!\left\\{\\!\\!\\!\small\begin{array}[]{l}0,\text{if}\sum_{n\in\mathcal{N}}(\widetilde{s}_{n}^{2}(\widetilde{\zeta}_{n}(0),0)z_{n}^{(k)}+\frac{a_{n}^{2}(\widetilde{\zeta}_{n}(0),0)}{4z_{n}^{(k)}})\leq
S_{e}^{\text{max}};\\\ \text{Solution to
}\sum_{n\in\mathcal{N}}(\widetilde{s}_{n}^{2}(\widetilde{\zeta}_{n}(\delta)\\!,\\!\delta)z_{n}^{(k)}\\!\\!\\!+\\!\frac{a_{n}^{2}(\widetilde{\zeta}_{n}(\delta),\delta)}{4z_{n}^{(k)}}\\!)\\!\\!=\\!\\!S_{e}^{\text{max}};\end{array}\right.$

8: update
$~{}a_{n}^{*}\\!\leftarrow\\!\max\\{\min\\{\widehat{a}_{n}(\widetilde{\zeta}_{n}(\delta^{*}),\delta^{*}),1\\},0\\}$,
and
$s_{n}^{*}\\!\\!\leftarrow\widetilde{s}_{n}(\widetilde{\zeta}_{n}(\delta^{*}),\delta^{*})$;

9: return The optimal variable value
$\mathcal{X}^{*}=[\bm{a}^{*},\bm{s}^{*}]$.

Theorem 1: The optimal solution of the proposed objective function can be
obtained by Algorithm 2 and is expressed as:

$\displaystyle\left\\{\\!\\!\begin{array}[]{l}a_{n}^{*}=\max\\{\min\\{\widehat{a}_{n}(\widetilde{\zeta}_{n}(\delta^{*}),\delta^{*}),1\\},0\\};\\\
s_{n}^{*}=\widetilde{s}_{n}(\widetilde{\zeta}_{n}(\delta^{*}),\delta^{*}).\end{array}\right.$
(23)

where $\widetilde{\zeta}_{n}(\delta)=\max\\{\widehat{\zeta}_{n}(\delta),0\\}$.
More specifically, $\widehat{a}_{n}(\zeta_{n},\delta)$ meets the condition
(18a) with
$D_{n}(a_{n},\zeta_{n},\delta)|_{a_{n}=\widehat{a}_{n}(\zeta_{n},\delta)}\\!=\\!0$
and $\widehat{\zeta}_{n}(\delta)$ satisfies condition (LABEL:kktCd) with
$I_{n}(\zeta_{n},\delta)|_{\zeta_{n}=\widehat{\zeta}_{n}(\delta)}\\!=\\!0$.

Proof: Please see Appendix References. $\hfill\blacksquare$

###  III-E Inter-Sub-Problem Algorithm

After completing the analysis of KKT conditions within each sub-problem, we
obtain the optimal solution for each sub-problem. In this section, we
introduce the inter-sub-problem algorithm, denoted as Algorithm 3, which
leverages iterative implementations of the SCA method. Its primary objective
is the determination of the globally optimal solution for the original
optimization problem $\mathbb{P}_{1}$.

Algorithm 3 Inter-Sub-Problem Programming

1: Initialize: $i=1$; Optimal solution space
$\bm{\mathcal{S}}^{(0)}\\!\\!\\!=\\!\\![\bm{a}^{(0)}\\!\\!,\bm{s}^{(0)}]$;

2: Adopt SCA method to obtain $\mathbb{P}_{2}$ from $\mathbb{P}_{1}$;

3: repeat

4: Obtain the $i$-th sub-problem by using $\bm{a}^{(i-1)}$;

5: Solve the $i$-th sub-problem by using Algorithm 1 and 2 to get the optimal
optimal variable value $\mathcal{X}^{*}\\!=\\![\bm{a}^{*}\\!,\bm{s}^{*}]$;

6: Let
$\bm{\mathcal{S}}^{i}[\bm{a}^{(i)},\bm{s}^{(i)}]\leftarrow\mathcal{X}^{*}\\!=\\![\bm{a}^{*}\\!,\bm{s}^{*}]$
of $i$-th sub-problem;

7: Set $i\leftarrow i+1$;

8: until
$|\bm{\mathcal{S}}^{(i)}-\bm{\mathcal{S}}^{(i-1)}|\leq\bar{\epsilon}_{0}$ or
reach the maximum iteration number $I$;

9: return $\bm{\mathcal{S}}^{(i)}[\bm{a}^{(i)},\bm{s}^{(i)}]$ as the optimal
solution of $\mathbb{P}_{1}$.

###  III-F Time complexity, Solution Quality and Convergence

Time Complexity. Based on the algorithms listed, the complexity of Algorithm 3
lies in the step 3-8. Assuming $K$ is the number of iterations in Algorithm 3,
then the complexity of adopting $\bm{a}^{i-1}$ to obtain the sub-problem in
step 4 is denoted as $\mathcal{O}(IN)$, where $N$ derives from computing
$a_{n}^{i-1}$ for each user in an iteration. In step 5, Algorithm 1 is mainly
solved by Algorithm 2 to obtain the optimization solution of the intra-sub-
problem, so the complexity of Algorithm 2 is first analyzed. The total
complexity of Algorithm 2 is $\mathcal{O}(2N)$ as the complexity of computing
$\bm{a}^{*}$ and $\bm{s}^{*}$ are both $\mathcal{O}(N)$ separately. Let $K$
denote the number of iterations in Algorithm 1. The complexity of the
Algorithm 1 is thus $\mathcal{O}(3KN)$ as step 5 of Algorithm 1 also costs
$\mathcal{O}(N)$. Therefore, the overall complexity can be derived as
$\mathcal{O}((I+3KI)N)$.

Solution quality and convergence. Algorithm 3 comprises SCA, Algorithm 1, and
Algorithm 2. Though the SCA method adopted to transform $\mathbb{P}_{1}$ to
$\mathbb{P}_{2}$ results in some loss of the optimality, the KKT analysis by
introducing the auxiliary variables in Algorithm 1 and Algorithm 3 both are
without loss of optimality. Thus Steps 3-8 of Algorithm 3 can guarantee to
find the global optimal solutions for $\mathbb{P}_{2}$. The convergence of
Algorithm 3 is also evident from the preceding analysis.

##  IV Experimental Results

###  IV-A Parameter Settings

In this experiment, we set the total number of users $N$ is 30. To generalize
the experimental result without any distortion and error, the original value
of the average error rate is set as 1, which implies the original content is
complete Gaussian noise. Based on the resource limitation, the fixed
discretization steps are $\Delta t_{n,0}=1/500$ and $\Delta t_{n,1}=1/1000$.
The computing capacity $g_{n}$ and frequency of mobile user $f_{n}$ are set as
$10$GHz and $1.5$GHz by default. The computation energy efficiency coefficient
$k_{n}$ and $k_{e}$ is $10^{-26}$. For the penalty parameter $\tau$, we set
the value as $10^{5}$.

###  IV-B Performance when Adapting Weight Parameters

(a)  (b)

(c)

Figure 1: Consumption under different maximum diffusion steps $S_{e}^{\max}$.

Figure 2: Joint optimization Performance.

As the weight factors $(c_{1},c_{2},c_{3})$ of cost functions have an effect
on the optimization process, we make modifications to these weights to tailor
the focus of the associated cost functions. To investigate the different costs
in the system, we augment $c_{1}$, $c_{2}$, and $c_{3}$ to respectively
heighten the sensitivity of computational time, average error rate, and energy
consumption. To further investigate the outcomes of our proposed algorithm, we
conduct experimental trials across the combinations of $(c_{1},c_{2},c_{3})$
and subsequently compare them to a baseline: random initialization. This
random initialization entails the arbitrary allocation of diffusion steps
following the stochastic selection of offloading decisions. For the baseline,
we assume that no preference exists in the algorithm, then the weight factor
is equal. Results of Fig. 2 show the computation $T$, average error of
computation results $\bar{\epsilon}$ and energy consumption $E$ under
different maximum assignable diffusion step $S_{e}^{\max}$. From Fig. 2(a), we
can see the total consumption time is decreased as $S_{e}^{\max}$ increase,
which means more complex computational tasks are processed on the edge server
and the fixed discretization step on the edge server is faster. As the energy
consumption is also larger than the local devices, the energy consumption also
increases with increasing reverse diffusion steps which is presented in Fig.
2(c). To provide a more intuitive representation of the quality of AIGC, we
adopt the average accuracy ($1-\bar{\epsilon}$) instead of average error as
shown in Fig. 2(b). For the different weights lines, the red line (proposed
equal weights) outperforms others comprehensively as no sacrifice exists in
this scenario compared with other lines. Furthermore, we can see that all
proposed methods outperform the baseline, demonstrating our method’s
superiority.

###  IV-C Performance when Setting Different Diffusion Steps

To conduct a comprehensive investigation into the performance of the proposed
method, we consider diverse joint optimization scenarios when the system
requires different $S_{e}^{\max}$. To emphasize the utility function, we set
the parameters $(w_{1},w_{2})$ as $(0.3,0.7)$. Three distinct experiments are
undertaken, each characterized by different local reverse diffusion steps
$S_{0}^{\max}$. As illustrated in Fig.2, we can see that as $S_{e}^{max}$
increases, all three lines exhibit an initial decline followed by a
stabilization phase upon reaching an optimal solution. This observation
underscores that as the system benefits from a higher resource allocation for
task processing, it tends to prioritize the augmentation of the utility
function to achieve superior joint optimization performance.

##  V Conclusion

In this paper, we propose one joint optimization algorithm that bridges the
gap between AIGC models and edge computing while mitigating the constraints
posed by resource limitations on devices. The presented algorithm offers an
enhanced approach by simultaneously optimizing offloading decisions and the
reverse diffusion steps of diffusion models, taking into account average error
and energy consumption. Our analysis of experimental results demonstrates the
effectiveness of the proposed algorithm in enhancing the system efficiency.

## References

  * [1] L. Floridi and M. Chiriatti, “GPT-3: Its Nature, Scope, Limits, and Consequences,” _Minds and Machines_ , vol. 30, pp. 681–694, 2020.
  * [2] X. Xu, Z. Wang, G. Zhang, K. Wang, and H. Shi, “Versatile Diffusion: Text, Images and Variations All in One Diffusion Model,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2023, pp. 7754–7765.
  * [3] Y. Mao, C. You, J. Zhang, K. Huang, and K. B. Letaief, “A Survey on Mobile Edge Computing: The Communication Perspective,” _IEEE communications surveys & tutorials_, no. 4, pp. 2322–2358, 2017.
  * [4] M. Xu, D. Niyato, H. Zhang, J. Kang, Z. Xiong, S. Mao, and Z. Han, “Sparks of GPTs in Edge Intelligence for Metaverse: Caching and Inference for Mobile AIGC Services,” _arXiv:2304.08782_ , 2023.
  * [5] S. Xie, Z. Zhang, Z. Lin, T. Hinz, and K. Zhang, “SmartBrush: Text and Shape Guided Object Inpainting With Diffusion Model,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2023, pp. 22 428–22 437.
  * [6] C. Saharia, W. Chan, S. Saxena, L. Li, J. Whang, E. L. Denton, K. Ghasemipour, R. Gontijo Lopes, B. Karagol Ayan, T. Salimans _et al._ , “Photorealistic Text-to-Image Diffusion Models with Deep Language Understanding ,” _Advances in Neural Information Processing Systems_ , vol. 35, pp. 36 479–36 494, 2022.
  * [7] M. Xu, H. Du, D. Niyato, J. Kang, Z. Xiong, S. Mao, Z. Han, A. Jamalipour, D. I. Kim, V. Leung _et al._ , “Unleashing the Power of Edge-Cloud Generative AI in Mobile Networks: A Survey of AIGC Services,” _arXiv preprint arXiv:2303.16129_ , 2023.
  * [8] R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer, “High-Resolution Image Synthesis with Latent Diffusion Models,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2022, pp. 10 684–10 695.
  * [9] F.-A. Croitoru, V. Hondru, R. T. Ionescu, and M. Shah, “Diffusion Models in Vision: A Survey,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2023.
  * [10] L. Yang, Z. Zhang, Y. Song, S. Hong, R. Xu, Y. Zhao, W. Zhang, B. Cui, and M.-H. Yang, “Diffusion Models: A Comprehensive Survey of Methods and Applications,” _ACM Computing Surveys_ , 2022.
  * [11] L. Zhang, A. Rao, and M. Agrawala, “Adding Conditional Control to Text-to-Image Diffusion Models,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2023, pp. 3836–3847.
  * [12] X. Li, Y. Liu, L. Lian, H. Yang, Z. Dong, D. Kang, S. Zhang, and K. Keutzer, “Q-Diffusion: Quantizing Diffusion Models,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2023, pp. 17 535–17 545.
  * [13] G. Franzese, S. Rossi, L. Yang, A. Finamore, D. Rossi, M. Filippone, and P. Michiardi, “How Much Is Enough? A Study on Diffusion Times in Score-Based Generative Models,” _Entropy_ , no. 4, p. 633, 2023.
  * [14] Z. Zhang, C. Li, W. Sun, X. Liu, X. Min, and G. Zhai, “A Perceptual Quality Assessment Exploration for AIGC Images,” _arXiv preprint arXiv:2303.12618_ , 2023.
  * [15] H. Chung, B. Sim, and J. C. Ye, “Come-Closer-Diffuse-Faster: Accelerating Conditional Diffusion Models for Inverse Problems Through Stochastic Contraction,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 12 413–12 422.
  * [16] H. Du, J. Wang, D. Niyato, J. Kang, Z. Xiong, and D. I. Kim, “AI-Generated Incentive Mechanism and Full-Duplex Semantic Communications for Information Sharing,” _IEEE Journal on Selected Areas in Communications_ , vol. 41, no. 9, pp. 2981–2997, 2023.
  * [17] H. Xiao, J. Zhao, J. Feng, L. Liu, Q. Pei, and W. Shi, “Joint Optimization of Security Strength and Resource Allocation for Computation Offloading in Vehicular Edge Computing,” _IEEE Transactions on Wireless Communications_ , 2023.
  * [18] Q. Liu, S. Huang, J. Opadere, and T. Han, “An Edge Network Orchestrator for Mobile Augmented Reality,” in _IEEE INFOCOM 2018 - IEEE Conference on Computer Communications_ , 2018, pp. 756–764.
  * [19] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-Minimization Algorithms in Signal Processing, Communications, and Machine Learning,” _IEEE Transactions on Signal Processing_ , no. 3, pp. 794–816, 2016.
  * [20] J. Zhao, L. Qian, and W. Yu, “Human-Centric Resource Allocation in the Metaverse over Wireless Communications,” _arXiv preprint arXiv:2304.00355_ , 2023.

Analysis of KKT conditions: In order to identify the optimal
$(\bm{a}^{*},\bm{s}^{*})$ that satisfies the KKT conditions, we present the
ensuing analysis rooted in principles of optimization theory.

Deriving from the functional expression $D_{n}(a_{n},\zeta_{n},\delta)$ as
delineated in condition (18a), two properties of condition (18a) can be
inferred contingent upon the values of auxiliary variables and Lagrange
multipliers: 1). $D_{n}(a_{n},\zeta_{n},\delta)$ is non-decreasing for $a_{n}$
and 2). Specifically, we obtain the explicit expression of $a_{n}$ by setting
$\beta_{n}=0$ and $\gamma_{n}=0$. Then we denote the explicit expression as
$\widehat{a}_{n}(\zeta_{n},\delta)$ and proceed to discuss the different cases
based on the condition (LABEL:kktCa) and (LABEL:kktCb) as outlined below:

  * •

Case 1: $\widehat{a}_{n}(\zeta_{n},\delta)\geq 1$. In this case, we can infer
that $D_{n}(1,\zeta_{n},\delta)\leq 0$ based on the given two properties.
Therefore, we can set $\widehat{a}_{n}(\zeta_{n},\delta)=1$ and $\beta_{n}=0$
to meet the conditions, and then the value of $\gamma_{n}$ is equal to
$-D_{n}(1,\zeta_{n},\delta)\geq 0$ which exactly meet the condition
(LABEL:kktCb).

  * •

Case 2: $0<\widehat{a}_{n}(\zeta_{n},\delta)<1$. We can simply set the value
of $a_{n}$ as $\widehat{a}_{n}(\zeta_{n},\delta)$ with $\beta_{n}=0$ and
$\gamma_{n}=0$ in this case.

  * •

Case 3: $\widehat{a}_{n}(\zeta_{n},\delta)\leq 0$. When
$\widehat{a}_{n}(\zeta_{n},\delta)\leq 0$, it means when $a_{n}=0$,
$D_{n}(0,\zeta_{n},\delta)$ is equal or better than 0. Therefore, we choose
the feasible solution as $\widehat{a}_{n}(\zeta_{n},\delta)=0$,
$\gamma_{n}=0$, and $\beta_{n}=D_{n}(0,\zeta_{n},\delta)$ which meet the
condition (LABEL:kktCa).

Thus, we summarise all of these cases as follows:

$\displaystyle\left\\{\\!\\!\begin{array}[]{l}\widetilde{a}_{n}(\zeta_{n},\delta)=\max\\{\min\\{\widehat{a}_{n}(\zeta_{n},\delta),1\\},0\\};\\\
\widetilde{\beta}_{n}(\zeta_{n},\delta)=\max\\{D_{n}(0,\zeta_{n},\delta),0\\};\\\
\widetilde{\gamma}_{n}(\zeta_{n},\delta)=-\min\\{D_{n}(1,\zeta_{n},\delta),0\\};\end{array}\right.$
(27)

Before we continue the analysis, one propriety of function
$s_{n}\\!\\!-\\!(\\!1\\!-\\!a_{n}\\!)S_{n,0}^{\text{max}}\\!-\\!a_{n}S_{n,1}^{\text{max}}$
which is denoted as $I_{n}(\zeta_{n},\delta)$ for simplify in the condition
(LABEL:kktCd) can be derived: $I_{n}(\zeta_{n},\delta)$ is non-increasing of
$\zeta_{n}$. Then we discuss the following cased by setting $\zeta_{n}=0$:

  * •

Case 1: $I_{n}(0,\delta)\leq 0$. Based on the propriety mentioned above, one
feasible solution meeting the condition (LABEL:kktCd) is $\delta=0$.

  * •

Case 2: $I_{n}(0,\delta)>0$. Conversely, if $I_{n}(0,\delta)>0$, the value of
$\zeta_{n}$ can be set as $\widehat{\zeta}_{n}(\delta)$ based on the
$I_{n}(\zeta_{n},\delta)|_{\zeta_{n}=\widehat{\zeta}_{n}(\delta)}=0$.

The solution can be summarized
as:$\widetilde{\zeta}_{n}(\delta)=\max\\{\widehat{\zeta}_{n}(\delta),0\\}$.

Drawing from the condition (18b), the expression of $s_{n}$ can be deduced and
we denote it as $s_{n}(\zeta_{n},\delta)$. Similarly, we discuss the cases
based on condition (LABEL:kktCc) by setting $\delta=0$:

  * •

Case
1:$\sum_{n\in\mathcal{N}}(\widetilde{s}_{n}^{2}(\widetilde{\zeta}_{n}(0),0)z_{n}^{(k)}+\frac{a_{n}^{2}(\widetilde{\zeta}_{n}(0),0)}{4z_{n}^{(k)}})\leq
S_{e}^{\text{max}}$. In this case, we can set $\delta=0$ to meet the
conditions (20d) and (LABEL:kktCc).

  * •

Case
2:$\sum_{n\in\mathcal{N}}(\widetilde{s}_{n}^{2}(\widetilde{\zeta}_{n}(0),0)z_{n}^{(k)}+\frac{a_{n}^{2}(\widetilde{\zeta}_{n}(0),0)}{4z_{n}^{(k)}})\textgreater
S_{e}^{\text{max}}$. Conversely, we need to find the optimal value of the
$\delta$ based on the condition (LABEL:kktCc) by using the bisection method.
Then we denote the solution as $\widetilde{\delta}$.

Consequently, the optimal value of $\delta$ is expressed as:

$\displaystyle\delta^{*}=\left\\{\\!\\!\begin{array}[]{ll}0,\\!\\!\\!&\sum_{n\in\mathcal{N}}(\widetilde{s}_{n}^{2}(\widetilde{\zeta}_{n}(0),0)z_{n}^{(k)}\\!\\!+\frac{a_{n}^{2}(\widetilde{\zeta}_{n}(0),0)}{4z_{n}^{(k)}})\\!\leq\\!S_{e}^{\text{max}},\\\
\widetilde{\delta},\\!\\!\\!&\sum_{n\in\mathcal{N}}(\widetilde{s}_{n}^{2}(\widetilde{\zeta}_{n}(0),0)z_{n}^{(k)}\\!\\!+\frac{a_{n}^{2}(\widetilde{\zeta}_{n}(0),0)}{4z_{n}^{(k)}})\\!>\\!S_{e}^{\text{max}}.\end{array}\right.$
(30)

Until now, we can get the optimal solution of variables
$[\bm{a}^{*},\bm{s}^{*}]$ and the Lagrange multipliers
$[\bm{\beta}^{*},\bm{\gamma}^{*},\bm{\zeta}^{*},\delta^{*}]$ by analyzing KKT
conditions. $\hfill\blacksquare$


DESY 12-091 arXiv:1205.7051 [math.NT]

On Multiple Zeta Values of Even Arguments
Michael E. Hoffman111Supported by a grant from the German Academic Exchange
Service (DAAD) during the preparation of this paper. The author also thanks
DESY for providing facilities and financial support for travel.
Dept. of Mathematics, U. S. Naval Academy
Annapolis, MD 21402 USA
and
Deutsches Elektronen-Synchrotron DESY
Platanenalle 6, D-15738 Zeuthen, Germany
<EMAIL_ADDRESS>
June 8, 2012

Keywords: multiple zeta values, symmetric functions, Bernoulli numbers
MR Classifications: Primary 11M32; Secondary 05E05, 11B68

For $k\leq n$, let $E(2n,k)$ be the sum of all multiple zeta values with even
arguments whose weight is $2n$ and whose depth is $k$. Of course $E(2n,1)$ is
the value $\zeta(2n)$ of the Riemann zeta function at $2n$, and it is well
known that $E(2n,2)=\frac{3}{4}\zeta(2n)$. Recently Z. Shen and T. Cai gave
formulas for $E(2n,3)$ and $E(2n,4)$ in terms $\zeta(2n)$ and
$\zeta(2)\zeta(2n-2)$. We give two formulas for $E(2n,k)$, both valid for
arbitrary $k\leq n$, one of which generalizes the Shen-Cai results; by
comparing the two we obtain a Bernoulli-number identity. We also give an
explicit generating function for the numbers $E(2n,k)$.

##  1 Introduction and Statement of Results

For positive integers $i_{1},\dots,i_{k}$ with $i_{1}>1$, we define the
multiple zeta value $\zeta(i_{1},\dots,i_{k})$ by

$\zeta(i_{1},\dots,i_{k})=\sum_{n_{1}>\dots>n_{k}\geq
1}\frac{1}{n_{1}^{i_{1}}\cdots n_{k}^{i_{k}}}.$ (1)

The multiple zeta value (1) is said to have weight $i_{1}+\dots+i_{k}$ and
depth $k$. Many remarkable identities have been proved about these numbers,
but in this note we will concentrate on the case where the $i_{j}$ are even
integers. Let $E(2n,k)$ be the sum of all the multiple zeta values of even-
integer arguments having weight $2n$ and depth $k$, i.e.,

$E(2n,k)=\sum_{\begin{subarray}{c}\text{$i_{1},\dots,i_{k}$ even}\\\
i_{1}+\dots+i_{k}=2n\end{subarray}}\zeta(i_{1},\dots,i_{k}).$

Of course

$E(2n,1)=\zeta(2n)=\frac{(-1)^{n-1}B_{2n}(2\pi)^{2n}}{2(2n)!},$ (2)

where $B_{2n}$ is the $2n$th Bernoulli number, by the classical formula of
Euler. Euler also studied double zeta values (i.e., multiple zeta values of
depth 2) and in his paper [2] gave two identities which read

$\displaystyle\sum_{i=2}^{2n-1}(-1)^{i}\zeta(i,2n-i)$
$\displaystyle=\frac{1}{2}\zeta(2n)$
$\displaystyle\sum_{i=2}^{2n-1}\zeta(i,2n-i)$ $\displaystyle=\zeta(2n)$

in modern notation. From these it follows that

$E(2n,2)=\frac{3}{4}\zeta(2n),$

though Gangl, Kaneko and Zagier [3] seem to be the first to have pointed it
out in print. Recently Shen and Cai [10] proved the formulas

$\displaystyle E(2n,3)$
$\displaystyle=\frac{5}{8}\zeta(2n)-\frac{1}{4}\zeta(2)\zeta(2n-2),\ n\geq 3$
(3) $\displaystyle E(2n,4)$
$\displaystyle=\frac{35}{64}\zeta(2n)-\frac{5}{16}\zeta(2)\zeta(2n-2),\ n\geq
4.$ (4)

Identity (3) was also proved by Machide [9] using a different method.

This begs the question whether there is a general formula of this type for
$E(2n,k)$. The pattern

$\frac{3}{4},\quad\frac{3}{4}\cdot\frac{5}{6}=\frac{5}{8},\quad\frac{3}{4}\cdot\frac{5}{6}\cdot\frac{7}{8}=\frac{35}{64}$

of the leading coefficients makes one curious. In fact, the general result is
as follows.

######  Theorem 1.

For $k\leq n$,

$E(2n,k)=\frac{1}{2^{2(k-1)}}\binom{2k-1}{k}\zeta(2n)\\\
-\sum_{j=1}^{\lfloor\frac{k-1}{2}\rfloor}\frac{1}{2^{2k-3}(2j+1)B_{2j}}\binom{2k-2j-1}{k}\zeta(2j)\zeta(2n-2j).$

The next two cases after (4) are

$\displaystyle E(2n,5)$
$\displaystyle=\frac{63}{128}\zeta(2n)-\frac{21}{64}\zeta(2)\zeta(2n-2)+\frac{3}{64}\zeta(4)\zeta(2n-4)$
$\displaystyle E(2n,6)$
$\displaystyle=\frac{231}{512}\zeta(2n)-\frac{21}{64}\zeta(2)\zeta(2n-2)+\frac{21}{256}\zeta(4)\zeta(2n-4).$

We prove Theorem 1 in §3 below, using the generating function

$F(t,s)=1+\sum_{n\geq k\geq 1}E(2n,k)t^{n}s^{k}.$

In §2 we establish the following explicit formula.

######  Theorem 2.

$F(t,s)=\frac{\sin(\pi\sqrt{1-s}\sqrt{t})}{\sqrt{1-s}\sin(\pi\sqrt{t})}.$

Our proof uses symmetric functions. We define a homomorphism
$\mathfrak{Z}:\operatorname{Sym}\to\mathbf{R}$, where $\operatorname{Sym}$ is
the algebra of symmetric functions, and a family
$N_{n,k}\in\operatorname{Sym}$ such that $\mathfrak{Z}$ sends $N_{n,k}$ to
$E(2n,k)$. We then obtain a formula for the generating functions

$\mathcal{F}(t,s)=1+\sum_{n\geq k\geq
1}N_{n,k}t^{n}s^{k}\in\operatorname{Sym}[[t,s]]$

and apply $\mathfrak{Z}$ to get Theorem 2.

From the form of $\mathcal{F}(t,s)$ we show that it satisfies a partial
differential equation (Proposition 1 below), which is equivalent to a
recurrence for the $N_{n,k}$. From the latter we obtain a formula for the
$N_{n,k}$ in terms of complete and elementary symmetric functions, to which
$\mathfrak{Z}$ can be applied to give the following alternative formula for
$E(2n,k)$.

######  Theorem 3.

For $k\leq n$,

$E(2n,k)=\frac{(-1)^{n-k-1}\pi^{2n}}{(2n+1)!}\sum_{i=0}^{n-k}\binom{n-i}{k}\binom{2n+1}{2i}2(2^{2i-1}-1)B_{2i}.$

Note that the sum given by Theorem 3 has $n-k+1$ terms, while that given by
Theorem 1 has $\lfloor\frac{k-1}{2}\rfloor+1$ terms. Yet another explicit
formula for $E(2n,k)$ can be obtained by setting $d=1$ in Theorem 7.1 of
Komori, Matsumoto and Tsumura [7]. That formula expresses $E(2n,k)$ as a sum
over partitions of $k$, and it is not immediately clear how it relates to our
two formulas.

Comparison of Theorems 1 and 3 establishes the following Bernoulli-number
identity.

######  Theorem 4.

For $k\leq n$,

$\sum_{i=0}^{\lfloor\frac{k-1}{2}\rfloor}\binom{2k-2i-1}{k}\binom{2n+1}{2i+1}B_{2n-2i}=\\\
(-1)^{k}2^{2k-2n}\sum_{i=0}^{n-k}\binom{n-i}{k}\binom{2n+1}{2i}(2^{2i-1}-1)B_{2i}.$

It is interesting to contrast this result with the Gessel-Viennot identity
(see [1, Theorem 4.2]) valid on the complementary range:

$\sum_{i=0}^{\lfloor\frac{k-1}{2}\rfloor}\binom{2k-2i-1}{k}\binom{2n+1}{2i+1}B_{2n-2i}=\frac{2n+1}{2}\binom{2k-2n}{k},\quad
k>n.$ (5)

Note that the right-hand side of equation (5) is zero unless $k\geq 2n$.

##  2 Symmetric Functions

We think of $\operatorname{Sym}$ as the subring of
$\mathbf{Q}[[x_{1},x_{2},\dots]]$ consisting of those formal power series of
bounded degree that are invariant under permutations of the $x_{i}$. A useful
reference is the first chapter of Macdonald [8]. We denote the elementary,
complete, and power-sum symmetric functions of degree $i$ by $e_{i}$, $h_{i}$,
and $p_{i}$ respectively. They have associated generating functions

$\displaystyle E(t)$
$\displaystyle=\sum_{j=0}^{\infty}e_{j}t^{j}=\prod_{i=1}^{\infty}(1+tx_{i})$
$\displaystyle H(t)$
$\displaystyle=\sum_{j=0}^{\infty}h_{j}t^{j}=\prod_{i=1}^{\infty}\frac{1}{1-tx_{i}}=E(-t)^{-1}$
$\displaystyle P(t)$
$\displaystyle=\sum_{j=1}^{\infty}p_{j}t^{j-1}=\sum_{i=1}^{\infty}\frac{x_{i}}{1-tx_{i}}=\frac{H^{\prime}(t)}{H(t)}.$

As explained in [5] and in greater detail in [6], there is a homomorphism
$\zeta:\operatorname{Sym}^{0}\to\mathbf{R}$, where $\operatorname{Sym}^{0}$ is
the subalgebra of $\operatorname{Sym}$ generated by $p_{2},p_{3},p_{4},\dots$,
such that $\zeta(p_{i})$ is the value $\zeta(i)$ of the Riemann zeta function
at $i$, for $i\geq 2$ (in [5, 6] this homomorphism is extended to all of
$\operatorname{Sym}$, but we do not need the extension here). Let
$\mathcal{D}:\operatorname{Sym}\to\operatorname{Sym}$ be the degree-doubling
map that sends $x_{i}$ to $x_{i}^{2}$. Then
$\mathcal{D}(\operatorname{Sym})\subset\operatorname{Sym}^{0}$, so the
composition $\mathfrak{Z}=\zeta\mathcal{D}$ is defined on all of
$\operatorname{Sym}$. (Alternatively, we can simply think of $\mathfrak{Z}$ as
sending $x_{i}$ to $1/i^{2}$: see [8, Ch. I, §2, ex. 21].) Note that
$\mathfrak{Z}(p_{i})=\zeta(2i)\in\mathbf{R}$. Further, $\mathfrak{Z}$ sends
the monomial symmetric function $m_{i_{1},\dots,i_{k}}$ to the symmetrized sum
of multiple zeta values

$\frac{1}{|\operatorname{Iso}(i_{1},\dots,i_{k})|}\sum_{\sigma\in
S_{k}}\zeta(2i_{\sigma(1)},2i_{\sigma(2)},\dots,2i_{\sigma(k)}),$

where $S_{k}$ is the symmetric group on $k$ letters and
$\operatorname{Iso}(i_{1},\dots,i_{k})$ is the subgroup of $S_{k}$ that fixes
$(i_{1},\dots,i_{k})$ under the obvious action.

Now let $N_{n,k}$ be the sum of all the monomial symmetric functions
corresponding to partitions of $n$ having length $k$. Of course $N_{n,k}=0$
unless $k\leq n$, and $N_{k,k}=e_{k}$. Then $\mathfrak{Z}$ sends $N_{n,k}$ to
$E(2n,k)$. Also, if we define (as in the introduction)

$\mathcal{F}(t,s)=1+\sum_{n\geq k\geq 1}N_{n,k}t^{n}s^{k},$

then $\mathfrak{Z}$ sends $\mathcal{F}(t,s)$ to the generating function
$F(t,s)$. We have the following simple description of $\mathcal{F}(t,s)$.

######  Lemma 1.

$\mathcal{F}(t,s)=E((s-1)t)H(t)$.

###### Proof.

Evidently $\mathcal{F}(t,s)$ has the formal factorization

$\prod_{i=1}^{\infty}(1+stx_{i}+st^{2}x_{i}^{2}+\cdots)=\prod_{i=1}^{\infty}\frac{1+(s-1)tx_{i}}{1-tx_{i}}=E((s-1)t)H(t).$

∎

###### Proof of Theorem 2.

Using the well-known formula for $\zeta(2,2,\dots,2)$ [4, Cor. 2.3],

$\mathfrak{Z}(e_{n})=\zeta(\underbrace{2,2,\dots,2}_{n})=\frac{\pi^{2n}}{(2n+1)!}.$
(6)

Hence

$\mathfrak{Z}(E(t))=\frac{\sinh(\pi\sqrt{t})}{\pi\sqrt{t}},$

and the image of $H(t)=E(-t)^{-1}$ is

$\mathfrak{Z}(H(t))=\frac{\pi\sqrt{-t}}{\sinh(\pi\sqrt{-t})}=\frac{\pi\sqrt{t}}{\sin(\pi\sqrt{t})}.$

Thus from Lemma 1 $F(t,s)=\mathfrak{Z}(\mathcal{F}(t,s))$ is

$\mathfrak{Z}(E((s-1)t)H(t))=\frac{\sinh(\pi\sqrt{(s-1)t})}{\pi\sqrt{(s-1)t}}\frac{\pi\sqrt{t}}{\sin(\pi\sqrt{t})}=\frac{\sin(\pi\sqrt{(1-s)t})}{\sqrt{1-s}\sin(\pi\sqrt{t})}.$

∎

Taking limits as $s\to 1$ in Theorem 2, we obtain

$F(t,1)=\frac{\pi\sqrt{t}}{\sin\pi\sqrt{t}}$

and so, taking the coefficient of $t^{n}$, the following result.

######  Corollary 1.

For all $n\geq 1$,

$\sum_{k=1}^{n}E(2n,k)=\frac{2(2^{2n-1}-1)(-1)^{n-1}B_{2n}\pi^{2n}}{(2n)!}.$

Another consequence of Lemma 1 is the following partial differential equation.

######  Proposition 1.

$t\frac{\partial{\mathcal{F}}}{\partial{t}}(t,s)+(1-s)\frac{\partial{\mathcal{F}}}{\partial{s}}(t,s)=tP(t)\mathcal{F}(t,s).$

###### Proof.

From Lemma 1 we have

$\displaystyle\frac{\partial{\mathcal{F}}}{\partial{t}}(t,s)$
$\displaystyle=(s-1)E^{\prime}((s-1)t)H(t)+E((s-1)t)H^{\prime}(t)$
$\displaystyle\frac{\partial{\mathcal{F}}}{\partial{s}}(t,s)$
$\displaystyle=tE^{\prime}((s-1)t)H(t)$

from which the conclusion follows. ∎

Now examine the coefficient of $t^{n}s^{k}$ in Proposition 1 to get the
following.

######  Proposition 2.

For $n\geq k+1$,

$p_{1}N_{n-1,k}+p_{2}N_{n-2,k}+\dots+p_{n-k}N_{k,k}=(n-k)N_{n,k}+(k+1)N_{n,k+1}.$

It is also possible to prove this result directly via a counting argument like
that used to prove the lemma of [6, p. 16].

The preceding result allows us to write $N_{n,k}$ explicitly in terms of
complete and elementary symmetric functions as follows.

######  Lemma 2.

For $r\geq 0$,

$N_{k+r,k}=\sum_{i=0}^{r}(-1)^{i}\binom{k+i}{i}h_{r-i}e_{k+i}.$

###### Proof.

We use induction on $r$, the result being evident for $r=0$. Proposition 2
gives

$\sum_{i=1}^{r+1}p_{i}N_{k+r+1-i,k}=(r+1)N_{k+r+1,k}+(k+1)N_{k+r+1,k+1},$

which after application of the induction hypothesis becomes

$\sum_{i=1}^{r+1}\sum_{j=0}^{r+1-j}(-1)^{j}p_{i}\binom{k+j}{j}h_{r+1-i-j}N_{k+j,k+j}=\\\
(r+1)N_{k+r+1,k}+(k+1)\sum_{j=0}^{r}\binom{k+1+j}{j}h_{r-j}N_{k+1+j,k+1+j}.$

The latter equation can be rewritten

$\sum_{j=0}^{r}(-1)^{j}\binom{k+j}{j}N_{k+j,k+j}\sum_{i=1}^{r+1-j}p_{i}h_{r+1-i-j}=\\\
(r+1)N_{k+r+1,k}-(k+1)\sum_{j=1}^{r+1}(-1)^{j}\binom{k+j}{j-1}h_{r+1-j}N_{k+j,k+j}.$

Now the inner sum on the left-hand side is $(r+1-j)h_{r+1-j}$ by the
recurrence relating the complete and power-sum symmetric functions, so we have

$(r+1)N_{k+r+1,k}-(r+1)N_{k,k}h_{r+1}=\\\
\sum_{j=1}^{r+1}(-1)^{j}h_{r+1-j}N_{k+j,k+j}\left((r+1-j)\binom{k+j}{j}+(k+1)\binom{k+j}{j-1}\right),$

and the conclusion follows after the observation that
$(k+1)\binom{k+j}{j-1}=j\binom{k+j}{j}$. ∎

###### Proof of Theorem 3.

Rewrite Lemma 2 in the form

$N_{n,k}=\sum_{i=0}^{n-k}\binom{n-i}{k}(-1)^{n-k-i}h_{i}e_{n-i}$

and apply the homomorphism $\mathfrak{Z}$, using equation (6) and

$\mathfrak{Z}(h_{i})=\frac{2(2^{2i-1}-1)(-1)^{i-1}B_{2i}\pi^{2i}}{(2i)!}.$

∎

##  3 Proof of Theorems 1 and 4

From the introduction we recall the statement of Theorem 1:

$E(2n,k)=\frac{1}{2^{2(k-1)}}\binom{2k-1}{k}\zeta(2n)\\\
-\sum_{j=1}^{\lfloor\frac{k-1}{2}\rfloor}\frac{1}{2^{2k-3}(2j+1)B_{2j}}\binom{2k-2j-1}{k}\zeta(2j)\zeta(2n-2j).$

We note that Euler’s formula (2) can be used to write the result in the
alternative form

$E(2n,k)=\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-1)^{j}\pi^{2j}\zeta(2n-2j)}{2^{2k-2j-2}(2j+1)!}\binom{2k-2j-1}{k}$
(7)

which avoids mention of Bernoulli numbers.

We now expand out the generating function $F(t,s)$. We have

$F(t,s)=\frac{1}{\sqrt{1-s}\sin\pi\sqrt{t}}\sin(\pi\sqrt{t}\sqrt{1-s})\\\
=\frac{\pi\sqrt{t}}{\sin\pi\sqrt{t}}\sum_{j=0}^{\infty}\frac{(-1)^{j}\pi^{2j}t^{j}(1-s)^{j}}{(2j+1)!}=\sum_{k=0}^{\infty}s^{k}G_{k}(t),$

where

$G_{k}(t)=(-1)^{k}\frac{\pi\sqrt{t}}{\sin\pi\sqrt{t}}\sum_{j\geq
k}\frac{(-1)^{j}\pi^{2j}t^{j}}{(2j+1)!}\binom{j}{k}.$ (8)

Then Theorem 1 is equivalent to the statement that

$G_{k}(t)=\sum_{n\geq
k}t^{n}\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-1)^{j}\pi^{2j}\zeta(2n-2j)}{2^{2k-2j-2}(2j+1)!}\binom{2k-2j-1}{k}$

for all $k$. We can write the latter sum as

$\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-2}(2j+1)!}\binom{2k-2j-1}{k}\sum_{n\geq
j+1}\zeta(2n-2j)t^{n-j}-\\\
\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-2}(2j+1)!}\binom{2k-2j-1}{k}\sum_{n=j+1}^{k-1}\zeta(2n-2j)t^{n-j}=\\\
\frac{1}{2}(1-\pi\sqrt{t}\cot\pi\sqrt{t})\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-2}(2j+1)!}\binom{2k-2j-1}{k}-\\\
\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-2}(2j+1)!}\binom{2k-2j-1}{k}\sum_{n=j+1}^{k-1}\zeta(2n-2j)t^{n-j},$
(9)

where we have used the generating function

$\frac{1}{2}(1-\pi\sqrt{t}\cot\pi\sqrt{t})=\sum_{i=1}^{\infty}\zeta(2i)t^{i}.$

Note that the last sum in (9) is a polynomial that cancels exactly those terms
in

$\frac{1}{2}(1-\pi\sqrt{t}\cot\pi\sqrt{t})\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-2}(2j+1)!}\binom{2k-2j-1}{k}$
(10)

of degree less than $k$. Thus, to prove Theorem 1 it suffices to show that

$G_{k}(t)=\text{terms of degree $\geq k$ in expression (\ref{expr}).}$

From equation (8) it is evident that

$G_{k}(t)=\frac{\pi\sqrt{t}}{\sin\pi\sqrt{t}}\cdot\frac{(-t)^{k}}{k!}\cdot\frac{d^{k}}{dt^{k}}\left(\frac{\sin\pi\sqrt{t}}{\pi\sqrt{t}}\right).$
(11)

We use this to obtain an explicit formula for $G_{k}(t)$.

######  Lemma 3.

For $k\geq 0$,

$G_{k}(t)=P_{k}(\pi^{2}t)\pi\sqrt{t}\cot\pi\sqrt{t}+Q_{k}(\pi^{2}t),$

where $P_{k},Q_{k}$ are polynomials defined by

$\displaystyle P_{k}(x)=$
$\displaystyle-\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4x)^{j}}{2^{2k-1}(2j+1)!}\binom{2k-2j-1}{k}$
$\displaystyle Q_{k}(x)=$
$\displaystyle\sum_{j=0}^{\lfloor\frac{k}{2}\rfloor}\frac{(-4x)^{j}}{2^{2k}(2j)!}\binom{2k-2j}{k}.$

###### Proof.

In view of equation (11), the conclusion is equivalent to

$f^{(k)}(t)=(-1)^{k}k!t^{-k}P_{k}(\pi^{2}t)\cos\pi\sqrt{t}+(-1)^{k}k!t^{-k}Q_{k}(\pi^{2}t)f(t),$

where $f(t)=\sin\pi\sqrt{t}/\pi\sqrt{t}$. Differentiating, one sees that the
polynomials $P_{k}$ and $Q_{k}$ are determined by the recurrence

$\displaystyle(k+1)P_{k+1}(x)$
$\displaystyle=kP_{k}(x)-xP_{k}^{\prime}(x)-\frac{1}{2}Q_{k}(x)$
$\displaystyle(k+1)Q_{k+1}(x)$
$\displaystyle=\frac{2k+1}{2}Q_{k}(x)-xQ_{k}^{\prime}(x)+\frac{x}{2}P_{k}(x)$

together with the initial conditions $P_{0}(x)=0$, $Q_{0}(x)=1$. The
recurrence and initial conditions are satisfied by the explicit formulas
above. ∎

###### Proof of Theorem 1.

Using Lemma 3, we have

$G_{k}(t)=-\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-1}(2j+1)!}\binom{2k-2j-1}{k}\pi\sqrt{t}\cot\pi\sqrt{t}\\\
+\sum_{j=0}^{\lfloor\frac{k}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k}(2j)!}\binom{2k-2j}{k}=\\\
\frac{1}{2}(1-\pi\sqrt{t}\cot\pi\sqrt{t})\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-4\pi^{2}t)^{j}}{2^{2k-2}(2j+1)!}\binom{2k-2j-1}{k}\\\
+\text{terms of degree $<k$},$

and this completes the proof. ∎

###### Proof of Theorem 4.

Using Theorem 1 in the form of equation (7), eliminate $\zeta(2n-2j)$ using
Euler’s formula (2) and then compare with Theorem 3 to get

$\sum_{j=0}^{\lfloor\frac{k-1}{2}\rfloor}\frac{(-1)^{n-1}\pi^{2n}B_{2n-2j}}{2^{2k-2n-1}(2n-2j)!(2j+1)!}\binom{2k-2j-1}{k}=\\\
\frac{(-1)^{n-k-1}\pi^{2n}}{(2n+1)!}\sum_{i=0}^{n-k}\binom{n-i}{k}\binom{2n+1}{2i}2(2^{2i-1}-1)B_{2i}.$

Now multiply both sides by $(-1)^{n-1}2^{2k-2n-1}\pi^{-2n}(2n+1)!$ and rewrite
the factorials on the left-hand side as a binomial coefficient. ∎

## References

  * [1] W. Y. C. Chen and L. H. Sun, Extended Zeilberger’s algorithm for identities on Bernoulli and Euler polynomials, _J. Number Theory_ 129 (2009), 2111-2132.
  * [2] L. Euler, Meditationes circa singulare serierum genus, _Novi Comm. Acad. Sci. Petropol._ 20 (1775), 140-186; reprinted in _Opera Omnia_ , ser. I, vol. 15, B. G. Teubner, Berlin, 1927, pp. 217-267.
  * [3] H. Gangl, M. Kaneko, and D. Zagier, Double zeta values and modular forms, in _Automorphic Forms and Zeta Functions_ , S. Böcherer et. al. (eds.), World Scientific, Singapore, 2006, pp. 71-106.
  * [4] M. E. Hoffman, Multiple harmonic series, _Pacific J. Math._ 152 (1992), 275-290.
  * [5] M. E. Hoffman, The algebra of multiple harmonic series, _J. Algebra_ 194 (1997), 477-495.
  * [6] M. E. Hoffman, A character on the quasi-symmetric functions coming from multiple zeta values, _Electron. J. Combin._ 15 (2008), res. art. 97.
  * [7] Y. Komori, K. Matsumoto and H. Tsumura, A study on multiple zeta values from the viewpoint of zeta-functions of root systems, preprint arXiv:1205.0182.
  * [8] I. G. Macdonald, _Symmetric Functions and Hall Polynomials_ , 2nd ed., Oxford Univ. Press, New York, 1995.
  * [9] T. Machide, Extended double shuffle relations and the generating function of triple zeta values of any fixed weight, preprint arXiv:1204.4085.
  * [10] Z. Shen and T. Cai, Some formulas for multiple zeta values, _J. Number Theory_ 132 (2012), 314-323.


###### Abstract

Rapid and luminous flares of non-thermal radiation observed in blazars require
an efficient mechanism of energy dissipation and particle acceleration in
relativistic active galactic nuclei (AGN) jets. Particle acceleration in
relativistic magnetic reconnection is being actively studied by kinetic
numerical simulations. Relativistic reconnection produces hard power-law
electron energy distributions
$N(\gamma)\propto\gamma^{-p}\exp(-\gamma/\gamma_{\rm max})$ with index $p\to
1$ and exponential cut-off Lorentz factor $\gamma_{\rm max}\sim\sigma$ in the
limit of magnetization $\sigma=B^{2}/(4\pi w)\gg 1$ (where $w$ is the
relativistic enthalpy density). Reconnection in electron-proton plasma can
additionally boost $\gamma_{\rm max}$ by the mass ratio $m_{\rm p}/m_{\rm e}$.
Hence, in order to accelerate particles to $\gamma_{\rm max}\sim 10^{6}$ in
the case of BL Lacs, reconnection should proceed in plasma of very high
magnetization $\sigma_{\rm max}\gtrsim 10^{3}$. On the other hand, moderate
mean jet magnetization values are required for magnetic bulk acceleration of
relativistic jets, $\sigma_{\rm mean}\sim\Gamma_{\rm j}\lesssim 20$ (where
$\Gamma_{\rm j}$ is the jet bulk Lorentz factor). I propose that the
systematic dependence of $\gamma_{\rm max}$ on blazar luminosity class—the
blazar sequence—may result from a systematic trend in $\sigma_{\rm max}$ due
to homogeneous loading of leptons by pair creation regulated by the energy
density of high-energy external radiation fields. At the same time,
relativistic AGN jets should be highly inhomogeneous due to filamentary
loading of protons, which should determine the value of $\sigma_{\rm mean}$
roughly independently of the blazar class.

###### keywords:

blazars; relativistic jets; magnetic reconnection

x 10.3390/—— 4 Academic Editors: Jose L. Gómez, Alan P. Marscher and Svetlana
G. Jorstad Received: 15 July 2016; Accepted: 3 September 2016 Applying
Relativistic Reconnection to Blazar Jets Krzysztof Nalewajko Krzysztof
Nalewajko

##  1 Introduction

Blazars are persistent extragalactic sources of non-thermal radiation
extending from the radio to the gamma-ray band and characterized by stochastic
variability over a wide range of time scales. Multiwavelength observations of
blazars typically reveal two major broad spectral components, with the low-
energy one (radio—UV/X-ray) interpreted universally as synchrotron radiation
of electrons. High-resolution interferometric radio/mm observations reveal a
core-jet structure with individual jet substructures propagating with
apparently superluminal velocity Jorstad et al. (2001); Lister et al. (2016).
Blazars are associated with active galactic nuclei (AGN) equipped with jets
(radio loud; Urry & Padovani (1995)), with one of the jets pointing closely at
the observer Blandford & Königl (1979), which leads to a dramatic relativistic
boost of apparent luminosity Rees (1966). Blazars are typically classified
into more luminous flat-spectrum radio quasars (FSRQs) and less luminous BL
Lac objects (BL Lacs). An anticorrelation between the radio luminosity and the
synchrotron peak frequency is known as the blazar sequence Fossati et al.
(1998).

High apparent luminosities of blazars, up to $L_{\gamma}\sim 10^{50}\;{\rm
erg\cdot s^{-1}}$ in the $\gamma$-ray band Abdo et al. (2011), require
efficient in situ dissipation of the jet power, efficient particle
acceleration and efficient radiative mechanisms. Modeling of the spectral
energy distributions of blazars can be performed in two basic scenarios
Böttcher et al. (2013): in the leptonic scenario, the high-energy spectral
component is interpreted as inverse Compton scattering of soft radiation
fields by energetic electrons; in the hadronic scenario, it is interpreted as
due to various mechanisms involving relativistic protons. In luminous blazars,
the requirement of high radiative efficiency favors the leptonic scenario
Sikora et al. (2009). However, in any case, the inferred characteristic
energies of electrons producing the low-energy synchrotron component are very
different in FSRQs ($\gamma_{\rm max}\sim 10^{3}$) and in BL Lacs
($\gamma_{\rm max}\sim 10^{6}$).

In this work, I attempt to address several fundamental questions about the
physics of relativistic AGN jets.

  1. (1)

What is the origin of the blazar sequence? What determines $\gamma_{\rm max}$?

  2. (2)

How is matter introduced into relativistic magnetized jets? What determines
$\Gamma_{\rm j}$?

  3. (3)

Can relativistic magnetic reconnection explain energy dissipation and particle
acceleration in blazar jets?

In Section 2, I summarize the current understanding of particle acceleration
in relativistic magnetic reconnection. In Section 3, I sketch two independent
mechanisms of mass loading of relativistic jets by leptons and protons. In
Section 4, I discuss why the blazar sequence is unlikely to be regulated by
radiative cooling rates. In Section 5, I propose a novel picture of the
composition of relativistic jets and an alternative explanation of the blazar
sequence.

##  2 Particle Acceleration in Relativistic Magnetic Reconnection

Rapid progress was achieved recently in understanding particle acceleration
during relativistic magnetic reconnection thanks to large-scale kinetic
particle-in-cell simulations. It is now established that relativistic
reconnection starting from the tearing instability of Harris-type current
layers produces power-law energy distributions with exponential cut-off
$N(\gamma)\propto\gamma^{-p}\exp(-\gamma/\gamma_{\rm max})$ with index $p$
decreasing with increasing background (upstream) plasma magnetization $\sigma$
and approaching $p\to 1$ in the limit of $\sigma\gg 1$ Sironi & Spitkovsky
(2014); Guo et al. (2014); Werner et al. (2016). In this limit, most energy is
contained in the highest-energy particles, and hence the maximum particle
energy $\gamma_{\rm max}\sim\sigma$. The high-energy cut-off in the particle
energy distribution is exponential for sufficiently large simulation domains
Werner et al. (2016), and no soft power-law tails that could extend
$\gamma_{\rm max}$ beyond the cut-off were observed.

The acceleration mechanism was described as first-order Fermi process, in
which more energetic particles achieve larger energy gains by interacting with
larger plasmoids Guo et al. (2014). A theoretical explanation of the $p=1$
electron distribution index was proposed Guo et al. (2014). Particles can gain
energy at three location types: at magnetic X-points (reconnection electric
field), between merging plasmoids (anti-reconnection electric field), and
within accelerating plasmoids (curvature drift) Nalewajko et al. (2015).
Energy gain within contracting plasmoids Drake et al. (2006) appears to be
limited to non-relativistic reconnection. In Harris-layer reconnection,
curvature drift has been shown to dominate global energy gain Guo et al.
(2014), although the highest energy particles were found to have passed
through a major X-point Sironi & Spitkovsky (2014). Under strong radiative
cooling, the relative importance of acceleration sites changes, with rapid
acceleration by strong reconnection electric fields at magnetic X-points
allowing the synchrotron photon energy to exceed the radiation reaction upper
limit Cerutti et al. (2013).

Relativistic reconnection in electron-ion plasma was investigated in two
regimes. In the case when electrons are relativistic and ions are non-
relativistic, ions gain up to 2/3 of the total dissipated energy, however
their energy distributions are much softer than those of the electrons Melzani
et al. (2014). In the case when both electrons and ions are relativistic, they
obtain similar energy distributions with hard power-law spectra Guo et al.
(2016), hence this case is qualitatively similar to reconnection in
relativistic pair plasma. In any case involving relativistic electrons, their
maximum Lorentz factor is given by $\gamma_{\rm max}\sim\sigma_{\rm
e}\sim(n_{\rm p}m_{\rm p}/n_{\rm e}m_{\rm e})\sigma_{\rm p}\lesssim
10^{3}\sigma_{\rm p}$ (Werner et al., in preparation).

An alternative numerical setup was investigated, starting from smooth
magnetostatic equilibria dubbed “ABC” fields Nalewajko et al. (2016); Yuan et
al. (2016); Lyutikov et al. (2016). Such an initial condition allows one to
study the formation and dynamical evolution of the current layer and
simultaneous particle acceleration and radiation. In ABC reconnection, energy
gain by the reconnection electric field appears to be more important.

##  3 Mass Loading in Relativistic Jets

Energy can be extracted from rotating black holes by purely electromagnetic
jets Blandford & Znajek (1977), hence, there is no fundamental principle for
the presence of matter at the base of jets. Without introducing matter from
outside, there is no fundamental upper limit on the value of magnetization
$\sigma$. Here I consider mass loading by pairs and protons as independent
mechanisms.

Mass loading by pairs can be due to photon-photon pair creation. This requires
the presence of roughly isotropic soft gamma-ray radiation or interaction of
roughly isotropic X-ray radiation with relativistically boosted X-ray
radiation of the jet Sikora & Madejski (2000). If this mechanism can operate
efficiently in the first place, it can be expected to operate uniformly across
the jet volume, providing a lower limit on plasma density, and hence an upper
limit on magnetization $\sigma$. Since this mechanism is very sensitive to the
energy density of external high-energy radiation fields, the upper limit on
$\sigma$ can be expected to depend systematically on blazar type. In the case
of FSRQs, we expect the pair creation to be more efficient, and the maximum
magnetization value to be lower, as compared with the BL Lacs.

On the other hand, mass loading by protons can be due to penetration of the
jet by gas clouds or even stars Barkov et al. (2010) or due to various plasma
instabilities, e.g., the magnetic Kelvin-Helmholtz instability or magnetic
Rayleigh-Taylor instability McKinney et al. (2012). These processes can be
expected to be highly non-uniform, leading to high density contrasts with most
of the protons concentrated in compact clouds/filaments/comae/tails. We may
also expect a significant radial stratification of the jets, with the inner
cores essentially free of protons, unless the entire jet is subject to a
global current driven instability. This will likely result in high contrasts
of magnetization $\sigma$. The proton loading mechanisms can be expected to be
independent of the radiative environment of the jet. Whether they should
depend on the jet power or the mechanical jet environment, is an interesting
open question.

##  4 On the Origin of Blazar Sequence

The blazar sequence is an observational correlation between the frequency of
the spectral peak of the synchrotron component $\nu_{\rm syn}$ and the
synchrotron luminosity measured in the radio band Fossati et al. (1998). It
should be noted that the observational evidence for the reality of blazar
sequence is being questioned as due to selection effects Giommi et al. (2012).
Nevertheless, for individual blazars it is quite striking that despite
dramatic variability in all bands their broad-band SEDs remain stable. Given
similar values of the co-moving magnetic field strength $B^{\prime}$ and the
jet bulk Lorentz factor $\Gamma_{\rm j}$, the synchrotron peak frequency is
determined by the characteristic energy of radiating electrons $\nu_{\rm
syn}\propto\gamma_{\rm peak}^{2}$. By modeling the spectral energy
distributions of blazars belonging to different classes, a theoretical
anticorrelation between $\gamma_{\rm peak}$ and the combined co-moving energy
density of magnetic fields and external radiation fields $u_{\rm
cool}^{\prime}=u_{\rm B}^{\prime}+u_{\rm ext}^{\prime}$ was demonstrated
Ghisellini et al. (1998); Celotti & Ghisellini (2008): $\gamma_{\rm peak}\sim
10^{3}(u_{\rm cool,0}^{\prime})^{-0.5}$, where $u_{\rm
cool,0}^{\prime}~{}=~{}u_{\rm cool}^{\prime}/(1\;{\rm erg\cdot cm^{-3}})$. It
was then proposed that the blazar sequence arises from acceleration of
electrons being controlled by their radiative cooling rate $t_{\rm
acc}^{\prime}\sim t_{\rm cool}^{\prime}$.

Acceleration of electrons under radiative cooling is governed by the following
equation:

$\frac{{\rm d}\gamma}{{\rm d}t^{\prime}}=\frac{eE^{\prime}}{m_{\rm
e}c}-\frac{4\sigma_{\rm T}u_{\rm cool}^{\prime}\gamma^{2}}{3m_{\rm e}c}\,,$
(1)

where $E^{\prime}$ is the effective electric field parallel to the
instantaneous particle velocity vector. In the particular case of relativistic
magnetic reconnection, but also in any other acceleration mechanism, we can
relate $E^{\prime}$ to the local magnetic field strength as $E^{\prime}=\beta
B^{\prime}$. In the case of reconnection, $\beta\equiv\beta_{\rm rec}$ is the
reconnection rate, and in other cases it can be interpreted as a
characteristic relative velocity of interacting magnetized flows. The
equilibrium electron energy is found simply from ${\rm d}\gamma/{\rm
d}t^{\prime}=0$ to be

$\gamma_{\rm eq}=\sqrt{\frac{3e\beta B^{\prime}}{4\sigma_{\rm T}u_{\rm
cool}^{\prime}}}\simeq 2\times 10^{7}\sqrt{\frac{\beta B_{0}^{\prime}}{u_{\rm
cool,0}^{\prime}}}\,,$ (2)

where $B_{0}^{\prime}=B^{\prime}/(1\;{\rm G})$. This relation predicts that
$\gamma_{\rm eq}\propto u_{\rm cool}^{\prime-0.5}$, but only when radiative
cooling is dominated by Comptonization of the external radiation and when
$B_{0}^{\prime}\sim{\rm const}$. However, in order to match the empirical
relation between $\gamma_{\rm peak}$ and $u_{\rm cool}^{\prime}$, with
$B_{0}^{\prime}\sim 1$ in parsec-scale jets Pushkarev et al. (2012), one would
require $\beta\sim 10^{-9}$. Even considering the maximum energy of the
particle distributions obtained from SED modeling $\gamma_{\rm max}\sim
10\gamma_{\rm peak}$, one would need $\beta\sim 10^{-7}$. In other words, _for
the maximum electron energies to be determined by radiative losses, the
acceleration mechanism should be extremely slow._ Relativistic magnetic
reconnection, with $\beta_{\rm rec}\sim 0.1$ Liu et al. (2015), should easily
accelerate electrons to the energies of order $\gamma_{\rm max}\sim 10^{7}$.
In such a case, the SEDs of blazars would be more similar to that of the Crab
Nebula with the synchrotron component extending up to $\sim$100 ${\rm MeV}$ de
Jager et al. (1996). An analogous problem in the context of particle
acceleration by relativistic shocks was identified by Inoue & Takahara (1996).

##  5 Proposition

Here I propose a qualitatively new picture of blazars and relativistic jets.
At given distance scale, relativistic jets are composed of roughly homogeneous
magnetic energy density, roughly homogeneous pair density, and highly
inhomogeneous proton density. The magnetic field and the proton components are
roughly independent of the blazar class, and the pair component scales with
the high-energy external radiation density $u_{\rm ext,he}$. The pair density
determines the local maxima of the jet magnetization $\sigma_{\rm max}$
decreasing with increasing $u_{\rm ext,he}$, and the proton density determines
the mean jet magnetization $\sigma_{\rm mean}$ roughly independent of $u_{\rm
ext,he}$. The maximum jet magnetization determines the maximum energy of
electrons $\gamma_{\rm max}$, potentially explaining the blazar sequence. The
mean jet magnetization determines the jet bulk Lorentz factor $\Gamma_{\rm
j}$, which depends only very weakly on the blazar class.

This model provides a suitable background for relativistic reconnection as the
main dissipation mechanism in blazars. Without the existence of very highly
magnetized regions, with $\sigma_{\rm max}\sim 10^{3}$, acceleration of
electrons at reconnection sites cannot produce the energy distributions
required for TeV BL Lacs. The existence of very highly magnetized regions in
relativistic jets is not inconsistent with the dominant emitting regions being
roughly in equipartition, since in relativistic reconnection the emitting
regions are in general expected to be in sharp contrast to the main
acceleration regions Sironi et al. (2015). In any case, blazar emission
produced in highly magnetized relativistic jets would be expected to be
characterized by very high synchrotron luminosity and low Compton dominance
Janiak et al. (2015).

##  6 Conclusions

Relativistic reconnection is a promising dissipation mechanism in relativistic
jets. Rapid progress in understanding particle acceleration during
relativistic reconnection has been achieved recently by means of kinetic
numerical simulations. We have learned that relativistic reconnection can
accelerate particles very efficiently, producing hard power-law distributions
with index $p=1$ and exponential cut-off at $\gamma_{\rm max}\sim\sigma_{\rm
e}$. Applying this mechanism to blazars requires very high jet magnetization
values, $\sigma_{\rm max}\sim 10^{3}$ in the case of TeV BL Lacs. I suggest
that such high magnetizations may be present locally in relativistic jets due
to highly inhomogeneous mass loading of protons by plasma instabilities. Fast
reconnection rates $\beta_{\rm rec}\sim 0.1$ can easily compete with radiative
cooling rates, allowing in principle acceleration of particles to much higher
energies than required by SED modeling, especially in the case of FSRQs.
Therefore, I argue that the blazar sequence can hardly be regulated by
radiative cooling rates. Instead, I suggest that the blazar sequence arises
due to homogeneous loading of leptons by pair creation regulated by external
radiation fields. Several aspects of the proposed picture of blazar physics
require detailed theoretical or numerical verification.

###### Acknowledgements.

I thank Marek Sikora and Mitch Begelman for stimulating discussions. This work
was supported by the Polish National Science Centre grant 2015/18/E/ST9/00580.
The author declares no conflict of interest.

## References

  * Jorstad et al. (2001) Jorstad, S.G.; Marscher, A.P.; Mattox, J.R.; Wehrle, A.E.; Bloom, S.D.; Yurchenko, A.V. Multiepoch very long baseline array observations of EGRET-detected quasars and BL lacertae objects: Superluminal motion of gamma-ray bright blazars. Astrophys. J. Suppl.  2001, 134, 181–240.
  * Lister et al. (2016) Lister, M.L.; Aller, M.F.; Aller, H.D.; Homan, D.C.; Kellermann, K.I.; Kovalev, Y.Y.; Pushkarev, A.B.; Richards, J.L.; Ros, E.; Savolainen, T. MOJAVE: XIII. Parsec-scale AGN jet kinematics analysis based on 19 years of VLBA observations at 15 GHz. Astron. J. 2016, 152, 12.
  * Urry & Padovani (1995) Urry, C.M.; Padovani, P. Unified schemes for radio-loud active galactic nuclei. Publ. Astron. Soc. Pac. 1995, 107, 803–845.
  * Blandford & Königl (1979) Blandford, R.D.; Königl, A. Relativistic jets as compact radio sources. Astrophys. J. 1979, 232, 34–48.
  * Rees (1966) Rees, M.J. Appearance of relativistically expanding radio sources. Nature 1966, 211, 468–470.
  * Fossati et al. (1998) Fossati, G.; Maraschi, L.; Celotti, A.; Comastri, A.; Ghisellini, G. A unifying view of the spectral energy distributions of blazars. Mon. Not. R. Astron. Soc. 1998, 299, 433–448.
  * Abdo et al. (2011) Abdo, A.A.; Ackermann, M.; Ajello, M.; Allafort, A.; Baldini, L.; Ballet, J.; Barbiellini, G.; Bastieri, D.; Bellazzini, R.; Berenji, B.; et al. Fermi gamma-ray space telescope observations of the gamma-ray outburst from 3C454.3 in November 2010. Astrophys. J. 2011, 733, L26.
  * Böttcher et al. (2013) Böttcher, M.; Reimer, A.; Sweeney, K.; Prakash, A. Leptonic and hadronic modeling of fermi-detected blazars. Astrophys. J. 2013, 768, 54.
  * Sikora et al. (2009) Sikora, M.; Stawarz, Ł.; Moderski, R.; Nalewajko, K.; Madejski, G.M. Constraining emission models of luminous blazar sources. Astrophys. J. 2009, 704, 38–50.
  * Guo et al. (2014) Guo, F.; Li, H.; Daughton, W.; Liu, Y.-H. Formation of hard power laws in the energetic particle spectra resulting from relativistic magnetic reconnection. Phys. Rev. Lett. 2014, 113, 155005.
  * Sironi & Spitkovsky (2014) Sironi, L.; Spitkovsky, A. Relativistic reconnection: An efficient source of non-thermal particles. Astrophys. J. 2014, 783, L21.
  * Werner et al. (2016) Werner, G.R.; Uzdensky, D.A.; Cerutti, B.; Nalewajko, K.; Begelman, M.C. The extent of power-law energy spectra in collisionless relativistic magnetic reconnection in pair plasmas. Astrophys. J. 2016, 816, L8.
  * Nalewajko et al. (2015) Nalewajko, K.; Uzdensky, D.A.; Cerutti, B.; Werner, G.R.; Begelman, M.C. On the distribution of particle acceleration sites in plasmoid-dominated relativistic magnetic reconnection. Astrophys. J. 2015, 815, 101.
  * Drake et al. (2006) Drake, J.F.; Swisdak, M.; Che, H.; Shay, M.A. Electron acceleration from contracting magnetic islands during reconnection. Nature 2006, 443, 553–556.
  * Cerutti et al. (2013) Cerutti, B.; Werner, G.R.; Uzdensky, D.A.; Begelman, M.C. Simulations of particle acceleration beyond the classical synchrotron burnoff limit in magnetic reconnection: An explanation of the crab flares. Astrophys. J. 2013, 770, 147.
  * Melzani et al. (2014) Melzani, M.; Walder, R.; Folini, D.; Winisdoerffer, C.; Favre, J.M. Relativistic magnetic reconnection in collisionless ion-electron plasmas explored with particle-in-cell simulations. Astron. Astrophys. 2014, 570, A111.
  * Guo et al. (2016) Guo, F.; Li, X.; Li, H.; Daughton, W.; Zhang, B.; Lloyd-Ronning, N.; Liu, Y.-H.; Zhang, H.; Deng, W. Efficient production of high-energy nonthermal particles during magnetic reconnection in a magnetically dominated ion-electron plasma. Astrophys. J. 2016, 818, L9.
  * Lyutikov et al. (2016) Lyutikov, M.; Sironi, L.; Komissarov, S.; Porth, O. Particle acceleration in explosive relativistic reconnection events and crab nebula gamma-ray flares. 2016, arXiv:1603.05731.
  * Nalewajko et al. (2016) Nalewajko, K.; Zrake, J.; Yuan, Y.; East, W.E.; Blandford, R.D. Kinetic simulations of the lowest-order unstable mode of relativistic magnetostatic equilibria. Astrophys. J. 2016, 826, 115.
  * Yuan et al. (2016) Yuan, Y.; Nalewajko, K.; Zrake, J.; East, W.E.; Blandford, R.D. Kinetic study of radiation-reaction-limited particle acceleration during the relaxation of unstable force-free equilibria. 2016, arXiv:1604.03179.
  * Blandford & Znajek (1977) Blandford, R.D.; Znajek, R.L. Electromagnetic extraction of energy from Kerr black holes. Mon. Not. R. Astron. Soc. 1977, 179, 433–456.
  * Sikora & Madejski (2000) Sikora, M.; Madejski, G. On pair content and variability of subparsec jets in quasars. Astrophys. J. 2000, 534, 109–113.
  * Barkov et al. (2010) Barkov, M.V.; Aharonian, F.A.; Bosch-Ramon, V. Gamma-ray flares from red giant/jet interactions in active galactic nuclei. Astrophys. J. 2010, 724, 1517–1523.
  * McKinney et al. (2012) McKinney, J.C.; Tchekhovskoy, A.; Blandford, R.D. General relativistic magnetohydrodynamic simulations of magnetically choked accretion flows around black holes. Mon. Not. R. Astron. Soc. 2012, 423, 3083–3117.
  * Giommi et al. (2012) Giommi, P.; Padovani, P.; Polenta, G.; Turriziani, S.; D’Elia, V.; Piranomonte, S. A simplified view of blazars: Clearing the fog around long-standing selection effects. Mon. Not. R. Astron. Soc. 2012, 420, 2899–2911.
  * Celotti & Ghisellini (2008) Celotti, A.; Ghisellini, G. The power of blazar jets. Mon. Not. R. Astron. Soc. 2008, 385, 283–300.
  * Ghisellini et al. (1998) Ghisellini, G.; Celotti, A.; Fossati, G.; Maraschi, L.; Comastri, A. A theoretical unifying scheme for gamma-ray bright blazars. Mon. Not. R. Astron. Soc. 1998, 301, 451–468.
  * Pushkarev et al. (2012) Pushkarev, A.B.; Hovatta, T.; Kovalev, Y.Y.; Lister, M.L.; Lobanov, A.P.; Savolainen, T.; Zensus, J.A. MOJAVE: Monitoring of jets in active galactic nuclei with VLBA experiments. IX. Nuclear opacity. Astron. Astrophys. 2012, 545, A113.
  * Liu et al. (2015) Liu, Y.-H.; Guo, F.; Daughton, W.; Li, H.; Hesse, M. Scaling of magnetic reconnection in relativistic collisionless pair plasmas. Phys. Rev. Lett. 2015, 114, 095002.
  * de Jager et al. (1996) De Jager, O.C.; Harding, A.K.; Michelson, P.F.; Nel, H.I.; Nolan, P.L.; Sreekumar, P.; Thompson, D.J. Gamma-ray observations of the crab nebula: A study of the synchro-compton spectrum. Astrophys. J. 1996, 457, 253–266.
  * Inoue & Takahara (1996) Inoue, S.; Takahara, F. Electron acceleration and gamma-ray emission from blazars. Astrophys. J. 1996, 463, 555–564.
  * Sironi et al. (2015) Sironi, L.; Petropoulou, M.; Giannios, D. Relativistic jets shine through shocks or magnetic reconnection? Mon. Not. R. Astron. Soc. 2015, 450, 183.
  * Janiak et al. (2015) Janiak, M.; Sikora, M.; Moderski, R. Magnetization of jets in luminous blazars. Mon. Not. R. Astron. Soc. 2015, 449, 431–439