Homodyne Detection Quadrature Phase Shift Keying Continuous-Variable Quantum Key Distribution with High Excess Noise Tolerance

Discrete-modulated continuous-variable quantum key distribution with homodyne detection is widely recognized for its ease of implementation, efficiency with respect to error correction, and its compatibility with modern optical communication devices. However, recent studies report that the application of homodyne detection obtains poor tolerance to excess noise and insufficient transmission distance, hence seriously restricting the large-scale deployment of quantum secure communication networks. In this paper, we propose a homodyne detection protocol using the quadrature phase shift keying technique. By limiting information leakage, our proposed protocol enhances excess noise tolerance to a high level. Furthermore, we demonstrate that homodyne detection performs better than heterodyne detection in quaternary-modulated continuous-variable quantum key distribution under the untrusted detector noise scenario. The security is analyzed using the tight numerical method against collective attacks in the asymptotic regime. Our results imply that the current protocol is able to distribute keys in nearly intercity area and thus paves the way for constructing low-cost quantum secure communication networks.

I. INTRODUCTION
In recent years, the rapid development of quantum computing [1,2] has signaled the coming quantum era, which threatens modern secure communications. For example, the most widely used public key cryptography, the Rivest-Shamir-Adleman cryptosystem [3], whose security relies on the complexity of factorizing a large number, can be cracked by the Shor algorithm [4] with quantum computers. Fortunately, the one-time pad algorithm [5] can ensure information-theoretically secure communication provided that two remote users share a large number of identical secret keys. Quantum key distribution (QKD) [6,7] is the only solution that enables two distant users to share secure keys against the most general attacks [8][9][10].
Since the first protocol was proposed by Bennett and Brassard in 1984 [6], a large number of discrete-variable quantum key distribution (DV-QKD) protocols, based on various discrete degrees of freedom of a single photon, have been developed, such as polarization [11][12][13] and time bin [14][15][16]. Many attempts have been made to realize the network deployment of DV-QKD [9,10,17,18]. Compared with DV-QKD, continuous-variable quantum key distribution (CV-QKD) [8,19] has the competitive advantages of lower-cost implementations and higher secret key rates over the metropolitan areas, and it has thus attracted widespread attention.
CV-QKD encodes keys on the quadratures of the quantized electromagnetic field by preparing and measuring coherent states [20][21][22] or squeezed states [23]. From * hlyin@nju.edu.cn † zbchen@nju.edu.cn an experimental perspective, CV-QKD based on Gaussian modulation has made remarkable achievements [24][25][26][27][28][29][30][31][32][33][34][35], such as the provision of security against all detection loopholes [33], integration of all optical components on chips [34], and record-breaking distance up to 200 km [35]. Because Gaussian modulation, which samples quadratures from Gaussian distribution, possesses U (n) symmetry, it has a relatively complete security analysis [36][37][38][39][40] based on the optimality of Gaussian attacks [41,42]. However, the continuous modulation of quadratures is not experimentally feasible [43]. Practical implementations sample from a discrete and finite distribution of coherent states to approximate the Gaussian distribution, resulting in an unclosed gap between security analysis and experimental implementation. Moreover, Gaussian modulation requires complex post-processing procedures mainly because of the low signal-to-noise ratio and complex data structure, which often consumes large amounts of computing resources [27,44]. Unlike Gaussian modulation, discrete modulation [45][46][47][48][49] can sample from a quite small set of coherent states with different phases and, thus, can perform an easier error correction. The security proof of discrete modulation carefully considers the effect of discrete distribution of prepared states. Despite that discrete modulation lacks U (n) symmetry, several effective numerical [50][51][52] and analytical methods [53,54] have emerged for security analysis, and robustness to collective attacks in the asymptotic regime has been verified. Recently, the security proof of a binary-modulated protocol has been proposed in the finite-key-size regime against general coherent attacks [55]. With security improvement and simplified implementations, discrete modulation becomes a trend of CV-QKD [9,10].
In this paper, we present a quaternary-modulated protocol with homodyne detection. Unlike in previous quaternary-modulated protocols [50,51], we prepare each signal with the phase chosen from { π 4 , 3π 4 , 5π 4 , 7π 4 }, which can be implemented using classical quadrature phase shift keying (QPSK) technology. We require the sender to generate raw keys according to quadrature choices of the receiver. According to this interesting operation, our protocol shows good tolerance to excess noise and promises significantly longer transmission distances than the previous homodyne detection protocol [51]. In addition, the secret key rate of the proposed protocol is comparable to that of the well-performed heterodyne detection protocol [51]. Considering the imperfections of detectors, our protocol even performs better than the heterodyne detection protocol.
We introduce our protocol in Section II, with an instruction of a possible experimental setup. In Section III, the formula of the secret key rate is given against collective attacks in the asymptotic regime. In Section IV, we numerically simulate the performance of our proposed protocol, which demonstrates an increase in transmission distance and secret key rate. Further discussions are presented in Section V.

II. PROTOCOL DESCRIPTION
We illustrate our protocol in this section using the preparing and mapping method in Fig. 1  tions show that no secret keys can be generated, they abort the protocol. Otherwise, they proceed. (4) Raw key generation.-Alice and Bob can obtain their raw keys using the remaining undisclosed rounds. Assuming that M rounds remain, they number these rounds according to the order of sending. For the kth round, Alice labels the corresponding state |φ k ∈ αe i π 4 , αe i 3π as a k ∈ {00, 10, 11, 01}. Then, the raw key bit x k is equal to the first bit of the label a k when the round is measured inq, and it is equal to the second bit of a k when the round is measured inp. Thus, she obtains her string X ′ = (x 1 , .., x k , ..., x M ). For the kth round, Bob maps the outcome greater than ∆ into the raw key bit z k = 0 and the outcome smaller than −∆ into z k = 1. The outcome with other values is mapped into z k = ⊥. Thus, he obtains his string Z ′ = (z 1 , .., z k , ..., z M ). ∆ is a non-negative parameter and related to post-selection. A protocol without post-selection can set ∆ = 0. If ∆ > 0, Bob communicates with Alice about positions of bits with value ⊥ and they get rid of these positions from their strings X ′ and Z ′ . Finally, they obtain their raw key strings X and Z. and Bob choose suitable methods to conduct error correction and privacy amplification. In the end, they generate secret keys. A realistic experimental setup of our protocol is proposed in Fig. 2. The required components are displayed with unnecessary details omitted. To prevent security loopholes of the LO, we adopt the local LO scheme with secure phase compensation [30][31][32]. Alice randomly [32] sends reference pulses and signal pulses. Signal pulses are modulated using QPSK technology, and are used to generate keys. Setting the phase of signal pulses before QPSK as the reference phase, Alice prepares reference pulses that are brighter than signal pulses and contribute nothing to the generation of keys. Reference and signal pulses are transmitted to Bob's side through quantum channels. The polarization of any pulse received by Bob is corrected by the electrical polarization controller. Bob generates a local LO with the phase randomly chosen from {0, π 2 } using the phase modulator, which will be sent to the beam splitter together with signal pulses or reference pulses for interference. The photodetector with a subtractor reveals the interference outcomes, according to which Bob evaluates the exact phase difference between the reference phase and local LO. Then, Bob uses this phase information to conduct phase compensation [32]. Thereafter, Alice and Bob conduct postprocessing described by steps (3)-(5), and we do not show this part in the figure.

III. SECRET KEY RATE
To evaluate the secret key rate, we analyze the security of an equivalent entanglement-based protocol of our prepare-and-measure protocol. Then, we calculate the secret key rate against collective attacks in the asymptotic regime by applying the numerical method [51,52].

A. Entanglement-based protocol
An equivalent entanglement-based protocol of our prepare-and-measure protocol can be described below. Alice prepares an entangled state where the subscripts A and A ′ represent two entangled systems.
x represents one of four labels in {00, 10, 11, 01}, and p x = 0.25 refers to the probability that the entangled state will collapse into the state labeled by x. {|x } is the set of orthogonal bases in system A, which can be measured by a set of positive operator valued mea- , respectively. The system A is kept by Alice and A ′ is sent to Bob via a quantum channel. Thus, the final state is given by where E A ′ →B is a completely positive and tracepreserving mapping. This mapping includes the influence of the environment and attacks by Eve. Alice randomly projects the system A to an eigenstate by {M x }. Then, she can write down the corresponding label x without performing any additional operations. Bob measures the system B as with the step (2) in our prepare-and-measure protocol and follows the remaining steps.

B. Key rate formula
Applying reverse reconciliation means that Alice corrects her string X according to Bob's string Z, which implies that adversary Eve should derive Bob's string unscrupulously. In the case where Eve performs collective attacks in the asymptotic regime, the secret key rate formula [51,56] is given by where p pass is the sifting probability of preserving a round in the post-selection step to generate a raw key. The conditional von Neumann entropy H(Z|E) describes the uncertainty of the string Z from Eve's perspective. Eve's maximal knowledge of Bob's string Z leads to the minimum uncertainty of Z under certain density matrix ρ AB . To evaluate Eve's maximal knowledge when ρ AB is uncertain for Alice and Bob, we must find the minimum conditional entropy H(Z|E) of all ρ AB that conform to constraints S. δ EC is the amount of information leakage of each round in the error correction step. Considering the reconciliation efficiency, the leakage is where H(Z) is the total information of Z related to the probability distribution of Bob's measurement outcomes. I(X; Z) is the classical mutual information between two strings. H(Z|X) is the conditional entropy describing the uncertainty of the string Z conditioned on knowing the string X. The parameter β represents the efficiency of error correction. Because Bob announces the measurement quadrature of each round, we can analyze the security of each quadrature separately. To this end, we rewrite the secret key rate formula as The coefficient 1 2 indicates that one half of states are measured inq and one half inp. The secret key rate comprises the secret key rate inq and the secret key rate inp. Note that Eve manipulates the density matrix ρ AB before Bob measures it, and she thus has no knowledge of Bob's choice of quadrature at that time. We can move the summation sign into the minimization problem because ρ AB of different quadratures are shared.
According to Refs. [57,58], the conditional entropy terms in formula (5) with their coefficients can be reformulated as where D(ρ||σ) = Tr(ρ log 2 ρ) − Tr(ρ log 2 σ) is the quantum relative entropy, G y describes the post-processing of different quadratures, and Z is a pinching quantum channel that reads out key information. The post-processing mapping G y (ρ) = K y ρK † y corresponds to the measurementŷ ∈ {q,p} and is given by Here |b R is the state of the key register R, which is determined by interval operators {I b y }. I b y projects the system B to one of two subspaces spanned by the eigenstates of operatorŷ in terms of the mapping rule: Here G y is in a simplified form introduced by Ref. [51]. Kraus operators {K y } describe the part of the measurement and announcement: Alice's operation is described by an identity matrix because Alice announces nothing and has no effect on the register R. Her measurement in the entanglement-based protocol can be moved after the announcement and never influences K y . Bob would announce the choice of quadrature, and his measurement is represented by interval operators. When considering the security, we can discuss two quadratures separately, denoted by the parameter y ∈ {q, p} in the above formulas.
Finally, the key rate formula reads

C. Numerical method
To realize the calculation of secret key rates, a photonnumber cutoff assumption [51] is adopted. Based on this assumption, we can describe operators and density matrices in the photon-number representation with finite dimensions N c . Because δ EC and p pass are only concerned with the information held by Alice and Bob, the term y∈{q,p} p y pass δ y EC in the key rate formula (9) can be easily calculated according to general definitions [51]. We calculate the rest of the formula numerically. The rest is a minimization problem that can be described by The variable is the density matrix ρ AB subject to constraints S [51]. The first four constraints come from the experimental outcomes, where x belongs to {00, 10, 11, 01} and q x , p x , n x , and d x are expectation values of operators when Bob measures states labeled by x. Homodyne detection directly outputs the outcomes of operatorsq andp, whilen = 1 2 (q 2 +p 2 − 1) andd =q 2 −p 2 correspond to the second moments ofq andp. The next constraint about partial trace of system B comes from the requirement of the completely positive and trace-preserving mapping, which implies that the quantum channel cannot influence the system A of Alice. The last two constraints are natural requirements because ρ AB is a density matrix.
This minimization problem can be solved using many methods. We adopt one using a linearization method, given in Appendix A. Briefly speaking, the numerical method involves two steps. First, we approach the optimal value of the minimization problem (10) by at most N i iterations. Second, the dual problem of the minimization problem is considered to guarantee that the result is less than or equal to the optimal value.

IV. PERFORMANCE
Here, we investigate the performance of our protocol by performing numerical simulations. The simulation method is described in Appendix B, where we set the reconciliation efficiency β = 0.95.
To simulate the performance, we adopt the phaseinvariant Gaussian channel model in the absence of Eve [51], which adds an effective excess noise to the input of the channel and transmits the input under the loss-only scenario. In this channel, the transmittance is η = 10 − aL 10 for distance L with a = 0.2 dB/km, and the effective excess noise is ξ = (∆q obs ) 2 (∆qvac) 2  theq quadrature. Both variances are in the natural unit. Because ξ is normalized by (∆q vac ) 2 , the excess noise is under the shot-noise unit without vacuum noise.
We first emphasize the significant performance improvement of our proposed protocol by comparing it with Protocol 1 in Ref. [51], which also prepares four types of coherent states and performs homodyne detection without post-selection.
As shown in Fig. 3, we investigate the excess noise tolerance at different distances. The excess noise tolerance is calculated by maximizing the excess noise ξ that enables the protocol to generate keys. The cutoff photon number N c is truncated at 12 with N i = 100 iterations in the first step of the numerical method. The excess noise tolerance of our protocol is approximately double that of Protocol 1.
As shown Fig. 4, we investigate the best secret key rates of our protocol for different levels of excess noise by optimizing the amplitude α of signal states in the interval [0.6, 1.1] with steps of 0.01. We also calculate the best secret key rates of Protocol 1 by optimizing the amplitude α of signal states in the interval [0.35, 0.6] with steps of 0.01. Different search ranges are used because the optimal range of α varies for different protocols. The cutoff photon number N c is 12 and the maximal iteration number N i of the first step in numerical method is 300. The solid lines show the performance of our proposed protocol, whereas the dashed lines show the performance of Protocol 1. When the excess noise is experimentally feasible, such as ξ = 0.02, our protocol enables the distri- bution of keys for around 200 km with meaningful 10 −6 bit secret keys per pulse. Moreover, when Protocol 1 can hardly generate keys for ξ = 0.04, our protocol can still distribute secret keys over 50 km. In Fig. 5, we show that there is an improvement in the transmission distance for a small excess noise with N c = 12 and N i = 300. When the excess noise is sufficiently small, but still possible, our protocol has the potential to distribute secret keys over 350 km. The transmission distance and secret key rates are increased significantly, which makes our protocol applicable and useful in the quantum secure communication network.
We further compare our protocol with a wellperformed quaternary-modulated heterodyne detection protocol called Protocol 2 in Ref. [51]. Owing to differences in the detection methods, the comparison is made under two scenarios: ideal detector scenario and untrusted detector noise scenario.
Under the ideal detector scenario, both excess noise and loss are from the channel. We set the detection (electronic) noise of detectors ξ det = 0 and the detection efficiency η det = 1. In Fig. 6, our protocol and Protocol 2 are simulated with parameters ξ = 0.04, N c = 12 and N i = 300. We search the best secret key rates of our protocol and Protocol 2 by optimizing the amplitude α of signal states in the interval [0.6, 1.1] with steps of 0.01. Protocol 2 can generate a higher secret key rate compared with our proposed protocol under the ideal detector scenario because heterodyne detection can accumulate twice the amount of raw key data than homo- The difference between the secret key rates of this work and Protocol 2 is normalized by the secret key rate of this work R het /Rour − 1 and is shown with the y axis on the right. dyne detection. We call the secret key rate of Protocol 2 R het and the secret key rate of our protocol R our . Then, we use R het /R our − 1 to represent the proportion of the secret key rate difference R het − R our in R our , as shown by the dashed line in Fig. 6. The secret key rate of Protocol 2 is no more than 20% higher than that of our proposed protocol within 30 km. The uncertainty relation limits the accuracy of measuring both quadratures simultaneously and increases the bit error rate of Protocol 2. Therefore, the improvement of the secret key rate by heterodyne detection is not significant. We also show the secret key rate of the Gaussian-modulated homodyne detection protocol in Fig. 6. The modulated variance in the Gaussian-modulated protocol is the parameter that has been optimized. In Fig. 7, we regard the detector as a noisy and lossy device, and Eve can control the imperfection of the device, which is a more practical condition. We should note that the heterodyne detector often comprises two homodyne detectors and one 50:50 beam splitter that splits signal pulses. Therefore, two points of imperfection are considered under the untrusted detector noise scenario that we used here. First, the detection noise of a homodyne detector is ξ hom = 0.002, and the detection noise of  [51]. The excess noise of our protocol is ξ = ξ oth + ξ hom η , and the excess noise of the heterodyne protocol is ξ = ξ oth + 2ξ hom a heterodyne detector is at least ξ het = 2ξ hom = 0.004. Second, the imperfection in the beam splitter usually causes an additional insertion loss of about 0.7 dB. By setting other excess noise as ξ oth = 0.01, the total excess noise of our protocol is ξ = ξ oth + ξ hom η , and the total excess noise of Protocol 2 with heterodyne detection is ξ = ξ oth + 2ξ hom η×10 − 0.7 10 . The cutoff photon number N c is 12 with N i = 100 iterations in the first step of the numerical method. Under this untrusted detector noise scenario, our protocol can perform better than the heterodyne detection protocol.
In Fig. 8, we compare the heterodyne protocol [48,51] under ideal, trusted, and untrusted detector noise scenarios with our protocol under ideal and untrusted detector noise scenarios. The data of the heterodyne protocol were reported by the authors of Ref. [48]. The imperfect detector has excess noise ξ hom = 0.01 and detection efficiency η d = 0.719. The other excess noise is ξ oth = 0.01 and the cutoff photon number N c is 15. We simulate our proposed protocol using the same parameters. Thus, the total excess noise of our protocol under the untrusted detector noise scenario is ξ = ξ oth + ξ hom η×η d . The iteration in the first step of the numerical method is N i = 100. Our protocol is a little worse than the heterodyne protocol under the ideal scenario, but it is better under the untrusted detector noise scenario. The heterodyne protocol   [48] is the blue line in the middle. The data of the heterodyne protocol has been reported by the authors of Ref. [48]. The parameters are the same as those used in Fig. 3 in Ref. [48]. Specifically, the imperfect detector has excess noise ξ hom = 0.01 and detection efficiency η d = 0.719. ξ oth = 0.01 is the excess noise and does not include detection noise. The total excess noise of our protocol under the untrusted detector noise scenario is ξ = ξ oth + ξ hom with trusted detector noise has a key rate that is comparable to the same protocol with the ideal detector, which may imply good performance when we also consider our protocol under the trusted detector noise scenario.
For the experiments, we show the optimal amplitudes that are searched for different transmission distances with parameters N c = 12 and N i = 300. As shown in Fig. 9, in the case of a long distance, the optimal amplitude decreases as the transmission distance increases, except for the case when there are some jitters, and a larger excess noise induces a lower optimal amplitude with little violation. The optimal amplitude for remote users is about 0.66. This conclusion also applies to the condition where N c = 12 and N i = 100.
The post-selection is a useful method to further increase the transmission distance and improve the secret key rate. In Fig. 10, we optimize the sifting parameter ∆ using the amplitude of signal states that have been optimized without post-selection. The cutoff photon number N c is truncated at 12 with N i = 100 iterations. The available ∆ has steps of 0.01. At distances where secret key rates are zero without post-selection, we search the best post-selection parameter ∆ by setting the amplitude to 0.66. This choice is reasonable considering that the optimal amplitude for the long distance is around 0.66 in average, as depicted in Fig. 9. The post-selection is more important for the scenario with the higher excess noise. For ξ = 0.04, the secret key rates are improved, and the largest transmission distance is extended. For smaller excess noise ξ = 0.03, the improvement of postselection is subtle. Over a long distance, there is no improvement in the secret key rates for ξ = 0.03. This is because as the excess noise increases, the overlap of outcome distributions obtained by measuring different kinds of signal states also increases, which causes higher bit error rates. The post-selection discards outcomes around zero because the overlap of outcome distributions most likely happens around zero. By post-selection, we discard many error bits, thus reducing the bit error rates. Protocols with heterodyne detection enjoy the same advantages of post-selection [59]. The optimal ∆ is larger under longer transmission distance. The sifting probability related to post-selection decreases to around 0.15 over the longest transmission distance for ξ = 0.04; the post-selection thus makes sense in our protocol. This operation can work well when the amplitude cannot be prepared optimally in the experiment.
Note that the curves in the above figures should have been smooth in theory. Although the secret key rate formula is a convex function with convex domain ρ AB , it is difficult to unify the imprecision of ρ AB under different amplitudes and transmittances. The inaccuracy also accumulates along with steps of calculations, resulting in some calculation noise. Finally, as a numerical method, limited by the calculation accuracy of computers and the stability of the convex optimization algorithm, unstable points and unsmooth curves appear. The quality of the algorithm largely determines the secret key rates. When the ideal secret key rate is quite small that it reaches the computational accuracy limit, numerical results sometimes show that no keys can be generated. Achieving improvements in the algorithm may help further enhance the transmission distance and stabilize the key rates. The calculation speed of the numerical method is often slow and is approximately 20 min per set of parameters. The speed largely depends on the number of iterations and the truncated photon number. Intriguingly, recent work [60] can utilize a neural network model to quickly predict the secure key rate of homodyne detection discrete-modulated CV-QKD with high accuracy.

V. DISCUSSION
Homodyne detection is the first adopted measurement method of CV-QKD with noticeable performance experimentally, and it provides excellent tolerance in terms of excess noise in Gaussian-modulated protocols. Our main goal is to inherit this mature technology and apply it to the discrete-modulated protocol. We propose a quaternary-modulated homodyne detection protocol with very high excess noise tolerance and long transmission distance that has never been reached by previous discrete-modulated homodyne detection protocols [51,61,62]. There is a quaternary-modulated heterodyne detection protocol [51] that shows high secret key rate and long transmission distance under large excess noise based on simulations. However, the heterodyne detector [63][64][65][66][67] is more sophisticated than the homodyne detector, and it introduces more noise and loss in experiments [68][69][70][71]. We show that our proposed protocol can offer a comparable performance with this heterodyne detection protocol under the ideal detector scenario. Our protocol can even perform better than this heterodyne detection protocol when considering the untrusted detector noise condition.
We determine why the previous quaternary-modulated homodyne detection protocol [51] cannot transmit over a long distance although they also send four kinds of states. In particular, from a practical point of view, our protocol and the previous protocol both generate one bit per pulse by homodyne detection. Compared with Eq. (17) in Ref. [51], the mapping K y of our protocol in Eq. (7) is an identity matrix at Alice's end, which means that Alice remains silent in our protocol. This mathematical difference implies that additional announcements cause more information leakage.
In addition, the adoption of QPSK is required during the preparation. Note that QPSK is a mature method in classical communication, shifting the phase of states by { π 4 , 3π 4 , 5π 4 , 7π 4 }. In previous quaternary-modulated protocols [50,51], they use states with phases {0, π 2 , π, 3π 2 }. The probability distribution of a state in quadraturep is symmetric about the original point if the state has phase 0 or π. Thus, states with phase 0 or π cannot represent keys when Bob appliesp measurement. Similarly, states with π 2 or 3π 2 cannot represent keys when Bob applieŝ q. The additional π 4 phase shift in our protocol offers chances to generate keys in both quadratures, and Alice can thus announce nothing.
Eventually, owing to the tolerance of excess noise, the transmission distance is improved, which removes the obstruction of the proceedings of the large-scale discretemodulated CV-QKD network. It is worth recalling that commercial instruments for QPSK have been updated to over 50 Gbps [72], which implies that the speed of preparation in our proposed protocol can also be increased to 50 Gbps. Our protocol will maximize the use of the mature homodyne detection technology to achieve the longdistance key distribution. The experimental demonstration of our protocol is of great practical value for network security and information security. Moreover, we can consider the trusted detector noise model [48]. With a simpler preparation, our proposed protocol has a chance to exceed the transmission record of 202.81 km set by the Gaussian modulation protocol [35] with an ultralow-loss fiber and small excess noise parameters. The experimental schematic and the search of optimal amplitudes in the text provide instructions for the experimental implementation.
Further studies of our protocol's security are also desirable. Note that the photon-number cutoff assumption can be removed [73] in the heterodyne detection protocol, which inspires us to present a more rigorous security analysis of our protocol in future work. Moreover, we can focus on determining a suitable adjustment of security theory to deal with the finite-size effects. The security against coherent attacks is another interesting and important direction for discrete-modulated CV-QKD. The quantum de Finetti representation theorem [74] is a possible path to equalize the information leaked by coherent attacks and collective attacks.

ACKNOWLEDGMENTS
We thank Jie Lin for valuable discussions on the numerical method in Ref. [51] and providing data points from Ref. [48] for the comparison in Fig. 8 The minimization problem in Eq. (10) is a convex optimal problem with one unique solution. It can be described as follows.
The independent variable is the density matrix ρ AB . The dimension number of Alice's system is 4, which is determined by the number of different states she prepared. |x x| A with x ∈ {00, 10, 11, 01} projects Alice's system into a one-dimensional subspace, which implies that Alice sends state |φ x ∈ αe i π 4 , αe i 3π to Bob. However, the state that Bob receives is infinitedimensional. The photon-number cutoff assumption is proposed to cut off the dimension of Bob's system. According to the assumption, ρ AB can be represented by a finite-dimensional matrix in Fock representation. Thus, the numerical method can be applied. The operatorŝ q,p,n,d correspond to the measurement on Bob's system. Considering the photon-number representation, they can be defined by the annihilation operatorâ: To calculate them numerically, the matrix ofâ in photonnumber representation is also cut off according to the photon-number cutoff assumption. We introduce notation q x , p x , n x , d x to represent the expectation values of these operators. The relative entropy D(ρ||σ) is the dependent variable that we need to minimize. G(ρ) = KρK † is the map in post-processing with Kraus operator K.
To implement the convex optimization, a two-step numerical method is introduced [58] below.
For convenience, we use f (ρ) to represent the function D(G(ρ)||Z[G(ρ)]). The constraints can be described by a set S: where H + is the set of positive semidefinite operators. Γ i is the Hermitian operator and γ i is the corresponding expectation value. To fit this form, one should rewrite the fifth constraint in Eq. (10). Specifically, one needs to transform the partial trace on Bob into some Hermitian operators.
Using the linearization method, we can obtain where ρ * is the optimal density matrix (i.e., the solution of minimization problem Eq. (10) and ǫ represents a perturbation parameter, by which the gradient of f (ρ) can exist. Owing to the existence of ǫ, the term is applied to correct the difference between f (ρ) and f ǫ (ρ) caused by perturbation. In our simulation, we set ǫ = 10 −12 to decrease the value of ζ ǫ . Using ǫ, the dependent variable becomes The gradient of f ǫ (ρ) is d ′ is the dimension of the matrix G(ρ), which is twice the dimension of density matrix ρ.
Using the linearization expressed in Eq. (A3), the optimization problem is transformed into a semi-definite program, such as the second or third term of the first or the second line, respectively, in Eq. (A3). Therefore, we only need to find ∆ρ or σ = ρ + ∆ρ that minimizes the related term. We take ∆ρ as an example, and the first step is to find a density matrix ρ that is quite close to the optimal one. To accomplish this step, we perform the algorithm as follows: 1. Begin with any ρ 0 ∈ S and set i = 0. 2. Compute min ρ+∆ρ∈S [Tr((∆ρ) T ∇f ǫ (ρ i ))]. 3. If Tr((∆ρ) T ∇f ǫ (ρ i )) is sufficiently small, STOP. 4. Find λ ∈ (0, 1) that minimizes f ǫ (ρ i + λ∆ρ). 5. Set ρ i+1 = ρ i + λ∆ρ, i ← i + 1 and go to 2.
When solving the minimization problem in step 2 of the algorithm, we find that the results are often sensitive to the setting of ρ 0 . Different convergence paths affect the convergence effect. In our simulation, we often choose the result of feasible verification under S or the result that maximizes the minimal eigenvalue of ρ under S as the primal density matrix ρ 0 . The iteration of 2-5 sometimes needs to be broken manually; otherwise, it may take a large amount of time to converge or even never stop. The largest iteration number in our simulation is 300 or 100. After 300 or 100 iterations, we choose ρ i with the smallest Tr((∆ρ) T ∇f ǫ (ρ i )) among all iterations as our result of the first step.
Reference [58] does not offer any method to deal with the minimization of f ǫ (ρ i + λ∆ρ) because its DV-QKD example is simple for calculation. However, it is not an easy task for CV-QKD. We provide our own method inspired by BinarySearch: 1. Set a starting point λ s = 0, an ending point λ e = 1, and a middle point λ m = 0.5.
2. Calculate f ǫ (ρ i + λ∆ρ) in three points, and compare them. 3. If λ s is the smallest one, update the ending point by λ e = λ m and the middle point by λ m = 0.5(λ s + λ e ). Then, go to 2. 4. If λ e is the smallest one, update the starting point by λ s = λ m and the middle point by λ m = 0.5(λ s + λ e ). Then, go to 2. 5. If λ m is the smallest one, two temporary points are generated by λ 1 = 0.5(λ s + λ m ) and λ 2 = 0.5(λ m + λ e ). 6. Calculate f ǫ (ρ i + λ∆ρ) of two temporary points and compare them with λ m . 7. If λ m is still the smallest one, set λ s = λ 1 and λ e = λ 2 , and go to 2. 8. If λ 1 is the smallest one, set λ e = λ m and λ m = λ 1 , and go to 2. 9. If λ 2 is the smallest one, set λ s = λ m and λ m = λ 2 , and go to 2.
The calculation results of some points can be temporarily stored for the next iteration. The termination of the iteration is decided by the calculation accuracy of computers.
After the first step, we seemingly find an appropriate density matrix ρ. However, considering the imprecision of numerical calculations, we cannot guarantee that the solution of the first step is the optimal one. In fact, we find an upper bound of the solution. To protect the security of the generated keys, we use the dual problem in the second step.
One can transform the primal problem min σ∈S [Tr(σ T ∇f ǫ (ρ))] into the dual problem max y∈S * (ρ) γ · y. The variable y obeys the constraints: In general, primal and dual problems satisfy the weak duality. The weak duality shows that the optimal value of the dual problem is always less than or equal to the optimal value of the primal problem. According to this principle, the final solution that we obtain from the dual problem must be a lower bound of the primal problem even if it is not optimal. The lower bound of the primal problem is also the lower bound of the secret key rates according to Eq. (A3).
In addition, we can consider the imprecision of the floating-point representations of {Γ i } and {γ i }.
where {Γ i } and {γ i } are the approximate representations used in the calculation. Furthermore, the solution ρ may not satisfy the constraints, especially the positive semidefinite constraint. Ref. [58] encourages one to transform ρ into positive semi-definite matrix ρ ′ by subtracting a term λ min I when the smallest eigenvalue λ min of ρ is negative. Then, we restrict the imprecision by After the first step, we extract the quantity ǫ ′ = max(ǫ rep , ǫ sol ) from the output ρ. This quantity influences the second step. The dual problem solved in the second step is adjusted to max ( y, z)∈S * (ρ) (A11) The description above is about one quadrature. Because f q (ρ) forq and f p (ρ) forp are added linearly, we can replace f ǫ (ρ) by f q,ǫ (ρ) + f p,ǫ (ρ). Then, ∇f ǫ (ρ) in primal and dual problems is replaced by ∇f q,ǫ (ρ) + ∇f p,ǫ (ρ). In addition, the perturbation compensation ζ ǫ is multiplied by 2.
Thus far, the lower bound of the minimization problem is solved numerically. Note that we solve primal and dual problems using the SDPT3 solver in CVX 1.22 package on MATLAB R2020b. This package is for specifying and solving convex programs, and SDPT3 is a free solver [75,76]. The code of our protocol is uploaded to the opensource code website [77].

Appendix B: Simulation Methods
We simulate the experiment's statistical results by assuming a phase-invariant Gaussian channel, as referred in Ref. [51]. When Alice sends a coherent state |α , Bob obtains a displaced thermal state centered at √ ηα with the variance 1 2 + ηξ 2 for each quadrature. The first term 1 2 in the variance is the vacuum noise. Thus, ρ x B can be given by with α x ∈ αe i π 4 , αe i 3π 4 , αe i 5π 4 , αe i 7π 4 and V A = ηξ 2 . ξ is the excess noise and η = 10 − aL 10 is the transmittance with a = 0.2 dB/km in the distance L. Thus, the expectation values of measurement operators are q x = 2ηRe(α x ), To simulate the error correction, the conditional probability between Alice's sending state and Bob's mapping result can be given by where y ∈ {q, p}. Then, we can find the sifting probability p y pass = x p x (P y (0|x) + P y (1|x)), and normalize the probability distribution after discarding bits ⊥. The conditional entropy for each quadrature is thus calculated by H(Z y |X y ) = x p x h P y (0|x) P y (0|x) + P y (1|x) ,


FPGA implementation of LDPC soft-decision decoders based DCSK for spread spectrum applications

Spread spectrum (SS) communications have attracted interest because of their channel attenuation immunity and low intercept potential. Apart from some extra features such as basic transceiver structures, chaotic communication would be the analog alternative to digital SS systems. Differential chaos shift keying (DCSK) systems, non-periodic and random characteristics among chaos carriers as well as their interaction with soft data are designed based on low-density parity-check (LDPC) codes in this brief. Because of simple structure, and glorious ability to correct errors. Using the Xilinx kintex7 FPGA development kit, we investigate the hardware performance and resource requirement tendencies of the DCSK communication system based on LDPC decoding algorithms (Prob. Domain, Log Domain and Min-Sum) over AWGN channel. The results indicate that the proposed system model has substantial improvements in the performance of the bit error rate (BER) and the real-time process. The Min-Sum decoder has relatively fewer FPGA resources than the other decoders. The implemented system will achieve 10-4 BER efficiency with 5 dB associate E b /N o as a coding gain.

INTRODUCTION
Being wide-band, chaotic signals are well suited for communication with the spread-spectrum [1]. Among the proposed digital schemes, the most comprehensive research was on chaos shift keying (CSK) and differential CSK (DCSK) [2]. The former was originally introduced for such a coherent implementation [3], requiring synchronized replicas of the chaotic base functions at the receiver. This requirement has still not been proved practical, though. On the other hand, the DCSK scheme represents a rather more robust noncoherent scheme [4] in which the receiver does not need the exact information of chaotic basis functions. Rather the two signal samples were correlated as well as the result of the correlation is compared with a threshold. The performance of digital communication systems based on chaos under an additive white gaussian noise (AWGN) environment has been studied thoroughly [5]- [8].
Low-density parity-check (LDPC) is one of the most effective techniques amongst the error correction codes. LDPC codes have gained a lot of attention recently since they can achieve exceptional performance close to the Shannon limit over the binary symmetric channel (BSC) and also the additive white gaussian noise (AWGN) channel [9]. Decoding an LDPC code allows for a high degree of parallelism, making it ideally suited for high data rate applications including wide-band wireless multimedia communications and magnetic storage systems. The parity check matrix's low-density nature thus contributes to both great distance properties and the relatively low complexity of its decoding algorithm [10], [11]. Moreover, the excessive noise derives from a wide scale of mobile communication, the bit error rate (BER) is acceptable in modern communication with an even more high-speed data rate. Hence, support for the DCSK communication system through LDPC codes becomes essential to mitigate the high error [12]. The field programmable gate array (FPGA) is used to evaluate a system that included each DCSK and LDPC codes in a real-time environment. Gallager 's suggested basic decoding algorithm for LDPC codes in 1962 was soft decisions such as the sum-product algorithm (SPA) and hard decision such as the bit flipping (BF) [13]. SPA requires more steps to addition and multiplication which enables implementation in the case of real-time implementation with FPGA due to the simplicity of this operation. A crucial trade-off between 'complexity' and 'efficiency' is needed in iterative decoding. In this paper, we propose a model of a communication system including a DCSK as a modulation technique supported by three SPA (Prob. Domain, Log Domain and Min-Sum which is designed to reduce hardware complexity) to improve the system performance. The proposed system will be implemented using an FPGA Kintex 7 development kit integrated with Vivado 2017.4 software.

DCSK WITH LDPC ENCODED COMMUNICATION SYSTEM
The proposed system diagram is shown in Figure 1 that included a binary data source that will be encoded through the LDPC block with a code rate of 0.5. The codeword is sent to a DCSK modulator that uses a chaotic carrier for spread the digital signal across a wide frequency band to achieve a modulated signal with such a spreading factor value equal to 16. AWGN noise is indicated with the modulated signal which is mostly used for the experimental applications for simplicity. Non-coherent demodulation will then be adapted on the receiver side to recover the received code-word, that will be decoded to regenerate the original information.

HARDWARE COMMUNICATION SYSTEM DESIGN
As mentioned in section 2 the system consists of three main parts: transmitter, channel and receiver. A transmitter section is established by the subsequent, Logistic-map as input chaotic data generator, LDPC encoder by using the systematic form of H matrix and computing the parity check equation for each row of this matrix and DCSK modulator. The receiver is established by the subsequent, AWGN noise, DCSK demodulator and LDPC decoder with three types of decoder algorithms such as; Prob. Domain, Log Domain and Min-Sum algorithms. The system model will be implemented using Xilinx SG as shown in Figure 2, each block is designed with specific parameters to match the overall system implementation as will be explained in detail individually. It is worth to mention that all the three decoder algorithms have the same design architecture except the part deal with decoder block.

LDPC encoder
The Bernoulli random binary number generated message word with a length of k=10 was fed to the LDPC block through gateway block to obtain a code word with a length of n=20 at the output of a certain block. There is a serial to parallel block within the LDPC encoder block shown in Figure 3 which will be used to convert a group of samples serially presented at the input to single samples present at the output is shown in Figure 4. In this paper, we have = 10 as mentioned previously represented by 1 − 10 in (1) while 11 − 20 represented the redundancy bits ( = − ) [14].
The parity bits are incorporated using (2) and (3) which reflect the first and second parity bits corresponds respectively to first and second rows of the (4) matrix. The remaining parity is generated in a  [15]. Figure 5 illustrates the hardware implementation of these equations that are implemented using 10 XOR gates to generate 10 parity check bits, based on the input message. 11 = 3 ⊕ 6 ⊕ 9 (2) 12 = 5 ⊕ 6 ⊕ 7 ⊕ 8 Now, the 20-bit code-word is supplied to concat block, concatenating two or more bits of inputs to get a symbol of these bits in the output. Eventually, the parallel to serial block induces the code word as output block. So that each sample provided in this block's input became multiple samples displayed serially in the output.

DCSK modulator
The code-word generated by LDPC was fed to a Xilinx SG DCSK modulator designed. Throughout the DCSK modulator, each bit Si is expressed by 2 sets of β-period chaotic signal samples where 2β refers to the spreading factor with such an integer number. The first and the second sets respectively represented the segment of reference and the segment of data. Samples from the data segment are equal to the reference segment for transmission +1, and an inverted variant of a reference segment is being used for transmission of bit -1. The transmitter ek 's output during the i th bit period is [16]- [18]. Figure 6 Shows the Xilinx SG block of DCSK modulation system with equal to 8.
Both Chaos generator and mapping blocks are shown respectively in Figures 7 and 8. Chebyshev polynomial function (CPF) of order two as in (6) is selected as a chaotic generator with such an initial condition value of 0.1 to be implemented simply using a multiplexer, multiplier and adder blocks:

Channel
If the modulated signal is transmitted over a channel, noise in the channel represented by additive white gaussian noise (AWGN), , can corrupt it. Such noise characterized by a wide frequency spectrum, which is statistically random radio noise. The hardware-implemented for these channels is illustrated with seed value 512 in Figure 9.

DCSK demodulation
To convert the 16 serial to parallel samples, the obtained signal passes to an S2P block to despread within a single bit. Figure 10 illustrates the details of the Xilinx SG S2P block consisting of 16 latches and its corresponding delay blocks. S2P 's 16 output samples are linked to a correlator, in which the reference samples as well as the corresponding information samples were correlated using multipliers. i th correlator output is the variable i th decision, Di [19]. The i th bit is demodulated using a zero threshold, measuring the sign of the final correlator output. The DCSK demodulator Xilinx SG is shown in Figure 11.

Decoders
Prob. Domain, Log Domain and Min-Sum algorithms represent a type of soft decision algorithms and it works by passing messages between the CNs and the variable nodes (VNs). It is used to correct the received bits to obtain 10-bit symbols which will represent the original signal after the decoding process. The demodulated bits are firstly fed to the S2P converter block to match each other within SG blocks. The S2P will be used to convert the 20-stream of bits to parallel bits to be initialized. In general, these algorithms perform in main three steps. − Initialization − Horizontal Step − Vertical Step and Decoding/Estimation The step details for each type of these algorithms will be described with a zoom-in view below.  (1) and P(0)= , (0) is obtained from (9) & (10) which represents the message that will be sent from the VN n to CN m. Each one of the messages has the pair , (0) as well as , (1) referring to the indication that rx has been one or zero [20].
Where (1)and (0) are the posterior probabilities and their values can be found based on the received signal, also N0 representing noise variance. Horizontal step: The horizontal step calculation is created depending on the number of 1s over all columns of ten rows throughout the matrix (11)-(13). , These variable node has been designed by using the basic add-sub and mult Xilinx SG blocks to evaluate the by taking the product of subtraction , ′ (1) also , ′ (0) excluding bit ′, , (0) as well as , (1) which representing the probability which check m can be designed easily by using add-sub to add or sub one from value then mux the results by 0.5 using Xilinx SG block. Figures 12 and 13 illustrate the first row computations according to the (10)-(12) of the horizontal step.  (1) Where ki n (0) and ki n (1) are scales factors which are implemented using Mult Xilinx SG blocks. Also, a decision for each bit is made according to (15).
The Xilinx SG block of column twelve is described in Figure 14.

Log domain algorithm
Initialization: In Log Domain decoder the log-likelihood ratio (LLR) of prior (which represents the receiving messages from the channel) and posterior (which represents the medial messages moved between CNs and VNs) probability is used. So, the prior messages sent from BN n to the CN m represent the LLR [21].
This process is applied to 20 received bits serially. Then, the absolute value and sign value (will calculate the sign value whether it is positive or negative or equal to zero) will be calculated for each bit. The implementation of the initialization block is done by using Mult, Mux and relational Xilinx SG block. The process of initialization block and the sign block will be declared in Figures 15 and 16 respectively. Horizontal step: The computations of the horizontal step are created depend upon the amount of 1s per column of ten rows throughout the matrix. To Computing the extrinsic messages for each set of bits connected to CN m by excluding the bit ′ the equation below will be used.
Where , represents the probability that parity-check m has been achieved in the case when bit n is supposed to be a 1 for the LLR. Compute the summation of , excluding the bit n.
Then compute , : Get products of ℎ , excluding the bit n: Finally, compute To reduce the complexity of Xilinx SG design related to the logarithm (log) process in (19), a lookup table is proposed by taking the values of rx and applying (19) as shown in Table 1. Figure 17 shows the details of the Xilinx SG design for the Pi calculation according to (19) for each bit equal to one.

4803
For example, the Xilinx SG design for the first row will be implemented by using Mult and Add-Sub block according to (20)- (23) as illustrated in Figure 18. The details for each one of the blocks in Figure 18 can be declared in Figures 19-21. Furthermore, to reduce the complexity of Xilinx SG design related to the logarithm (log) process according to (21), a look-up table is proposed as illustrated in Figure 20 and Table 2.
Where represent the collective log-likelihood ratio for n' th digit. Also, a decision is made according to (25) for each bit. The design model is implemented using Xilinx SG blocks (Add-Sub, Mux, and Relational Xilinx blocks) to manage the value of the variable node. The details for the Xilinx SG block of the twelve columns are described in Figure 22.

Min-sum algorithm
The receiver side of the Min-Sum algorithm is very similar to the Log Domain algorithm except for the initialization process and horizontal step. It's implemented by replacing the implementation of (19)-(23) that belongs to the Log Domain initialization and horizontal step block by (26)-(29) respectively [22][23][24].
This process is implemented using Xilinx SG blocks (Multi, Mux, and Relational Xilinx blocks) to manage the value of the variable node as illustrated. The details are illustrated in Figures 23-28.

Down sample
Finally, the fully flexible LDPC decoder design output bits were connected to down-sample block from Xilinx SG blocks with a sampling rate equal to 320 which is used to control the rate of the output signal at the receiver. These bits from the output of the down-sample blocks will be converted from parallel bit to stream bits after passing through concat and P2S blocks respectively to get the information data as seen in Figure 29. Once the Simulink models are developed for the proposed system, the VHDL code can be generated automatically using the SG block.

SIMULATION RESULTS
The performances of codes obtained from the proposed system were measured by BER as a function of SNR. To complete the simulation with an optimum parameter was the No. of one in each column (w c )=3, No. of iteration (iter)=6, frame length (F)=100, code rate=1/2 and spreading factor (2 ) value used are equal to 16. The results illustrated in Figure 30 are clear-out the comparisons between these communication systems with and without LDPC code.
From previous results, it is clear that the performance of the hard decision BF decoders will be worse than the soft decision decoder [25]. Also, the results show that close performances are illustrated when using Prob. Domain, Log Domain and Min-Sum, but all give approximately the same amount of gain. The system needs more than 13 dB SNR to approach 10 -3 BER without code. However, using soft-decision decoder algorithms (Prob. Domain, Log Domain and Min-Sum) will improve the performance of the communication system at low and high SNR. All results are summarized in Table 3 for the spreading factor 16 at BER of 10 -3 .

FPGA SYNTHESIS RESULTS
Tables 4 show the hardware resources devices utilization required for (Prob. Domain, Log Domain and Min-Sum) algorithms respectively. To measure the amount of complexity for the overall system in each type of soft decision decoder, Table 5 summarized the resource device utilization for each one. From all the comparisons it is easy to conclude that the optimum decoding algorithms in the term of complexity and performance are the Min-Sum algorithms due to their BER performance and suitable resource utilization.  Figure 31 shows the Xilinx SG simulation test result among the Bernoulli binary transmitted signal and the received signal where all type of decoders is to record the original signal with 960 ns delay. The hardware Co-Simulation results of the user data information and recovered information for the proposed algorithms are shown in Figure 32 respectively where there is a delay due to the operation of the extraction process and there are some Xilinx blocks causes delay.

CONCLUSION
In this work, the DCSK communication system has been developed with LDPC codes by using Prob. Domain, Log Domain and Min-Sum decoding algorithms and implementation have been done on the Kintex 7 FPGA development kit using Xilinx SG tools. The simulation results show that Min-Sum is outperformed from other decoders were achieved 8.8 dB improvement gain while when using Prob. Domain and Log Domain the improvement gain is approximately the same amount as 9 dB for the 16 Spreading factors. With FPGA the results show that the Prob. Domain algorithm has the highest complexity in term of the resources utilization where its consumed 1% of the slice registers, 64% of DSP and 4% of LUT while Log Domain has consumed 3% of the slice registers, 27% of DSP and 41% of LUT. The consumption of the Min-Sum algorithm was 1% of the slice registers, 23% of DSP and 2% of LUT which means that the optimum decoding algorithms are the Min-Sum algorithm due to its BER performance and suitable resource utilization. The implementation of hardware confirms that the proposed system is suitable for future communications, especially in real-time applications.


Fractional model for the spread of COVID-19 subject to government intervention and public perception

COVID-19 pandemic has impacted people all across the world. As a result, there has been a collective effort to monitor, predict, and control the spread of this disease. Among this effort is the development of mathematical models that could capture accurately the available data and simulate closely the futuristic scenarios. In this paper, a fractional-order memory-dependent model for simulating the spread of COVID-19 is proposed. In this model, the impact of governmental interventions and public perception are incorporated as part of the nonlinear time-varying transmission rate. In addition, an algorithm for approximating the optimal values of the fractional order and strength of governmental interventions is provided. This approach makes our model suitable for capturing the given data set and consequently reliable for future predictions. The model simulation is performed using the two-step generalized exponential time-differencing method and tested for data from Mainland China, Italy, Saudi Arabia and Brazil. The simulation results demonstrate that the fractional order model calibrates to the data better than its integer order counterpart. This observation is further endorsed by the calculated error metrics.

Introduction
The novel coronavirus pandemic, COVID-19, has impacted all aspects of our life and the disease continues to threaten global economy, public health, human interactions, and safety. Since the outbreak of the disease, there has been a collective effort to monitor, predict, and control its spread. For this purpose, different biological and non-biological techniques, tools, interventions, analysis, and experiments have been employed to understand the dynamics of disease spread.
Epidemiological mathematical models continue to play a pivotal role in understanding, capturing and predicting the spread of infectious diseases. In particular, first-order models such as SIR, SEIR, and many others, are intensively explored to mimic the collected data and make effective predictions. In the recent months, several such models have been developed to understand and predict the transmission of COVID-19 pandemic. See for example [1,2] .
Fractional-order models have been recognized as a powerful mathematical tool to study anomalous behaviors observed in many physical processes with prominent memory and hereditary properties. In general, such properties might not be satisfactorily captured by the integer-order counterpart models. This motivated researchers to consider fractional-order compartmental models for better understanding of dynamics of infectious diseases.
Since the outbreak of COVID-19 pandemic, a variety of fractional-order models for the spread of this disease are formulated [11][12][13][14] These models differ mainly in the considered compartments, mitigation effects, multi-fractional orders, and parameter forms.
Fractional-order epidemiological models are often obtained from their integer-order counterparts by replacing the firstorder derivative with a fractional one. However, this leads to a mismatch between the dimensions if the same parameters are used. This balance of dimensionality has been highlighted in some models but without much epidemiological interpretation of the new parameters. A common interpretation is that the parameters are observed in inhomogeneous time scale [15] corresponding to the fractional power of time unit. These parameters in the fractional model are commonly approximated by taking the fractional power of the classical parameters or by introducing an auxiliary parameter [16][17][18][19][20] .
In this article, we propose and simulate a system of fractional order differential equations for the dynamics of COVID-19 pandemic. The proposed model uses the nonlinear time varying transmission rate given in [1,21] which incorporates the effects of governmental interventions and public perception. In addition, an algorithm is developed so that for a given data set, the optimal fractional order of the model together with empirical values imitating the strength of governmental interventions are approximated. Consequently, reliable predictions for the future scenarios could be made.
The proposed model and algorithms are applied to the data from Mainland China, Saudi Arabia, Italy and Brazil. The simulation results suggest that the fractional order model has better agreement with the data than its integer order counterpart. This outcome is also supported by the error metrics considered in this paper. In addition, prediction trajectories for Saudi Arabia (a case with several changes in governmental interventions) and Brazil (a case with no governmental interventions) are included.
Since the focus is to provide adequate study and simulation of the fractional order model, we adopt the two-step generalized time differencing GETD2 scheme developed by Garrappa and Popolizio [22] to obtain the numerical results. However, to achieve a cost-effective implementation of this scheme, we use the partial fraction decomposition of the global Padé approximants of the two-parameter Mittag-Leffler function [23] .
The rest of this paper is arranged as follows. Next, the memory-dependent SEIR model is introduced. Section 3 is devoted to the discussion of reproduction number of the disease and local stability analysis. The computational algorithm is described in Section 4 . The simulation results are presented and discussed in Section 5 .

Memory-dependent COVID-19 Model
Following [1] , we assume the total population N at any given time to be distributed among the four compartments: susceptible, exposed, infectious, and recovered, labeled by S, E, I and R, respectively. That is, The assumed dynamics among these populations is depicted in Fig. 1 and can be describes as follows: • Susceptible individuals enter the exposed compartment at the rate βI/N, where β is a nonlinear transmission rate. • Exposed individuals progress to the infectious compartment at the rate σ, while the infectious individuals recover from the disease at a rate γ . • The population is assumed to have a natural death rate μ, which is given by μ = 1 ×365 , where denotes the lifeexpectancy.
• The birthrate of the population is assumed to be equal to the death rate.
Furthermore, we consider a variable P (t) to represent the public perception of the risk of the pandemic. This quantity increases with the number of infections at a rate φγ , and diminishes over time at a rate ω. The parameter φ is the proportion of severe cases.
The values σ −1 , γ −1 , and ω −1 can be interpreted as the mean latent period, mean infectious period, and mean duration of public reaction, respectively. These quantities together with μ, and β all have units of (time) −1 .
Based on the above assumptions, we propose the following model for the dynamics of COVID-19 pandemic: subject to the non-negative initial conditions The operator c D q denotes the Caputo fractional derivative of order q ∈ (0 , 1) defined by This derivative has been used successfully in the literature for compartmental models of different inf ectious diseases [3,4,6] .
In addition, the Caputo derivative allows to specify the initial data in a natural way.
In the special case q = 1 , the model in (1) reduces to a modification of COVID-19 model introduced in [1] . As observed in the right-hand side of (1) , the q th power of the parameters is used as suggested by Diethelm [3] and Angstmann et al. [18] . This is necessary to ensure the consistency in units since the left-hand side has the units of (time) −q .
For the transmission rate β in (1) , as in [1] , we use the formula In the integer order SEIR models, q = 1 , the parameters γ , μ, σ and β 0 describe the transfer rates between the different compartments and are measured per unit time, i.e. (time ) −1 . These parameters are the intuitive standard epidemiological rates such as, recovery, death, incubation, and transmission rates, which can be identified empirically or statistically. As for the fractional order model, q ∈ (0 , 1) , these transfer rates retain similar epidemiological meanings. However, they are to be measured on time scale (time ) −q which might not be feasible. Raising these parameters to the power q is an approach to approximate the corresponding transfer rates measured on time scale (time ) −q . An advantage of this approach is that the model reduces to the classical counterpart when q = 1 .

Basic reproduction number and equilibria stability
The basic reproduction number, R 0 , is the expected number of secondary infections that a single infectious individual will generate on average within a susceptible population. This number is arguably the most important quantity in determining the nature of an epidemic. In practice, a disease spreads whenever the reproduction number exceeds 1 ( R 0 > 1 ). On the other hand, an outbreak with a reproduction number below 1 will gradually disappear.
An equilibrium point is defined as follows.
Definition 3.1. [24] A point u * is an equilibrium point of the Caputo fractional dynamical system c D α Using this definition, we see that is an equilibrium point of (1) . This point is called a disease-free equilibrium. By applying the next-generation operator approach (see [25,26] ), the reproduction number R 0 is obtained as the spectral radius of the next generation matrix , which is the transmission rate β evaluated at u * . Thus, R 0 of (1) is given by Next, we analyze the local asymptotic stability of the disease-free equilibrium. For this purpose, we state the following result which follows by linearizing the system about the equilibrium point (as in Ahmed et al. [27] ) and then applying Theorem 2 of Matignon [28] .

Lemma 3.2.
Consider the following fractional-order system: then u * is locally asymptotically stable.
Based on Lemma 3.2 , we have the following stability result for system (1) .
Proof. The Jacobian of the right hand side of (1) at u * is given by The eigenvalues of the matrix J u * are given by and thus λ 6 is also negative. Consequently, if R 0 < 1 , then condition (3) is satisfied and Lemma 3.2 implies that u * is locally asymptotically stable.

Computational algorithm
Here, we describe briefly the numerical algorithm used to simulate model (1) . For this purpose, we consider the matrix form of the system, Table 1 Model parameters with possible range of values together with references.

days [1]
Proportion of severe cases, φ 0.2 [1] Mean duration of public reaction, ω −1 11.2 days [1,21] To solve this system numerically we use the GETD2 scheme developed in [22] . It has been shown that this scheme when applied to (4) , has a (1 + q ) −rate of convergence. In addition, it has elegant stability properties that makes it robust for stiff systems without need to impose severe restrictions on the step size.
where g j = g(v j ) . To initialize the scheme we take v 0 = u 0 and compute v 1 as For efficient implementation of the GETD2 scheme, we use a technique based on the partial fraction decomposition of the global Padé approximant of type (7,2) developed in [23] for the Mittag-Leffler function.

Model simulation
In this section, we perform numerical simulations of the fractional order COVID-19 model (1) using the GETD2 scheme (5) . We consider different case studies in which we calibrate the model to the data on confirmed COVID-19 infections in Mainland China, Italy, Saudi Arabia and Brazil. In addition, we present some predictions regarding possible spike and/or decay in the number of infections based on the intensity of behavioral response of population. For these purposes, we begin with a discussion about the parameters and their values.
The range of values of the parameters of the COVID-19 model (1) are listed in Table 1 . For our simulations, we use the following values.
On the other hand, as outlined in Algorithm 1 , the optimal values of the fractional order q and governmental intervention strength ν will be identified using the given data sets.
Based on the timeline of governmental interventions, the numerical values assumed by the parameter ν may vary on different time intervals. For the implementation of Algorithm 1 , feasible values of this parameter should reflect the intensity of the interventions. For example, periods with strict regulations should have no smaller values of ν than those with relaxed interventions. As such, at the initial phase of the pandemic and before any governmental intervention, we set ν to 0. On the other hand, based on the studies in [29,31] , substantial values of ν can be considered when accounting for the effect of mandatory use of face mask and physical distancing in public places.
The intensity κ was estimated in the study of Influenza pandemic by He et al. [21] as κ = 1117 . 3 . However, by studying and matching different data sets of COVID-19 infections, we observed that to better reflect the human behavioral response, piecewise adjustments to this value are to be considered to reflect the interventions timeline. In particular, we use κ = 0 at the early stage of the outbreak since it is assumed that adequate awareness and precautions are not yet in place. For approximating q and ν in Algorithm 1 , we use the value κ = 1117 . 3 for later time intervals. Nonetheless, when most of the governmental lockdown measures are being lifted, it is noted that the intensity of human behavioral response plays a major role in the decline of the number of infections. Thus, once the optimal values of ν and q are approximated, we fine-tune κ in order to achieve better matching.
Algorithm 1 Approximating the optimal q and ν pair Step 1 Create a list { q j } l j=1 ∈ [ a, 1] , a > 0 , for values of q .
Step 3 On each sub-interval S i , create a list V i of feasible values of ν for time-interval S i .
Step 4 Construct a collection C of all ordered combinations c Step 5 Define an empty list M . For q j , j = 1 , 2 , . . . , l Define an empty list E.
For c ∈ C Solve system (4) Use an error metric (say mean square error) to quantify the error betweenthe computed infection numbers and the reported ones. Store the triple ( q j , c, error) in E. Select from E the triple with minimum error value and store it in M . Select from M the triple with minimum error value and return the corresponding c and q j .
In all our simulations, the values of the RMSE between the approximated infected population and the actual data is provided for optimal q and q = 1 . Moreover, to give a perspective of the corresponding relative error, we also present the values of the mean absolute percentage error (MAPE).

Case study: outbreak in Mainland China
In Table 2 , we consider the reported data of cumulative COVID-19 infections in Mainland China from January 19, 2020 to March 6, 2020, as provided in [32] . To simulate the outbreak of the pandemic, we use the β 0 = 0 . 52 , μ = 1 / (76 . 91 × 365) , together with those given in (7) . Furthermore, the initial populations are As per the lockdown records of Chinese Administrative Division, we consider the governmental intervention in two phases. The First phase corresponds to the period Jan 20 -Feb 02 and the second one to the period Feb 03 -March 06.
Applying Algorithm 1 , we obtain q = 0 . 93 as an approximation of the optimal fractional order. The corresponding values of governmental intervention strength ν are given in Table 3 . For comparison purposes, we also provide the optimal values of ν for the integer order problem, q = 1 . Furthermore, the values of RMSE and MAPE in approximating the infected population are also appended. The error metrics in Table 3 suggest that the fractional order model performs better in terms of the RMSE. However, the MAPE seems to support the integer order model. This blended outcome can be understood in view of Fig. 2 . As can be observed, the absolute errors for the fractional order model are smaller throughout most of the period and thus smaller RMSE. On the other hand, although the relative errors are smaller in general for the fractional model, they attained some large values towards the end of the period, which in turn impacted the MAPE.   5 provides an insight on the effect of governmental interventions. As observed, the change in governmental intervention strength ν on Feb 3rd yields an instantaneous change in the rate of change of exposed population. This phenomenon is expected, because governmental interventions such as inter-city lockdown, curfew, school closure among others, are expected to mitigate the contacts between susceptible and infectious individuals. However, some increase in the number of infections is still observed. This accounts for the individuals who were incubating the virus before implementation of the governmental interventions.
Our numerical experiments also demonstrate the robustness of the fractional-order model as a generalization of the integer-order model. It is observed that the fractional order model gives extra degree of freedom, via the order of differentiation q, for data fitting.

Case study: outbreak in Saudi Arabia
Since the outbreak of COVID-19 in Saudi Arabia, the Kingdom has witnessed several phases of governmental interventions aimed at mitigating spread of the diseases. Due to seasonality, special breaks and festivals, there have been a number of relaxations and tightening of these interventions. To simulate the outbreak of the Pandemic within the Kingdom, we only adopt those interventions that were implemented at a national scale, as provided in https://www.garda.com/crisis24/news-alerts/334921/saudi-arabia-curfewhours-modified-for-ramadan-april-21-update-29 . We consider the parameters in (7) together with β 0 = 0 . 54 , μ = 1 / (75 . 13 × 365) .

The initial conditions are
Based on the data provided in https://ourworldindata.org/coronavirus-source-data , our simulation is performed for the period of March 03 to July 17 excluding the days 4th, 5th and 7th of March, where there was no reported data.
Using Algorithm 1 , we obtained an approximation of the optimal q and ν pair, with q = 0 . 87 . For comparison purposes, we present in Table 4 the ν values corresponding to both q = 0 . 87 and q = 1 . In addition, the values of the RMSE and MAPE in approximating the infected population for these values of q are presented. Brief remarks on the governmental actions and the values of intensity κ are also provided.
We observe that the values of ν presented in Table 4 reflect the intensity of the control measures in the sense that periods with more restrictive interventions have higher values of ν. From an epidemiology perspective, these interventions such as imposing a curfew, restricting the gatherings size, among others, reduce interactions and consequently mitigate the spread of the disease. As a result, the fitted values of ν tend to reflect the epidemiological reality since bigger values of ν will imply smaller values of the transmission rate β in (2) . Fig. 6 includes plots of the daily infections obtained from the fractional model (1) and plots of the parameters ν, κ and β. It is observed that the polynomial behavior of the data at the early stage of the outbreak is better captured by the power-law corresponding to q = 0 . 87 than by the exponential or almost exponential law, q = 1 or q = 0 . 98 . However, the final stage of simulation is better captured with q = 0 . 98 . This Further illustrates the effectiveness of fractional order models  for making predictions. The outcome of our simulation for the cumulative infections is presented in Fig. 7 . We note that q = 0 . 87 yields a satisfactory matching. Using the calibrated model (1) , predictions for the spread of the disease subject to different values of κ are shown in Fig. 8 . We assume that the governmental intervention continues to be the same as recent period. It is shown that, a drop in the intensity κ below the value between 21-June and 17-July will yield a spike in the infection numbers. On the other hand, maintaining same intensity will help obtain a slow decline of the infection numbers. However, a notable decay in the number of infections could be achieved if there is a significant increase in the intensity of the behavioral response. This corresponds to improved compliance with proper hand-washing, face masking, physical distancing, among others.

Case study: outbreak in Italy
In the early outbreak of COVID-19 in Italy, the country was hit badly and a number of governmental interventions were initiated. To simulate this outbreak using (1) , we consider   together with those provided in (7) . The initial conditions are The timeline of governmental interventions considered for this simulation are based on the report by Briscese et al. [33] , up to 3-April, together with other reopening information available on https://www.bbc.com/news/world-europe-52435273 . Our simulation covers the period of 31-January 2020 to 17-July 2020, based on the data provided in https://ourworldindata.org/ coronavirus-source-data . However for the period February 01-21, 2020, no new infections were reported. Consequently, we excluded these dates from our simulation.
Algorithm 1 applied to the fractional order model (1) yields an approximation for the optimal q and ν pair, with q = 0 . 88 . The corresponding values of ν are provided in Table 5 . In addition, brief descriptions of the governmental actions together with the adopted values of κ are also provided. Similar to Saudi Arabia case, we observe that the fitted values of ν reflect the intensity of governmental actions and are consistent from an epidemiological point of view.
Plots of the daily and cumulative infections obtained from the model (1) are presented in Figs. 9 and 10 , respectively. We observe from these plots that the fractional order model has a better matching. This observation is also supported by the values of the error metrics included in Table 5 .

Case study: outbreak in Brazil
Information from various sources indicate that there were no governmental control measures to mitigate the COVID-19 outbreak in Brazil. Some of the events in the early days of the outbreak can be found in an article by Reuters https://www.reuters.com/article/us-health-coronavirus-brazil-timeline/timeline-key-moments-inbolsonaros-handling-of-covid-19-crisis-idINKBN2482UL . Thus, for simulating the outbreak we set ν = 0 . The simulation covers the period from 26th February to 17th July using the data provided in https://ourworldindata.org/ coronavirus-source-data . We exclude the days February 27-29 and March 02-04, 08 and 10, since no infections were reported. The other parameters used in simulating (1) for the outbreak in Brazil are β 0 = 0 . 75 , μ = 1 / (75 . 88 × 365) , together with those provided in (7) . The initial conditions are With ν = 0 for the whole simulation time, implementation of Algorithm 1 reduces to approximating only the optimal fractional order q, which is approximated as q = 0 . 70 .
Our examination of the infection data in Fig. 11 shows that there is a significant change in the reported numbers around day sixty. Under no governmental interventions, the parameter in our model that can be used to account for this dramatic change is the human behavioral response κ. For this purpose, we split our time interval into [0, 60] and [61,134]. Then, we fit appropriate values of κ for both q = 0 . 70 (the optimal value of q ) and q = 1 . The calculated values are shown in Table 6 together with the RMSE and MAPE in approximating the daily number of infections.
Plots of the daily infections obtained from (1) are presented in Fig. 12 . As can be observed, the integer order model could not provide a reasonable matching of the data. On the other hand, the fractional order model with q = 0 . 70 is able to capture well the initial spread of the disease. Plots of the cumulative infections are also presented in 13 .
Overall, we remark that modeling the COVID-19 outbreak in Brazil is not a trivial task due to the erratic nature of the reported data. However, using the extra degree of freedom offered by the fractional order model through the fractional order q, we are able to find a possible fit for daily infection data.  Table 6 for q = 0 . 70 . Then we predict possible  infection numbers for the period from July 18 to August 16. The plots suggest that for a significant reduction in infections, a substantial increase in κ is required. This corresponds to improved compliance with proper hand-washing, face masking, physical distancing, among other public health measures.

Discussion
The proposed model aims for improved fitting of the real data and thus more realistic predictions. It provides extra degrees of freedom through the fractional order q and the transmission rate β(t) that includes the government intervention and intensity of human behavioral response.   From the outcomes of our numerical simulations, we note the following observations which can help in policy making for public health interventions.
• The data matching reveals that strict governmental lockdown measures result in large values of the parameter ν, and thus rapid decrease in the transmission rate. This is notably observed in the case of China.
• The intensity of human behavioral response, such as adherence to masking, hand washing, physical distancing, among others, compliments the governmental interventions in mitigating the transmission of the disease. This is demonstrated by the predictions in Fig. 8 , where we see that for fixed governmental interventions, the number of infections is inversely proportional to the behavioral response κ.
• Complete absence of governmental interventions results in growing transmission of the disease even with strong human behavioral response. This scenario is demonstrated in the Brazil model, where we set ν to zero since there is no record of governmental interventions. As can be observed, the human behavioral response alone is not adequate to mitigate the spread of the disease. Rather, it complements the governmental interventions. • As a consequence of the above observations, gradual lifting of governmental measures is advisable. For example, interventions such as mandatory physical distancing and masking in public places should be in place. These less strict interventions can still be well complemented by the behavioral responses and consequently help flatten the curve. The prediction for Saudi Arabia in Fig. 8 is a typical reference.

Conclusion
In this paper, we proposed a nonlinear fractional order model to study the dynamics of COVID-19 pandemic. Using the data of the outbreak in Mainland China, Saudi Arabia, Italy, and Brazil, the effectiveness of the proposed model in simulating the spread of the disease is demonstrated. Predictions based on varying intensity of human behavioral responses are provided. By examining the mean square error of the predicted number of infections and the reported data, we found that the proposed fractional model yields better matching than the integer order model. From the perspective of relative error, the values of mean absolute percentage error confirm this conclusion for the simulation of outbreak in Saudi Arabia, Italy and Brazil. On the other hand, for China, these two error metrics have mixed outcomes.