# On the design of multiplex control to reject disturbances in nonlinear
network systems affected by heterogeneous delays

Shihao Xie1 and Giovanni Russo2 1Shihao Xie is with School of Electrical and
Electronic Engineering, University College Dublin, Ireland
<EMAIL_ADDRESS>Russo is with the Department of Information
and Electrical Engineering and Applied Mathematics, University of Salerno,
Italy<EMAIL_ADDRESS>
###### Abstract

We consider the problem of designing control protocols for nonlinear network
systems affected by heterogeneous, time-varying delays and disturbances. For
these networks, the goal is to reject polynomial disturbances affecting the
agents and to guarantee the fulfilment of some desired network behaviour. To
satisfy these requirements, we propose an integral control design implemented
via a multiplex architecture. We give sufficient conditions for the desired
disturbance rejection and stability properties by leveraging tools from
contraction theory. We illustrate the effectiveness of the results via a
numerical example that involves the control of a multi-terminal high-voltage
DC grid.

##  I Introduction

Over the past decades, the size and complexity of network systems have
considerably evolved thanks to the rapid development of computing and
communication technologies. Much research efforts have been devoted to study
collective behaviours such as consensus, synchronisation, formation control. A
key challenge when designing the control protocols is to achieve desired
behaviours despite imperfect communications, exogenous disturbances and
delays.

In this context, we study the problem of designing distributed integral
control protocols that guarantee the fulfilment of the desired network
behaviour, while rejecting certain classes of disturbances. These requirements
are captured via an Input-to-State Stability (ISS) property and we give
sufficient conditions for this property based on non-Euclidean contraction
theory.

#### Related works

the design of integral control protocols for network systems that are able to
reject constant disturbances has been investigated in e.g. [1, 2]. In [3], a
PI controller is delivered via multiplex architecture to achieve consensus.
Recently, in [4], integral actions delivered by multiplex architecture with
multiple layers are shown to be effective in rejecting higher order polynomial
disturbances while guaranteeing a scalability property. The results in this
paper are based on ideas from contraction theory [5], particularly leveraging
the use of non-Euclidean norms [6, 7, 8]. We refer to [9, 10, 11] and
references therein for details. In the context of delay-free networks,
leveraging contraction theory, conditions for the synthesis of distributed
controls using separable control contraction metrics are given in [12];
contracting recurrent network is introduced in [13] with guarantees of
stability and robustness. For network systems affected by delays, [14] shows
the preservation of contraction for a time-delayed network using Euclidean
contraction metric and, in [4], conditions are given for networks with
homogeneous delays.

#### Statement of Contributions

we present a distributed multiplex integral control design for nonlinear
network systems affected by both heterogeneous time-varying delays and
disturbances (possibly with polynomial components). The goal of the control
protocol is to guarantee, for the network: (i) rejection of polynomial
disturbances; (ii) the fulfilment of some desired behaviours. These properties
are rigorously formalised in Section III. Specifically, our technical
contributions can be summarised as follows: (i) we formalise the control
problem as an Input-Output Stability problem and give sufficient conditions to
assess this property. While the results of this paper leverage some of the
tools from [4], here, differently from [4], we consider a weaker stability
property for networks with heterogeneous delays. In [4], although a stronger
scalability property was considered, results were obtained by assuming that
the network has homogeneous delays; (ii) we show that our results can serve as
design guidelines for the control protocol; (iii) the results are validated on
the problem of designing a control protocol for multi-terminal high-voltage DC
(MTDC) grid and we show how the conditions can be fulfilled by solving a
convex optimisation problem. Simulations confirm the effectiveness of the
results.

##  II Mathematical preliminaries

Given a norm $\|\cdot\|$, we denote by $\|A\|$ its induced matrix norm with
respect to a $m\times m$ real matrix $A$ and the corresponding matrix measure
$\mu(A)=\lim_{h\rightarrow 0^{+}}\frac{\|I+hA\|-1}{h}$. The symmetric part of
a matrix $A$ is denoted as $[A]_{s}:=\frac{A+A^{\mathsf{T}}}{2}$. For
$\eta\in\mathbb{R}_{>0}^{n}$, we define a diagonal matrix by
$[\eta]\in\mathbb{R}^{n\times n}$ with $[\eta]_{ii}=\eta_{i}$,
$i\in\\{1,\ldots,n\\}$. The diagonally weighted $\ell_{\infty}$-norm of
$x\in\mathbb{R}^{n}$ is defined as
$\|x\|_{\infty,[\eta]^{-1}}:=\max_{i}\\{|x_{i}|/\eta_{i}\\}$ with the induced
matrix norm
$\|A\|_{\infty,[\eta]^{-1}}:=\max_{i}\\{\sum_{j}\frac{\eta_{j}}{\eta_{i}}|A_{ij}|\\}$
and matrix measure $\mu_{\infty,[\eta]^{-1}}(A):=\max_{i}\\{A_{ii}+\sum_{j\neq
i}\frac{\eta_{j}}{\eta_{i}}|A_{ij}|\\}$. We denote by $\mathbb{I}_{n}$ the
$n\times n$ identity matrix, by $0_{m\times n}$ the $m\times n$ zero matrix
(if $m=n$ we simply write $0_{n}$) and by $\mathbf{1}_{n}\in\mathbb{R}^{n}$
the one vector. Let $f$ be a sufficiently smooth function, we denote by
$f^{(n)}$ the $n$-th derivative of $f$. We recall that a continuous function
$\alpha:[0,a)\rightarrow[0,\infty)$ is said to belong to class $\mathcal{K}$
if it is strictly increasing and $\alpha(0)=0$. It is said to belong to class
$\mathcal{K}_{\infty}$ if $a=\infty$ and $\alpha(r)\rightarrow\infty$ as
$r\rightarrow\infty$. A continuous function
$\beta:[0,a)\times[0,\infty)\rightarrow[0,\infty)$ is said to belong to class
$\mathcal{KL}$ if, for each fixed $s$, the mapping $\beta(r,s)$ belongs to
class $\mathcal{K}$ with respect to $r$ and, for each fixed $r$, the mapping
$\beta(r,s)$ is decreasing with respect to $s$ and $\beta(r,s)\rightarrow 0$
as $s\rightarrow\infty$.

The following results, originally introduced in [15], can be found in its
current form in [16].

######  Lemma 1.

Given $r$ positive integers $n_{1},\ldots,n_{r}$ such that
$n_{1}+\cdots+n_{r}=n$. Consider the vector
$x:=[x_{1}^{\mathsf{T}},\ldots,x_{r}^{\mathsf{T}}]^{\mathsf{T}}\in\mathbb{R}^{n}$,
$x_{i}\in\mathbb{R}^{n_{i}}$. We let the composite norm
$\|x\|_{\text{cmpst}}:=\|\left[\|x_{1}\|_{1},\ldots,\|x_{r}\|_{r}\right]\|_{\text{agg}}$,
with $\|\cdot\|_{i}$($\|\cdot\|_{\text{agg}}$) being local(aggregating) norms
on $\mathbb{R}^{n_{i}}$($\mathbb{R}^{n}$), and induced matrix norm
$\|\cdot\|_{i}$($\|\cdot\|_{\text{agg}}$) and matrix measure
$\mu_{i}(\cdot)$($\mu_{\text{agg}}(\cdot)$). Finally, given an $n\times n$
block matrix $A$ with $A_{ij}\in\mathbb{R}^{n_{i}\times n_{j}}$, define:

  1. (i)

the aggregate majorant $\lceil A\rceil$ with $(\lceil
A\rceil)_{ij}=\|A_{ij}\|_{ij}$;

  2. (ii)

the aggregate Metzler majorant $\lceil A\rceil_{\text{M}}$ with

$(\lceil A\rceil_{\text{M}})_{ij}=\begin{cases}\|A_{ij}\|_{ij},\text{if }j\neq
i,\\\ \mu_{i}(A_{ii}),\text{if }j=i,\end{cases}$

where $\|A_{ij}\|_{{ij}}:=\sup_{\|x_{j}\|_{j}=1}\|A_{ij}x_{j}\|_{i}$.
Then,

  * (i)

the composite norm $\|\cdot\|_{\text{cmpst}}$ is well-defined, i.e. satisfying
the norm properties;

  * (ii)

If the aggregating norm $\|\cdot\|_{\text{agg}}$ is monotonic, then:

    1. 1.

$\mu_{\text{cmpst}}(A)\leq\mu_{\text{agg}}(\lceil A\rceil_{\text{M}})$;

    2. 2.

$\|A\|_{\text{cmpst}}\leq\|\lceil A\rceil\|_{\text{agg}}$.

If the norms $\|\cdot\|_{i}$, $\|\cdot\|_{\text{agg}}$ in Lemma 1 are
$p$-norms (with the same $p$) then $\|\cdot\|_{\text{cmpst}}$ is again a
$p$-norm defined on a larger space. The next lemma follows from [17, Theorem
$2.4$].

######  Lemma 2.

Let $u:[-\tau_{0},+\infty)\rightarrow\mathbb{R}_{\geq 0}$ , $\tau_{0}<+\infty$
and assume that

$D^{+}u(t)\leq au(t)+b\sup_{t-\tau(t)\leq s\leq t}u(s)+c,\ \ t\geq t_{0}$

with: (i) $\tau(t)$ being bounded and non-negative, i.e.
$0{\leq}\tau(t)\leq\tau_{0}$, $\forall t$; (ii) $u(t)=|\varphi(t)|$, $\forall
t\in[t_{0}-\tau_{0},t_{0}]$ where $\varphi(t)$ is bounded and continuous in
$[t_{0}-\tau_{0},t_{0}]$; (iii) $a<0$, $b\geq 0$ and $c\geq 0$. Assume that
there exists some $\sigma>0$ such that $a+b\leq-\sigma<0,\forall t\geq t_{0}$.
Then:

$u(t)\leq\sup_{t_{0}-\tau_{0}\leq s\leq
t_{0}}u(s)e^{-\hat{\lambda}t}+\frac{c}{\sigma}$

where $\hat{\lambda}:=\inf_{t\geq
t_{0}}\\{\lambda|\lambda(t)+a+be^{\lambda(t)\tau(t)}=0\\}$ is positive.

##  III System description and problem formulation

Consider a network system comprised of $N$ agents with the dynamics of the
$i$-th agent described by

$\displaystyle\begin{split}&\dot{x}_{i}(t)=f_{i}(x_{i}(t),t)+u_{i}(t)+d_{i}(t)\\\
&{y_{i}(t)=g_{i}(x_{i}(t))}\end{split}$ (1)

$t\geq t_{0}$, $i\in\\{1,\ldots,N\\}$. In the above expression,
$x_{i}(t)\in\mathbb{R}^{n}$ denotes the agent state,
$u_{i}(t)\in\mathbb{R}^{n}$ denotes the control input,
$f_{i}:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{n}$ is a smooth
function and $g_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{r}$ is the output
function, which we assume to be Lipschitz. The term
$d_{i}(t)\in\mathbb{R}^{n}$ models the external disturbance of the form:

$\displaystyle
d_{i}(t)=w_{i}(t)+\bar{d}_{i}(t):=w_{i}(t)+\sum_{k=0}^{m-1}\bar{d}_{i,k}\cdot
t^{k}$ (2)

where $\bar{d}_{i}(t)$ represents the polynomial component of the disturbance
of the order $m-1$ ($m\in\mathbb{Z}_{>0}$) with $\bar{d}_{i,k}$ being constant
vectors and $w_{i}(t)$ is a piecewise continuous signal capturing _residual_
terms in the disturbance that are not modelled with the polynomial. In order
to reject polynomials up to order $m-1$, we design the control protocol
$u_{i}(t)$ following [4] with integral actions delivered by $m$ multiplex
layers as illustrated in Figure 1.

Figure 1: An example of the multiplex architecture (layer $0$ to layer $m$).
The agents are represented by orange nodes and the leader is represented by
the black node. Layers are not required to have the same topology. Each layer
outputs a signal to the layer immediately below (if any), i.e. layer $k$
outputs $[r_{1,k}^{\mathsf{T}},\ldots,r_{N,k}^{\mathsf{T}}]^{\mathsf{T}}$ to
layer $k-1$. Only one disturbance, e.g. $d_{1}$ on agent $1$, is shown.

$\displaystyle\begin{split}u_{i}(t)&=h_{i,0}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)+h_{i,0}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\\\
&+r_{i,1}(t)\\\
\dot{r}_{i,1}(t)&=h_{i,1}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)+h_{i,1}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\\\
&+r_{i,2}(t)\\\ \vdots\\\
\dot{r}_{i,m}(t)&=h_{i,m}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)+h_{i,m}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\end{split}$
(3)

where $\\{x_{j}\\}_{j\in\mathcal{N}_{i}}$ denotes the state of the neighbours
of agent $i$, $x_{l}(t)$ is the exhogenous reference signal from e.g., a
leader. The functions $h_{i,k},h_{i,k}^{(\tau)},\forall k$ are smooth coupling
functions on $k$-th multiplex layer modelling delay-free and delayed
communications from neighbours of agent $i$ and the leader. For delayed
communications, we consider time-varying delays
$\tau_{ij}(t),\tau_{il}(t)\in(0,\tau_{\max}]$ when information is transmitted
to agent $i$ from agent $j$ and from the leader, respectively. Note that in
general, $\tau_{ij}(t)\neq\tau_{ji}(t)$. In what follows, time dependence
inside these coupling functions are omitted for notational convenience.

######  Remark 1.

Protocol in (3) can be used for both leaderless and leader-follower networks.
It arises in a wide range of applications. For example, the classic diffusive-
type protocol can be written as in (3) with
$u_{i}(t)=h_{i,0}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)=\sum_{j\in\mathcal{N}_{i}}a_{ij}(x_{j}(t-\tau_{ij})-x_{i}(t-\tau_{ij}))$,
i.e. when all the communications are delayed and no integral actions are
applied in a leaderless network.

Let $x(t)=[x_{1}^{\mathsf{T}}(t),\ldots,x_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$,
$u(t)=[u_{1}^{\mathsf{T}}(t),\ldots,u_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$,
$y(t)=[y_{1}^{\mathsf{T}}(t),\ldots,y_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$ and
$d(t)=w(t)+\bar{d}(t)$ where
$w(t)=[w_{1}^{\mathsf{T}}(t),\ldots,w_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$,
$\bar{d}(t)=[\bar{d}_{1}^{\mathsf{T}}(t),\ldots,\bar{d}_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$,
the interconnected system (1) can then be written in a compact form as

$\displaystyle\begin{split}&\dot{x}(t)=f(x(t),t)+u(t)+d(t)\\\
&{y(t)=g(x(t))}\end{split}$ (4)

$t\geq t_{0}$, where
$f(x(t),t)=[f_{1}^{\mathsf{T}}(x_{1}(t),t),\ldots,f_{N}^{\mathsf{T}}(x_{N}(t),t)]^{\mathsf{T}}$
and
$g(x(t))=[g_{1}^{\mathsf{T}}(x_{1}(t)),\ldots,g_{N}^{\mathsf{T}}(x_{N}(t))]^{\mathsf{T}}$.
We also let
$r_{k}(t)=[r_{1,k}^{\mathsf{T}}(t),\ldots,r_{N,k}^{\mathsf{T}}(t)]^{\mathsf{T}}$,
$k\in\\{1,\ldots,m\\}$. Next, we state the control goal in terms of the
desired solution of the network when there are no disturbances and delays. The
desired solution is the solution of the unperturbed system satisfying: (i)
$\dot{x}_{i}^{\ast}=f_{i}(x_{i}^{\ast},t)$, $\forall i$, $\forall t\geq t_{0}$
and (ii) $r_{i,k}^{\ast}(t)=0$, $\forall i,\forall k$. The desired output is
then
$y^{\ast}(t):=[y_{1}^{*\mathsf{T}}(t),\dots,y_{N}^{\ast\mathsf{T}}(t)]^{\mathsf{T}}$,
with ${y}^{\ast}_{i}(t)=g_{i}(x_{i}^{\ast}(t)),\forall i$. In what follows, we
simply term
$x^{\ast}(t):=[x_{1}^{\ast\mathsf{T}}(t),\ldots,x_{N}^{\ast\mathsf{T}}(t)]^{\mathsf{T}}$
as desired solution. We are ready to give the following:

######  Definition 1.

The closed-loop system (4) affected by disturbance $d(t)=w(t)+\bar{d}(t)$ is
Input-Output Stable with respect to $w(t)$ if there exists class
$\mathcal{KL}$ functions $\alpha(\cdot,\cdot)$, $\beta(\cdot,\cdot)$, a class
$\mathcal{K}$ function $\gamma(\cdot)$, such that

$\displaystyle\begin{split}&{\|y(t)-y^{*}(t)\|}\leq\alpha\left(\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}\|x(s)-x^{*}(s)\|,t-t_{0}\right)\\\
&+\beta\left(\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}\sum_{k=1}^{m}\|r_{k}(s)+\bar{d}^{(k-1)}(s)\|,t-t_{0}\right)\\\
&+\gamma\left(\sup_{t}\|w(t)\|\right)\end{split}$

holds $\forall t\geq t_{0}$, $\forall x(s)$, $\forall
r_{k}(s),k\in\\{1,\ldots,m\\}$, where $x(s)=\varphi(s),r_{k}(s)=\phi_{k}(s)$
with $\varphi(s),\phi_{k}(s)$ being continuous and bounded functions in
$[t_{0}-\tau_{\max},t_{0}]$.

##  IV Technical results

In this section, we give a sufficient condition guaranteeing that the closed-
loop system (4) affected by disturbances of the form (2) is Input-Output
Stable. The results are stated in terms of:

  1. 1.

a composite norm ([16, Section 2.4.4])
$\|x\|_{\text{cmpst}}=\|\left[\|x_{1}\|_{p},\ldots,\|x_{N}\|_{p}\right]\|_{\infty,[\eta]^{-1}}$
where $[\eta]^{-1}\in\mathbb{R}^{N\times N}$ is a diagonal weighting matrix
with $\eta\in\mathbb{R}^{N}_{>0}$;

  2. 2.

a block diagonal coordinate transformation matrix
$T=\text{diag}\\{T_{1},\ldots,T_{N}\\}$.

For the statement of our result, it is also useful to relabel the delays
affecting the network. Specifically, we define $\tau_{k}(t),k=1,\ldots,q$ with
$q\leq N^{2}$, as an element of the set

$\\{\tau_{ij}(t):i,j=1,\ldots,N,i\neq
j\\}\bigcup\\{\tau_{il}(t):i=1,\ldots,N\\}$

######  Proposition 1.

Consider the closed-loop system (1) affected by external disturbance (2).
Assume that, $\forall t\geq t_{0}$ and for some matrices
$T_{i},T_{j},i,j=1,\ldots,N$, the following conditions are satisfied for some
$0<\underline{\sigma}<\bar{\sigma}<\infty$:

  * C1

$h_{i,k}(x_{i}^{*},\\{x_{j}^{\ast}\\}_{j\in\mathcal{N}_{i}},x_{l},t)=h_{i,k}^{(\tau)}(x_{i}^{*},\\{x_{j}^{\ast}\\}_{j\in\mathcal{N}_{i}},x_{l},t)=0$,
$\forall i,k$;

  * C2

$\mu_{p}(T_{i}\tilde{A}_{ii}(t)T_{i}^{-1})+\sum_{j}\frac{\eta_{j}}{\eta_{i}}\|T_{i}\tilde{A}_{ij}(t)T_{j}^{-1}\|_{p}\leq-\bar{\sigma}$,
$\forall i$, $\forall x$, $\forall x_{l}$;

  * C3

$\sum_{k=1}^{q}\sum_{j}\frac{\eta_{j}}{\eta_{i}}\|T_{i}(\tilde{B}_{k}(t))_{ij}T_{j}^{-1}\|_{p}\leq\underline{\sigma}$,
$\forall i$, $\forall x$, $\forall x_{l}$.

In the above expression,

$\displaystyle\begin{split}\tilde{A}_{ii}(t)&=\left[\begin{matrix}\frac{\partial
f_{i}}{\partial x_{i}}+\frac{\partial h_{i,0}}{\partial
x_{i}}&\mathbb{I}_{n}&0_{n}&\cdots&0_{n}\\\ \frac{\partial h_{i,1}}{\partial
x_{i}}&0_{n}&\mathbb{I}_{n}&\cdots&0_{n}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ \frac{\partial h_{i,m-1}}{\partial
x_{i}}&0_{n}&0_{n}&\cdots&\mathbb{I}_{n}\\\ \frac{\partial h_{i,m}}{\partial
x_{i}}&0_{n}&0_{n}&\cdots&0_{n}\end{matrix}\right]\\\
\tilde{A}_{ij}(t)&=\left[\begin{matrix}\frac{\partial h_{i,0}}{\partial
x_{j}}&0_{n}&\cdots&0_{n}\\\ \vdots&\vdots&\ddots&\vdots\\\ \frac{\partial
h_{i,m}}{\partial x_{j}}&0_{n}&\cdots&0_{n}\end{matrix}\right]\\\
(\tilde{B}_{k})_{ij}(t)&=\left[\begin{matrix}\frac{\partial
h_{i,0}^{(\tau)}}{\partial x_{j}}&0_{n}&\cdots&0_{n}\\\
\vdots&\vdots&\ddots&\vdots\\\ \frac{\partial h_{i,m}^{(\tau)}}{\partial
x_{j}}&0_{n}&\cdots&0_{n}\end{matrix}\right]\end{split}$ (5)

Then, the system is Input-Output Stable. In particular,

$\displaystyle\begin{split}&{\|y(t)-y^{*}(t)\|_{\text{cmpst}}}\leq\|T^{-1}\|_{\text{cmpst}}\|T\|_{\text{cmpst}}\cdot\\\
&\cdot e^{-\lambda(t-t_{0})}\biggl{(}\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}\|x(s)-x^{*}(s)\|_{\text{cmpst}}\\\ &+\sup_{t_{0}-\tau_{\max}\leq s\leq
t_{0}}\sum_{k=1}^{m}\|\sum_{b=0}^{m-k}\frac{(m-1-b)!}{(m-k-b)!}\cdot\bar{d}_{m-1-b}\cdot
s^{m-k-b}\\\
&+r_{k}(s)\|_{\text{cmpst}}\biggr{)}+\frac{\|T^{-1}\|_{\text{cmpst}}\|T\|_{\text{cmpst}}}{\bar{\sigma}-\underline{\sigma}}\sup_{t}\|w(t)\|_{\text{cmpst}}\end{split}$

where $\lambda>0$ is the solution of the following equation:

$\displaystyle\lambda-\bar{\sigma}+\underline{\sigma}e^{\lambda\tau_{\max}}=0$
(6)

######  Remark 2.

Condition $C1$ implies that $u_{i}(t)=0$ at the desired solution which
guarantees $x^{\ast}(t)$ is a solution of the unperturbed dynamics. $C2$ and
$C3$ give conditions on the Jacobian of delay-free and delayed part of agent
dynamics, respectively. These conditions are concerned with the couplings
between agents and depend on the number of neighbours of the agents and on the
coupling strength.

######  Remark 3.

If $C2$ and $C3$ are satisfied, then the network is connective stable [18,
Chapter 2.1]. In the special case where the agents are contracting and there
are no delays, the conditions yield a known property of contracting systems:
convergence of the network is achieved even without coupling (however, in this
case, there is no guarantee that the solutions towards which the agents
converge is the desired solution).

######  Remark 4.

The diagonal weighting matrix $[\eta]^{-1}$ in $\|x\|_{\text{cmpst}}$ can be
carefully designed to achieve sharper bounds for the induced matrix norm,
according to [16]. Hence, we could find such matrix $[\eta]^{-1}$ to minimise
the upper bound in $C2$ or $C3$. On the other hand, the selection of the
coordinate transformation matrix $T$ is also of great significance as in many
cases it turns out to be difficult to design controllers fulfilling conditions
proposed for the original system. The conditions in Proposition 1 not only
provides guideline for the design of the controller, but also for the
computation of matrix $T$, see also [4].

######  Remark 5.

Following Proposition 1, the closed-loop network has a convergence rate
$\lambda$ depending on $\tau_{\max}$. As is highlighted in (6), larger
$\tau_{\max}$ yields lower $\lambda$ which means slower convergence towards
the system desired solution. Hence (6) gives an implicit condition on
communication delays for networks with desired convergence rates.

The proof is inspired by the derivations of [4, Proposition $1$]. Hence, some
technical steps are omitted here for brevity.

###### Proof.

Consider the augmented state
$z_{i}(t)=[x_{i}^{\mathsf{T}}(t),\zeta_{i,1}^{\mathsf{T}}(t),\ldots,\zeta_{i,m}^{\mathsf{T}}(t)]^{\mathsf{T}}$
where

$\displaystyle\zeta_{i,k}(t)$
$\displaystyle=r_{i,k}(t)+\sum_{b=0}^{m-k}\frac{(m-1-b)!}{(m-k-b)!}\cdot\bar{d}_{i,m-1-b}\cdot
t^{m-k-b}$

$k\in\\{1,\cdots,m\\}$. The desired solution is
$z_{i}^{\ast}(t)=[x_{i}^{\ast\mathsf{T}}(t),0_{1\times n},\ldots,0_{1\times
n}]^{\mathsf{T}}$. Hence the error dynamics
$e_{i}(t):=z_{i}(t)-z^{\ast}_{i}(t)$ follows

$\displaystyle\begin{split}\dot{e}_{i}(t)&=\tilde{f}_{i}(z_{i},t)-\tilde{f}_{i}(z_{i}^{\ast},t)+\tilde{h}_{i}(z,t)+\tilde{h}_{i}^{(\tau)}(z,t)+\tilde{w}_{i}(t)\end{split}$

where $\tilde{f}_{i}(z_{i},t)=[f_{i}^{\mathsf{T}}(x_{i},t),0_{1\times
n},\ldots,0_{1\times n}]^{\mathsf{T}}$ and

$\displaystyle\begin{split}\tilde{h}_{i}(z,t)=\left[\begin{array}[]{*{20}c}h_{i,0}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)+\zeta_{i,1}(t)\\\
h_{i,1}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)+\zeta_{i,2}(t)\\\
\vdots\\\
h_{i,m-1}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)+\zeta_{i,m}(t)\\\
h_{i,m}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\end{array}\right]\end{split}$

and

$\displaystyle\begin{split}\tilde{h}_{i}^{(\tau)}(z,t)=\left[\begin{array}[]{*{20}c}h_{i,0}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\\\
h_{i,1}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\\\ \vdots\\\
h_{i,m}^{(\tau)}(x_{i},\\{x_{j}\\}_{j\in\mathcal{N}_{i}},x_{l},t)\end{array}\right]\end{split}$

and $\tilde{w}_{i}(t)=[w_{i}^{\mathsf{T}}(t),0_{1\times n},\ldots,0_{1\times
n}]^{\mathsf{T}}$. Then following [19], let $\eta_{i}(\rho)=\rho
z_{i}+(1-\rho)z_{i}^{\ast}$,
$\eta(\rho)=[\eta_{1}^{\mathsf{T}}(\rho),\ldots,\eta_{N}^{\mathsf{T}}(\rho)]^{\mathsf{T}}$
and we can rewrite the error dynamics in a compact form as

$\displaystyle\dot{e}(t)=A(t)e(t)+\sum_{k=1}^{q}B_{k}(t)e(t-\tau_{k}(t))+\tilde{w}(t)$
(7)

where
$\tilde{w}=[\tilde{w}_{1}^{\mathsf{T}}(t),\ldots,\tilde{w}_{N}^{\mathsf{T}}(t)]^{\mathsf{T}}$.
The Jacobian matrix $A(t)$ has entries
$A_{ij}(t)=\int_{0}^{1}\tilde{A}_{ij}(t)d\rho$ where
$\tilde{A}_{ij}(t)=J_{\tilde{f}_{i}}(\eta_{j}(\rho),t)+J_{\tilde{h}_{i}}(\eta_{j}(\rho),t)$
and $B_{k}(t)$ with entries
$(B_{k})_{ij}(t)=\int_{0}^{1}(\tilde{B}_{k})_{ij}(t)d\rho$ where
$(\tilde{B}_{k})_{ij}(t)=J_{\tilde{h}_{i}^{(\tau)}}(\eta_{j}(\rho),t)$. Then
let $\tilde{e}(t)=Te(t)$ where $T=\text{diag}\\{T_{1},\ldots,T_{N}\\}$, we
have

$\displaystyle\dot{\tilde{e}}(t)=TA(t)T^{-1}\tilde{e}(t)+\sum_{k=1}^{q}TB_{k}(t)T^{-1}\tilde{e}(t-\tau_{k}(t))+T\tilde{w}(t)$

Taking the Dini derivative of $\|\tilde{e}(t)\|_{\text{cmpst}}$, we obtain

$\displaystyle\begin{split}&D^{+}\|\tilde{e}(t)\|_{\text{cmpst}}=\limsup_{h\rightarrow
0^{+}}\frac{1}{h}\left(\|\tilde{e}(t+h)\|_{\text{cmpst}}-\|\tilde{e}(t)\|_{\text{cmpst}}\right)\\\
&\leq\mu_{\text{cmpst}}(TA(t)T^{-1})\|\tilde{e}(t)\|_{\text{cmpst}}+\sup_{t}\|T\tilde{w}(t)\|_{\text{cmpst}}\\\
&\phantom{=}+\sum_{k=1}^{q}\|TB_{k}(t)T^{-1}\|_{\text{cmpst}}\sup_{t-\tau_{\max}\leq
s\leq t}\|\tilde{e}(s)\|_{\text{cmpst}}\end{split}$

Condition $C2$ and $C3$ implies that
$\mu_{\text{cmpst}}(TA(t)T^{-1})\leq-\bar{\sigma}$ and
$\sum_{k=1}^{q}\|TB_{k}(t)T^{-1}\|_{\text{cmpst}}\leq\underline{\sigma}$.
Then, following similar steps to the ones of the proof of Proposition $1$ in
[4] and leveraging the Lipschitz assumption of the output function, we get the
upper bound of the output deviation of (1) from its desired output:

$\displaystyle\begin{split}&{\|y(t)-y^{*}(t)\|_{\text{cmpst}}}\leq\\\
&\|T^{-1}\|_{\text{cmpst}}\|T\|_{\text{cmpst}}e^{-\lambda(t-t_{0})}\biggl{(}\sup_{t_{0}-\tau_{\max}\leq
s\leq t_{0}}\|x(s)-x^{*}(s)\|_{\text{cmpst}}\\\ &+\sup_{t_{0}-\tau_{\max}\leq
s\leq
t_{0}}\sum_{k=1}^{m}\|\sum_{b=0}^{m-k}\frac{(m-1-b)!}{(m-k-b)!}\cdot\bar{d}_{m-1-b}\cdot
s^{m-k-b}\\\
&+r_{k}(s)\|_{\text{cmpst}}\biggr{)}+\frac{\|T^{-1}\|_{\text{cmpst}}\|T\|_{\text{cmpst}}}{\bar{\sigma}-\underline{\sigma}}\sup_{t}\|w(t)\|_{\text{cmpst}}\end{split}$

∎

##  V Application Example

We consider an MTDC grid model from [20]. In our example, the grid has $30$
terminals arranged in a ring (an example illustration of a ring consisting of
$5$ terminals is given in Figure 2). The dynamics of the terminals are
described by

$\displaystyle\begin{split}c_{i}\dot{v}_{i}(t)&=-\sum_{j\in\mathcal{N}_{i}}I_{ij}(t)+u_{i}(t)+d_{i}(t),\
i\in\\{1,\ldots,5\\}\\\ y_{i}(t)&=v_{i}(t)\end{split}$ (8)

where the current $I_{ij}(t)=\frac{1}{R_{ij}}(v_{i}(t)-v_{j}(t))$ according to
Ohm’s law and note that $v_{j}(t)$ is not delayed in this equation as
communication is not required. In the above expression, $v_{i}$ denotes the
voltage deviation of terminal $i$ from the nominal voltage $v^{\text{nom}}$
which is assumed to be identical for all terminals, $c_{i}$ is the capacitance
of terminal $i$, $R_{ij}$ is the line resistance between terminal $i$ and
terminal $j$ satisfying $R_{ij}=R_{ji}$, $u_{i}$ is the injected control
current, $d_{i}$ is the disturbance due to e.g. load changes at the terminal
and $y_{i}$ is the output. The design of control protocols to reject constant
disturbances for (8) was considered in [21] but without delays and in [22]
with constant delays. In this example, the capacitance $c_{i}=1$mF and the
resistance $R_{ij}=20\Omega$. We now leverage Proposition 1 to consider the
case where delays are heterogeneous and the disturbances are polynomial.

Figure 2: Illustration of the MTDC composed of $5$ terminals. Terminal $1$ is
affected by a disturbance $d_{1}$.

###  V-A Controller design

In this example, we consider a first order disturbance, which can model e.g.
the rapid increase of the current in the terminal caused by fault [23]. To
reject such disturbance before being diagnosed, we design the controller with
$2$ multiplex layers following:

$\displaystyle\begin{split}u_{i}(t)&=-k_{0}v_{i}(t)-\sum_{j\in\mathcal{N}_{i}}k_{0}^{(\tau)}(v_{i}(t-\tau_{ij}(t))-v_{j}(t-\tau_{ij}(t)))\\\
&+r_{i,1}(t)\\\
\dot{r}_{i,1}(t)&=-k_{1}v_{i}(t)-\sum_{j\in\mathcal{N}_{i}}k_{1}^{(\tau)}(v_{i}(t-\tau_{ij}(t))-v_{j}(t-\tau_{ij}(t)))\\\
&+r_{i,2}(t)\\\
\dot{r}_{i,2}(t)&=-k_{2}v_{i}(t)-\sum_{j\in\mathcal{N}_{i}}k_{2}^{(\tau)}(v_{i}(t-\tau_{ij}(t))-v_{j}(t-\tau_{ij}(t)))\end{split}$
(9)

In the above expression, the delay occurs because the voltage information
needs to be transmitted from terminal $j$ to terminal $i$ via communication
(see also [22]). For our design, we consider the composite norm
$\|x\|_{\text{cmpst}}=\|\left[\|x_{1}\|_{2},\ldots,\|x_{N}\|_{2}\right]\|_{\infty,{[\eta]^{-1}}},{\eta\in\mathbb{R}_{>0}^{N}}$
and the coordinate transformation matrix $T$ with identical diagonal blocks

$\displaystyle T_{i}={\tilde{T}=}\left[\begin{matrix}1&\alpha&0\\\
0&1&\beta\\\ 0&0&1\end{matrix}\right],\ \alpha,\beta\in\mathbb{R},\ \forall i$

The desired solution of (8) is $v_{i}^{\ast}(t)=0$, $\forall i$. Hence, $C1$
is guaranteed by design of control (9). We recast the fulfilment of $C2$ and
$C3$ as an optimisation problem following

$\displaystyle\min_{\xi}\ \mathcal{J}$ (10a) $\displaystyle s.t.\quad
k_{0}\geq 0,k_{1}\geq 0,k_{2}\geq 0,k_{0}^{(\tau)}\geq 0,k_{1}^{(\tau)}\geq
0,k_{2}^{(\tau)}\geq 0$ (10b)
$\displaystyle\phantom{s.t.\quad}k_{0}+k_{0}^{(\tau)}>0,k_{1}+k_{1}^{(\tau)}>0,k_{2}+k_{2}^{(\tau)}>0$
(10c) $\displaystyle\phantom{s.t.\quad}\bar{\sigma}>0,\underline{\sigma}\geq
0,\bar{\sigma}-\underline{\sigma}>0$ (10d)
$\displaystyle\phantom{s.t.\quad}\mu_{2}(\tilde{T}\tilde{A}_{ii}\tilde{T}^{-1})+\sum_{j}\frac{\eta_{j}}{\eta_{i}}\|\tilde{T}\tilde{A}_{ij}\tilde{T}^{-1}\|_{2}\leq-\bar{\sigma},\
\ \forall i,j$ (10e)
$\displaystyle\phantom{s.t.\quad}\sum_{k=1}^{q}\sum_{j}\frac{\eta_{j}}{\eta_{i}}\|\tilde{T}(\tilde{B}_{k})_{ij}\tilde{T}^{-1}\|_{2}\leq\underline{\sigma},\
\ \forall i,j$ (10f)

where
$\xi:=[k_{0},k_{1},k_{2},k_{0}^{(\tau)},k_{1}^{(\tau)},k_{2}^{(\tau)},\bar{\sigma},\underline{\sigma}]$
are the decision variables and
$\mathcal{J}=-k_{0}^{(\tau)}-k_{1}^{(\tau)}-k_{2}^{(\tau)}$ is the cost
function chosen, in this case, to maximise the coupling between terminals. In
the above expression, the constraints on matrix measure and matrix norm
($\tilde{A}_{ii},\tilde{A}_{ij},(\tilde{B}_{k})_{ij}$ are computed from (5)
and are constant matrices) can be recast as LMIs, for fixed $\alpha$, $\beta$
and $\eta$. Indeed, (10e) is satisfied if

$\displaystyle\begin{split}&\mu_{2}(\tilde{T}\tilde{A}_{ii}\tilde{T}^{-1})\leq-
b_{1},\ \ \forall i\\\
&\sum_{j}\frac{\eta_{j}}{\eta_{i}}\|\tilde{T}\tilde{A}_{ij}\tilde{T}^{-1}\|_{2}\leq
b_{2},\ \ \forall i,j\\\ &b_{1}>0,b_{2}\geq
0,-b_{1}+b_{2}\leq-\bar{\sigma}\end{split}$ (11)

where $b_{1},b_{2}$ are auxiliary variables. Also, note that due to the ring
topology, for terminal $i$ there are two neighbours. It is then useful to
recast
$\sum_{j}\frac{\eta_{j}}{\eta_{i}}\|\tilde{T}\tilde{A}_{ij}\tilde{T}^{-1}\|_{2}\leq
b_{2}$ as

$\displaystyle\begin{split}&\frac{\eta_{j}}{\eta_{i}}\|\tilde{T}\tilde{A}_{i{j}}\tilde{T}^{-1}\|_{2}\leq
b_{3},\ \ \forall i,j\\\ &b_{3}\geq 0,2b_{3}\leq b_{2}\end{split}$ (12)

with auxiliary variable $b_{3}$. Analogously, we recast (10f) as

$\displaystyle\begin{split}&\|\tilde{T}(\tilde{B}_{k})_{ii}\tilde{T}^{-1}\|_{2}\leq
b_{4},\ \ \forall i\\\
&\frac{\eta_{j}}{\eta_{i}}\|\tilde{T}(\tilde{B}_{k})_{i{j}}\tilde{T}^{-1}\|_{2}\leq
b_{5},\ \ \forall i,j\\\ &b_{4}\geq 0,b_{5}\geq
0,q(b_{4}+b_{5})\leq\underline{\sigma}\end{split}$ (13)

where $b_{4},b_{5}$ are auxiliary variables and where $q=60$ is the
cardinality of the set of delays. Now following the steps in [4, Appendix B],
(11)-(13) can be recast as the following

$\displaystyle\begin{split}&[\tilde{T}\tilde{A}_{ii}\tilde{T}^{-1}]_{s}\leq-
b_{1}\mathbb{I}_{3},\ \ \forall i\\\
&\left[\begin{matrix}b_{3}\mathbb{I}_{3}&\frac{\eta_{j}}{\eta_{i}}(\tilde{T}\tilde{A}_{i{j}}\tilde{T}^{-1})^{\mathsf{T}}\\\
\frac{\eta_{j}}{\eta_{i}}\tilde{T}\tilde{A}_{i{j}}\tilde{T}^{-1}&b_{3}\mathbb{I}_{3}\end{matrix}\right]\succeq
0,\ \ \forall i,j\\\
&\left[\begin{matrix}b_{4}\mathbb{I}_{3}&(\tilde{T}(\tilde{B}_{k})_{ii}\tilde{T}^{-1})^{\mathsf{T}}\\\
\tilde{T}(\tilde{B}_{k})_{ii}\tilde{T}^{-1}&b_{4}\mathbb{I}_{3}\end{matrix}\right]\succeq
0,\ \ \forall i\\\
&\left[\begin{matrix}b_{5}\mathbb{I}_{3}&\frac{\eta_{j}}{\eta_{i}}(\tilde{T}(\tilde{B}_{k})_{i{j}}\tilde{T}^{-1})^{\mathsf{T}}\\\
\frac{\eta_{j}}{\eta_{i}}\tilde{T}(\tilde{B}_{k})_{i{j}}\tilde{T}^{-1}&b_{5}\mathbb{I}_{3}\end{matrix}\right]\succeq
0,\ \ \forall i,j\\\ &b_{1}>0,b_{2}\geq 0,b_{3}\geq 0,b_{4}\geq 0,b_{5}\geq
0,2b_{3}\leq b_{2}\\\
&-b_{1}+b_{2}\leq-\bar{\sigma},60(b_{4}+b_{5})\leq\underline{\sigma}\end{split}$
(14)

We then solve the optimisation problem (10) replacing (10e)-(10f) with (14)
for a grid of $\alpha$ and $\beta$ and for the fixed $\eta$111The code solving
the optimisation problem and the data, i.e. $\eta$, can be found in
https://tinyurl.com/yc3frafb.. The optimal solution for this problem is:
$k_{0}=0.7445,k_{1}=1.3399,k_{2}=0.5052,k_{1}^{(\tau)}=0.00057,k_{1}^{(\tau)}=0.00076,k_{2}^{(\tau)}=0.00048$,
with corresponding $\alpha=-0.5,\beta=-1$.

###  V-B Simulation

We consider a disturbance acting on a random terminal, say terminal $1$, which
is $d_{1}(t)=3+t+e^{-0.2t}\sin t$. The heterogeneous delays are selected as
$\tau_{ij}(t)=\tau_{i}(t)=0.1+0.1\sin(t+i)$ seconds. The grid is assumed to be
initiated at $v_{i}(0)\sim\mathcal{N}(0,1),\forall i$. Figure 3 (top panel)
illustrates the voltage deviation of all the terminals from the nominal value.
It shows that all the deviations finally reduce to $0$ including the perturbed
terminal $1$. This is in accordance with the theoretical prediction as the
designed controller injected a ramp current to compensate for the ramp
component in the disturbance, as illustrated in the bottom panel.

Figure 3: Top panel: Voltage deviation from nominal voltage; Bottom panel:
Injected control current.

##  VI Conclusions and future work

We considered the problem of designing distributed multiplex integral control
protocols for nonlinear networks affected by delays and disturbances. The
designed control protocol, delivered via multiplex architecture and fulfilling
certain conditions leveraging non-Euclidean contraction theory, is able to:
(i) reject polynomial disturbances; (ii) achieve Input-to-State Stability for
nonlinear networks affected by heterogeneous time-varying delays. We validated
the results by considering the problem of controlling an MTDC grid and
simulations confirmed the effectiveness of the results. Apart from devising
approaches to systematically determine the right coordinate transformation
matrices $T$ and the weighting matrix $[\eta]^{-1}$, see e.g. [24], we would
like to tackle the challenge of designing protocols that guarantee the
stronger scalability property studied in [25, 4] for the network systems with
heterogeneous delays considered in this paper.

## References

  * [1] G. F. Silva, A. Donaire, A. McFadyen, and J. J. Ford, “String stable integral control design for vehicle platoons with disturbances,” Automatica, vol. 127, p. 109542, 2021.
  * [2] S. Knorn, A. Donaire, J. C. Agüero, and R. H. Middleton, “Passivity-based control for multi-vehicle systems subject to string constraints,” Automatica, vol. 50, no. 12, pp. 3224–3230, 2014.
  * [3] D. A. B. Lombana and M. Di Bernardo, “Multiplex PI control for consensus in networks of heterogeneous linear agents,” Automatica, vol. 67, pp. 310–320, 2016.
  * [4] S. Xie and G. Russo, “On the design of integral multiplex control protocols for nonlinear network systems with delays,” arXiv preprint arXiv:2206.03535 (submitted to Automatica), 2022.
  * [5] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,” Automatica, vol. 34, no. 6, pp. 683–696, 1998.
  * [6] S. Jafarpour, A. Davydov, M. Abate, F. Bullo, and S. Coogan, “Robust training and verification of implicit neural networks: A non-Euclidean contractive approach,” arXiv preprint arXiv:2208.03889, 2022.
  * [7] A. Davydov, S. Jafarpour, and F. Bullo, “Non-Euclidean contraction theory for robust nonlinear stability,” IEEE Transactions on Automatic Control, 2022.
  * [8] G. Russo, M. Di Bernardo, and E. D. Sontag, “Global entrainment of transcriptional systems to periodic inputs,” PLoS computational biology, vol. 6, no. 4, p. e1000739, 2010.
  * [9] Z. Aminzare and E. D. Sontag, “Contraction methods for nonlinear systems: A brief introduction and some open problems,” in 53rd IEEE Conference on Decision and Control, pp. 3835–3847, IEEE, 2014.
  * [10] M. di Bernardo, D. Fiore, G. Russo, and F. Scafuti, “Convergence, consensus and synchronization of complex networks via contraction theory,” Complex Systems and Networks, pp. 313–339, 2016.
  * [11] H. Tsukamoto, S.-J. Chung, and J.-J. E. Slotine, “Contraction theory for nonlinear stability analysis and learning-based control: A tutorial overview,” Annual Reviews in Control, vol. 52, pp. 135–169, 2021.
  * [12] H. S. Shiromoto, M. Revay, and I. R. Manchester, “Distributed nonlinear control design using separable control contraction metrics,” IEEE Transactions on Control of Network Systems, vol. 6, no. 4, pp. 1281–1290, 2018\.
  * [13] M. Revay, R. Wang, and I. R. Manchester, “Recurrent equilibrium networks: Flexible dynamic models with guaranteed stability and robustness,” arXiv preprint arXiv:2104.05942, 2021.
  * [14] W. Wang and J.-J. Slotine, “Contraction analysis of time-delayed communications and group cooperation,” IEEE Transactions on Automatic Control, vol. 51, no. 4, pp. 712–717, 2006.
  * [15] G. Russo, M. di Bernardo, and E. D. Sontag, “Stability of networked systems: A multi-scale approach using contraction,” in 49th IEEE Conference on Decision and Control (CDC), pp. 6559–6564, 2010.
  * [16] F. Bullo, Contraction Theory for Dynamical Systems.  Kindle Direct Publishing, 1.0 ed., 2022.
  * [17] L. Wen, Y. Yu, and W. Wang, “Generalized Halanay inequalities for dissipativity of Volterra functional differential equations,” Journal of Mathematical Analysis and Applications, vol. 347, no. 1, pp. 169–178, 2008\.
  * [18] D. D. Siljak, Decentralized control of complex systems.  Courier Corporation, 2011.
  * [19] C. Desoer and H. Haneda, “The measure of a matrix as a tool to analyze computer algorithms for circuit analysis,” IEEE Transactions on Circuit Theory, vol. 19, no. 5, pp. 480–486, 1972.
  * [20] M. Axelson-Fisk and S. Knorn, “Online distributed design for control cost reduction,” Automatica, vol. 141, p. 110312, 2022.
  * [21] M. Axelson-Fisk and S. Knorn, “A graph theoretic approach to ensure a scalable performance measure in multi-agent systems with MIMO agents,” IFAC-PapersOnLine, vol. 55, no. 13, pp. 103–108, 2022.
  * [22] M. Andreasson, M. Nazari, D. V. Dimarogonas, H. Sandberg, K. H. Johansson, and M. Ghandhari, “Distributed voltage and current control of multi-terminal high-voltage direct current transmission systems,” IFAC Proceedings Volumes, vol. 47, no. 3, pp. 11910–11916, 2014.
  * [23] J. Sneath and A. D. Rajapakse, “Fault detection and interruption in an earthed HVDC grid using ROCOV and hybrid DC breakers,” IEEE Transactions on Power Delivery, vol. 31, no. 3, pp. 973–981, 2014.
  * [24] V. Centorrino, A. Gokhale, A. Davydov, G. Russo, and F. Bullo, “Euclidean contractivity of neural networks with symmetric weights,” arXiv preprint arXiv:2302.13452, 2023.
  * [25] S. Xie, G. Russo, and R. H. Middleton, “Scalability in nonlinear network systems affected by delays and disturbances,” IEEE Transactions on Control of Network Systems, vol. 8, no. 3, pp. 1128–1138, 2021.


# Evolution of Non-parametric Morphology of Galaxies in the JWST CEERS Field
at $z\simeq$0.8-3.0

Yao Yao Deep Space Exploration Laboratory / Department of Astronomy,
University of Science and Technology of China, Hefei 230026, China  Jie Song
Deep Space Exploration Laboratory / Department of Astronomy, University of
Science and Technology of China, Hefei 230026, China  Xu Kong  Deep Space
Exploration Laboratory / Department of Astronomy, University of Science and
Technology of China, Hefei 230026, China  School of Astronomy and Space
Sciences, University of Science and Technology of China, Hefei, 230026, China
Guanwen Fang  Institute of Astronomy and Astrophysics, Anqing Normal
University, Anqing, 246133, China  Hong-Xin Zhang  Deep Space Exploration
Laboratory / Department of Astronomy, University of Science and Technology of
China, Hefei 230026, China  Xinkai Chen  Deep Space Exploration Laboratory /
Department of Astronomy, University of Science and Technology of China, Hefei
230026, China  Kavli Institute for Astronomy and Astrophysics, Peking
University, Beijing, 100871, China

###### Abstract

Galaxy morphology is one of the most fundamental ways to describe galaxy
properties, but the morphology we observe may be affected by wavelength and
spatial resolution, which may introduce systematic bias when comparing
galaxies at different redshift. Taking advantage of the broad wavelength
coverage from optical to near-IR and high resolution NIRCam instrument of
JWST, we measure the non-parametric morphological parameters of a total of
1376 galaxies at $z\simeq$0.8-3.0 in the CEERS field through an optimized code
called statmorph_csst. We divide our sample into three redshift intervals and
investigate the wavelength- and redshift-dependence of the morphological
parameters. We also explore how the widely-used galaxy type classification
methods based on the morphological parameters depend on wavelength and spatial
resolution. We find that there are variations in all morphological parameters
with rest-frame wavelength ($\lambda_{\rm rf}$), especially at the short
wavelength end, and the $\lambda_{\rm rf}$ mainly affects the classification
between late-type and early-type galaxy. As the $\lambda_{\rm rf}$ increases,
the galaxies on the $G-M_{20}$ diagram move to the upper left with a slope of
-0.23$\pm$0.03 on average. We find that spatial resolution mainly affects the
merger identification. The merger fraction in F200W resolution can be
$\gtrsim$2 times larger than that in F444W resolution. Furthermore, We compare
the morphological parameter evolution of galaxies with different stellar
masses. We find that there are differences in the morphological evolution of
high- and low-mass (log$M_{*}\geqslant$10 and 9$<$log$M_{*}<$10) galaxies in
the studied redshift range, which may be caused by their different evolution
paths.

Galaxy structure — Galaxy morphology — High redshift — JWST

††journal: ApJ††facilities: JWST††software: astropy (Collaboration et al.,
2022), Numpy (Harris et al., 2020), Matplotlib (Hunter, 2007), statmorph
(Rodriguez-Gomez et al., 2019), reproject (Robitaille et al., 2020), photutils
(Bradley et al., 2022), Cython (Behnel et al., 2011)

##  1 Introduction

Morphology of galaxies is a fundamental aspect that provides crucial insights
into the structure and formation of galaxies. It allows a direct
representation of the spatial distribution of stars and dust in galaxies,
which holds the key to understanding their evolution and assembly history. The
galaxy morphology of the local universe has rich diversity. Hubble (1926) has
classified galaxies into the now well-known tuning fork through their visual
structures. His classification includes three distinct groups: ellipticals,
spirals, and irregulars, where ellipticals are also usually called Early-Type
Galaxies (ETGs), and spirals and irregulars are also usually called Late-Type
Galaxies (LTGs). There are not only significant differences in morphology
between ETGs and LTGs but also significant differences in star formation rate
and other physical properties.

In the last decades, thanks to the sky surveys of Hubble Space Telescope (HST)
and ground-based telescopes, the relation between galaxy morphology and galaxy
physical evolution have been researched deeper (e.g., Kauffmann et al., 2003;
Baldry et al., 2004; Conselice, 2006; Mortlock et al., 2013) and more methods
to classify galaxy morphology have been developed, including Sérsic model
(Sérsic, 1963), bulge/disk decomposition (Freeman, 1970; Kent, 1985), and non-
parametric measurements (e.g., Conselice, 2003; Lotz et al., 2004; Freeman et
al., 2013; Ferrari et al., 2015). These methods are more objective than manual
visual inspection and enable comparison between results from different
researchers. In recent years, with the rapid evolution of computer technology,
machine-learning algorithms for galaxy morphological classification and merger
identification have been widely developed and utilized. (e.g., Domínguez
Sánchez et al., 2018; Tang et al., 2020; Yao et al., 2022; Fang et al., 2023).
At the same time, the powerful computing power has also facilitated a lot of
galaxy morphological studies based on cosmological simulation (e.g., Snyder et
al., 2015; Bignone et al., 2017; Rodriguez-Gomez et al., 2019; Bignone et al.,
2020), which provide a direct way to link observed morphology with physical
processes.

However, limited by the spatial resolution, sensitivity, and red wavelength
range of HST (reddest filter F160W), the morphological study of high-redshift
galaxies in rest-frame optical band stops at $z\simeq$3\. The new James Webb
Space Telescope (JWST) solves these problems well. Its superior spatial
resolution and longer wavelength filter set (reddest filter F444W of NIRCam)
not only allow us to explore galaxy structure at $z>$6 (e.g., Ferreira et al.,
2022a; Kartaltepe et al., 2023; Treu et al., 2023) but also make the images of
galaxies at lower redshifts have higher quality and wider wavelength range,
which can reduce the measurement uncertainties caused by the resolution,
signal-to-noise ratio (S/N) and morphological K-correction.

In this work, we make use of the JWST’s observation data and study the non-
parametric morphology of galaxies at $z\simeq$0.8-3.0. This paper is organized
as follows. In Section 2, we describe our selection criteria of the sample and
the non-parametric measurements from the data. In Section 3, we present our
results of non-parametric morphology measurement and discuss the effects of
the rest-frame wavelength ($\lambda_{\rm rf}$) and the point spread function
(PSF) on morphology and the differences between galaxies with different
stellar masses ($M_{*}$) at different redshifts. In Section 4, we summarize
the main results of this work. Throughout this paper, we adopt a $\Lambda$-CDM
cosmology with $\Omega_{\Lambda}$= 0.7, $\Omega_{m}$=0.3, and H0=70 km s-1
Mpc-1.

##  2 Data and Measurements

###  2.1 CEERS Images and Catalogs

CEERS111https://ceers.github.io/index.html (Finkelstein et al., in prep.) is
one of 13 early release science (ERS) surveys of JWST. It is designed to
obtain data covering the EGS (Extended Groth Strip, Davis et al., 2007)
extragalactic deep field early in Cycle 1. The observations, tests, and data
products will pave the way for Cycle 2 observations. CEERS is based around a
mosaic of 10 NIRCam (Rieke et al., 2005) pointings, with six obtaining NIRSpec
(Jakobsen et al., 2022) in parallel and four with MIRI (Rieke et al., 2015) in
parallel. The NIRCam images of CEERS includes 7 bands, from the short
wavelength (SW) detector (F115W, F150W, and F200W) and the long wavelength
(LW) detector (F277W, F356W, F410M, and F444W).

Here we make use of the first four CEERS NIRCam pointings, obtained on 2022
June 21, known as CEERS1, CEERS2, CEERS3, and CEERS6. We directly use the
final mosaic images and catalogs provided by Valentino et al.
(2023)222https://s3.amazonaws.com/grizli-v2/JwstMosaics/v4/index.html#ceers,
which is processed with the grizli (Brammer et al., 2022) software pipeline.
Besides zeropoint correction333https://github.com/gbrammer/grizli/pull/107
relative to jwst_0942.pmap, the 1/f, snowballs, and wisp correction have also
been applied by grizli. The pipeline also uses SEP (Barbary, 2016), a Python
library implementing Source-Extractor (Bertin & Arnouts, 1996) to extract
sources from the detection image combining of all the “wide” NIRCam LW filters
available (typically F277W+F356W+F444W) optimally weighted by their noise
maps. It also outputs JWST+HST photometric catalogs and EAZY-
py444https://github.com/gbrammer/eazy-py (Brammer et al., 2008) photometric
redshift catalogs for all pointings. Spectroscopic redshifts have also been
added to these catalogs, if available. The sampling rate of the final SW
mosaic images is 0″.02 pixel-1, and that of LW images is 0″.04 pixel-1.

###  2.2 Sample Selection

We select a sample to analyze morphology from the photometric and redshift
catalogs of Valentino et al. (2023), which contain 26438 sources. We select
galaxies that have mag_auto$<$27.0 in detection image, and cover at least one
pixel on the segmentation maps. A total of 2969 sources meet the above
conditions, of which 75 sources do not have effective photometric or
spectroscopic redshift and $M_{*}$ (in unit of $M_{\odot}$) in the above
catalogs. We divided the remaining 2894 sources into 3 redshift bins: 0.8-1.3,
1.3-2.0, and 2.0-3.0. The reason why we choose these values is that we can
match the NIRCam filters of different bins at a similar $\lambda_{\rm rf}$. In
addition, the reason why we do not choose a higher redshift is that the sample
size at the higher redshift is limited, the strong cosmic dimming effect, and
few filters can be matched with the low-redshift galaxies. Thus, in these
three redshift bins, we have four bands that can be matched at a similar
$\lambda_{\rm rf}$, at $\sim$570 nm (close to $V$), $\sim$760 nm (close to
$i$), $\sim$1000 nm, and $\sim$1340 nm respectively.

We use the catalog containing 2894 sources as the input of our process
pipeline for morphological parameter measurement (see Section 2.3). According
to the result of our measurement, we exclude sources that have an invalid
morphological parameter in any filters, mainly contributed by the diffraction
spikes of starlight. Lotz et al. (2004, 2006) and Treu et al. (2023) have
pointed out that the morphological parameter measurement accuracy is affected
when the average S/N per pixel ($\langle$S/N$\rangle$) goes below 2. So we
only select sources whose $\langle$S/N$\rangle$ are higher than 2 in all bands
we analyzed. As a result, the number of sources decreased from 2894 to 2068.
Section 4.3.2 of Rodriguez-Gomez et al. (2019) describes the method of
obtaining $\langle$S/N$\rangle$.

Figure 1 shows the distribution of the absolute magnitude of F444W(MF444W) and
the $M_{*}$ of our galaxies with redshift, and Table 1 summaries the details
of three redshift bins which we select for our main research. We only analyze
galaxies whose log$M_{*}>$9 to balance the completeness and range of $M_{*}$,
ensuring that the sample size is greater than 1000 and the lower redshift bin
is mass-completed. Finally, there are 504, 530, and 342 galaxies in each bin,
resulting in 1376 galaxies in total for following analysis, in which 299
galaxies have spectroscopic redshifts and the remaining 1077 galaxies have a
mean uncertainty of photometric redshift of $\sim$0.097. The mean uncertainty
of log$M_{*}$ of all galaxies in our three bins is $\sim$0.145 and the mean
uncertainty of MF444W is $\sim$0.007.

Figure 1: The distribution of the MF444W and the $M_{*}$ of our galaxies with
redshift. The blue dots represent sources with high $\langle$S/N$\rangle$s and
valid morphological parameters (2068). The gray dots represent other sources
removed due to not meeting the conditions above (826). The horizontal dashed
black line represents our cut of $M_{*}$. The three different colored
rectangular regions divided by vertical dashed black lines represent our three
redshift bins, in which the numbers represent the number of galaxies.  Table
1: Summary of redshift bins

Redshift bin | Mean redshift | Total | Filter | Resolution | $\lambda_{\rm rf}$ |  $G-M_{20}$ classification
---|---|---|---|---|---|---
$z$ | $\left\langle z\right\rangle$ |  |  |  | (nm) | Merger | LTG | ETG
|  |  | F115W | F200W |  $\sim$570 | 72 (14.3%) | 336 (66.7%) | 96 (19.0%)
|  |  | F115W | F444W |  $\sim$570 | 33 (6.5%) | 426 (84.5%) | 45 (8.9%)
0.8-1.3 | 1.05 | 504 | F150W | F444W |  $\sim$760 | 25 (5.0%) | 435 (86.3%) | 44 (8.7%)
|  |  | F200W | F444W |  $\sim$1000 | 23 (4.6%) | 436 (86.5%) | 45 (8.9%)
|  |  | F277W | F444W |  $\sim$1340 | 22 (4.4%) | 428 (84.9%) | 54 (10.7%)
|  |  | F150W | F200W |  $\sim$570 | 111 (20.9%) | 329 (62.1%) | 90 (17.0%)
|  |  | F150W | F444W |  $\sim$570 | 36 (6.8%) | 460 (86.8%) | 34 (6.4%)
1.3-2.0 | 1.65 | 530 | F200W | F444W |  $\sim$760 | 26 (4.9%) | 469 (88.5%) | 35 (6.6%)
|  |  | F277W | F444W |  $\sim$1000 | 31 (5.8%) | 443 (83.6%) | 56 (10.6%)
|  |  | F356W | F444W |  $\sim$1340 | 37 (7.0%) | 438 (82.6%) | 55 (10.4%)
|  |  | F200W | F200W |  $\sim$570 | 92 (26.9%) | 201 (58.8%) | 49 (14.3%)
|  |  | F200W | F444W |  $\sim$570 | 19 (5.6%) | 309 (90.4%) | 14 (4.1%)
2.0-3.0 | 2.36 | 342 | F277W | F444W |  $\sim$760 | 17 (5.0%) | 307 (89.8%) | 18 (5.3%)
|  |  | F356W | F444W |  $\sim$1000 | 37 (7.0%) | 438 (82.6%) | 55 (10.4%)
|  |  | F444W | F444W |  $\sim$1340 | 28 (8.2%) | 293 (85.7%) | 21 (6.1%)

Note. — From left to right, the columns (1), (2), and (3) correspond to the
redshift range of each bin, mean redshift of all galaxies in the bin, and the
quantity of galaxy with log$M_{*}>9$. The columns (4), (5), and (6) correspond
to the NIRCam filter used for the measurement, the target filter used to match
PSF, and the rest-frame pivot wavelength (approximation) of the filter at the
mean redshift. The columns (7)-(9) are the quantity and fraction of different
types of galaxy classified by the $G-M_{20}$ diagram.

###  2.3 Non-parametric Morphology Measurement

Non-parametric morphology measurement is designed to capture the major
features of the underlying structures of galaxies, without any assumed light
distribution model (e.g. Sérsic profile). At high redshift universe, galaxies
are usually irregular and cannot be fitted by a simple model (Conselice,
2014). Non-parametric measurement is ideal for characterizing galaxy evolution
over many epochs and has been applied in many high-redshift works (e.g., Law
et al., 2012; Ferreira et al., 2022b; Treu et al., 2023).

We analyze all wide bands of the CEERS NIRCam images (F115W, F150W, F200W,
F277W, F356W, and F444W) by the statmorph code (Rodriguez-Gomez et al., 2019),
which calculates non-parametric morphological diagnostics such as the
Concentration-Asymmetry-Smoothness system (CAS, Conselice, 2003) and
Gini-$M_{20}$ system ($G$–$M_{20}$, Lotz et al., 2004). We have modified the
code of statmorph to improve performance, resulting in a 3-5 times improvement
in computation speed without sacrificing measurement accuracy. This
optimization was mainly accomplished through acceleration by Cython (Behnel et
al., 2011). Appendix A shows our optimization performance and a comparison
between the outputs of the original version of statmorph and our optimized
code. In addition, the modified version also includes the call to Source-
Extractor. Users can enter parameters of Source-Extractor to do photometry and
generate a segmentation map, or skip this step and provide their own
segmentation map. Furthermore, we also add some extended optional
morphological parameters (e.g., $\xi$, Papovich et al. 2003; $\Psi$, Law et
al. 2007; $G_{2}$, Rosa et al. 2018), which can be selectively calculated
according to the user’s specific needs. Finally, we change the flag used for
warnings to a bitmap format, which can store different types of warnings
rather than just indicating their presence or absence. All the above
optimizations and improvements are made to meet the requirements of the China
Space Station Telescope (CSST, Zhan, 2011; Gong et al., 2019) survey in the
future, so we rename the modified code to statmorph_csst555Optimized code is
available at https://github.com/pentyum/statmorph_csst and
https://gitee.com/pentyum/statmorph_csst. CSST is a 2-meter space telescope
planned to be launched in the early 2020s. It will operate in the same orbit
as the China Manned Space Station and survey 17,500 deg2 sky area over 10
years of operation. In this work, we directly use the segmentation map
provided by grizli as the input of statmorph_csst.

####  2.3.1 CAS Statistics

The most-widely used non-parametric morphological indicator is the CAS
statistics. It includes three parameters, concentration ($C$), asymmetry
($A$), and smoothness ($S$, also clumpiness). The definitions and the
corrections of these parameters have changed and developed over the past
decades (e.g., Kent, 1985; Abraham et al., 1996; Bershady et al., 2000;
Conselice, 2003; Conselice et al., 2003; Lotz et al., 2004; Conselice et al.,
2008; Law et al., 2012), and the definition given in Conselice (2014) is
generally adopted in recent literature.

####  2.3.2 $G-M_{20}$ Statistics

Another set of widely used non-parametric morphological indicators is the
$G-M_{20}$ classification system (Abraham et al., 2003; Lotz et al., 2004,
2008b). The measurement of these parameters does not need to provide a
predefined location of the galactic center, which is particularly desirable
for irregular and merging galaxies at high redshift.

To calculate the $G-M_{20}$ statistics, statmorph creates a new $G-M_{20}$
segmentation map that only selects pixels whose brightness are brighter than
the average brightness at major-axis radius of $r_{p}$. This procedure is the
same as Lotz et al. (2004), and its detailed implements by statmorph are
described in Section 4.4.1 of Rodriguez-Gomez et al. (2019).

###  2.4 PSF Matching

To accurately measure the morphological parameters of galaxies in each NIRCam
filter, it is essential to account for the difference of PSFs across filters.
Since JWST is a space telescope, its PSF is diffraction-limited PSF. The
longer the wavelength is, the larger the FWHM of the PSF is. According to the
simulation by WebbPSF (Perrin et al., 2014), FWHMs of the NIRCam bands vary
from 0″.037 (F115W) to 0″.140 (F444W), which may introduce systematic
wavelength-dependent bias in the measurement of morphological parameters
(e.g., Bershady et al., 2000; Lotz et al., 2004; Law et al., 2012), especially
in $A$ (Bignone et al., 2020). Therefore, it is necessary to perform PSF
matching, which creates convolution kernels from PSFs with small FWHM to PSFs
with large FWHM. The following in this subsection describes our PSF matching
procedure. The overall process is similar to Yao et al. (2022), but some
optimizations have been made.

The PSFs we use are from grizli-psf-
library666https://github.com/gbrammer/grizli-psf-library/tree/main/ceers.
These PSFs are sampled at a pixel scale of 0″.04. When we apply these PSFs to
the SW images of NIRCam, whose pixel scale is 0″.02. We resample the PSFs to
match the pixel scale.

We make use of the photutils777https://github.com/astropy/photutils (Bradley
et al., 2022) package to create a convolve kernel from short-wavelength image
to long-wavelength. We convolve all images with wavelength shorter than F444W
to the F444W, whose PSF is the largest among all filters. In order to study
the influence of different PSFs, besides the convolution to F444W, we also
convolve the F115W and F150W image to the F200W whose PSF is the largest among
the SW detectors.

We use an optimized method to create convolve kernels. Aniano et al. (2011)
provide two parameters to measure kernel performance. One is $D$, describing
the accuracy in redistribution of PSF power. The other is $W_{-}$, describing
the negative values of the kernel. The lower the two values, the better. For a
good PSF matching performance, $D\textless$0.1 and $W_{-}\textless$1 are
required. We take this as our optimization goal and build a simple loss
function with $D$ and $W_{-}$ as parameters. We find the parameter of each
window function888https://photutils.readthedocs.io/en/stable/psf_matching.html
provided by photutils that minimizes the loss function, and then compare the
minimum value of the loss function generated by different types of window
function, and then select the lowest one. Table 2 lists the types of window
functions we finally adopt and the corresponding $D$ and $W_{-}$ during the
PSF matching process.

Table 2: Window type and the kernel performance of PSF matching

Original filter | Target filter | Window type | $D$ | $W_{-}$
---|---|---|---|---
F115W | F444W | Tukey | 0.032 | 0.069
F150W | F444W | CosineBell | 0.047 | 0.084
F200W | F444W | CosineBell | 0.074 | 0.068
F277W | F444W | TopHat | 0.017 | 0.345
F356W | F444W | TopHat | 0.019 | 0.629
F115W | F200W | CosineBell | 0.080 | 0.341
F150W | F200W | CosineBell | 0.088 | 0.386

Note. — Photutils provides the following window function types: HanningWindow,
TukeyWindow, CosineBellWindow, SplitCosineBellWindow, and TopHatWindow. Here
we have tried all the window function type, except the SplitCosineBellWindow,
because it contains multiple parameters, which is not convenient to use
scipy.optimize.minimize_scalar, and other types have been enough to meet our
optimization goals.

Then we use the reproject999https://github.com/astropy/reproject package to
project the convolved SW images (0″.02 pixel-1) to the same sampling rate as
the LW images (0″.04 pixel-1). In this way, we can get multi-band images with
the same spatial resolution and sampling rate.

Given the redshift range ($z$=0.8-3.0) of the galaxies measured by us, the
physical scale of given angular size does not vary much (7.5-8.5 kpc
arcsec-1), so we do not attempt to make further redshift-dependent adjustment
of the PSF. The physical scale corresponding to the FWHM of the PSF in F444W
band is 1.05-1.19 kpc, and the physical scale of the pixel in our final images
is 0.30-0.34 kpc.

##  3 Results and Discussion

###  3.1 Morphology Evolution with Rest-frame Wavelength

It is well known that the color of different age of stellar population is
different. At shorter wavelengths (e.g. $U$-band), the morphology and
quantitative structures trace the distribution of younger stellar populations,
whereas longer wavelengths (e.g. $R$-band) tend to trace older stellar
populations that usually dominate the total stellar mass. Furthermore, the
distribution of dust and emission nebula can also affect the morphology at
different $\lambda_{\rm rf}$ (Baes et al., 2020; Nersesian et al., 2023).
Therefore, to ensure a fair comparison of galaxy morphologies at different
redshift, it is necessary to make a “K-correction” (Kuchinski et al., 2001;
Conselice, 2003) when analyzing them. Wuyts et al. (2012) have analyzed some
morphological parameters ($R_{e}$, $C$, $G$, and, $M_{20}$) of massive SFGs in
the HST GOODS-South field ($z$=0.5-2.5), and finds that SFGs to be more
concentrated, smaller, and smoother in red bands or their stellar mass
distribution than in blue bands. These differences can be explained by the
inside-out growth scenario in combination with a ”Christmas tree” model
(Lowenthal et al., 1997).

Figure 2: The four morphological parameters, $C$, $A$, $G$, and $M_{20}$ vary
with $\lambda_{\rm rf}$ in the three redshift bins colored by red, green, and
blue individually. The lines and dots represent the mean value of all galaxies
in each bin. The semi-transparent filled areas indicate their 32% and 68%
(1$\sigma$) quantiles. The thick solid lines represent the values that have
been matched to F444W PSF. The thin dashed lines represent the values that
have been matched to F200W PSF, while the thin solid lines represent values
without any PSF matching. The arrows at long wavelength ends are the rough
corrections of the underestimated $A$ in F444W.

Thanks to the wide NIR wavelength range of JWST instruments, we are able to
investigate the morphology of galaxies and study their evolution at different
redshifts at similar $\lambda_{\rm rf}$. Figure 2 illustrates how the average
of the four morphological parameters, $C$, $A$, $G$, and $M_{20}$, vary with
$\lambda_{\rm rf}$. Our analysis reveals that all morphological parameters of
all galaxies, regardless of their redshift, are influenced by the
$\lambda_{\rm rf}$. Specifically, the morphological parameters have a stronger
wavelength dependence at shorter wavelengths ($<$1000nm), while at longer
wavelengths ($>$1000nm), the parameters become nearly independent of the
$\lambda_{\rm rf}$.

At the same resolution, we observe that the $C$ values noticeably increase as
the wavelength increases. Furthermore, the $C$ values exhibit little
difference among galaxies at different redshifts at short $\lambda_{\rm rf}$.
However, at longer $\lambda_{\rm rf}$ ($\sim$1340nm), the $C$ values of
galaxies at $z$=2.0-3.0 are 0.076 lower than those of galaxies at $z$=0.8-1.3.
For each redshift bin we examined ($z$=0.8-1.3, $z$=1.3-2.0, $z$=2.0-3.0), the
$C$ values increase by $\sim$0.18 ($\sim$6.8%), $\sim$0.18 ($\sim$6.5%), and
$\sim$0.12 ($\sim$4.5%), respectively, from rest-frame $\sim$570nm to
$\sim$1340nm. Wuyts et al. (2012) have also found that the wavelength
dependence of $C$ at high redshift is weaker than low redshift, but the
amplitude of the change in our work is lower than that obtained by Wuyts et
al. (2012). This may be because we calculate the difference between the
optical and NIR bands, while Wuyts et al. (2012) calculate the difference
between the UV and optical bands, but $C$ is more sensitive to $\lambda_{\rm
rf}$ at the short wavelength end.

The $A$ values also vary with the $\lambda_{\rm rf}$. We notice that the $A$
values all “drop” suddenly at the longest $\lambda_{\rm rf}$ in all redshift
bins. Bignone et al. (2020) have pointed out that the $A$ value may be
dominated by the $A_{bgr}$ when the $\langle$S/N$\rangle$ is low. The
excessive correction can result in artificially low or even negative asymmetry
values. We have compared the $A_{bgr}$ of different bands, and find that the
$A_{bgr}$ of F444W corresponding to the longest $\lambda_{\rm rf}$, is
systematically $\sim$0.02 higher than that of other bands indeed, most likely
due to the relatively lower $\langle$S/N$\rangle$ of this band than other
bands. This may result in an underestimation of $A$ at the F444W. This trend
is consistent with the previous quantitative studies on the robustness of
morphological measurement (e.g., Lotz et al., 2004; Thorp et al., 2021; Treu
et al., 2023). If we ignore F444W bands, for each redshift bin we examined,
($z$=0.8-1.3, $z$=1.3-2.0, $z$=2.0-3.0), the $A$ values decrease by
$\sim$0.022 ($\sim$13.3%), $\sim$0.009 ($\sim$5.8%), and $\sim$0.014
($\sim$8.4%) respectively, from rest-frame $\sim$570nm to $\sim$1000nm.

The $G$ parameter is similar to $C$, in the sense that higher $G$ means more
concentrated light distribution. Therefore, $G$ is expected to have a similar
wavelength-dependent trend with $C$. For each redshift bin, the $G$ values
only increase by $\sim$0.016 ($\sim$3.3%), $\sim$0.016 ($\sim$3.2%), and
$\sim$0.014 ($\sim$2.7%), respectively, from rest-frame $\sim$570nm to
$\sim$1340nm. The $M_{20}$, which is also an indicator of concentration, shows
an anti-correlation with $\lambda_{\rm rf}$. For each redshift bin, the
$M_{20}$ values decrease by $\sim$0.072 ($\sim$4.4%), $\sim$0.085
($\sim$5.3%), and $\sim$0.034 ($\sim$2.2%), respectively, from rest-frame
$\sim$570nm to $\sim$1340nm. The trends of $G$ and $M_{20}$ are similar to
Wuyts et al. (2012), but the amplitude is lower. The reason is the same as
$C$. The wavelength dependence of $G$ and $M_{20}$ at high redshift is also
weaker than low redshift.

Figure 3: The $C-A$ and $G-M_{20}$ diagrams at rest-frame $\sim$570nm and
$\sim$1000nm. The first column represents the values that have been matched to
F200W PSF, and other columns represent the values that have been matched to
F444W PSF. The red, green, and blue crosses are present the mean values and
1$\sigma$ quantiles of $z$=0.8-1.3, $z$=1.3-2.0, and $z$=2.0-3.0 bins,
individually. The sloped lines on the $C-A$ diagrams are the boundaries
between LTG and intermediate galaxies defined by Bershady et al. (2000), and
the horizontal lines are used to divide mergers and LTG at $A$=0.35 defined by
Conselice (2003). The classification of galaxies on the $G-M_{20}$ diagrams is
based on Lotz et al. (2008b). The solid orange lines ($G$=0.14$M_{20}$+0.80)
show the boundaries between LTG and ETG, and the blue lines
($G$=$-$0.14$M_{20}$+0.33) divide mergers and non-mergers. The dashed orange
lines are the boundaries between LTG and ETG described by Equation 1, which
have corrected the $\lambda_{\rm rf}$-dependent morphology biases.

There is currently ongoing debate about the evolution of morphological
parameters at $\lambda_{\rm rf}$. Some studies at $z<$6 have reported similar
trends to our results (e.g., Conselice et al., 2008; Wuyts et al., 2012),
while others at higher redshifts ($z>$6) have found no obvious wavelength
dependence (e.g., Yang et al., 2022; Treu et al., 2023). Treu et al. (2023)
suggests that the lack of dependence may be due to the shorter time available
for star formation in high-redshift galaxies. The star formation of extremely
high redshift galaxies selected by Lyman-break seems to be global. The old and
young stellar populations are mixed together spatially, and cannot be
separated by any filters. In contrast, galaxies at lower redshifts especially
in local universe have more complex star formation histories, resulting in
differences in the spatial distributions of old and young stellar populations.
Our findings of the weaker wavelength dependence in the high-redshift bins are
in line with this trend.

###  3.2 $C-A$ and $G-M_{20}$ Diagrams

Combinations of the morphological parameters can be used to classify galaxies
into different types. The $C-A$ (Bershady et al., 2000; Conselice, 2003) and
$G-M_{20}$ (Lotz et al., 2004, 2008b) diagrams are the most common non-
parametric methods for classification. In this section, we utilize these
diagrams to examine how the redshift, the $\lambda_{\rm rf}$, and the spatial
resolution impact the distribution of galaxies in morphological parameter
space.

####  3.2.1 $C-A$ Diagram

The first row of Figure 3 shows the $C-A$ diagrams of our galaxy sample at
different $\lambda_{\rm rf}$ in different resolution. The diagram features the
boundary that divides LTG and intermediate galaxy as proposed by Bershady et
al. (2000) along with the threshold separating mergers from normal galaxies at
$A$=0.35, defined by (Conselice, 2003). From these diagrams we observe that
the average position of the galaxies (colored crosses) changes minimally with
redshift but is affected by the $\lambda_{\rm rf}$ and resolution. When the
$\lambda_{\rm rf}$ increases, C increases while A decreases, resulting in the
galaxies shifting towards the lower right of the $C-A$ diagram. When the
spatial resolution decreases, both $C$ and $A$ decrease, causing the galaxies
to move in the lower left direction of the $C-A$ diagram. The impact of PSF
will be discussed in more details in Section 3.3.

####  3.2.2 $G-M_{20}$ Diagram

Lotz et al. (2008b) have defined a simple partitioning of the $G-M_{20}$
diagram that can be used to roughly separate ETG, LTG, and merging galaxies at
$z$=0.2-0.4. The second row of Figure 3 shows the $G-M_{20}$ diagrams of our
galaxy sample at different $\lambda_{\rm rf}$ in different resolution. We can
see that, unlike Lotz et al. (2008b), our sample does not exhibit a clear
bimodality between LTG and ETG, which is similar to what has been seen at
other redshifts and in simulations (e.g., Lotz et al., 2008a; Bignone et al.,
2020; Kartaltepe et al., 2023). Instead, most of our galaxies are located in
the LTG area in any $\lambda_{\rm rf}$, just like what is classified by the
$C-A$ diagram in the first row of Figure 3. The 7th-9th columns of Table 1
list the quantity and fraction of each classification of the $G-M_{20}$
diagram.

The relationship between redshift and $G-M_{20}$ classification is also
examined in conjunction with the results presented in Figure 2. As redshift
increases, $M_{20}$ and $G$ exhibit an upward trend, leading to a shift of
galaxies towards the upper right region of the $G-M_{20}$ diagram.
Additionally, the fraction of mergers also increases. We note that the
fraction of ETG decreases with the increase of redshift in all resolutions and
$\lambda_{\rm rf}$, but the fraction of LTG decreases significantly in high
spatial resolution (F200W) but increases in low resolution (F444W) in
contradictory. This can be attributed to the large number of LTGs being
reclassified as mergers in high resolution. This will be discussed more in
Section 3.3.

Figure 4: The movements of the mean value of galaxy morphological parameters
on the $G-M_{20}$ diagram with the increase of $\lambda_{\rm rf}$. Different
colors represent different redshift bins, and black dots represent the overall
sample.

Furthermore, the $\lambda_{\rm rf}$ is an additional factor that influences
galaxy distribution in the $G-M_{20}$ diagram. Figure 2 and 4 show that an
increase in $\lambda_{\rm rf}$ results in a decrease in $M_{20}$ and an
increase in $G$, causing galaxies to move towards the upper left region of the
$G-M_{20}$ diagram with a slope of -0.23$\pm$0.03. The path of the movement is
almost parallel to the boundary between merger and non-merger, with an angle
of less than 5$\arcdeg$. In terms of the proportion of different types shown
in Table 1, the fraction of LTGs decreases while the fraction of ETGs
increases, with little change observed in the fraction of mergers. This
reflects that the long rest-frame wave is more able to track the distribution
of the old stellar population.

Because the division of Lotz et al. (2008b) is based on the rest-frame $B$
($\sim$450nm). To reduce potential wavelength-dependent biases in the
morphological classification of galaxies, we suggest that the boundaries
between LTG and ETG (the orange line in Figure 3) should rise with the
increase of $\lambda_{\rm rf}$ to do the “K-correction”. A simple modified
form is

$G=0.14M_{20}+0.8+\alpha\left(\frac{\lambda_{\rm rf}~{}\rm nm}{450~{}\rm
nm}-1\right).$ (1)

The last item is the “K-correction” item and $\alpha$ is the correction
factor. If the classification of a galaxy is not affected by $\lambda_{\rm
rf}$, the ratio of LTG to ETG of the total sample should be stable at all
$\lambda_{\rm rf}$. We obtain a correction factor of $\alpha=0.006$, which
minimize the relative standard deviation of the ratios of LTG to ETG at
different wavelengths for all redshift bins of our sample. When using the new
corrected boundary, the ratio of LTG to ETG for our all galaxies is about
13.8:1. The ratio of our three redshift bins are about 11.0:1, 13.4:1, and
24.6:1, respectively. The corrected boundaries are shown as dashed orange
lines in Figure 3. It should be noted that the above corrections are all
carried out in the resolution of F444W (FWHM$\simeq$1.1kpc physically). This
correction may not be applicable to other redshift ranges due to the potential
different wavelength dependence of morphological parameters at different
redshift. Another additional note is that the wavelength dependence of
morphological parameters is not linear for an expanded wavelength range (see
Figure 3 in Baes et al., 2020), meaning that the corrected form of Equation 1
may not be applicable to a wide wavelength range.

###  3.3 PSF Influences on Morphological Measurements

In this section, we aim to investigate the impact of spatial resolution on the
non-parametric morphology of galaxies. We not only calculate the morphological
parameters convolved to F444W, but also calculate the parameters that have not
been convolved and only convolved to F200W. The thin lines in Figure 2 show
the effect of different spatial resolutions on morphological parameters, and
the first two rows of each bin in Table 1 show how the spatial resolution
affects the classification of galaxies. We find that the smoothing effect of
PSF can reduce $C$, $A$, and $G$ significantly ($\gtrsim$10%). The
classification of galaxies is obviously affected by the size of PSF. The
merger fraction in F200W resolution can be $\gtrsim$2 times larger than that
in F444W resolution, regardless of the classification method of $C-A$ or
$G-M_{20}$. This is because if the spatial resolution decrease (i.e. the PSF
FWHM increases), $C$, $A$, $M_{20}$, and $G$ would both decrease, the galaxies
would move to the lower left direction in the $C-A$ and $G-M_{20}$ diagram.
So, the fraction of merger would decrease significantly. This result can also
explain why the line divided merger and non-merger galaxies defined by Lotz et
al. (2008b) is moved down compared with Lotz et al. (2004). Because the sample
of Lotz et al. (2004) is closer than Lotz et al. (2008b), and has a better
spatial resolution. In addition, at the same time, the fraction of LTG has
also increased with the decrement of resolution, because many of the merging
galaxies in F200W resolution are classified as LTG in F444W resolution.

Bignone et al. (2020) have used simulated images from the EAGLE simulation
(Schaye et al., 2015; Trayford et al., 2017) to study the effect of spatial
resolution on non-parametric morphology. They have got the same trend as us,
but they have found that $C$ only shows a small systematic deviation of less
than 2% in even for the worst FWHM of 1.5 kpc, which is much smaller than the
variation we find. The average value of $C$ at $z$=2.0-3.0 in F115W resolution
(FWHM$\simeq$0.33kpc) and F115W band (i.e. $\lambda_{\rm rf}\simeq$340nm) is
$\sim$3% and $\sim$10% higher than that of F200W (FWHM$\simeq$0.54kpc) and
F444W (FWHM$\simeq$1.18kpc) resolution respectively. This may be due to the
differences in galaxy sample, where our sample at higher redshift consists of
more irregular and clumpier galaxies that are more affected by PSF than their
sample at $z$=0.1 which contains more elliptical ETGs that are less affected.
This morphology-dependent statistics can also be found in our sample. For each
redshift bin we examined ($z$=0.8-1.3, $z$=1.3-2.0, $z$=2.0-3.0), the $C$
values at rest-frame $\sim$570nm increase by $\sim$6%, $\sim$7%, and $\sim$8%,
respectively, from the PSF of F444W to F200W. Therefore, it is particularly
important to perform PSF matching for morphological measurements of high
redshift galaxies.

###  3.4 Uncertainties of Morphological Measurements

It is difficult to obtain the uncertainties of morphological parameters
directly, so we adopt the empirical relationship between $\langle$S/N$\rangle$
and error ($\Delta$) to estimate the uncertainties of our morphological
parameters. Lotz et al. (2004) have found that $G$, $M_{20}$, and $C$ are
reliable to within $\sim$10% ($\Delta\leq$ 0.05, 0.2, and 0.3, respectively)
for galaxy images with $\langle$S/N$\rangle\geq$2, and $A$ shows offsets less
than 0.1 at $\langle$S/N$\rangle\geq$5.

During our sample selection process (Section 2.2), we only include sources
with $\langle$S/N$\rangle\geq$2 in all bands we analyzed, so the values of
$\Delta$ given by Lotz et al. (2004) above can be considered as the upper
limits of uncertainties. However, the average $\langle$S/N$\rangle$ of our all
galaxies is $\gtrsim$10 in all six bands, indicating that the actual typical
uncertainties should be much smaller than the upper limits. Additionally,
increasing the $\langle$S/N$\rangle$ cutoff to 5 during the sample selection
process would not significantly affect the results (seeing Appendix B.1). We
also conduct a simulation moving low-redshift galaxies to higher redshift in
Appendix B.2. Therefore, we can conclude that the wavelength evolution we
observed above in Figure 2 and 4 is significantly existed.

###  3.5 Morphology Evolution for Different $M_{*}$

Figure 5: The four morphological parameters, $C$, $A$, $G$, and $M_{20}$
evolve with redshift at four matched $\lambda_{\rm rf}$ for the two stellar
mass bins, 9$<$logM${}_{*}<$10 and logM${}_{*}\geqslant$10, colored by blue
and red individually. The lines and dots represent the mean value in each bin.
The semi-transparent filled areas indicate their 32% and 68% (1$\sigma$)
quantiles. The first column represents the values that have been matched to
F200W PSF, and the following columns represent the values that have been
matched to F444W PSF. The arrows at the longest wavelength diagrams are the
rough corrections of the underestimated $A$ in F444W.

Galaxies with different $M_{*}$ usually have different evolutionary histories
(Madau & Dickinson, 2014). We divide our sample into two $M_{*}$ bins
(9$<$log$M_{*}<$10 and log$M_{*}\geqslant$10) to investigate the difference in
morphological evolution between high-mass and low-mass galaxies. Figure 5
shows how the four morphological parameters, $C$, $A$, $G$, and $M_{20}$
evolve with redshift for different $M_{*}$. We can see that at
$z\simeq$0.8-3.0, galaxies with log$M_{*}\geqslant$10 have higher $C$, $A$,
and $G$, and lower $M_{20}$ than galaxies with log$M_{*}<$10 at any redshift
bins, $\lambda_{\rm rf}$, and spatial resolutions. These differences are more
significant at lower redshift bins, signifying different evolutionary paths
for high- and low-mass galaxies. This evolutionary trend could not be caused
by incompleteness of $M_{*}$ in the highest redshift bin, because the low-mass
galaxies own lower values of $C$ and $G$, and higher values of $M_{20}$, while
we find that the high-redshift bin which lacks low-mass galaxies has already
owned these features. If we make the completeness of $M_{*}$ higher, this
trend will only become stronger, rather than weaken.

In particular, for the evolution of $C$, our results are different from
Whitney et al. (2021). They have found that both high-mass and low-mass bins
show a decrease in $C$ with cosmic time. This difference may be due to
different observation bands and different data processing procedures. We
discover that in blue $\lambda_{\rm rf}$ bands ($\lambda_{\rm rf}<$1000nm),
the $C$ of galaxies with log$M_{*}\geqslant$10 increases with cosmic time,
while the $C$ of galaxies with log$M_{*}<$10 decreases with cosmic time. In
red rest-frame bands ($\lambda_{\rm rf}\geqslant$1000nm), the $C$ of both
high-mass and low-mass galaxies increases with cosmic time. The $G$ of low-
mass galaxies also demonstrates an evolutionary difference, albeit less
pronounced than $C$’s. This evolutionary differences at different
$\lambda_{\rm rf}$ can be explained as a scenario that the low-mass galaxies
is still undergoing continuous star formation in the outside of their bulges,
while the high-mass galaxies have already been assembled, and are experiencing
merging mainly during this period (Conselice, 2014; Madau & Dickinson, 2014)
with higher $A$. The inside-out formation process leads to different
evolutionary trends in different bands.

##  4 Summary and Conclusions

Benefiting from the multi-wavelength broadband images from JWST/NIRCam, we
have made use of the data processed by grizli to investigate the non-
parametric morphology of galaxies across the redshift ranges of
$z\simeq$0.8-3.0 in the CEERS field. Our analysis have taken into account the
PSF difference among different wave bands, and we have measured $C$, $A$, $G$,
and $M_{20}$ of a total of 1376 galaxies divided into three redshift bins. The
main results are as follows:

  1. 1.

Non-parametric morphological parameters change with the $\lambda_{\rm rf}$ at
$z\simeq$0.8-3.0, especially at the short wavelength end. Different parameters
and different redshift bins have different sensitivities on $\lambda_{\rm
rf}$. The wavelength dependence of high redshift is weaker than that at low
redshift. $C$ and $A$ are more sensitive than $G$ and $M_{20}$.

  2. 2.

The wavelength dependence of the morphological parameters may affect their
usage in the classification of galaxy types, especially the distinction
between ETG and LTG. We suggest that the boundaries between LTG and ETG in
$G-M_{20}$ diagram should shift upwards as the rest-wavelength increases. We
give a form of correction (Equation 1, $\alpha$=0.006) that can keep the ratio
of LTG to ETG stable at different $\lambda_{\rm rf}$ for our current sample.

  3. 3.

Spatial resolution also significantly affects the morphological parameters and
classification of galaxies, especially the fraction of merger. The merger
fraction in F200W resolution can be $\gtrsim$2 times larger than that in F444W
resolution, regardless of the diagnostic parameter combinations we adopt. The
influences of PSF become stronger as the redshift increases.

  4. 4.

The morphological parameters (especially the concentration parameter $C$) of
high-mass (log$M_{*}\geqslant$10) and low-mass (9$<$log$M_{*}<$10) galaxies
have a different evolution with redshift. At high redshift ($z\gtrsim$2), the
difference between high-mass and low-mass galaxies is small, while at low
redshift ($z\lesssim$1.5), the difference appears. In short $\lambda_{\rm rf}$
bands ($\lambda_{\rm rf}<$1000nm), the $C$ of high-mass and low-mass galaxies
have opposite evolutionary trend, while in longer rest-frame bands
($\lambda_{\rm rf}\geqslant$1000nm), the trend becomes similar for both high-
and low-mass galaxies. This suggests different evolution processes between
high-mass and low-mass galaxies. The star formation of low-mass galaxies still
continues in their disk, while galaxy mergers have played an important role in
the evolution of high-mass galaxies during the period of $z\simeq$0.8-3.0.

This work is supported by the Strategic Priority Research Program of Chinese
Academy of Sciences (Grant No. XDB 41000000), the National Science Foundation
of China (NSFC, Grant No. 12233008, 11973038, 11973039), the China Manned
Space Project (No. CMS-CSST-2021-A07), the Cyrus Chun Ying Tang Foundations
and the Frontier Scientific Research Program of Deep Space Exploration
Laboratory.

##  Appendix A Performance evaluation of starmorph_csst

We make use of COSMOS/HST image to evaluate the multiplier of performance
improvement achieved by statmorph_csst. We measured the CAS and $G-M_{20}$
morphological parameters of galaxies in the F814W band 064 image using both
statmorph and statmorph_csst, respectively, and record the time consumed for
each galaxy. Finally, we get the morphological parameters of 1889 galaxies
without any error flags. Figure 6 shows the distribution histogram of the
multipliers of performance, with an average optimization multiplier of
$\sim$5.57, and a median of $\sim$3.34.

Figure 6: The distribution of the optimized multipliers, which is the ratio of
the time consumed by the original version (statmorph) to the time consumed by
the optimized version (statmorph_csst).

Although we use Cython to rewrite the code of statmorph for optimization,
considering that we reduce the max times of iteration during solving
equations, and Cython and Python have different precision in calculating
floating numbers, it is also necessary to verify whether the results of the
optimized version differ significantly from the original version. Here we use
the same sample above in the COSMOS field for testing. The results shown in
Figure 7 indicate that our optimization has almost no effect on the
calculation of morphological parameters. Additionally, some other
optimizations are also included in statmorph_csst, but it will not be enabled
by default due to their significant increase in $\sigma$.

Figure 7: The comparison of the results of five main morphological parameters
between statmorph_csst and the original statmorph. $\sigma$ on the upper left
is the standard deviation.

##  Appendix B The Influence of $\langle$S/N$\rangle$ on the Robustness of
Measurement

We employ two methods to assess the extent to which noise introduces
significant biases in our morphological measurements. These methods include
increasing the cutoff of $\langle$S/N$\rangle$ and conducting simulations for
observations at higher redshifts.

###  B.1 Increasing the Cutoff of $\langle$S/N$\rangle$

To reduce the uncertainties caused by noise and evaluate the influence of
uncertainties on morphological parameter measurement. Here we increase the
cutoff of $\langle$S/N$\rangle$ increases to 5 during the sample selection
(Section 2.2). The sample size is reduced to 770, and there are 378, 273, and
119 galaxies in each redshift bin. Figure 2 shows how the average of the four
morphological parameters vary with $\lambda_{\rm rf}$. Due to the discarding
of low $\langle$S/N$\rangle$ sources, the underestimation of $A$ in F444W is
weakened. So the arrows at the longest wavelength diagrams which present the
rough corrections of the underestimation are shortened to 0.015. The
wavelength evolution shown in Figure 8 is not significantly different from
Figure 2.

Figure 8: Same as Figure 2, but the cutoff of $\langle$S/N$\rangle$ increases
to 5.

###  B.2 Moving the Low-redshift Images to Higher Redshift

In addition, We conduct simulations to move the images of galaxies from the
low-redshift bin ($z$=0.8-1.3) to the mid-redshift ($z$=1.3-2.0) and high-
redshift bin ($z$=2.0-3.0). The steps are as follows:

  1. 1.

We cut out the F115W stamp images of the 504 galaxies in the low-redshift bin.
The mid-redshift images will be observed in F150W, and the high-redshift
images will be observed in F200W.

  2. 2.

We resample the images and weight images to match the angular diameters at the
target redshift. This relative variation of angular diameter is very weak,
less than $\sim 10\%$.

  3. 3.

We darken the resampled images to match the surface brightness at the target
redshift. The dimming coefficient is

$k=\left[\frac{D_{A}(z_{\rm target})}{D_{L}(z_{\rm
target})}\times\frac{D_{L}(z_{\rm origin})}{D_{A}(z_{\rm
origin})}\right]^{2},$ (B1)

where $D_{L}$ is the luminosity distance, $D_{A}$ is the angular diameter
distance, $z_{\rm origin}$ is the original redshift before moving, and $z_{\rm
target}$ is the target redshift. At high redshift, the cosmological dimming
can have a significant impact on the surface brightness. For the high-redshift
bin, the $k$ is about 0.1.

  4. 4.

Due to the fact that the step above also weaken the noise of F115W images, It
is necessary to add additional noise to the darkened images. The additional
noise is generated from normal distribution whose $\sigma=\sqrt{\sigma_{\rm
targetband}^{2}-k^{2}\sigma_{\rm F115W}^{2}}$. The $\sigma_{\rm band}$ is
obtained from the weight images of bands. Thus, the simulated mid-redshift and
high-redshift images from F115W will own the same noise level as the F150W and
F200W respectively.

  5. 5.

We convolve the original F115W images and the simulated higher-redshift images
to the F444W PSF with the kernel created in Section 2.4, and re-measure their
morphological parameters by statmorph_csst. We subtract the parameters before
moving from the parameters after moving to obtain the $\Delta$ (i.e.,
$\Delta$=simulated $-$ original).

The results shown in Figure 9 and Table 3 indicate that when moving to high
redshift, the $\langle$S/N$\rangle$ deteriorates, and random errors increase
significantly, but the systematic errors show minor changes. Since our results
here represent the average values of numerous galaxies’ parameters, the main
consideration is the impact caused by systematic errors. Therefore, it can be
proven that the evolution of morphological parameter with $\lambda_{\rm rf}$
and redshift is significant and reliable.

Figure 9: The systematic error (mean of $\Delta$) of $C$,$A$,$G$, and $M_{20}$
of the F115W images of galaxies in the low-redshift bin. The error bars are
the random error (standard deviation of $\Delta$). Detailed data is listed in
Table 3.  Table 3: Systematic and random error (mean and standard deviation of
$\Delta$)

Parameter | Systematic error | Random error
---|---|---
Move to F150W | Move to F200W | Move to F150W | Move to F200W
$C$ | -0.0320 | -0.0331 | 0.1332 | 0.2425
$A$ | -0.0448 | -0.0741 | 0.0816 | 0.1575
$G$ | -0.0006 | 0.0050 | 0.0259 | 0.0388
$M_{20}$ | 0.0139 | 0.0318 | 0.0936 | 0.1591

Note. — Move to F150W/F200W: move the low-redshift F115W images ($\lambda_{\rm
rf}\sim 570nm$) to the mid-redshift/high-redshift bin, which will be able to
observed in the F150W/F200W filters.

## References

  * Abraham et al. (1996) Abraham, R. G., van den Bergh, S., Glazebrook, K., et al. 1996, ApJS, 107, 1, doi: 10.1086/192352
  * Abraham et al. (2003) Abraham, R. G., van den Bergh, S., & Nair, P. 2003, ApJ, 588, 218, doi: 10.1086/373919
  * Aniano et al. (2011) Aniano, G., Draine, B. T., Gordon, K. D., & Sandstrom, K. 2011, PASP, 123, 1218, doi: 10.1086/662219
  * Baes et al. (2020) Baes, M., Nersesian, A., Casasola, V., et al. 2020, A&A, 641, A119, doi: 10.1051/0004-6361/202038470
  * Baldry et al. (2004) Baldry, I. K., Glazebrook, K., Brinkmann, J., et al. 2004, ApJ, 600, 681, doi: 10.1086/380092
  * Barbary (2016) Barbary, K. 2016, The Journal of Open Source Software, 1, 58, doi: 10.21105/joss.00058
  * Behnel et al. (2011) Behnel, S., Bradshaw, R., Citro, C., et al. 2011, Computing in Science and Engineering, 13, 31, doi: 10.1109/MCSE.2010.118
  * Bershady et al. (2000) Bershady, M. A., Jangren, A., & Conselice, C. J. 2000, AJ, 119, 2645, doi: 10.1086/301386
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Bignone et al. (2020) Bignone, L. A., Pedrosa, S. E., Trayford, J. W., Tissera, P. B., & Pellizza, L. J. 2020, MNRAS, 491, 3624, doi: 10.1093/mnras/stz3014
  * Bignone et al. (2017) Bignone, L. A., Tissera, P. B., Sillero, E., et al. 2017, MNRAS, 465, 1106, doi: 10.1093/mnras/stw2788
  * Bradley et al. (2022) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2022, astropy/photutils: 1.5.0, 1.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.6825092
  * Brammer et al. (2022) Brammer, G., Strait, V., Matharu, J., & Momcheva, I. 2022, grizli, 1.5.0, Zenodo, Zenodo, doi: 10.5281/zenodo.6672538
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503, doi: 10.1086/591786
  * Collaboration et al. (2022) Collaboration, A., Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Conselice (2003) Conselice, C. J. 2003, ApJS, 147, 1, doi: 10.1086/375001
  * Conselice (2006) —. 2006, MNRAS, 373, 1389, doi: 10.1111/j.1365-2966.2006.11114.x
  * Conselice (2014) —. 2014, ARA&A, 52, 291, doi: 10.1146/annurev-astro-081913-040037
  * Conselice et al. (2003) Conselice, C. J., Bershady, M. A., Dickinson, M., & Papovich, C. 2003, AJ, 126, 1183, doi: 10.1086/377318
  * Conselice et al. (2008) Conselice, C. J., Rajgor, S., & Myers, R. 2008, MNRAS, 386, 909, doi: 10.1111/j.1365-2966.2008.13069.x
  * Davis et al. (2007) Davis, M., Guhathakurta, P., Konidaris, N. P., et al. 2007, ApJ, 660, L1, doi: 10.1086/517931
  * Domínguez Sánchez et al. (2018) Domínguez Sánchez, H., Huertas-Company, M., Bernardi, M., Tuccillo, D., & Fischer, J. L. 2018, MNRAS, 476, 3661, doi: 10.1093/mnras/sty338
  * Fang et al. (2023) Fang, G., Ba, S., Gu, Y., et al. 2023, AJ, 165, 35, doi: 10.3847/1538-3881/aca1a6
  * Ferrari et al. (2015) Ferrari, F., de Carvalho, R. R., & Trevisan, M. 2015, ApJ, 814, 55, doi: 10.1088/0004-637X/814/1/55
  * Ferreira et al. (2022a) Ferreira, L., Adams, N., Conselice, C. J., et al. 2022a, ApJ, 938, L2, doi: 10.3847/2041-8213/ac947c
  * Ferreira et al. (2022b) Ferreira, L., Conselice, C. J., Sazonova, E., et al. 2022b, arXiv e-prints, arXiv:2210.01110, doi: 10.48550/arXiv.2210.01110
  * Freeman (1970) Freeman, K. C. 1970, ApJ, 160, 811, doi: 10.1086/150474
  * Freeman et al. (2013) Freeman, P. E., Izbicki, R., Lee, A. B., et al. 2013, MNRAS, 434, 282, doi: 10.1093/mnras/stt1016
  * Gong et al. (2019) Gong, Y., Liu, X., Cao, Y., et al. 2019, ApJ, 883, 203, doi: 10.3847/1538-4357/ab391e
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hubble (1926) Hubble, E. P. 1926, ApJ, 64, 321, doi: 10.1086/143018
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80, doi: 10.1051/0004-6361/202142663
  * Kartaltepe et al. (2023) Kartaltepe, J. S., Rose, C., Vanderhoof, B. N., et al. 2023, ApJ, 946, L15, doi: 10.3847/2041-8213/acad01
  * Kauffmann et al. (2003) Kauffmann, G., Heckman, T. M., White, S. D. M., et al. 2003, MNRAS, 341, 54, doi: 10.1046/j.1365-8711.2003.06292.x
  * Kent (1985) Kent, S. M. 1985, ApJS, 59, 115, doi: 10.1086/191066
  * Kuchinski et al. (2001) Kuchinski, L. E., Madore, B. F., Freedman, W. L., & Trewhella, M. 2001, AJ, 122, 729, doi: 10.1086/321181
  * Law et al. (2007) Law, D. R., Steidel, C. C., Erb, D. K., et al. 2007, ApJ, 656, 1, doi: 10.1086/510357
  * Law et al. (2012) Law, D. R., Steidel, C. C., Shapley, A. E., et al. 2012, ApJ, 745, 85, doi: 10.1088/0004-637X/745/1/85
  * Lotz et al. (2008a) Lotz, J. M., Jonsson, P., Cox, T. J., & Primack, J. R. 2008a, MNRAS, 391, 1137, doi: 10.1111/j.1365-2966.2008.14004.x
  * Lotz et al. (2006) Lotz, J. M., Madau, P., Giavalisco, M., Primack, J., & Ferguson, H. C. 2006, ApJ, 636, 592, doi: 10.1086/497950
  * Lotz et al. (2004) Lotz, J. M., Primack, J., & Madau, P. 2004, AJ, 128, 163, doi: 10.1086/421849
  * Lotz et al. (2008b) Lotz, J. M., Davis, M., Faber, S. M., et al. 2008b, ApJ, 672, 177, doi: 10.1086/523659
  * Lowenthal et al. (1997) Lowenthal, J. D., Koo, D. C., Guzmán, R., et al. 1997, ApJ, 481, 673, doi: 10.1086/304092
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615
  * Mortlock et al. (2013) Mortlock, A., Conselice, C. J., Hartley, W. G., et al. 2013, MNRAS, 433, 1185, doi: 10.1093/mnras/stt793
  * Nersesian et al. (2023) Nersesian, A., Zibetti, S., D’Eugenio, F., & Baes, M. 2023, arXiv e-prints, arXiv:2303.01907, doi: 10.48550/arXiv.2303.01907
  * Papovich et al. (2003) Papovich, C., Giavalisco, M., Dickinson, M., Conselice, C. J., & Ferguson, H. C. 2003, ApJ, 598, 827, doi: 10.1086/378976
  * Perrin et al. (2014) Perrin, M. D., Sivaramakrishnan, A., Lajoie, C.-P., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9143, Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, ed. J. Oschmann, Jacobus M., M. Clampin, G. G. Fazio, & H. A. MacEwen, 91433X, doi: 10.1117/12.2056689
  * Rieke et al. (2015) Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584, doi: 10.1086/682252
  * Rieke et al. (2005) Rieke, M. J., Kelly, D., & Horner, S. 2005, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 5904, Cryogenic Optical Systems and Instruments XI, ed. J. B. Heaney & L. G. Burriesci, 1–8, doi: 10.1117/12.615554
  * Robitaille et al. (2020) Robitaille, T., Deil, C., & Ginsburg, A. 2020, reproject: Python-based astronomical image reprojection, Astrophysics Source Code Library, record ascl:2011.023.  http://ascl.net/2011.023
  * Rodriguez-Gomez et al. (2019) Rodriguez-Gomez, V., Snyder, G. F., Lotz, J. M., et al. 2019, MNRAS, 483, 4140, doi: 10.1093/mnras/sty3345
  * Rosa et al. (2018) Rosa, R. R., de Carvalho, R. R., Sautter, R. A., et al. 2018, MNRAS, 477, L101, doi: 10.1093/mnrasl/sly054
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41.  https://ui.adsabs.harvard.edu/abs/1963BAAA....6...41S
  * Snyder et al. (2015) Snyder, G. F., Lotz, J., Moody, C., et al. 2015, MNRAS, 451, 4290, doi: 10.1093/mnras/stv1231
  * Tang et al. (2020) Tang, Y., Chen, Q., Zhang, H.-X., et al. 2020, ApJ, 897, 79, doi: 10.3847/1538-4357/ab98fd
  * Thorp et al. (2021) Thorp, M. D., Bluck, A. F. L., Ellison, S. L., et al. 2021, MNRAS, 507, 886, doi: 10.1093/mnras/stab2201
  * Trayford et al. (2017) Trayford, J. W., Camps, P., Theuns, T., et al. 2017, MNRAS, 470, 771, doi: 10.1093/mnras/stx1051
  * Treu et al. (2023) Treu, T., Calabrò, A., Castellano, M., et al. 2023, ApJ, 942, L28, doi: 10.3847/2041-8213/ac9283
  * Valentino et al. (2023) Valentino, F., Brammer, G., Gould, K. M. L., et al. 2023, ApJ, 947, 20, doi: 10.3847/1538-4357/acbefa
  * Whitney et al. (2021) Whitney, A., Ferreira, L., Conselice, C. J., & Duncan, K. 2021, ApJ, 919, 139, doi: 10.3847/1538-4357/ac1422
  * Wuyts et al. (2012) Wuyts, S., Förster Schreiber, N. M., Genzel, R., et al. 2012, ApJ, 753, 114, doi: 10.1088/0004-637X/753/2/114
  * Yang et al. (2022) Yang, L., Morishita, T., Leethochawalit, N., et al. 2022, ApJ, 938, L17, doi: 10.3847/2041-8213/ac8803
  * Yao et al. (2022) Yao, Y., Chen, G., Liu, H., et al. 2022, A&A, 661, A112, doi: 10.1051/0004-6361/202243104
  * Yao et al. (2022) Yao, Y., Liu, H., Kong, X., et al. 2022, ApJ, 926, 57, doi: 10.3847/1538-4357/ac3b55
  * Zhan (2011) Zhan, H. 2011, Scientia Sinica Physica, Mechanica & Astronomica, 41, 1441, doi: 10.1360/132011-961