Subsystem complexity after a local quantum quench

We study the temporal evolution of the circuit complexity after the local quench where two harmonic chains are suddenly joined, choosing the initial state as the reference state. We discuss numerical results for the complexity for the entire chain and the subsystem complexity for a block of consecutive sites, obtained by exploiting the Fisher information geometry of the covariance matrices. The qualitative behaviour of the temporal evolutions of the subsystem complexity depends on whether the joining point is inside the subsystem. The revivals and a logarithmic growth observed during these temporal evolutions are discussed. When the joining point is outside the subsystem, the temporal evolutions of the subsystem complexity and of the corresponding entanglement entropy are qualitatively similar.


Introduction
A quantum circuit constructs a target state from a given reference state through a sequence of gates chosen within a set of allowed gates. The circuit complexity has been introduced in quantum information theory as the minimum number of allowed gates employed to construct the circuit [1][2][3][4][5]. During the past few years the circuit complexity has been investigated in the context of quantum gravity through the gauge/gravity (holographic) correspondence [6][7][8][9][10][11][12][13][14][15]. It is worth exploring the circuit complexity also in quantum field theory and in quantum many-body systems.
Within the class given by the quantum many-body systems, it is natural to start with quantum circuits made by Gaussian states in free systems [16,17]. The complexity of these circuits when only pure states are allowed has been investigated, obtaining explicit expressions both for bosonic and for fermionic lattices [18][19][20][21][22][23][24][25][26].

JHEP08(2021)135
It is important to quantify the complexity also for quantum circuits made by mixed states, which are characterised by density matrices [27][28][29][30][31][32]. The reduced density matrices provide an important class of mixed states that are crucial to study the bipartite entanglement. Given the spatial bipartition A ∪ B of a quantum system in a state characterised by the density matrix ρ whose Hilbert space can be factorised accordingly as H = H A ⊗ H B , the reduced density matrix ρ A ≡ Tr H B ρ of the spatial subsystem A (normalised by Tr H A ρ A = 1) characterises a mixed state. An important quantity to consider is the entanglement entropy S A = −Tr(ρ A log ρ A ). When the entire system is in a pure state, S A = S B (see [33][34][35][36][37][38][39] for reviews). The subsystem complexity C A is defined as the complexity of a circuit where both the reference state and the target state are reduced density matrices associated to the same spatial subsystem A.
It is insightful to compare C A and S A . At equilibrium, this comparison has been discussed in lattice models [29,30] and in various gravitational backgrounds within the holographic correspondence [11,15,28,[40][41][42], finding that a major distinction occurs in the leading divergence: while for S A it is determined by the volume of the boundary of A (this law is violated in a 2D conformal field theory in its ground state, e.g. when A is an interval and therefore the boundary of A is made by two points [43][44][45]), for C A it grows like the volume of A, with a power that depends on the choice of the cost function [29,30].
Quantum quenches are interesting protocols to study the dynamics of isolated quantum systems out of equilibrium (see [46,47] for reviews). Given a system prepared in a state |ψ 0 , consider a sudden change at t = 0 that provides the time-evolved state |ψ(t) = e −iHt |ψ 0 for t > 0. Since typically |ψ 0 is not an eigenstate of the evolution Hamiltonian H, this timeevolved pure state is highly non trivial. The kind of sudden change leads to identify two main classes of quantum quenches. Global quenches are characterised by sudden changes that involve the entire system (e.g. a modification of a parameter in the Hamiltonian) [48][49][50]. Instead, in local quenches the sudden change occurs only at a point. For instance, local quenches where either two systems are joined together [51,52] or a local operator is inserted at some point [53,54] have been explored. The temporal evolutions of the entanglement entropy S A after various quantum quenches (either global or local) have been widely studied during the past few years [55][56][57][58][59][60][61]. For systems in a finite volume, revivals occur in some temporal evolutions [58,62,63].
It is worth investigating the temporal evolutions of the circuit complexity for the entire system and of the subsystem complexity after different quantum quenches. For some global quenches, holographic prescriptions have been employed to determine numerically the temporal evolutions of the complexity for the entire system [64][65][66] and of the subsystem complexity [67][68][69][70]. In free lattice models, the temporal evolutions of the complexity after some global quenches have been studied, both for the entire system [71][72][73][74][75][76] and for subsystems [72,76,77]. The temporal evolution of the holographic entanglement entropy after a local quench has been explored in [78][79][80][81][82][83][84]. For the local quench corresponding to an operator insertion the temporal evolution of the holographic subsystem complexity has been considered [85,86], while for the local quench where two systems are joined together only the holographic entanglement entropy has been studied [79,83].

JHEP08(2021)135
< l a t e x i t s h a 1 _ b a s e 6 4 = " P 4 Q 3 X N 9 u Q v 9 Y B I 2 x 1 w U k U 6 b x q u 0 = " > A A A B 6 H i c b Z C 7 S g N B F I b P x l t c b 1 F L m 8 E g W I X d F G o j B m 2 s J A F z g W Q J s 5 O z y Z j Z C z O z Q g h 5 A h s L R W z 1 Y e x t x L d x c i k 0 8 Y e B j / 8 / h z n n + I n g S j v O t 5 V Z W l 5 Z X c u u 2 x u b W 9 s 7 u d 2 9 m o p T y b D K Y h H L h k 8 V C h 5 h V X M t s J F I p K E v s O 7 3 r 8 Z 5 / R 6 l 4 n F 0 q w c J e i H t R j z g j G p j V W 7 a u b x T c C Y i i + D O I H / x Y Z 8 n 7 1 9 2 u Z 3 7

b H V i l o Y Y a S a o U k 3 X S b Q 3 p F J z J n B k t 1 K F C W V 9 2 s W m w Y i G q L z h Z N A R O T J O h w S x N C / S Z O L + 7 h j S U K l B 6 J v K k O q e m s / G 5 n 9 Z M 9 X B m T f k U Z J q j N j 0 o y A V R M d k v D X p c I l M i 4 E B y i Q 3 s x L W o 5 I y b W 5 j m y O 4 8 y s v Q q 1 Y c E 8 K x Y q T L 1 3 C V F k 4 g E M 4 B h d O o Q T X U I Y q M E B 4 g C d 4 t u 6
s R + v F e p 2 W Z q x Z z z 7 8 k f X 2 A w h x k B Y = < / l a t e x i t > N < l a t e x i t s h a 1 _ b a s e 6 4 = " I z p x t 4 c l E C 9 1 r A v / Z H b u z S l N Q u M = " > A A A B 6 n i c b Z D L S g M x F I b P e K 3 j r e r S T b A I r s p M F + p G L L p x J R X t B d q h Z N J M G 5 r J h C Q j l K G P 4 M a F I i 7 1 X d y 7 E d / G 9 L L Q 1 h 8 C H / 9 / D j n n h J I z b T z v 2 1 l Y X F p e W c 2 t u e s b m 1 v b + Z 3 d m k 5 S R W i V J D x R j R B r y p m g V c M M p w 2 p K I 5 D T u t h / 3 K U 1 + + p 0 i w R d 2 Y g a R D j r m A R I 9 h Y 6 1 a h 6

a + 4 B W 9 s d A 8 + F M o n H + 4 Z / L t y 6 2 0 8 5 + t T k L S m A p D O N a 6 6 X v S B B l W h h F O h 2 4 r 1 V R i 0 s d d 2 r Q o c E x 1 k I 1 H H a J D 6 R Q l C j 7 h E F j 9 d H h m O t B F o K 2 N s e n o 2 G 5 n / Z c U R K d B x o R M D R V k 8 l G U c m Q S N N o b d Z i i x P C B B U w U s 7 M i 0 s M K E 2 O v 4 9 o j + L M r z 0 O t V P S P i 6 U b r 1 C + g I l y s A 8 H c A Q + n E A Z r q A C V S D Q h Q d 4 g m e H O 4 / O i / M 6 K V 1 w p j 1 7 8 E f O + w 8 w l p C 8 < / l a t e x i t >
rN < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 c O v f N b T p x 2 4 x w V n 9 W j 9 g U u b E L E = " > A A A B 6 H i c b Z C 7 S g N B F I b P x l t c b 1 F L m 8 E g W I X d F G o j B m 0 s L B I w F 0 i W M D s 5 m 4 y Z v T A z K 4 S Q J 7 C x U M R W H 8 b e R n w b J 5 d C E 3 8 Y + P j / c 5 h z j p 8 I r r T j f F u Z p e W V 1 b X s u r 2 x u b W 9 k 9 v d q 6 k 4 l Q y r L B a x b P h U o e A R V j X X A h u J R B r 6 A u t + / 2 q c 1 + 9 R K h 5 H t 3 q Q o B f S b s Q D z q g 2 V u W m n c s 7 B W c i s g j u D P I X H / Z 5 8 v 5 l l 9 u 5 z 1 Y n Z m m I k W a C K t V 0 n U R 7 Q y o 1 Z w J H d i t V m F D W p 1 1 s G o x o i M o b T g Y d k S P j d E g Q S / M i T S b u 7 4 4 h D Z U a h L 6 p D K n u q f l s b P 6 X N V M d n H l D H i W p x o h N P w p S Q X R M x l u T D p f I t B g Y o E x y M y t h P S o p 0 + Y 2 t j m C O 7 / y I t S K B f e k U K w 4 + d I l T J W F A z i E Y 3 D h F E p w D W W o A g O E B 3 i C Z + v O e r R e r N d p a c a a 9 e z D H 1 l v P w V p k B Q = < / l a t e x i t > L < l a t e x i t s h a 1 _ b a s e 6 4 = " 4 c O v f N b T p x 2 4 x w V n 9 W j 9 g U u b E L E = " > A A A B 6 H i c b Z C 7 S g N B F I b P x l t c b 1 F L m 8 E g W I X d F G o j B m 0 s L B I w F 0 i W M D s 5 m 4 y Z v T A z K 4 S Q J 7 C x U M R W H 8 b e R n w b J 5 d C E 3 8 Y + P j / c 5 h z j p 8 I r r T j f F u Z p e W V 1 b X s u r 2 x u b W 9 k 9 v d q 6 k 4 l Q y r L B a x b P h U o e A R V j X X A h u J R B r 6 A u t + / 2 q c 1 + 9 R K h 5 H t 3 q Q o B f S b s Q D z q g 2 V u W m n c s 7 B W c i s g j u D P I X H / Z 5 8 v 5 l l 9 u 5 z 1 Y n Z m m I k W a C K t V 0 n U R In this manuscript, we are interested in the temporal evolution of the circuit complexity after a local quench. We consider the local quench described by Eisler and Zimborás in [87], where two harmonic chains containing rN and (1 − r)N sites are joined at t = 0 (here r is a rational number 0 r 1), as shown in the top panel in figure 1. We focus on circuits made only by Gaussian states. First we study the temporal evolution of the complexity for the entire chain; then we investigate the temporal evolution of the subsystem complexity C A for the subsystem A given by a block of L consecutive sites (the spatial bipartitions considered in the manuscript are shown in the bottom panels of figure 1), by employing the complexity for mixed states based on the Fisher information geometry [30].
The outline of this manuscript is as follows. In section 2 we describe the local quench protocol, introducing the covariance matrices characterising the states involved in the construction of the optimal circuit. In section 3 we evaluate numerically the temporal evolution of the circuit complexity for the entire chain, choosing the (pure) initial state at t = 0 as the reference state and the (pure) state at generic time t > 0 after the local quench as the target state. In section 4 we discuss our numerical analysis of the temporal evolutions of the subsystem complexity for a block of consecutive sites. Some conclusions and open questions are drawn in section 5. The appendices A, B and C contain further technical details and supplementary results.

Covariance matrix after the quench
The Hamiltonian of the harmonic chain made by N sites (we set = 1) reads where the position and the momentum operatorsq i andp i are hermitean operators satisfying the canonical commutation relations [q i ,q j ] = [p i ,p j ] = 0 and [q i ,p j ] = iδ i,j . The Dirichlet boundary conditions (DBC)q 0 =q N +1 = 0 andp 0 = 0 are imposed at the endpoints.

JHEP08(2021)135
The initial state is given by the following pure state |Ψ (l,r) 0 ≡ |ψ l ⊗ |ψ r (2.2) where |ψ l is the ground state of the Hamiltonian H l , defined by (2.1) for the sites labelled by 0 i N l , with the physical parameters m l , ω l and κ l and DBC imposed at i = 0 and i = N l + 1. Similarly, |ψ r is the ground state of the Hamiltonian H r in (2.1) for the chain made by the sites labelled by N l i N l + N r + 1, with parameters m r , ω r and κ r and DBC imposed at i = N l and i = N l + N r + 1. Thus, the initial state (2.2) depends on N l , m l , ω l , κ l , N r , m r , ω r and κ r . The total number of sites is N ≡ N l + N r . Equivalently, we can describe the initial state in terms of N and of the position parameter 0 r 1 (in unit of N ), that determines the separation between the left and the right chain; indeed N l = N r and N r = N (1 − r).
Given the state (2.2) at t = 0, the time evolved state at t > 0 through (2.1) is This setup describes different quantum quenches. A global quench can be obtained by setting N l = N and N r = 0 (or viceversa, equivalently), κ l = κ r ≡ κ 0 , m l = m r ≡ m 0 and ω l = ω r ≡ ω 0 . In this case the initial state is the ground state of a single chain made by N sites. If κ 0 = κ, m 0 = m and ω 0 = ω, the global quench involves all the parameters occurring in the Hamiltonian (2.1). An important special case is the global quench of the frequency parameter [49] discussed in appendix A.1. In appendix A.2 we consider the global quench of the spring constant and of the frequency, which corresponds to m 0 = m, κ 0 = κ and ω 0 = ω.
In this manuscript we consider the local quench described in [87], where two disconnected harmonic chains, containing N l = rN and N r = (1 − r)N sites, are joined at t = 0, as represented pictorially in the top panel of figure 1. This quench protocol corresponds to m r = m l = m, κr = κl = κ and ωr = ω l = ω.
The bosonic Gaussian states in harmonic lattices are fully characterised by their covariance matrix [16,17]. In the quench protocol that we are considering [87], the initial state (2.2) is Gaussian and the time evolution generated by (2.1) preserves its Gaussian nature; hence (2.3) is Gaussian too, for any t > 0. The time evolved state (2.3) is completely characterised by the 2N × 2N covariance matrix γ (l,r) (t), whose generic element is defined as The covariance matrix of the initial state (2.2) at t = 0 reads [87] γ (l,r) where the superscript indicates that this covariance matrix corresponds to an initial configuration made by two disjoint chains, containing N l and N r sites respectively, and

JHEP08(2021)135
with T 0 being the following diagonal matrix written in terms of the dispersion relations The matrix V 0 in (2.5) is block diagonal too and it can be written as where the elements of the V N l and V Nr are given by Notice that, since the matrices V Ns defined by (2.10) are orthogonal, V 0 is symplectic and orthogonal. Thus, also V 0 in (2.9) is symplectic and orthogonal.
We find it worth remarking that, in the expression (2.5) for γ (l,r) 0 , the parameters m l , ω l , κ l , m r , ω r and κ r occur in Q 0 and P 0 , while V 0 depends only on N l and N r .
The covariance matrix (2.4) of the pure state (2.3) at any t > 0 after the quench is written in terms of the covariance matrix of the initial state (2.5) as [87] where the time dependence occurs only through the matrix E(t).
A convenient block decomposition of (2.11), which is employed in section 4.1, reads where Q(t), P (t) and M (t) are the N ×N correlation matrices whose elements are |ψ 0 respectively. All these three matrices provide a non trivial temporal dependence. An insightful decomposition for the matrix E(t) in (2.11) is [87] in terms of the matrix V N , whose generic element is given by (2.10) with N s replaced by N , and of the matrix E(t), whose block decomposition reads with D, A and B being diagonal matrices whose elements are respectively

JHEP08(2021)135
where Ω k is the dispersion relation given by (2. 16) Since V N is orthogonal, V is symplectic and orthogonal. This observation and the fact that E(t = 0) = 1 lead to E(t = 0) = 1; hence, from (2.11), we have that γ (l,r) (t = 0) = γ (l,r) 0 , as expected. Notice that, by using (2.15), one finds that E(t) in (2.14) is symplectic; hence, since V is symplectic too, we conclude that E(t) in (2.13) is symplectic. Thus, E(t) implements on the initial covariance matrix the unitary transformation on the initial state given in (2.3).
In order to investigate the circuit complexity, we find it worth employing also the Williamson's decompositions [16,88] of the covariance matrices of the reference and of the target states.
The Williamson's decomposition of the initial covariance matrix γ (l,r) 0 in (2.5) reads where the symplectic matrix V 0 has been defined in (2.9). Since the initial state (2.2) is pure, all the symplectic eigenvalues of its covariance matrix γ (l,r) 0 are identical and equal to 1/2. Notice that X 2 0 = Γ 0 , where Γ 0 has been introduced in (2.5). Plugging the Williamson's decomposition (2.17) into (2.11), it is straightforward to obtain the Williamson's decomposition of the covariance matrix γ (l,r) (t) at any t > 0, which characterises the pure state (2.3). It reads where the matrices E(t) and W 0 have been defined in (2.13) and (2.17) respectively. Since both W 0 and E(t) are symplectic matrices, the matrix W (t) is symplectic too.

Complexity for the harmonic chain
In this section we discuss the temporal evolution of the complexity for the entire chain after the local quench defined in section 2; hence both the reference and the target states are pure. We focus on the simplified setup where the quantum circuits are made only by bosonic Gaussian states with vanishing first moments.

Optimal circuit and complexity
The reference and the target states are fully characterised by their covariance matrices, which are γ R and γ T respectively. The circuit complexity obtained from the Fisher-Rao distance between γ R and γ T reads [89,90] C ≡

JHEP08(2021)135
This complexity, which corresponds to the F 2 cost function, has been studied for both pure states [18,19] and mixed states [30]. The optimal circuit that allows to construct γ T from γ R is made by the following sequence of covariance matrices [90] Denoting by t R and t T the values of time t corresponding to the reference and to the target states respectively, for their covariance matrices we have In the most general setup, these matrices depend on the sets of parameters given by Y S ≡ {m l,S , κ l,S , ω l,S , m r,S , κ r,S , ω r,S , m S , κ S , ω S }, with S = R and S = T for the reference and the target state respectively. The corresponding states can be interpreted as the states obtained through the time evolutions at t = t R 0 and t = t T t R respectively, through two different quenches determined by the parameters Y R and Y T respectively, as described in section 2. The circuit complexity (3.1) can be evaluated by finding the eigenvalues of γ T γ −1 R . From (2.5), (2.11) and (2.13) for the reference and the target states (where V is or- This expression is difficult to deal with mainly because of V V t 0 , which encodes the spatial geometries before and after the local quench. Notice that, when t R = t T = 0, we have E R = E T = 1; hence (3.4) simplifies to Since V 0 defined in (2.9) is orthogonal and the matrices Γ 0,T and Γ 0,R in (2.6) are diagonal, the eigenvalues of γ T γ −1 R in (3.5) are the ratios of the entries of Γ 0,T and Γ 0,R . By employing this observation and (3.1), we obtain the following expression for the circuit complexity in terms of the dispersion relations (2.8). In the special case where m l,R = m l,T , m r,R = m r,T and either N l = 0 or N r = 0, the expression (3.6) becomes the result obtained in [18]. The above results can be employed to study the temporal evolution of the complexity after a global quench, as already mentioned in section 2. In appendix A.1 we discuss the case of the quench of the mass parameter, showing that the analysis of [87] allows to recover the correlators obtained in [50], which have been employed to evaluate the temporal evolutions JHEP08(2021)135 both of the complexity of the entire chain [73,77] and of the subsystem complexity [77] after this kind of global quench. Instead, in appendix A.2 the temporal evolution of the complexity after the global quench of the spring constant is mainly considered, with the initial state (which is also the reference state) given by the unentangled product state (i.e. the ground state of the hamiltonian (2.1) with a certain frequency and vanishing spring constant), that has been adopted as the reference state in various studies about complexity [18-20, 29, 71].
In this manuscript we are interested in circuits whose reference and the target states are pure states along the time evolution of a given local quench at different times t R and t T . This can be done by choosing the parameters introduced in section 2 as follows From (2.6), we have that Γ 0,T and Γ 0,R do not depend on time and the setting given by (3.7) leads to Γ 0,T = Γ 0,R ≡ Γ 0 . Thus, (3.4) simplifies to In our analysis we mainly consider the initial state (2.2) as the reference state. This choice corresponds to set t R = 0 in (3.3) and γ R = γ (l,r) 0 given by (2.5). In this case, from (2.14) and (2.15), one finds that E R = E(t = 0) = 1 and that (3.8) simplifies to where the notation E T = E has been introduced. Finding the eigenvalues of (3.9) analytically is complicated; hence we study them numerically. When both the reference and the target states are pure, the Williamson's decompositions of their covariance matrices (3.3) read respectively where W R and W T are symplectic matrices. By introducing the following symplectic matrix it is straightforward to realise that the complexity (3.1) becomes [19] In the case where the reference and the target states are pure states along the time evolution of a given local quench and, furthermore, the initial state is chosen as the reference state, we have that W R = W 0 and W T = W (t), where W 0 and W (t) are defined in (2.17) JHEP08 (2021)135 and (2.18) respectively. From the expressions of W 0 and W (t) and the fact that V 0 in (2.5) is orthogonal, for (3.11) one obtains By using that X 2 0 = Γ 0 , where Γ 0 is given in (2.5), we find that 14) The diagonalisation of this matrix is as difficult as the one of (3.9). However, this form could be helpful in future attempts to obtain analytic results for the complexity (3.12). The Euler decomposition (also known as Bloch-Messiah decomposition) of a symplectic matrix S reads [91] where the diagonal matrix Λ = diag(Λ 1 , . . . , Λ N ) contains the squeezing parameters Λ j 0. From (3.12), it is straightforward to realise that the complexity of circuits made by pure states can be written in terms of the squeezing parameters (Λ TR ) j corresponding to the symplectic matrix W TR as follows [19] (3.16) The symplectic matrix E(t) t can be decomposed into four N × N blocks which are diagonal matrices (see (2.14)); hence we can find its Euler decomposition E(t) t = L E X E R E (where all the three matrices can depend on t) by following the procedure discussed in appendix B. Plugging this decomposition into (3.13), one obtains This expression does not provide the Euler decomposition of W TR because of the occurrence of the diagonal matrix X 0 , which is not orthogonal. By contradiction, if X 0 were orthogonal, (3.17) would be the Euler decomposition of W TR with the squeezing parameters given by X E because L E , R E , V and V 0 are symplectic and orthogonal matrices. This would lead to a complexity (3.16) independent of the position of the joining point because X E depends only on E in (2.14), which is determined by the parameters characterising the evolution Hamiltonian. The numerical analysis performed in section 3.3 shows that this is not the case.

Initial growth
It is worth exploring the leading term of the initial growth of the temporal evolution of the complexity (3.1) for the entire chain when the reference state is the initial state (i.e. γ R = γ (l,r) 0 in (2.5)) and the target state is the state at time t after the local quench that we are exploring (i.e. γ T = γ (l,r) (t) in (2.11)).

JHEP08(2021)135
By expanding E(t) in (2.14) as t → 0, one finds where Ω k is given in (2.16). By employing (3.18) in (2.13) and the fact that V N is orthogonal, we obtain By using the expansion (3.19), for the covariance matrix (2.11) we find where the O(t) term is symmetric, as expected. This straightforwardly leads to This expansion provides the following linear growth for the complexity (3.1) where for the coefficient c 1 we find (see the appendix C for its derivation) Simplifying further the last term in this expression is complicated, hence we evaluate it numerically, as done in the bottom panel of figure 4 to determine the dashed straight line. As a consistency check for (3.23), let us consider the trivial case where the quench does not occur, which corresponds to set N l = N and N r = 0 (or viceversa), implying that V 0 = V N . By using that V N is orthogonal, (2.6) and (2.7), one finds that the last term in (3.23) simplifies to N k=1 Ω 2 k . Then, since N r = 0, the second sum in (3.23) does not occur and therefore c 1 = 0, as expected, consistently with the fact that the initial state does not evolve.

Numerical results
In this section we discuss some temporal evolutions of the complexity for the entire chain after a local quench where two chains are joined (see section 2), evaluated numerically through (3.1). The reference and the target states are respectively the initial state (t R = 0) and the pure state corresponding to a generic value of t T ≡ t 0 along the evolution after the quench. The parameters of this quench protocol are set as in (3.7). The data points reported in all the figures shown in the main text have been obtained for m = 1 and κ = 1.
In figure 2 and figure 3, the temporal evolutions corresponding to critical Hamiltonians are considered; i.e. ω = 0. Since the volume is kept finite, revivals are observed, as already discussed for the temporal evolutions of other quantities [62]. The different cycles correspond to p < t/(2N +2) < p+1, with p being a non-negative integer. This approximate periodic behaviour is observed also in the correlators providing the covariance matrix. For instance, in figure 2 the cycles corresponding to p = 0 and p = 1 are displayed.
Within each cycle we can identify three temporal regimes: (I) p < t/(2N + 2) < p + r, characterised by an initial growth, a local maximum and a subsequent decrease; (II) p+r < t/(2N + 2) < p + 1 − r, where the evolution is almost stationary (a slight convexity of the curves is observed by zooming in); (III) p + 1 − r < t/(2N + 2) < p + 1, characterised by a growth until a local maximum is reached and a subsequent decrease. The last regime is very similar the first one, after a time reversal; indeed, the curve of C(t) within each cycle remains roughly invariant after a reflection with respect to the value of t corresponding to the center of the cycle.
In the special case of r = 1/2 (see the top left panel of figure 2 and the black symbols in figure 3), the second regime does not occur; hence the cycles correspond to p < t/(N +1) < p + 1, with p being a non-negative integer. In the top panels of figure 2, the temporal evolutions of C − 1 5 log(N + 1) are displayed for r = 1/2 (left panel) and r = 1/4 (right panel). When N is large enough, the data for different values of N nicely collapse, except for the beginning and the end of each cycle, as discussed below. In the top left panel of figure 2, where r = 1/2, also the following curve is shown which nicely agrees with the data points in the middle of each cycle p < t/(N + 1) < p + 1, when N is large enough.
In the bottom panels of figure 2, we consider the temporal regime of initial growth for C subsequent to the early linear growth (3.22). We find that the data corresponding to different values of N nicely collapse on the curve given by with const 0.5346 within a temporal regime whose width increases with N . Notice that (3.25) does not correspond to the leading term of (3.24) when t/(N + 1) → 0 because the coefficients multiplying the logarithms are different. This is consistent with the fact that the data in the top panels of figure 2 do not collapse at the beginning and at the end of each cycle. Taking t/(2N + 2) instead of t as the independent variable in the bottom panels JHEP08(2021)135 of figure 2, the data collapse is observed for C − 1 4 log(N + 1) and not for C − 1 5 log(N + 1), which is plotted in the top panels of the same figure. This is consistent with the data corresponding to the black symbols in the right panels figures 6 and 8 and in the left panels of figure 10, which describe the complexity of the entire chain. Let us anticipate that also for the subsystem complexity C A different temporal regimes occur where the data points for C A − α log(N + 1) corresponding to increasing values of N collapse, with different values of α in the different regimes (see section 4). By comparing the two bottom panels in figure 2, we observe that the initial growth of the complexity is independent of the value of r.
In figure 3 we consider a longer range of t, in order to include more cycles and to highlight the fact that the approximate periodicity persists, for various values of r. Notice that the values of the local maximum within each cycle increases with r until r = 1/2. Furthermore, the height of the plateaux characterising the second temporal regime within each cycle grows with r until certain value r * (from figure 3, we have 0.25 < r * < 0.35), then it decreases. Instead, the duration of this plateaux is always decreasing for 0 < r 1/2 and vanishes at r = 1/2. The symmetry of the problem straightforwardly leads to realise that the temporal evolution of the complexity for a given r is equal to the one corresponding to 1 − r, for the same choice of all the other parameters. We have obtained numerical data for the temporal evolutions of the complexity displayed in figure 3 also for N = 200, finding that the data points of C − 1 5 log(N + 1) for N = 100 and N = 200 approximatively collapse (see also the top panels of figure 2). In figure 3 we have reported only the numerical curves for N = 100 in order to display in a clear way the qualitative changes in the temporal evolutions corresponding to different r.
Some temporal evolutions of the complexity determined by gapped Hamiltonians after the local quench are shown in figure 4, where the different coloured curves correspond to different values of ωN 50.
In the top panels of figure 4, we show that the curves for C − 1 4 log(N + 1) corresponding to different values of N collapse. We remind that this collapse has been observed for C − 1 5 log(N + 1) when ω = 0 (see the top panels of figure 2). It would be interesting to understand this numerical observation. Furthermore, by comparing the temporal evolutions in the top panels of figure 4 with the ones in figure 3, we notice that the periodicity highlighted for ω = 0 does not occur when ωN > 0 in general.
When ωN 1, the initial part of the temporal evolution is similar to the one observed in the case of ω = 0 (see figure 2), as one realises from the curves corresponding to ωN = 1 in the top panels of figure 4. For large values of ωN 10, the temporal evolution of C − 1 4 log(N + 1) is roughly described by a complicated oscillation about a constant value. This constant value decreases with ωN and, when ωN is large enough, is independent of r. Also the amplitude of the oscillations about this constant value decreases as ωN increases.
The bottom panel of figure 4 focuses on the initial growth of C; hence it is instructive to compare it against the bottom panels of figure 2 where ω = 0. In the temporal regime considered in the bottom panel of figure  -0.5 We find it worth mentioning some results about the temporal evolution of the complexity evaluated within the gauge/gravity correspondence.
The temporal evolution of the holographic complexity in the Vaidya gravitational spacetimes, which model the formation of a black hole through the collapse of a shell and have been exploited to study the gravitational duals of global quenches [92,93], has been studied in [7-9, 12, 13, 64-66]. Qualitative comparisons between these results and the temporal evolution of the complexity in harmonic chains have been discussed in [71,72,77].
A gravitational background dual to the local quench obtained through the insertion of a local operator [53,54] has been proposed in [78]. The temporal evolution of the holographic complexity in this spacetime has been studied in [85,86]. However, this local quench is very different from the one considered in this manuscript, where two systems initially disconnected are glued together at some point. A gravitational dual for this local quench has been studied e.g. in [79,80,83] by employing the AdS/BCFT setup discussed in [94,95]. It would be interesting to investigate the temporal evolution of the holographic complexity for the entire system in this spacetime.

Evolution Hamiltonians made by two sites
In the simplest case, two separate systems containing only one site are joined at t = 0; hence N l = N r = 1 and N = 2, i.e. the evolution Hamiltonian describes two sites.
In order to specialise (3.9) to this case, one first observes that (2.10) gives which (from (2.9) and (2.13)) provide respectively V 0 = 1 and that is symmetric and orthogonal. Since V 0 = 1, in this case (3.9) simplifies to and, by using (2.14), we have being the diagonal matrices whose elements are given by (2.15) with N = 2. From (3.28) and (3.30) one observes that the structure of γ T γ −1 R is not very easy already in this simple case of N = 2. From (3.29) and (3.30), one obtains an explicit expression for γ T γ −1 R in (3.28), which is not reported here because we find it not very insightful. Its eigenvalues are g TR,1 , g −1 TR,1 , g TR,2 and g −1 TR,2 , in terms of where Ω k ≡ ω 2 + (2k − 1) κ m , with k = 1, 2. Notice that (3.31) depends only on the two dimensionless parameters ω ≡ ω/ κ/m and ωt. For any given k, the oscillatory behaviour is governed by the frequency π/ Ω k . This result is qualitatively similar from the one obtained for the temporal evolution of the complexity after the global quench of the mass; indeed, also in that case the temporal evolution associated to each mode is determined by a single JHEP08(2021)135 frequency (see eq. (3.7) of [77]). When the evolution is critical, (3.31) is written through (3.32) By employing the above expressions for the eigenvalues of γ T γ −1 R , we find that the complexity (3.1) in this case can be written as follows in terms of the expressions in (3.31). We find it worth investigating the asymptotic regime given byω ≡ ω/ κ/m → ∞, while ωt is kept fixed and finite; hence the condition ωt ω is imposed. In this regime, the expansion of (3.31) reads (3.35) and we find it worth remarking that b 1 (t) is independent of k. By employing (3.34), it is straightforward to obtain the first terms in the expansion of the complexity (3.33) in this asymptotic regime. The leading term reads In figure 5 we show some temporal evolutions for the complexity (3.33), where the evolution Hamiltonian is made by two sites.

JHEP08(2021)135
The expression (3.31) tells us that g TR,k with k = 1, 2, which provide the complexity (3.33), are oscillating functions whose periods are π/ Ω k . A straightforward numerical inspection shows that g TR,1 > g TR,2 in the whole ranges ofω and of ωt. This leads us to plot the complexity in terms of Ω 1 t in the left panel of figure 5, whereω is not too large (in this case this argument does not apply because of (3.34)). Interestingly, we observe that the local extrema of the curves having different ω occur approximatively at the same values of Ω 1 t.
Notice that the temporal evolution corresponding to ω = 0 in the left panel of figure 5 is very different from the one displayed in the top left panel of figure 2, obtained for N 1. In the right panel of figure 5, the data points, obtained numerically for N = 10 and N = 20 and three values of ω, are compared against the leading term given in (3.36). We find it worth remarking that (3.36) nicely agrees with the temporal evolution of the complexity for small values of t, even when N > 2. The agreement between the analytic curve and the data improves as ω grows, as expected.

Subsystem complexity
In this section we investigate the temporal evolution of the subsystem complexity C A after the local quench introduced in section 2, when the reference and the target states are the reduced density matrices of the block A in the configurations shown in the bottom panels of figure 1.

Optimal circuit and subsystem complexity
In the harmonic lattices in the pure states that we are considering, the reduced density matrix associated to a spatial subsystem A characterises a mixed Gaussian state which can be fully described through its reduced covariance matrix γ A [16,33,96], defined as the 2L × 2L real, symmetric and positive definite matrix (L denotes the number of sites in A) where Q A (t), P A (t) and M A (t) are the reduced correlation matrices, obtained by selecting the rows and the columns corresponding to A in (2.12), namely The reduced correlation matrices usually depend on the time t after the quench. In this section we study the circuit complexity when both the reference and the target states are mixed states corresponding to a subsystem A. In particular, we apply to the local quench that we are investigating the results for the circuit complexity of mixed states based on the Fisher information geometry [30], as done in [77] for a global quench.
We consider the reference state given by the reduced density matrix for the subsystem A at time t R 0 obtained through the local quench protocol characterised by {m R , κ R , ω R } and the target state given by the reduced density matrix of the same subsystem at time JHEP08(2021)135 t T t R , constructed through the quench protocol described by {m T , κ T , ω T } (see section 3.1). The corresponding reduced covariance matrices, denoted by γ R,A (t R ) and γ T,A (t T ) respectively, can be decomposed as done in (4.1).
The approach to the circuit complexity of mixed states based on the Fisher information geometry [30] provides also the optimal circuit connecting γ R, where 0 s 1, which is a covariance matrix for any s [97]. The length of the optimal circuit (4.2) is proportional to its complexity Both the reduced covariance matrices in (4.3) have the form (4.1), obtained by restricting to A the covariance matrix γ(t) in (2.11), as discussed above.
In our analysis we consider the simplest setup where the reference state is the initial state (i.e. t R = 0) and the target state corresponds to a generic value of t T = t 0 after the local quench. The remaining parameters are fixed to In this case the subsystem complexity (4.3) reads It is instructive to compare the temporal evolution of C A against the temporal evolution of the entanglement entropy S A after the same local quench, which can be evaluated from the symplectic spectrum of γ A (t) in the standard way [34,36,98,99]. The considerations above can be easily adapted to harmonic lattices in any number of spatial dimensions.

Numerical results
In the following we discuss some numerical results for the temporal evolution after a local quench of the subsystem complexity (4.4) in the case where the subsystem A is a block made by L consecutive sites in harmonic chains made by N sites. Let us remind that the reference state is the initial state (t R = 0) and the target state corresponds to the state at the generic value t T ≡ t 0 after the local quench, whose protocol is specified by the values of the parameters in (3.7) with m = 1 and κ = 1. For a given local quench, we display both the temporal evolution of the subsystem complexity C A and of the entanglement entropy S A .
The temporal evolutions in figures 6, 7, 8, 9, 11, 12 and 13 correspond to blocks A adjacent to the left boundary of the chain (as shown pictorially in the bottom left panel of figure 1) and for this bipartition the joining point is outside the subsystem whenever L < rN . The temporal evolutions in figure 10 correspond to blocks A whose midpoint coincides with the joining point (see figure 1, bottom right panel). While the temporal evolutions in figures 6, 7, 8, 9 and 10 are determined by the critical evolution Hamiltonian, for the ones in figures 11, 12 and 13 the evolution Hamiltonian is gapped.
Both the temporal evolutions of C A and S A exhibit revivals because our system has finite volume. For a generic values of r, the cycles correspond to p < t/(2N +2) < p+1, with p non negative integer (see figure 8 and figure 9), while only for r = 1/2 they correspond to p < t/(N + 1) < p + 1 (see figure 6 and figure 7) because of the symmetry provided by the fact that the joining point coincides with the midpoint of the chain [58,62].
Focussing on the temporal evolution during a single cycle, as N and L increase with L/N kept fixed, two different scalings are observed: one at the beginning and at the end of the cycle and another one in its central part. In these two temporal regimes, the curves obtained for different values of N collapse when the time independent quantity α log(N +1) is subtracted, with different values of α.
In figure 6 and figure 8 we show some temporal evolutions of C A − α log(N + 1) when r = 1/2 and r = 1/4 respectively. We find that α depends on (a) whether the joining point is outside (L < rN ) or inside (L > rN ) the subsystem; (b) the temporal regime within the JHEP08(2021)135 cycle where the collapse of the data is observed (either the central part of the cycle or its extremal parts). When the entangling point coincides with the joining point, i.e. L = rN , the collapses of the data in the different temporal regimes is observed for values of α that are slightly different from the ones adopted in the vertical axes of the panels in figure 6 and figure 8. In particular, when C A is not constant, the black curves in the left panels of these figures collapse with α 1/7, otherwise the data collapse is observed with α 1/10 (see figure 8, left panels).
The different scalings in the diverse temporal regimes within each cycle pointed out in (b) occur also for the temporal evolution of S A after a local quench [51,52,59]. Numerical results for the temporal evolution of S A after the same quench and for the same bipartition considered above (see the bottom left panel of figure 1) are reported in figure 7 and figure 9 for r = 1/2 and r = 1/4 respectively. In the case of r = 1/2, these numerical outcomes for S A are well described by the analytic curve discussed in [59], namely 1 1 See eq. (39) of [59] with c = 1 and v F = 1.

JHEP08(2021)135
(see the black dashed curves in figure 7) within the first cycle (then extended periodically to the subsequent cycles), where d ≡ 1 2 − L N +1 parameterises the distance between the entangling point and the joining point and we have introduced the temporal regimes T 0 ≡ (0, d) ∪ (1 − d , 1) and T 1 ≡ (d , 1 − d). The expression (4.5) holds only when r = 1/2 and the interval A is adjacent to one of the boundaries of the segment. The different scalings corresponding to the two different regimes within the cycle, which lead to subtract α log(N +1) with either α = 1/3 (top panels) or α = 1/6 (bottom panels), agree with (4.5). We remark that, since r = 1/2, the numerical curves in the left panels of figure 7 are identical to the ones in the right panels characterised by the same coloured marker: this is because the entanglement entropy of a subsystem is equal to the entanglement entropy of its complement when the entire state is in a pure state (this is the case for any t > 0 after the local quench that we are exploring).
As for the temporal evolution of the subsystem complexity C A , when r = 1/2 and L < N/2, hence the joining point is outside the subsystem (see the left panels of figure 6), we find that it is qualitatively similar to the temporal evolution of S A . Combining this observation with the different scalings obtained numerically, we are led to consider the following ansatz within the first cycle (the parameter d and the temporal regimes are introduced in (4.5)), which is then extended periodically to any value of t > 0. In the top left panel of figure 6, a remarkable agreement is observed between the numerical data and the ansatz (4.6), which corresponds to the black dashed curves.
Considering also the right panels of figure 6, where r = 1/2 again but L > N/2, we find that the temporal evolutions of C A for blocks that include the joining point are qualitatively different from the ones corresponding to blocks that do not contain the joining point. Indeed, in the right panels of figure 6, focussing e.g. on the first cycle and considering t/(N + 1) < 1/2 (the regime t/(N + 1) > 1/2 is obtained straightforwardly through a time reversal), we observe three regimes: an initial growth until a local maximum, followed by a fast decrease and then another growth, milder than the previous one (it becomes almost flat as L/N increases).
When r = 1/2, the symmetry under a spatial reflection with respect to the midpoint of the chain does not occur and more regimes are observed within a cycle, for both C A and S A .
The same quantities considered in figure 6 and figure 7, where r = 1/2, are shown in figure 8 and figure 9 for r = 1/4. Notice the different periodicity with respect to the case of r = 1/2, as already mentioned above. The main feature to highlight is the qualitative difference between the temporal evolutions of C A when the joining point lies outside A (see figure 8, left panels) and the ones corresponding to blocks that include the joining point (see figure 8, right panels). Furthermore, the temporal evolution of C A when the joining JHEP08(2021)135   point lies outside A is also qualitatively similar to the one of the corresponding S A (see figure 9, left panels).
When L < rN , focussing on the temporal evolutions in the first half of the first cycle (i.e. 0 < t 2N +2 < 1/2), for both C A and S A we observe three regimes (left panels of figure 8 and figure 9): first a flat curve, then a growth followed by a decrease and finally another regime where the evolution is almost constant. This means that, when L < rN , for the temporal evolutions within the first cycle we identify five regimes. The values of t

2N +2
at which the changes of regime occur are given by rN − L, L + rN , 2N − L − rN and 2N − rN + L, whose time ordering depends on the explicit values of N , r and L. In the special case of r = 1/2, we have only three regimes within the first cycle (first a flat regime, then a growth/decrease regime and finally another flat regime), as one can observe from figure 6 and figure 7, but also from the analytic expressions in (4.5) and (4.6).
When L > rN and therefore the joining point is inside the subsystem, by comparing the right panels of figure 8 against the right panels of figure 9, it is straightforward to realise that the temporal evolutions of C A and S A are qualitatively very different. In particular, while the temporal evolutions of S A in the right panels of figure 9 are similar to the ones displayed in the left panels of the same figure (e.g. the same five regimes mentioned above), as expected from the fact that S A = S B for the spatial bipartition A ∪ B of the system in a pure state (the qualitative difference is only due to the asymmetric position of the joining point), the temporal evolutions of C A in the right panels of figure 8 are more complicated than the ones in the left panels of the same figure, which correspond to blocks that do not include the joining point. For instance, considering the temporal evolution of C A immediately after the quench, a rapid initial growth is observed when L > rN (highlighted in the insets in the right panels of figure 8), while it remains stationary when L < rN . We remind that, whenever L = rN , also the temporal evolution of S A right after the quench remains stationary (see figure 7 and figure 9). As for initial growth of C A when L > rN , an interesting numerical observation that we find it worth remarking is the fact that the logarithmic curve (3.25), which has been first employed in the bottom panels of figure 2 to describe the logarithmic growth for the complexity of the entire chain, occurs also in the temporal evolution of the subsystem complexity; indeed it corresponds also to the dashed lines displayed in the bottom right panels of figure 6 and figure 8. Notice that, when the joining point is outside the subsystem (see the left panels of figure 6 and figure 8), a logarithmic growth right after the stationary regime is observed, but in this case the coefficient of the logarithm is different from the one in (3.25), as one can infer from the second line of (4.6) when r = 1/2 and t/(N + 1) d.
While in figures 6, 7, 8 and 9 the block A is adjacent to a boundary (see figure 1, bottom left panel) and therefore only one entangling point occurs, in figure 10 we consider some temporal evolutions of C A and S A when r = 1/2 and the joining point coincides with the midpoint of A (see figure 1, bottom right panel), hence two entangling points separate A from its complement B, which is made by two disjoint intervals adjacent to different boundaries. By construction, for this configuration the joining point is always inside the subsystem. The blocks providing the reduced covariance matrix (4.1) for this bipartition are obtained by restricting the indices of the matrices Q, P and M in (2.12) to i, j ∈ N 2 − L 2 + 1, . . . , N 2 + L 2 .

JHEP08(2021)135
The numerical results shown in figure 10 for some temporal evolutions of C A (left panels) and S A (right panels) after local quenches correspond to critical evolution Hamiltonians, i.e. with ω = 0. Also in this numerical analysis we subtract α log(N + 1) with the proper value of α, in order to observe collapses of data sets corresponding to the same L/N when N is large enough, finding that α depends both on the quantity (either C A or S A ) and on the temporal regime within the cycle where the data collapses are observed (either the central regime or the initial and final regimes). Interestingly, by comparing the left panels of figure 10 against the right panels of figure 6, we observe that, when the joining point is inside the subsystem A, the temporal evolutions of C A are qualitatively very similar, despite the fact that the number of entangling points is different in the two figures. Moreover, the values of α employed are the same, which are therefore independent of the number of the entangling points. Instead, let us remind that the values of α to employ for S A depend on the number of the entangling points, as one realises by comparing the right panels of figure 10 against the right panels of figure 7.
Focussing on bipartitions where the joining point lies inside the subsystem, by comparing the left and the right panels of figure 10, one notices that, while S A is constant at the beginning of its evolution, C A increases immediately. This feature has been highlighted also during the comparison of the right panels of figure 6 and figure 8 against the right panels of figure 7 and figure 9, where only one entangling point occurs.
Another interesting difference between the temporal evolutions corresponding to the two bottom panels in figure 1 is that the first local minimum occurs at t N +1 We find it worth mentioning some intriguing similarities between the temporal evolution of C A after the local quench discussed above and the one after the global quench studied in [77]. Let us consider the block made by L consecutive sites in the infinite chain and compare the temporal evolution of C A after the global quench of the mass parameter, as done in [77], against the one after the local quench where two half-lines are joined at the midpoint of A. The latter temporal evolution can be inferred by taking e.g. the red curves in the left panels of figure 10 for t N +1 < 1 2 , while the former one corresponds e.g. to the black data points in the top panel of figure 14 of [77]. These temporal evolutions are qualitatively very similar. However, important differences occur when these curves are studied quantitatively. For instance, while at the beginning a logarithmic growth is observed in the case of the local quench, as remarked above, a power law behaviour occurs in the case of the global quench [77]. Notice that, by performing the same comparison for the temporal evolutions of the corresponding S A , qualitatively different behaviours are observed (see the red data points in the right panels of figure 10 against the black data points in the bottom panel of figure 14 in [77]). It would be interesting to explore further these comparisons by considering different kind of quenches and performing a quantitative analysis.
In the final part of this discussion we consider temporal evolutions of C A and S A after local quenches characterised by gapped evolution Hamiltonians, for some fixed values of ωN > 0. Since for ω = 0 the qualitative behaviour of the temporal evolution of C A depends on whether the joining point is located inside or outside the block A, let us explore these two cases also when ω > 0. Considering the bipartition shown in the bottom left panel of figure 1, where r = 1/2, in figure 11 we display some temporal evolutions of C A for two fixed values of L/N such that the joining point is either outside (left panels) or inside (right panels) the block A. In figure 12 the same analysis is performed in the case of r = 1/4. These numerical results show that the temporal evolution of C A depends on whether the joining point is inside or outside the subsystem. The temporal evolutions of S A for these quenches are reported in figure 13 and, since S A = S B for any t > 0 (the entire chain is A ∪ B), whether the joining point is inside or outside the block does not influence the qualitative temporal evolution of S A , as already remarked above (once the eventual asymmetric position of the joining point is taken into account).

JHEP08(2021)135
The approximate periodicity highlighted in the evolutions corresponding to ω = 0 is not observed in general when ω > 0. For small values of ωN an approximate periodicity can be identified for a temporal regime whose duration decreases as ωN increases.
When the block A contains the joining point, C A has a non-trivial initial growth, while the evolution of the corresponding S A is constant at the beginning. This is the same feature highlighted for the critical evolution through the comparison of figure   Let us conclude our discussion by mentioning some results about the temporal evolutions of the subsystem complexity obtained within the gauge/gravity correspondence [67][68][69][70]85].
In the Vaidya gravitational spacetimes, the temporal evolution of the holographic subsystem complexity has been studied through the prescription based on the volume of a particular spacetime slice [67,68], finding curves that qualitatively agree with the temporal evolution of the subsystem complexity after a global quench of the mass parameter in the harmonic chains discussed in [77].
It would be interesting to perform a comparison between the qualitative behaviour of the temporal evolutions of the subsystem complexity discussed in this manuscript and the one of the temporal evolutions of the holographic subsystem complexity in the spacetime describing the gravitational dual of the joining local quench [79,83].

Single site in the chain made by two sites
In the following we discuss the temporal evolution of C A after the local quench that we are exploring for the chain made by two sites, described in section 3.4, and the subsystem made by a single site. From (3.29) and (3.30) we have that the covariance matrices of the reference and the target states are respectively The blocks occurring in (3.30) are 2 × 2 matrices that can be written as where the expressions for A j , B j and D j are obtained by specifying (2.15) to N = 2 and V 2 has been defined in (3.26).
Considering the subsystem A made by the first oscillator of the chain, for the 2 × 2 reduced covariance matrices of the reference and of the target states we find respectively 1 has been defined in (3.29) and we have introduced A ± ≡ A 1 ±A 2 , B ± ≡ B 1 ±B 2 and D ± ≡ D 1 ± D 2 , which allow to construct The quantities A j , B j and D j , with j ∈ {1, 2}, depend on t and on the parameters of the local quench m, κ and ω as reported in (2.15). The eigenvalues of the matrix γ T,A γ −1 R,A can be written in terms of these quantities as follows By employing these results into (4.4), we obtain the subsystem complexity whose explicit expression in terms of m, κ, ω and t is quite cumbersome; hence we have not reported it here. In the top left panel of figure 14, the subsystem complexity (4.13) is shown for various ω 2 and κ = 1 and m = 1. The local maxima of the curves corresponding to different ω occur at the same values of Ω 2 t given by multiple integers of π (see the vertical lines), where Ω k with k ∈ {1, 2} is defined below (3.31). The same feature is observed also for the local minima, if Ω 1 t is employed as the independent variable on the horizontal axis instead of Ω 2 t.
In the asymptotic regime given byω → ∞, whereω has been introduced in (3.34), while ωt is kept fixed and finite (introduced in section 3.4), the expansion of (4.12) reads g TR,± = 1 + 1 2ω 4 sin(ωt) 2 ± sin(ωt) 2 − ωt sin(2ωt) + ω 2 t 2 + O 1/ω 6 . (4.14) By employing this result, for the expansion of (4.13) one finds In the top right panel of figure 14, this expression corresponds to the black solid line, while the other curves have been drawn through the exact formula (4.13). In the same panel, we have also reported C A of half chains with N = 2L > 2 (coloured symbols). We find remarkable their agreement with (4.13) for early times, which improves as ω increases.  Comparing figure 14 and figure 5, we notice that, while for the complexity of the entire chain made by two sites only a main oscillatory behaviour is observed, for the subsystem complexity we can identify two kinds of oscillations: one has a larger amplitude and period πω 2 /ω and another one is characterised by a smaller amplitude and period π/ω. Whenω 1, the amplitude of the latter oscillation becomes negligible and we find that the temporal evolution of the single site subsystem complexity is nicely described by the following ansatz which is compared against the exact result (4.13) in the bottom panel of figure 14.

Conclusions
We studied the temporal evolutions of the circuit complexity and of the subsystem complexity after a local quench by considering harmonic chains in a segment with Dirichlet boundary conditions and the local quench where two finite chains made by rN and (1−r)N sites (with 0 < r < 1) are joined at t = 0 [87]. The subsystem complexity C A in (4.3) has JHEP08(2021)135 been evaluated by employing the complexity of mixed bosonic Gaussian states based on the Fisher information geometry [30], which provides also the optimal circuit (4.2). For the sake of simplicity, we considered only the case where the subsystem is a block of L consecutive sites and we mainly studied the complexity of circuits whose reference state is the initial state at t = 0.
The covariance matrices of the reference state and of the target state at time t > 0 along the temporal evolution have been introduced in section 2. Then, in section 3 they have been employed to evaluate some temporal evolutions of the circuit complexity for the entire harmonic chain.
For any value ω 0 of the mass occurring in the evolution Hamiltonian, we found that the initial growth of the complexity immediately after the quench is linear (see (3.22)) with a slope given by (3.23), which can be evaluated numerically, as done in the bottom panel of figure 4. When the evolution Hamiltonian is critical (i.e. ω = 0), after the above mentioned initial growth we observe a logarithmic growth independent of r (see (3.25) and the bottom panels of figure 2). We expect to observe this feature also when the system is infinite. In our numerical analysis we have considered only finite systems. The temporal evolutions of the complexity for finite systems and ω = 0 display revivals, independently of r. Three temporal regimes are observed within the first half of the temporal interval containing a single revival: a growth followed by a decrease and finally a regime where the complexity does not evolve (see the top right panel of figure 2 and figure 3). In the case of r = 1/2, the latter regime does not occur; hence this choice halves the duration of a revival (see the top left panel of figure 2). When ω > 0, the temporal evolutions of the complexity are more complicated; indeed, for instance, an approximate periodicity is not observed (see figure 4). When ωN is large, the complexity rapidly changes through small variations about a constant value that is independent of r. Importantly, we have identified different temporal regimes where different scaling behaviours are observed as N increases. It would be interesting to explain these scalings through quantum field theory methods.
In section 4 we have explored the temporal evolutions of C A for the bipartitions shown in the bottom panels of figure 1, where either one entangling point or two entangling points occur. One of our main results is given by the numerical evidences that the qualitative behaviour of the temporal evolutions of C A depends on whether the block A contains the joining point. In the case of the spatial bipartition shown in the bottom left panel of figure 1, where one entangling point occurs, this qualitative difference is evident once the left panels are compared against the corresponding right panels both in figure 6 and figure 8 (where r = 1/2 and r = 1/4 respectively) when ω = 0 and both in figure 11 and figure 12 (where, again, r = 1/2 and r = 1/4 respectively) when ω > 0. When the evolution Hamiltonian is critical and the joining point is inside the block, during the initial regime of the temporal evolution of C A we observe the same logarithmic growth (3.25) occurring in the temporal evolution of the complexity of the entire chain (compare the insets of figure 6 and figure 8 against the bottom panels of figure 2). Furthermore, in the case of r = 1/2 and when the joining point lies outside the block, we find that the analytic expression (4.6) for the temporal evolution of C A nicely reproduces the behaviour of the numerical data.

JHEP08(2021)135
It is very instructive to compare a temporal evolution of C A (in this manuscript we have considered only circuits where the reference state is the initial state) against the corresponding temporal evolution of S A , obtained for the same bipartition and the same quench protocol. The temporal evolutions of S A for the bipartition shown in the bottom left panel of figure 1 have been reported in figure 7 and figure 9 (where r = 1/2 and r = 1/4 respectively) for ω = 0 and in figure 13 for ω > 0. We remark that, whenever the block A does not contain the joining point, the temporal evolutions of C A and of S A are qualitatively similar (for instance, both these quantities do not evolve immediately after the quench whenever L = rN ). Instead, they are qualitatively very different when the joining point is inside the subsystem; indeed, for instance, when L = rN we find that C A rapidly grows immediately after the quench, while S A remains constant for a while.
In this manuscript we have considered both the spatial bipartitions shown in the bottom panels of figure 1, in order to investigate the influence of the number of entangling points on the temporal evolution of C A . By comparing the right panels of figure 6 against the left panels of figure 10, where r = 1/2 and ω = 0, we observed that, when the joining point is located inside the block, both the qualitative behaviour of the temporal evolution of C A and the values of the scaling parameter α are not influenced by the number of entangling points. Instead, the values of the scaling parameter α for S A do depend on the number of entangling points (see the right panels of figure 7 and figure 10).
We find it worth remarking that the logarithmic growth of C A highlighted in the inset of the bottom right panels of figure 6 and of figure 8 for one entangling point and of the bottom left panel of figure 10 for two entangling points is described by the same curve (3.25), including the additive constant, which has been found for the temporal evolution of the complexity of the entire chain (see the bottom panels of figure 2).
We have also explored the temporal evolutions of the complexity for the entire chain and of the subsystem complexity after the local quench in the minimal setup where the chain is made by two sites, and therefore the subsystem A contains only one site (see section 3.4 and section 4.3). In this simple setup, we have obtained the analytic expressions given by (3.31) and (3.33) (see also figure 5) for the complexity of the chain and by (4.12) and (4.13) (see also figure 14) for the subsystem complexity. While these analyses are useful to get some insights about some regimes of the parameters (e.g. small t and largẽ ω), they do not capture many important features observed for large values of N .
As for the local quench considered in this manuscript, it would be interesting to explore more systematically the temporal evolutions of the subsystem complexity when ω > 0 or for asymmetric bipartitions involving two or even more entangling points, to obtain analytic results in the thermodynamic limit, to find bounds that still describe some essential features of the temporal evolution of the subsystem complexity and also to study the thermalisation of the subsystem complexity, as done in [77] for a global quench of the mass parameter.
The Gaussian states can be employed to investigate the temporal evolution of the subsystem complexity after some local quenches in higher dimensions and also in free fermionic lattice models. It could be instructive to explore these temporal evolutions by employing the entanglement spectrum or the entanglement Hamiltonians [34,36,96,[100][101][102][103][104][105][106][107][108][109][110][111], as done in [30] at equilibrium. It would be insightful to explore also the dependence JHEP08(2021)135 of the temporal evolution of C A on the reference state, by adopting a state different from the initial one as the reference state (e.g. the unentangled product state).
In this manuscript we have compared the temporal evolutions of the subsystem complexity against the ones of the corresponding entanglement entropy, but it could be interesting to perform analogue comparisons against the temporal evolutions of other entanglement quantifiers like the entanglement negativity [87,[112][113][114][115][116][117][118][119], the entanglement contours [120,121] and the relative entropies [122].
We remark that investigating the temporal evolutions of the subsystem complexity after various quantum quenches through lattice methods and quantum field theory techniques in interacting models is a very challenging task that deserves future studies. Holography can provide important benchmarks. Interesting analyses have been performed [123][124][125][126][127][128][129][130][131][132][133][134][135][136] and it would be interesting to employ these methods to explore also the out-of-equilibrium dynamics of the circuit complexity.

JHEP08(2021)135
By employing (2.11) and (2.13) (where V N is orthogonal), we find that the temporal evolution of (A.6) reads From (A.6) and the block decomposition in (2.14), we get where Q(t), P(t) and M(t) are diagonal matrices whose diagonal elements are given respectively by (A.5) with Ω 0,k replaced by µ.
We consider the temporal evolution of the complexity where the reference and the target states are the states at the values of time t R and t T respectively along a given quench, namely γ R = γ(t R ) and γ T = γ(t T ), with γ(t) given in (A.7).
Following the analysis reported in [77] (see also [72,73]), we compute the eigenvalues of γ T γ −1 R which provide the complexity (3.1), finding with C TR,k = 1 + 1 2 where Ω k is given in (2.16). Notice that, when µ = ω, we can set ω = 0 and obtain a finite result for the complexity in (A.9). Since C TR,k in (A.10) is an oscillating function of |t T − t R | for any k, the complexity C is finite, also for large values of |t T − t R |.
In figure 15 we show some temporal evolutions of the complexity (A.9) when t R = 0 and t T = t. In the top panels we keep µ = ω, setting either ω = 0 (left panel) and ω > 0 (right panel). These temporal evolutions are qualitatively similar to the ones observed for the global quench discussed in appendix A.1, as expected from the fact that the role of Ω 0,k is played by µ in this case.
In the case where only the quench of the spring constant is performed, i.e. when µ = ω (see the bottom panels of figure 15), by using the explicit expression of Ω k in (2.16), one obtains C TR,k = 1 + 1 2 which provides the complexity through (A.9). In this case the critical evolution cannot be explored because (A.11) diverges for any k as ω → 0. For a given ω > 0, the coefficient of (sin[Ω k (t T − t R )]) 2 in the r.h.s. of (A.11) becomes negligible when k N +1 1, while it reaches its maximum when k N +1 1. Thus, the main contributions to the complexity given by (A.9) and (A.11) come from the modes such that with frequency approximately equal to 2, independently of ω; instead, when ω 2 4, the frequency of the oscillations is √ ω 2 + 4. This behaviour can be observed in the bottom panels of figure 15, where some temporal evolutions of the complexity given by (A.9) and (A.11) are shown, in the case of t R = 0 and t T = t, i.e. when the reference state is the (initial) unentangled product state. When ω ∈ {0.05, 0.2, 0.5}, in the initial part of the evolution the curves for C/ √ N collapse displaying the same oscillatory behaviour with frequency independent of ω and approximately equal to 2. Instead, when ω ∈ {1.5, 4} and therefore ω 2 ∼ 4, the frequency of the oscillations in the initial part of the curves depends also on ω and it is given by √ ω 2 + 4.

B Euler decomposition for a class of symplectic matrices
The Euler decomposition is a powerful tool to evaluate the circuit complexity for pure states. Indeed, (3.1) can be written also in terms of the squeezing parameters of the symplectic matrix W TR (see (3.16)), as discussed in section 3.1 [19,72,77]. In this appendix we derive the analytical expressions for the matrices in the Euler decomposition (3.15) for a specific class of symplectic matrices. This provides the Euler decomposition of the matrix E defined in (2.14), which is exploited to get (3.17). Consider a 2N × 2N matrix M partitioned into the four N × N blocks S, U , Y and Z which are diagonal matrices (we denote respectively by s k , u k , y k , z k their k-th element JHEP08(2021)135 The sign in (B.6) has to be fixed case by case, checking that the correct matrix M k is obtained through the decomposition (B.

C Derivation of the initial growth
In this appendix we report the derivation of (3.22) and (3.23), which provide the initial linear growth of the complexity (3.1) after a local quench.
As t → 0, from (3.21) we find the following expansion where the last step has been obtained by using that whose trace can be easily computed by using that V N is orthogonal and the matrix N defined in (3.18) is diagonal. From (C.3) we have that the first two terms within the square root in (C.2) are negative; thus, in order to have c 1 0, the term containing the anticommutator under the square root must be positive. For this term we find Tr E (1) , γ (l,r)

JHEP08(2021)135
From (3.19), (2.5) and (2.6), we get where V 0 has been defined in (2.9). Then, by exploiting (2.6) and (2.7), we obtain Tr E (1) where Ω (l) k and Ω (r) k are given by (2.8) and m l = m r = m for the local quench that we are considering. Since we are not able to simplify V N V t 0 , we cannot write an analytic expression for the last term in the r.h.s. of (C.7). Finally, by using (C.3), (C.5) and (C.7), we obtain that the slope c 1 in (C.2) becomes the expression (3.23) in the main text.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.