Non-Markovianity in a collision model with environmental block

We present an extended collision model to simulate the dynamics of an open quantum system. In our model, the unit to represent the environment is, instead of a single particle, a block which consists of a number of environment particles. The introduced blocks enable us to study the effects of different strategies of system-environment interactions and states of the blocks on the non-Markovianities. We demonstrate our idea in the Gaussian channels of an all-optical system and derive a necessary and sufficient condition of non-Markovianity for such channels. Moreover, we show the equivalence of our criterion to the non-Markovian quantum jump in the simulation of the pure damping process of a single-mode field. We also show that the non-Markovianity of the channel working in the strategy that the system collides with environmental particles in each block in a certain order will be affected by the size of the block and the embedded entanglement and the effects of heating and squeezing the vacuum environmental state will quantitatively enhance the non-Markovianity.


I. INTRODUCTION
The interaction between a quantum system and an environment leads to a non-unitary time-evolution of the state of system. Such irreversible dynamics can be described by the theory of open systems [1]. Understanding the dynamics of an open quantum system is an essential question in quantum information processing [2][3][4][5][6], quantum biology [7,8] and quantum optics [9][10][11]. Usually the dynamics can be classified as Markovian or non-Markovian cases. In the Markovian case, the dynamics is typically characterized by the master equation in the so-called Lindblad form which corresponding to a completely positive and trace-preserving (CPT) map. In particular, if the Lindbladian of the master equation is time-independent, it gives rise to a dynamical semigroup of maps [1]. However the Markovian demonstration is not always adequate. For instance, the CPT condition is violated if a quantum system is strongly coupled to the environment, so that the dynamics becomes non-Markovian [12]. In recent years a number of criteria characterizing non-Markovianity have been proposed, from different perspectives, basing on the dynamical divisibility [13][14][15][16], back-flow of information characterized by trace distance [4,17], Fisher information [18], mutual information [19], relative entropy [20,21], accessible information [22], Gaussian interferometric power [23] and response functions [24]. For recent reviews, see [25,26]. These criteria or measures help us to distinguish whether a dynamics is Markovian or not; however, in general they do not agree with each other in detecting the emergence and quantifying the degree of non-Markovianity.
An alternative approach to studying the dynamics of an open system is the so-called collision model (CM). In the CM based scheme, the continuous time-evolution of a system is simulated by a sequence of system-environment collisions representing the interactions between system and the environment. The environment is represented by an ensemble of uncorrelated identical particles. If the system collides with each environmental particle in a sequential way, the dynamics of the system is Markovian since the environmental particle in the upcoming collision is fresh and thus contains no information of the history. By introducing the intra-collision between environmental particles [27][28][29][30][31][32], long-range systemenvironment collisions [33], the correlations among environmental particles [34][35][36], and a composite structure of the system [37], the non-colliding environmental particle could have the possibility to carry the information of the history, thus the dynamics of the system may become non-Markovian. Very recently, the idea of CM is also adopted in the content of thermodynamics [38]. We note that the unit representing the environment is a single particle in the aforementioned modified CMs. However, this is abridged in simulating the details of system-environment interactions as well as the diverse states of the environment: on the one hand, the approaches of the memory recover, such as recovering from the latest to the earliest time or the reverse, cannot be manifested through a collision between the system and the single-particle environmental unit; on the other hand, the single-particle unit excludes the nonlocal many-body correlations of the environment. Therefore, a CM with more complicated environmental unit would be of interest not only in simulating the dynamics of open system but also in exploring the essence of the non-Markovian process.
In this work, we will consider an extended CM that the environment is represented by an ensemble of identical blocks. Each block consists of a number of particles. The systemenvironmental interactions are simulated by the collisions between the system and the environmental blocks. The internal structures of the environmental blocks enable us to explore how the system-environment interactions and the environmental states affect the non-Markovianity of the dynamics. Here, we will mainly discuss our CM in the realization of all-optical system. We will derive a necessary and sufficient condition of the non-Markovianity in our CM basing on the indivisibility of the dynamical maps and show the evidence of its equivalence to the non-Markovian quantum jump through the pure damping process of a single-mode field. Thanks to the internal structure of the environments introduced by the environmental blocks, we will investigate the effects the strategies of the system-environment collisions which are related to the approaches of memory recover in a realistic dynamical process. We will also study how the properties of the environmental state, such as temperature, squeezing and entanglement, affect the non-Markovianity. Because the collisions between different modes are taken place at the beam-splitters (BSs), we can use the Hamiltonian of two-mode linear mixing,Ĥ ∝â †b +âb † (â andb denote different modes), to create such interactions. The usage of the two-mode mixing Hamiltonian in describing the interactions between a bosonic mode and a reservoir implies the potential of our CM in simulating the dynamics of the open quantum optical system. This paper is organized as follows. In Sec. II, we introduce the idea of our model and apply this idea to the Gaussian channel in an all-optical system. In Sec. III A, we review the measure of non-Markovianity of the Gaussian channel recently proposed in Ref. [16] and derive an explicit expression of the measure of non-Markovianity for our model. In Sec. III B, we discuss the non-Markovianity of the Gaussian channel with vacuum environmental state through two strategies of system-environment collisions and simulate the pure damping process of a single-mode field. The necessary and sufficient condition of the non-Markovian channel in the vacuum environment is given as well. In Sec. III C, we investigate the effects of the temperature and squeezing on the non-Markovianity of the channels with generic Gaussian state. We compare the non-Markovianities of the channels with product and entangled states in in Sec. III D. The conclusion is drawn in the last section.

II. SIMULATING COLLISION MODEL IN ALL-OPTICAL SYSTEM
A. The CM with environmental block In our model the unit to represent the environment is a block rather than a single particle as that in the standard CM. An environmental block is consisted of a number of particles and all the blocks are supposed to be identical. For the explicitness of explanation, we label the l-th block as B l with l = 1, 2, ..., L and the particle in B l as E l, j with j = 1, 2, ..., L B . We have set the number of blocks to be L and the number of particles in a block to be L B . We will denote L B to the size of the block hereinafter. As shown in Fig. 1, our model works via the following steps.
Step 1. As the start, the system collides with the environmental block B l with l = 1. The S -B l collision is accomplished by the system sequentially colliding with the particles E l, j in a certain order, for example, the ascending order of j. After the S -B l collision, each particle in B l carries part of the information of system at discrete time points during the past. More precisely, the earliest information of system is stored in E l,1 while the latest information is stored in E l,L B .
Step 2. In this step the intra-collision of the environment takes place. The block B l collides with a fresh block B l+1 . The B l -B l+1 collision is accomplished by respectively colliding each pair of E l, j and E l+1, j particles. After the B l -B l+1 collision, part of the lost information has possibility to be transferred to the block B l+1 .
Step 3. The collision between the system S and the block B l+1 takes place. The S -B l+1 collision enables the system to Step 1 Step 2 Step 3  The discrete dynamics of an open system is simulated by a series of collisions between system and environmental blocks. In step 1, the system particle (circle labeled by "S ") collides with block B l (squares labeled by "B"). The S -B l collision is implemented by collides "S " with each particle E l, j (circles labeled by "E") in B l in ascending order of j. In step 2, block B l collides with B l+1 . The B l -B l+1 collision is implemented by colliding each pair of E l, j and E l+1, j particles, see (b). In step 3, the system collides with block B l+1 and then go to step 2 with l → l + 1 for iteration. There are two strategies of the S -B l+1 collision. In strategy 1, the S -E l+1, j collisions take place in the same order of j as that in S -B l collision, i.e. the information firstly input to the environment firstly outputs, see (c); in strategy 2, the sequence of S -E l+1, j collisions take place in the order of j which is opposite to that in the S -B l collision, i.e., the information firstly input to the environment lastly outputs, see (d). In Boxes (b)-(d), the system circles with dashed border denote the initial position and the ones with solid border denote the final position. The blue circles E l, j in the dashed rectangle belong to block B l .
get the lost information back. In this step, we may implement the S -B l+1 collision in two different strategies. In strategy 1, the system S always sequentially collides with E l+1, j s in the same order of j as in step 1, i.e., the information first input to block B l is first output. In strategy 2, the system collides with E l+1, j s in the reverse order of j as that in S -B l collision, i.e. the information first input to the environmental block is last output. Once the S -B l+1 collision is completed, we go to step 2 with l → l + 1 to iterate.
We would like to emphasize that in our CM the systemenvironment collision is represented by the S -B l collision. This implies a coarse-graining of the system evolution. We are only interested in the evolved system state after each complete S -B l collision, and the intermediate system states after each S -E l, j collisions are considered to be the details hiding in the system-environment collision.

B. A scheme in the all-optical system
Our CM can be implemented in the all-optical system which is composed of an array of BSs. The realistic optical system can be perfectly controlled, integrated and scaled up. The system and environmental particles are represented by the independent optical modes propagating along different paths. The collisions between any two particles can be realized by mixing two corresponding input modes at the BS. We recall that a BS transfers two input modesâ in = [â in 1 ,â in 2 ] T into output modesâ out = [â out 1 ,â out 2 ] T = Sâ in , whereâ is the annihilation operator (the superscript T denotes the transpose), and S is the scattering matrix with r = sin θ and t = cos θ being the reflectivity and transmissivity of the BS and θ ∈ [0, π/2]. In Eq. (1), we have set the reflected mode to be the output mode. For θ = π/2, both the input modes will be completely reflected thus indicating the strength of interaction between two modes is zero, while for θ = 0, both the input modes will completely transmit thus indicating a swap operation of two modes. In our model we denote the BS that mixes the system and environmental modes as BS1 with the reflectivity and transmissivity being r 1 = sin θ 1 and t 1 = cos θ 1 , and the BS that mixes two environmental modes as BS2 with the reflectivity and transmissivity being r 2 = sin θ 2 and t 2 = cos θ 2 . Therefore, the strengths of system-environment and environmentenvironment interactions can be tuned by varying θ 1 and θ 2 , respectively. In the following, we will describe how to realize the S -B l , B l -B l+1 and S -B l+1 collisions in the all-optical system.
S -B l collision. The collision between the system S and the environmental block B l can be simulated by sequentially mixing the system mode,â S , and each environmental mode,â l, j , at a series of identical BS1s. For instance, in the case of S interacts with each E l, j in the ascending order of j, theâ S mode firstly interacts withâ l,1 at the first BS1s, and then the output (reflected mode) ofâ S will interact withâ l,2 at the second BS1 and so on so forth. The S -B l collision is completed until thê a S mode has interacted with theâ l,L B mode.
B l -B l+1 collision. In order to simulate the B l -B l+1 collision, each of the output (reflected) mode ofâ l, j is guided to interact with the correspondingâ l+1, j mode of B l+1 individually. It should be noted that, different from the S -B l collision, the interactions betweenâ l, j andâ l+1, j take place at BS2s.
S -B l+1 collision. As mentioned before, there are two strategies for the S -B l+1 collisions. For strategies 1 and 2, the output mode ofâ S is guided to interact with the output ofâ l+1, j modes in the same and reverse orders of j to that in the previous S -B l collision, respectively. Again the S -B l+1 collisions take place at BS1s.
In Fig. 2, we show the setups working in strategies 1 and 2, respectively. Both setups are composed of the building blocks (dashed boxes) that realize the S -B l , B l -B l+1 and S -B l+1 collisions. The operators with primes denote the intermediate evolved state and the operators with superscripts 'in' In the S -B l collision, the system mode sequentially interacts with the (l, j)-th environmental modes at BS1s. After the S -B l collision, the output of the (l, j)-th environmental modes, i.e. the reflected modes, are guided to interact with the fresh (l + 1, j)-th modes of B l+1 individually to simulate the B l -B l+1 collisions. In the S -B l+1 collision, the output of (l, j)-th modes are discarded and the output of (l + 1, j)-th modes are guided to interact with the system mode. The difference between strategies 1 and 2 are presented in the dotted boxes (shaded gray). In strategy 1, the interactions between system and the (l + 1, j)-th modes are implemented with the same order of j to that in the previous S -B l collision, while in strategy 2, the interactions between system and the (l + 1, j)-th modes are implemented with the reverse order of j to that in the previous S -B l collision. By concatenating the B l -B l+1 and S -B l+1 building blocks, the stroboscopic evolution of the system mode can be simulated.
and 'out' denote the initial and final state, respectively. The dotted boxes (shaded gray) show the distinction of strategies 1 and 2, i.e. the system mode interacts with environment modes in different orders of j.
We can concatenate such processes to simulate the discrete dynamical evolution of the system mode.
Sup-pose the system will interact with L environmental blocks, then the channel can be described, with the help of scattering matrix S(L), by [â out S ,â out 1,1 , ...,â out l, j , ...,â out L, Hereinafter, we use the subscript S to denote the system mode and the pairwise (l, j) to denote the j-th mode in the l-th environmental block.
The (L B L + 1)-dimensional scattering matrix S(L), for L ≥ 2,has the form where the matrices in the cumulative product are in the descending order of l from right to left. In Eq. (2), the matrix S B l B l+1 describing B l -B l+1 collision is given by where I n is the n × n identity matrix. The matrixS S B l describing the S -B l collision has two different forms with respect to strategies 1 and 2. In strategy 1, we haveS S B l = S S B l for all l with S S B l = L B j=1 S l, j , while in strategy 2, we havẽ The matrix S l, j describing the interaction between system and the (l, j)-th environmental mode is given by Note that S S B l and R S B l are obtained by multiplying S l, j in the ascending and descending orders of j, respectively. From Eqs.
(2)-(4), we see that the property of the 'bare' channel (i.e. with vacuum environment) is determined by the reflectivities and transmissivities of BS1 and BS2.
C. The dynamics of the system mode It is convenient to study the dynamics in the channel with the characteristic function formalism [39]. Actually, the density operatorρ of a quantum state is equivalent to the characteristic function in presenting the probability distribution. The symmetrically ordered characteristic function is defined by χ(ν) = Tr[D(ν)ρ] with the Weyl displacement operator D(ν) = exp (νâ † − ν * â ). Thus we can represent the density operator of a bosonic system which is defined in an infinitedimensional Hilbert space with a complex function. In particular, in terms of the characteristic function, the first and second moments are sufficient to characterize the Gaussian state which is widely used in the quantum information processing with continuous variables system [40]. Reversely, the density operatorρ can be represented in the Weyl expansion with the χ(ν) acting as the weight function, i.e.,ρ = d 2 νχ(ν)D(ν)/π. In our model, the joint characteristic function of the multimode input stateρ in J is given by χ in j=1D l, j (ν l, j ). The subscript 'J' denotes the joint modes.
Initially the modes of the system and the environment are uncorrelated, the joint input characteristic function is thus calculated by where ν l = [ν l,1 , ν l,2 , ..., ν l,L B ] T and χ in l ( ν l ) is the characteristic function of the l-th block.
The input-output relation of the joint characteristic function after L times system-environment collisions is just determined by the scattering matrix S(L) through the following formula, as detailed in Appendix, where χ out,L J ( ν) is the joint characteristic function of the output modes and S −1 (L) = S † (L) is the inverse of the scattering matrix S(L). Since we are interested in the evolution of the system mode, we need to trace out all the environmental modes in Eq.
where c 1 = [c 1,1 , c 1,2 , ..., c 1,L B L+1 ] T is a column vector that equals to the transpose of the first row of S(L).
Here we concentrate on the cases that all the input modes are initially in Gaussian states. A state of continuous variable system is Gaussian if its characteristic function is Gaussian. The Gaussian state are the resources for a plethora of quantum information and communication protocols with continuous variables [40]. In our model, it is easy to prove that the channel with Gaussian environmental state will always keep the Gaussianity of the system state, therefore we may regard the channel described by Eq. (2) to be a Gaussian channel. We recall that the characteristic function of a generic Gaussian states is given by The real parameter A and complex parameters B and C are related with the properties of the Gaussian state as where n is the thermal mean photon number, r is the squeezing strength, φ is the rotating angle and α is the complex displacement [42]. So far, we are able to describe the dynamics of the system mode with the help of the scattering matrix in the characteristic function formalism. In our CM, the correlations between the system and each environmental blocks are built after each S -B l collisions and present during the whole evolution of the system mode, because all the environmental modes are traced out after the S -B L collision. This is different to the existing CMs in which the system-environment correlations are erased, before or after the B l -B l+1 collision, in each step [29,33]. It is shown that the system-environment correlations play important role in establishing the non-Markovianity [29], thus our CM has the potential in studying the role of systemenvironment correlations in a rather flexible way.

III. NON-MARKOVIANITY OF THE GAUSSIAN CHANNEL
A. Measure of non-Markovianity In Ref. [16], a measure of non-Markovianity of the Gaussian channel by quantifying the degree of the violation of dynamical divisibility is presented. We will employ this measure in our model. According to Eq. (7), we can represent the evolved system mode after l times system-environment collisions with the following dynamical map on the input characteristic function, or, in terms of the covariance matrix, The covariance matrix is the second moment of the characteristic function and its elements are defined by where {·, ·} is the anticommutator, · is the expected value and 2i. The symmetrically ordered moments can be computed can be computed as where the subscript "symm" denotes the symmetrical order. The dynamical map E l is always CPT and can be always formally split as the following, where Φ l,l−1 is an intermediate process that maps the χ out,l−1 S to χ out,l S , and the "•" represents the composition of the maps. The divisibility of the Gaussian channel can be determined by Φ l,l−1 . If Φ l,l−1 is CPT for all l, then the dynamics is divisible and hence Markovian. Otherwise, if Φ l,l−1 is non-CPT for some values of l, then the dynamics is indivisible and hence non-Markovian.
For a generic Gaussian channel, E l has the following form, where X l and Y l are 2 × 2 real matrices. The necessary and sufficient conditions of the CPT property of Φ l,l−1 is the semipositive definiteness of the following 2 × 2 matrix [43], with X l,l−1 = X l X −1 l−1 , Y l,l−1 = Y l − X l,l−1 Y l−1 X T l,l−1 , and Ω = [0, 1; −1, 0] being the single mode symplectic matrix. The negative eigenvalue of Λ l contributes to the non-CPT of Φ l,l−1 and, as a consequence, the non-Markovianity of the Gaussian channel. Thus the non-Markovianity of the Gaussian channel can be measured by the sum of the negative eigenvalues of all the Λ l , where λ l,k are the eigenvalues of Λ l . Eq. (17) is the expression of the non-Markovianity measure for our CM and will be used in the analysis hereinafter. We would like to point out that although Eq. (17) is sufficient and necessary in characterizing and quantifying the non-Markovianity, it is computable only when the channel can be completely characterized. Fortunately, it is possible to completely demonstrate the Gaussian channel of our CM in the all-optical system.
We restrict the initial system state to be a Gaussian state with the characteristic function as expressed in Eq. (8). The corresponding covariance matrix is given by Once the scattering matrix S(l) is constructed and the environmental Gaussian state is specified to A l, j = A E , B l, j = B E , and C l, j = C E , we can compute the evolved characteristic function of the system mode with the help of Eq. (6) and then obtain the corresponding covariance matrix as the following, Accordingly, we have the explicit forms of X l and Y l in Eq. (15) as X l = Re(c 1,1 (l)) −Im(c 1,1 (l)) Im(c 1,1 (l)) Re(c 1,1 (l)) , and The matrices X l and Y l can fully demonstrate the Gaussian channel. One can see that the properties of the Gaussian channel is determined by the r 1 = sin θ 1 and r 2 = sin θ 2 in terms of the c 1,k (l) as well as the properties of Gaussian environmental state in terms of A E and B E . It is straightforward to obtain the eigenvalues of Λ l as a function of A E , B E , and c 1,1 (l),

B. Vacuum environmental state
We start with the vacuum environmental state, i.e. A E and B E are both zero. For vanishing A E and B E , the eigenvalues of Λ l , Eq. (16), are λ l,+ = 1 − c 2 1,1 (l)/c 2 1,1 (l − 1) and λ l,− = 0. Using Eq. (17), we can obtain the non-Markovianity of the channel with vacuum environment, N vac (L), as The above expression indicates a necessary and sufficient condition of the non-Markovian Gaussian channel with vacuum environmental state, i.e., In Fig. 3, we show the Markovian and non-Markovian regions in the plane expanded by θ 1 and θ 2 in strategy 1 for different sizes of the environmental block. For L B = 1, our model is reduced to the standard CM and thus the values of θ 1 and θ 2 characterize directly the strengths of the systemenvironment and environment-environment interactions. For the case of θ 1 /π = 0.5, the system mode is complete reflected after each S -B l collision and thus isolated from the environment. As θ 1 decreasing, the system-environment interaction is activated. For the limit case of θ 2 /π = 0.5, the dynamics of system is Markovian since the strength of B l -B l+1 collision is zero. In the opposite side, θ 2 /π = 0, the dynamics of system is strongly non-Markovian since the B l -B l+1 collision is a perfect swap operation. As a consequence, for a fixed θ 1 , we can switch the channel from Markovian to non-Markovian cases by tuning θ 2 . There are critical θ 2 s that separate the Markovian and non-Markovian regions.
We remind that the CM simulates the dynamics of an open quantum system in a stroboscopic way, i.e., the time interval between two successive system-environment collisions is The stroboscopic evolution is cut off after the system colliding with L = 50 blocks. For θ 1 /π = 0.5, the system mode is isolated from the environments. For the limit of θ 2 = 0 the evolution of the system mode is unitary and θ 2 /π = 1 the dynamics of the system is Markovian. With the size of block increasing, the non-Markovian range shrinks and converges for L B > 8. The inset is a zoom in of the range of small θ 1 /π, i.e., the strong coupling of the system and environment.
τ. The overall effect of one collision between system and an environmental block is mapping the system state at time t to t + τ, regardless of the microscopic details in the collision. Namely, we may regard two S -B l collisions with different L B to be equivalent if they map the same initial state to the same final state. Basing on this idea, we are able to study the cases of L B > 1 in a unified frame. This can be realized in the following approach: if the reflectivity of BS1 is r 1 for L B = 1, then the reflectivity of BS1 is set to be r 1/L B 1 for L B > 1. This guarantees the identity of the effective strengths (or θ 1,eff ) of the system-environment interactions with different L B , because the successive S -E l, j collisions in block B l is Markovian. From Fig. 3 we see that the non-Markovian region shrinks with the size of block increasing. However the boundary of non-Markvoian region converges for large L B . In the plot we numerically compute the critical θ 2 as a function of θ 1 with the size of block up to L B = 16. For the case of weak coupling of system and environment the boundary converges fast, while for strong coupling of system and environment the boundary converges slow.
We note that for L B = 1 the channel with strategy 2 is equivalent to that with strategy 1. However, for strategy 2, the non-Markovian region in θ 1 -θ 2 plane does not affected by size of the environmental block. We will quantitatively investigate the non-Markovianties in both strategies.

Non-Markovianities in strategies 1 and 2
In this subsection, we will compare the non-Markovianties of both strategies 1 and 2 for, L B > 1, with the vacuum environmental state being vacuum. With the help of Eq. (23), we could compute the non-Markovianities for the channels with both strategies. In Fig. 4 we show the non-Markovianity as a function of L B with θ 1 /π = θ 2 /π = 1/6. The degrees of non-Markovianities of strategy 2 are always stronger than those of strategy 1. Moreover, in strategy 2, the non-Markovianity remains the same as that in the case of L B = 1, while, in strategy 1, the non-Markovianity decreases with the size of the block increasing and converges for L ≥ 8. Note that, in strategy 2, the successive S -B l−1 , B l−1 -B l and S -B l collisions construct an L B -level nested Mach-Zehnder interferometer of the system mode and L B (dissipative) environmental modes. Considering the normalization on the reflectivity of BS 1, i.e. the effective strength of the S -B l interactions for different L B are equal, we can conclude that the non-Markovianity in strategy 2 is independent of the block size.
In order to show the differences between the two strategies, we show the stroboscopic evolutions of the matrix element |c 1,1 (l)| and the nonzero eigenvalue of Λ l in Fig. 5. We see that in strategy 2 the revival of |c 1,1 (l)| is stronger than that in strategy 1. As a consequence, the negative eigenvalues contribute more to the indivisibility of the channel in strategy 2. It is easy to understand the advantages of non-Markovianity in strategy 2 in the limit of θ 2 = 0. In such a case, the time evolution of the system is unitary. Moreover the output ofâ l, j after S -B l collision are the input, with an additional π phase, ofâ l+1, j in S -B l+1 collision. This guarantees the time-reversal symmetry of the input and output of system states in two consecutive system-block collisions for strategy 2 and leads a strong non-Markovianity.

Pure damping process of a single-mode field
The CM with vacuum environmental state can be used to simulate the pure damping process of a single-mode field. The damping process of a single-mode field can be described by the Lindblad-type master equation, in the weak-coupling The contribution of the negative eigenvalues in strategy 2 is larger than that in strategy 1. This indicates that the Gaussian channel with strategy 2 violates the divisibility stronger than strategy 1. The parameters of BSs are θ 1 /π = 1/6 and θ 2 /π = 1/6. limit, whereâ is the annihilation operator,ρ(t) is the density operator of the field, g 1 is the coupling strength and γ(t) is the damping rate. We note that the evolution of a generic Gaussian state, governed by Eq. (25), can be described in terms of the covariance matrix, as shown in Eq. (15), with the matrices X(t) = exp [−Γ(t)/2]I 2 and Y(t) = {1 − exp [−Γ(t)]}I 2 /2 where Γ(t) = 2g t 0 γ(s)ds. The matrices X(t) and Y(t) coincide with X l and Y l in Eqs. (20) and (21), for the vacuum environment, through Above, we have set the elapsed time t = lτ where τ is the time interval between two successive system-environment collisions as mentioned before.
The non-Markovianity of the damping master equation, N PD , is measured, basing on the indivisibility of the dynamical map, through the time-dependent damping rate γ(t) [16,44] as where I are the intervals in which γ(t) < 0. It has been shown that N PD is proportional to the degree proposed by Rivas et al. which is measured by the increases in entanglement [13]. Eq. (27) indicates that the nonzero non-Markovianity originates from the negative γ(t) during the evolution. The correspondence between the damping rate in Eq. (25) and c 1,1 (l) in the stroboscopic CM is obtained as, via Eq. (26), .
Apparently, the necessary and sufficient condition of the non-Markovianity of the pure damping process, i.e. γ(t) < 0, is consistent with ours in the CM, i.e. |c 1,1 (l)| > |c 1,1 (l − 1)|. The contribution of the negative damping rate to the non-Markovian dynamics can be interpreted by the reverse quantum jump in the theory of non-Markovian quantum jump [45,46]. A quantum jump, occurring at positive γ(t), always interrupts the deterministic evolution, while the reverse jump, occurring at negative γ(t), will recover the coherence of the system of interest. In our CM, there is a similar process to the reverse jump in the non-Markovian evolution. Remind that the physical interpretation of |c 1,1 (l)| is the contribution of the input system mode to the output of the system mode. In the Markovian evolution, |c 1,1 (l)| decreases monotonically since the photons are always leaking. Contrastively, the nonmonotonic behavior of |c 1,1 (l)| means a photon reabsorption at some intermediate steps reminiscing the reverse jump.

C. Generic Gaussian environmental state
We now consider the case that the environmental state is a generic Gaussian state. By substituting Eq. (9) into Eq. (22), we obtain the non-Markovianity, N G (L), as the following, where n E is the thermal photon number and r E is the squeezing strength of the environmental states. One can see a generic Gaussian environmental states will enhance the non-Markovianity of vacuum environment and will not modify the boundary between Markovian and non-Markovian regions.

D. Entangled environmental state
In this subsection we will investigate effects of the entanglement embedded in the block on the non-Markovianity. We restrict our investigation to the case of L B = 2. The entanglement of the two-mode Gaussian state can be well characterized with the logarithmic negativity [47], which measures the entanglement by quantifying the violation of positive partial transpose separability criterion and has been proved to be a full entanglement monotone [48].
Let us consider that the two environment modes in the block are in a two-mode squeezed vacuum (TMSV) state with squeezing parameter ξ. A TMSV state |TMSV(ξ) l is generated from the vacuum via a two-mode squeezing operator, where the subscript l denotes the l-th block and |vac l stands for the vacuum state. Without loss of generality, we set ξ to be real, the characteristic function of Eq. (30) can be expressed as χ in l (ν l,1 , ν l,2 ) = exp − |ν l,1 | 2 + |ν l,2 | 2 2 cosh ξ The entanglement of the TMSV state measured by the logarithmic negativity is 2ξ [40]. Substituting Eq. (31) into Eq. (5) and following the procedures of computing Λ l , we can obtain the eigenvalues of Λ l as, where γ(l−1) = sinh (2ξ)c 2 1,1 (l−1) l l =1 c 1,2l (l)c 1,2l +1 (l) and γ(l) = sinh (2ξ)c 2 1,1 (l) l−1 l =1 c 1,2l (l − 1)c 1,2l +1 (l − 1) . The reduced state of each mode in a TMSV state is a thermal state ρ th with an effective particle number n th = sinh 2 (ξ). In order to investigate the effect of the entanglement embedded in the block, we compare the non-Markovianities of the channel with the states of the l-th block being the entangled state ρ l = |TMSV(ξ) l TMSV(ξ)| and the product state ρ l = ρ th ⊗ ρ th . We denote the non-Markovianities of each case as N TMSV and N prod , respectively, and the discrepancy δN = N TMSV − N prod .
In Fig. 6, we show δN for both strategies as functions of ξ with θ 1 /π = 0.2 and θ 2 /π = 0.1. For strategy 1, the non-Markovianity increases with ξ increasing. This indicates that the entanglement of the environmental particles in a block may enhance the non-Markovianity with the chosen parameters. In the inset of Fig. 6, we show δN as functions of θ 1 for different θ 2 with ξ = 1. Although the entanglement enhances, even maximally at some optimal θ 1 , the non-Markovianity for θ 2 /π = 0.1, 0.2 and 0.3, it does not affect the degree of non-Markovianity for θ 2 /π = 0.4. Whether the entanglement will affect the non-Markovianity depends on the intrinsic properties of the BSs. For strategy 2, the value of δN is always zero and irrelevant to ξ indicating that the entanglement of environment particle does not affect the degree of non-Markovianity.

IV. CONCLUSIONS
We have presented an extended CM to simulate the non-Markovian dynamics of a quantum system. In such a CM, the unit to represent the environment is a block consisted of a number of particles. The introduced environmental block enables us to study the non-Markovianity of a quantum channel through different strategies of the system-environmental interactions and states of the environmental units.
In our CM, the system-environment (S -B l ) collisions are implemented in two strategies: in strategy 1, the system mode S sequentially interacts with the environmental modes E l, j in the ascending order of j for all l; in strategy 2, the system mode interacts with the environmental modes E l, j in the ascending order of j for odd l and in the descending order of j for even l. We have adopted an all-optical system to implement the modified CM. By restricting the input modes to be Gaussian and the interactions to be linear, the dynamics of the system can be described via a Gaussian channel. With the help of the measure of non-Markovianity based on the indivisibility of dynamical maps, we have studied the effects of both strategies on the dynamics of the system mode. In strategy 1, it is shown that the non-Markoviantiy will be suppressed and converge with the size of block increasing. While in strategy 2, the non-Markovianity is independent on the size of the block.
We have also presented a necessary and sufficient condition of the non-Markovianity of the Gaussian channel. The physics behind the condition is that the contribution of the input system mode to the output of the system is nonmonotonic during the stroboscopic evolution, i.e. |c 1,1 (l)| > |c 1,1 (l − 1)| for some intermediate l. Such a process is similar to the reverse jump in the theory of non-Markovian quantum jump. Our measure of non-Markovianity is based on quantifying the extent by which the intermediate process fails to be CP. This corresponds to the quantification of the negative eigenvalues of the symmetric matrix associated with the intermediate process Φ l,l−1 . This measure coincides with other existing criterions, e.g. the one based on the quantifying the negative decoherence rate of the master equation in the canonical form [44], in detecting the non-Markovian features. However, since based on different point of views the existing measures may not agree with each other in quantifying the non-Markovianity of some specific channels, for instance the Gaussian channel with thermal environment [24]. It would be interesting to investigate the connections of our measure to other ones in the future work.
We have found that the generic Gaussian environment states with nonzero temperature and squeezing will quantitatively enhance the non-Markovianity of the channel with vacuum state. We have also investigated the effects of the entanglement embedded in the environmental block on the non-Markovianity. By comparing non-Markovianity in the cases of the environmental block being in TMSV state and the product state of the corresponding reduced (thermal) states, we found that, in strategy 1, if the entanglement will enhance the non-Markovianity depends on the intrinsic properties of the channel, i.e., the reflectivity and transmissivity of the BSs. However, in strategy 2, the entanglement does not play roles in the non-Markovianity.
We emphasize that, the environment, which is in permanent contact with the system in a realistic process, is modeled by an ensemble of identical blocks in the CM. Thus we can simulate various dynamics of the open system subjected to different reservoirs by specifying the states of system and environmental blocks. For instance, apart from the Gaussian channel, we can set the environment to be vacuum and at most one excitation in the system mode to simulate the qubit amplitudedamping channel [49].
Finally we would like to briefly discuss the possible experimental realization of our model. It could be implemented in the advanced integrated photonic quantum simulator [50][51][52][53][54]. Such a platform has the advantages of intrinsic phase stability, arbitrary control of the reflectivity (and transmissivity), and flexible scalability. The integrated photonic simulator has been used to observe the Anderson localization in disordered quantum walk composed of eight steps [55]. Such a scale of the concatenated interferometers is capable to witness the effects of the interaction strategies and entanglement on the non-Markovianities with L B = 2 and L = 4 in our model. (A2) The output joint characteristic function is calculated by