Circuit complexity for generalised coherent states in thermal field dynamics

In this work, we study the circuit complexity for generalized coherent states in thermal systems by adopting the covariance matrix approach. We focus on the coherent thermal (CT) state, which is non-Gaussian and has a nonvanishing one-point function. We find that even though the CT state cannot be fully determined by the symmetric two-point function, the circuit complexity can still be computed in the framework of the covariance matrix formalism by properly enlarging the covariance matrix. Now the group generated by the unitary is the semiproduct of translation and the symplectic group. If the reference state is Gaussian, the optimal geodesic is still be generated by a horizontal generator such that the circuit complexity can be read from the generalized covariance matrix associated to the target state by taking the cost function to be $F_2$. For a single harmonic oscillator, we discuss carefully the complexity and its formation in the cases that the reference states are Gaussian and the target space is excited by a single mode or double modes. We show that the study can be extended to the free scalar field theory.


Introduction
Complexity has been a focus in the recent study of the AdS/CFT correspondence [1] and blackhole physics. As firstly pointed out by L. Susskind [2], "entanglement is not enough" to describe the dynamics of the blackhole, especially the growth of the Einstein-Rosen-Bridge(ERB). Instead, he proposed [3] that the growth of the ERB should be dual to the growth of the quantum complexity of the evolving state, the thermofield double (TFD) state [5]. There are two proposals put forward by Susskind and his collaborators to quantify the size of the ERB: one is the "complexity=volume"(CV) conjecture [4] which states that the holographic complexity is given by the volume of the codimension-1 maximal spacelike surface in the bulk connecting the left and right sides; the other is the "complex-ity=action"(CA) conjecture [6,7] which states that the holographic complexity is captured by the gravitational action of the bulk region known as the Wheeler-de Witt (WdW) patch bounded by light sheets. Both conjectures introduce the gravitational observables which probe the spacetime region deep behind the black hole horizon, and therefore they have been intensely discussed since their birth [8,9,10,11,12,13,14,15,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. Nevertheless, the development is hindered by our poor understanding the quantum complexity in the dual field theory.
Originally the complexity is a concept in theoretical computer science [38,39], characterizing the difficulty in carrying out a task. In quantum computing, one may find a unitary operationÛ which maps an input quantum state for some number of qubits to an output quantum state with the same number of qubits [40,41,42]. In a circuit model,Û could be constructed from some elementary gates. There could be many ways in constructingÛ to some accuracy > 0. The circuit complexity of the unitaryÛ is given by the minimal number of elementary gates required to construct the desired unitaryÛ , up to some tolerance . However, to generalize the above definition to quantum field theory, even a free field theory, is highly nontrivial, due to the fact that there are infinite number of degrees of freedom in a field theory. In order to define the circuit complexity, one first needs to identify the reference state and the target state, and then identify the optimal circuit out of the infinite number of possible circuits connecting the reference state and the final target state.
There have been some initial steps in studying the complexity in quantum field theory. In [43], the circuit complexity of the ground state of a free scalar field theory was investigated.
The optimal circuit was determined geometrically by the minimal geodesic in the space of unitariesÛ with a suitable metric, as developed by Nielsen and his collaborators [44]. This approach has been applied to free fermionic theories in [45,46]. Another similar geometric definition of the complexity 1 based on the Fubini-Study metric has been explored in the free scalar field theory in [47]. The circuit complexity in interacting field theories has been discussed in [53].
The study of the complexity has been generalized to the TFD state in free scalar field theory [35,48,49]. In this case, the target state is the TFD state, while the reference state has different choices: in [35] the reference state was chosen to be composed of two copies of the reference state used in [43,47]; in [48,49] the reference state was two unentangled copies of the vacuum state. Due to the difference in reference state and other points, the complexity for the TFD state in two approaches differs in many ways, especially the one for the time-dependent TFD state.
In this work, we would like to study the circuit complexity of general coherent state, extending the study of complexity of coherent state in [37]. We will focus on the coherent states in thermal systems. There are two kinds of coherent states in a thermal system: the coherent thermal (CT) state and thermal coherent state. As these two kinds of states are somehow equivalent, we will consider the circuit complexity for coherent thermal state by applying the covariance matrix approach developed in [35]. Since the one-point function of the CT state is not vanishing, the two-point function is not enough to characterize the state. Nevertheless we show that the covariance matrix formalism is still applicable after some improvement. The essential point is that if the reference state remains Gaussian, the optimal geodesic will still be generated by a horizontal generator such that the circuit complexity can be read from the norm of the generator.
For a simple harmonic oscillator system, we compute the complexity and discuss the formation of the complexity by choosing various Gaussian reference states and the target states with single mode and double modes. We extend our study on the complexity of the CT state to the free scalar field theory by fixing the reference state to be the Dirac vacuum state.
The remaining parts of the paper are organized as follows. In section 2, we give a brief introduction to the coherent state and the thermal vacuum state in the harmonic oscillator system. In section 3, we introduce the general coherent states in thermal systems. In section 4, we study the circuit complexity for the coherent thermal state. Due to the loss of Gaussianity of the CT state, we need to generalize the covariance matrix approach and furthermore compute the circuit complexities for the Gaussian reference states. In section 5, we study the complexity of CT state in a free scalar field theory by choosing the Dirac vacuum state to be the reference state. We end with conclusions and discussions in section 6.

Preliminaries: coherent state and thermal vacuum state
In this section, we shall briefly review the construction of the Glauber coherent state and the thermal vacuum state. It will be shown in the next section that proper considerations from these two states lead to several different generalizations of the coherent state to the thermal field dynamics. Besides, to study the circuit complexity of the coherent thermal state for a free field theory, we will begin our story with a toy model: the harmonic oscillators. The complexity for the field theory will be discussed at last in Sec. 5.

Coherent state
For a single harmonic oscillator with a Hamilton H = p 2 /2m + 1 2 mω 2 q 2 , the annihilation and creation operator can be defined by with p = −i∂ q . They satisfy the commutation relation [a , a † ] = 1 .
The Hamilton can be expressed into the form of The vacuum state is defined by The energy eigenstates are defined by the creation operators acting on the vacuum It is known that these states form a complete basis in the Hilbert space, namely ∞ n=0 |n n| = 1 .
The Glauber coherent state is another kind of interesting excited state. It is defined by the eigenstate of the annihilation operator where in general α is a complex c-number since the operator a is not Hermitian. In fact, the coherent state can also be obtained by a particular operator D(α) acting on the vacuum where D(α) is called displacement, defined as Note that D(α) is anti-Hermitian because of D † (α) = D −1 (α) and it obeys By simple calculations, one finds that the coherent state is a superposition of the energy To end this subsection, we would like to introduce a new set of states {|n , α } that is of great importance in the construction of generalised coherent states in thermal field dynamics. The states are produced by the displacement operator acting upon the energy eigenstates |n , α ≡ D(α)|n , where n = 0 , 1 , 2 , · · · . These states are complete as well ∞ n=0 |n , α n , α| = 1 .
In fact, this set of states can be created and annihilated by the following operators In this case, the Glauber coherent state can be viewed as a ground state because of Moreover, one has so that These relations are similar to those between the energy eigenstates and the operators a , a † .

Thermal vacuum state
In thermal field dynamics, the finite temperature problems are treated by using the techniques developed for zero temperature quantum field theories. The price for this convenience is that one needs to deal with an enlarged Hilbert space, which is a direct product of two copies of the ordinary zero temperature Hilbert space. We denote where H L and H R stand for the ordinary Hilbert space for the zero temperature theory on the left hand side and the right hand side, respectively. In the following, all the operators and the state vectors for the left/right hand side will be assigned with a subscript "L/R".
The creation and annihilation operators obey the commutation relations It is worth emphasizing that the thermal vacuum state is a kind of excited states, rather than the true vacuum state of the Hilbert space H L ⊗ H R . This will be clear from the relation (26). In recent literatures, it is usually called the Thermal Field Double (TFD) state. We will use this for a shorthand notation throughout this paper.
To build the TFD state, we first introduce an anti-Hermitian operator where β = 1/T is the inverse of temperature and θ is related to the temperature by Note that tanh θ = e −βω/2 . Under the Bogoliubov transformation, the creation and annihilation operators are transformed as 2 The new operators obey the following commutation relations The TFD state is defined by or by the operator U (β) acting on the vacua The two definitions are equivalent, as one can check using the Bogoliubov transformation.
Furthermore, the TFD state can be written explicitly as where in the second equality we have adopted the operator identity However, since the state is introduced to deal with the finite temperature problems, it is interesting to see how the state become for a local observer, for example on the left hand side. Tracing over the degrees of freedom on the right hand side, one finds the reduced density matrix Clearly, this is a thermal density matrix, describing an ordinary thermal equilibrium state.
In fact, this should be a priori for the construction of TFD state in thermal field dynamics.
Likewise, it is a priori rule to test the generalisations of the coherent state: a suitable generalisation should not only be coherent in thermal field dynamics but also be thermal for the one-sided theory.

Time dependent TFD state
Since the theory under consideration is a direct product of two copies of ordinary zero temperature theory, the time evolution operator is given by e −i(H L t L +H R t R ) . The time dependent TFD state is obtained by In principle, the evolution on one side is independent from the one on the other. In this paper, we choose the symmetric case t L = t R = t/2 for the sake of convenience. It was established [35] that the time dependent TFD state can be written into a nice form similar to (26), by using the explicit expression (27) for the TFD state. However, in the following, we would like to provide a different derivation by using the operator algebra directly. The new approach is more neat and more suitable for our later purpose.
By settingH = (H L + H R )/2, one finds where U (β , t) ≡ e −iHt U (β)e iHt . After introducing the operators which obey the commutation relation one obtains where z = θe −iωt and in the last line we have adopted the Baker-Campbell-Hausdorff (BCH) formula to deduce the relation Moreover, it turns out that the operator U (β , t) can be recast into a more compact Gaussian form in terms of the canonical variables ξ a = (q L , q R , p L , p R ), where In other words, the TFD state is a Gaussian state up to an unimportant phase factor because Expressing the state into this form is particularly useful when computing the complexity in the covariance matrix approach [35]. We will turn to this point in Sec. 4.

Generalised coherent states in thermal field dynamics
Now we are ready to introduce the generalised coherent states in thermal field dynamics.
Interestingly, we find that there are two types of generalisations in literature [54,55,56,57,58,59,60,61]. In [61], the generalized states are called Coherent Thermal (CT) state and Thermal Coherent (TC) state, respectively. However, we will show that though the two states are indeed defined differently, they can be related to each other via a parameter transformation, see (55) and (57). As a consequence, the two definitions are equivalent in the sense that they just scan the eigenvalue spaces for the same set of states in different ways, leading to the apparent differences.
Before introducing these states, we explain some of our notations below at first. The ordinary Glauber coherent state will be represented as |α L |γ R , where α , γ characterize the eigenvalues of the left-and right-hand-side annihilation operators respectively. The state is obtained by acting the displacement on the vacua where D(α , γ) is the product of the displacements on two sides

Coherent thermal state
A CT state is defined by an anti-Hermitian operator U (β) acting upon the Glauber coherent state [61] |CT where it was understood that The operator U (β) is related to U (β) via a unitary transformation Therefore, one has In other words, the state is first thermalized and then displaced. Moreover, by using Eq.
(27) for the TFD state, the CT state can be expressed into a similar nice form It is clear that the CT state is a two-parameter generalisation of the thermal vacuum state.
On one hand, when α = γ = 0, it reduces to the latter. On the other hand, the state for a local observer becomes thermal because of Note that it describes the thermal equilibrium in terms of the set of states {|n , α }.
The time-dependent coherent thermal state is produced by where Using the relations e −iHt a L ,R e iHt = e iωt/2 a L ,R , e −iHt a † L ,R e iHt = e −iωt/2 a † L ,R , we deduce Furthermore, using Eq.(1), it can be expressed as an exponential of the superposition of canonical variables in the phase space where where f and f denote the real and the imaginary part of f , respectively. From this, it is clear that the CT state is non-Gaussian. In sec. 4, we will adopt the above result to compute the complexity of a CT state using the covariance matrix approach.

Thermal coherent state
Compared to the CT state, the thermal coherent (TC) state is defined by thermalizing a Glauber coherent state [61] Note the different order of the operators U (β) , D(α , γ) acting on the vacuum state in (43).
According to the thermal Bogliubov transformation (23), the TC state turns out to be the eigenstate of thermal annihilation operators, namely In addition, from (52), the TC state is related to the thermal vacuum state via a thermal displacement operator D(α , γ; β) where D(α , γ; β) is defined by It is interesting to note that the thermal displacement D(α , γ; β) is related to ordinary displacement D(α , γ) via the parameter transformation (55). One has Therefore, a TC state is related to a CT state by It strongly implies that the two set of states are equivalent. Of course, a particular TC state with given parameters α , γ is distinguished from the CT state with the same parameters since their dependence on the temperature are very different. However, the relation (57) tells us that the two sets of states just scan the eigenvalue space for the same set of states in different ways, leading to the apparent differences.
In fact, the equivalence between the two sets of states becomes even more clear for several special cases. For example, when γ = ±α * , it was argued [61] that both the TC state and the CT state belong to the generalized coherent states with respect to the Lie group E(1 , 1). We define It follows that the generators Thus, in this case, the two states can be collectively represented as From this, it is obvious that the two states must be related via a transformation of parameters since they both give a representation of the same group. Another interesting case is when α, γ are real. The parameter transformation (55) reduces to a boost. The Lie algebra The algebra turns out to be a subalgebra of two commuting sl(2 , R) algebras.
The time-dependent TC state is produced by where Thus, a time-dependent TC state is again related to a time-dependent CT state by which generalizes the static counterpart (57). According to this relation, the complexity for a TC state can always be obtained from that of a corresponding CT state To avoid redundancy, we will focus on the CT state in the computation of circuit complexity.

Circuit complexity for generalised coherent states
In this section, we will adopt the covariance matrix approach to calculate the circuit complexity for the CT state. However, the original approach developed for the TFD state [35] cannot be applied to our case directly since now our target state is non-Gaussian. One has from (31) and (46) Note that an overall phase factor e −iωt/2 has been dropped since it does not play any role in our discussions. The non-Gaussianity of the state implies a non-vanishing one-point function of the state. To deal with this and derive the complexity, we will generalize the covariance matrix approach appropriately.

Covariance matrix approach: generalisations and complexity
Considering a quantum system with canonical coordinates ξ a ≡ (q 1 , · · · , q N , p 1 , · · · , p N ), one has the commutation relations where Ω ab is an anti-symmetric tensor, given by In the following, we will use Ω ab to raise indices for tensors, for example It was established in [35] that for a pure Gaussian state which has a vanishing one-point function, it is completely characterized by its symmetric two-point function, usually referred to as the covariance matrix where the 1/2 factor is introduced for our later purpose. However, for the coherent state and its thermal generalisations, the one-point function is non-vanishing and hence the covariance matrix is not enough to characterize the state.
From (66), a CT state is connected to the reference state |ψ R , either Gaussian or non-Gaussian, by two kinds of unitary transformations The first unitary e −iλaξ a produces a translation for the canonical variables ξ a and hence breaks the Gaussianity of the target state. It turns out that to compute the circuit complexity for such a non-Gaussian state, we need to enlarge the covariance matrix (70) properly.
For this purpose, we first introduce the one-point function and then study how the one-point function ϕ a and the symmetric two-point function G ab transform under the unitaryÛ .
By simple calculations, we find the following relations where Making use of the BCH formula, we get where M ≡ e K . Combining the above results, we are led tô It follows that under the unitaryÛ the one-point function transforms as whilst the symmetric two-point function behaves as In matrix language, the above results can be expressed simply as This is a general transformation that connects two generally non-Gaussian states. For two Gaussian states, it reduces to ϕ = ϕ = 0 , G = M GM T , consistent with the result in [35].
Based on the transformation (80), we introduce an extended covariance matrix so that the transformation (80) can be expressed in a more compact form where forms a matrix representation for our unitary group.
The quantum circuit that connects the reference state |ψ R to the target state |ψ T will be constructed by a series of elements U. By definition, the complexity is given by the length of the minimal geodesic in the group manifold. However, the geodesic may not be unique and each of them may have a different length, as shown in [43,35]. The reason is that there exists a stabilizer subgroup V for the reference state, i.e., ∀ U V ∈ V , which leads to where B V is the generator of the subalgebra V. This implies that if a unitary e A (or a geodesic y(t) = e tA in the group manifold) connects the reference state to the target state, there will be a lot of unitaries (or geodesics) achieving the same goal because of Of course, the unitaries U V are redundant in the construction of the target state. In other words, if a geodesic that connects the reference state to the target state has unitaries U V , it will not have the minimal length. This in turn tells us that the optimal geodesic definitely should not have any unitary belonging to the stabilizer subgroup.
To find the optimal geodesic, we first define the inner product on the Lie algebra of the transformation group where A, B are two infinitesimal generators,Ḡ R is the extended covariant metric associated to the reference state andḡ R is its inverse matrix. Next we define the horizontal subspace that is transverse to the stabilizer subalgebra Clearly, a horizontal generator A ∈ H does not belong to the stabilizer subalgebra and vice versa. However, it does not mean that the Lie algebra G associated to the transformation group must be split as G = H ⊕ V. It is still possible that some of the generators do not belong to the product space H ⊕ V, depending on the group structure as well as how we choose the reference state and the target state. It has been shown in [35] that when both |ψ R and |ψ T are Gaussian, the transformation group is Sp(2N , R) and the Lie algebra can indeed split into the product form sp(2N , R) = H ⊕ V. In this case, the optimal geodesic will be generated by a horizontal generator. However, in our case, the target state is non-Gaussian and the group structure has been enlarged as well as the Lie algebra. It turns out that the above decomposition for our algebra becomes invalid. Nevertheless, as we will show soon, the extra generator in the algebra is trivial as long as the reference state is Gaussian. It does not connect the reference state to the target state so that the optimal geodesic in our case will still be generated by a horizontal generator, i.e., y(t) = e tA , where A ∈ H. The complexity of the target state will be given by where we have taken the cost function to be F 2 . Therefore, evaluating the complexity is equivalent to finding the horizonal generator A such that The uniqueness of the generator will be guaranteed by our derivations below.
In the remaining of this paper, to derive the complexity, we choose the reference state to be a Gaussian state, which has the extended covariance matrix In this case, the stabilizer subgroup can be defined by Consequently, its Lie algebra satisfies According to (87), this leads to It implies that the horizontal generators obey From Eq.(92), the Lie algebra of the stabilizer subgroup can be expressed as whilst Eq.(94) implies that the horizontal generators satisfy Notice that though O A ⊕ O B gives rise to the full subalgebra sp (2N , R), the decomposition G = H ⊕ V does not hold because an extra generator remains. One finds However, according to the transformation (82), the last generator just transforms a unity to a unity and hence is trivial. Therefore, in our case the optimal geodesic will still be generated by a horizontal generator.
For later convenience, we parameterize the horizontal generator that connects the reference state and the target state to be where ν is a vector which will be determined later and K ∈ sp(2N , R) obeying Evaluating the exponential of the horizonal generator yields Comparing to (83), one immediately finds χ = λ. The vector ν is determined by It is worth emphasizing that a nontrivial vector ν is essentially determined by a nonvanishing one-point function for the target state. On the contrary, when the target state is Gaussian, one has ϕ = 0 = λ, leading to ν = 0.
Now it is straightforward to derive the complexity However, we have not solved the generator K in terms of the information for |ψ R and |ψ T .
According to the transformation (80), the one-point function and the covariance matrix of the target state are given by (recall that the reference state is Gaussian) On the other hand, using (99), one has Then from (103), one finds where ∆ T ≡ G T g R is called the relative covariance matrix [35]. This solves the generator K and its matrix exponential. However, it is worth emphasizing that for any given Gaussian reference state, the quantity ∆ T as well as the generator K does not really depend on the non-Gaussianity of the target states (this is not hard to understand since K generates rotations only for the canonical variables ξ a , instead of translation). In sec.4.3, we will explicitly show that ∆ (0) T is nothing else but simply the relative covariance matrix for the ground thermofield double state, namely The complexity turns out to be In general, the result gives rise to the complexity between a non-Gaussian target state and a Gaussian reference state. However, it also establishes the relation between the complexity for non-Gaussian target states and that of the ground state. In fact, for any Gaussian reference state, one has according to (127) where the equality is taken when the target state becomes purely Gaussian. The result implies that for a same Gaussian reference state, the complexity for an excited state is always larger than that of the ground state.
Last but not least, despite that our target state has been specified in (66), we can still choose different Gaussian reference states. With different choices, the difficulty for the calculations and the results will be significantly different, as will be shown below.

Dirac vacuum as the reference state
In this subsection, we will calculate the complexity analytically by choosing the Dirac vacuum as the reference state, namely |ψ R = |0 L |0 R .
By definition, the covariance matrix for this reference state can be evaluated as It turns out that for this particular reference state, the unitary connecting |ψ R to |ψ T has already been given in (66) and hence the generator K can be read off directly. One has cb . The matrix form is given by Next we derive its matrix exponential M = e K . We find that the generator K can be diagonalized by a matrix S, i.e., K = SK 0 S −1 with K 0 = diag{−θ, −θ, θ, θ} and It is easy to see that Tr K 2 = Tr K 2 0 = 4θ 2 . The transformation matrix M can be evaluated as where e K 0 is a diagonalized matrix e K 0 = diag{e −θ , e −θ , e θ , e θ }. We deduce On the other hand, we have where λ a 's have been specified in (51). With these results in hand, it is straightforward to evaluate the vector ν using the formula (101). We obtain Finally, using (102) we read the complexity Remarkably, the result is independent of time! Notice that the complexity only depends on the quadratic polynomials of the real and imaginary parts of α and γ. Thus, we have In fact, for some special cases, the complexity enjoys even more symmetries.
a) Single excitation (γ = 0): At first, we would like to consider single excitations.
Without loss of generality, we set γ = 0. Then, we have It is immediately seen that the complexity depends only on |α|, implying that for any α = |α|e iϕ , C a (α , 0) = C a (α , 0). This is much stronger than the relation (117).
b) Two equal excitations (γ = α): Secondly, we move to the two equal excitation case. The Eq. (116) gives rise to One easily finds It is interesting to compare C b with C a and 2C a for a same eigenvalue α. We have Recall tanh θ = e −βω/2 > 0 which implies θ > 0. However, since α and α are arbitrary, the sign of the above two equations are not fixed. It implies that the two modes are highly entangled and have a strong temperature dependence.
Furthermore, we have which is positive for | α| > | α| and negative for | α| < | α|, independent of the temperature. d) Two conjugated excitations (γ = α * ): At last, let's turn to the two conjugated excitation case. The complexity is given by It is interesting to note that the result depends only on |α|, similar to the single excitation case. Moreover, comparing C d with C c and C b for a same α, we find which are both negative for any non-vanishing α as well as the temperature.
However, it is worth emphasizing that the above relations between the complexity for different types of excitations strongly depend on the reference state. We will turn to this point again in sec.4.3.2 (see the last paragraph for two conjugated excitations).

More general reference state
While we have chosen the reference state to be a Gaussian state, it is not necessarily to be the Dirac vacuum. A more general Gaussian state could be characterized by a covariance matrix where η > 0 and the state is no longer the Dirac vacua when η = 1. As a consequence, the generator K cannot be read off from (66) It is immediately seen that ∆ T does not rely on the translation generator λ. It simply gives the relative covariance matrix for the ground thermodfield double state.
However, now it is of great difficult to solve the generator K and the vector ν analytically.
Nevertheless, since we have extracted the covariance matrix G T for the target state, it is straightforward to solve them numerically. This is a purely algebraic problem: the generator K can be solved from (105) and then the vector ν will be determined by (101). Finally, it is straightforward to obtain the complexity through (102) or (107).
According to these considerations as well as the expression (51) for the translation generator λ, we can reasonably speculate that in general the complexity for the CT states is a time periodic function and may take the form of where N = 0 , 1 , 2 , · · · is some nonnegative integer. The specific value of N as well as the amplitudes a n and the initial phases ϕ n for above each term are all functions of the parameters (α , γ , β , η). Indeed, by carefully scanning the parameters space, we find that all of our numerical results can be perfectly fitted by the above formula. This in turn helps us to understand the numerical results better. For example, we can use the fitting formula to precisely test some of our relations for complexity that are read off from numerical results.
In the following, we will show how to extract some important information about the evolution of complexity via numerical analysis. However, before moving to this, we shall explain some of our notations at first. Since we are particularly interested in the dependence on the eigenvalues, the complexity for the CT states will be denoted by C (α ;γ) (t). It was understood that the states have eigenvalues α(γ) on the left (right) hand side. However, sometimes it is more convenient for us to label the complexity by using the real and imaginary parts of the eigenvalues directly. We denote C (α ;γ) (t) = C (P,Q ;U,V ) (t), where P , Q , U , V take real values and are related to the eigenvalues α and γ by For example, C (1,2;0,0) (t) denotes the time-dependent complexity of the target state with α = 1 + 2i , γ = 0, and C (+,−;0,0) (t) denotes C (P,Q ;U,V ) (t) with P > 0, Q < 0, U = V = 0. We will frequently switch between these two notations. The readers should not be confused.
In addition, we choose the initial time of the evolution to be t = 0. In our numerical results, we will simply show the evolution for the time regime t ≥ 0. However, it should be emphasized that there is no difficulty to analytically continue our results to the t < 0 regime. In fact, some of our formulas that are obtained from the numerical results should be understood in this way, for example C (α * ;α * ) (t) = C (α ;α) (2π/ω −t). To avoid redundancy, we will not emphasize this again in the following.

Single excitation
Let us start with a simple case in which the target state is only excited by a single mode.
Without loss of generality, we assume γ = 0. Here, we find that it is more convenient for us to study a new quantity which describes the growth of complexity (GC) from the value at the initial time. It should not be confused with the complexity itself.
In Fig. 1, we show the evolution of the GC for single-mode excitations (in the first period) with various η and α. From the figure, we can read off a lot of interesting features.
1. First, the (minimal) time period of the GC (or the complexity) is given by 2π/ω rather than π/ω. The latter is known to be the time period for the ground thermal field double state [35]. The difference can be attributed to the non-Gaussianity of our target state, which now has a non-vanishing vector ν. By carefully examining Eq. (102), we find that the first term Tr(K 2 ) gives lowest order terms as cos 2 (ωt) or sin 2 (ωt) while the second term νg R ν T results in terms like cos(ωt) or sin(ωt). Thus, in our framework, one can generally distinguish a non-Gaussian target state from a Gaussian state using the minimal evolving time period of the complexity.
2. Secondly, the GC is sensitive to the parameters. For example, when only the real part of α mode is excited, C (+,0 ;0,0) (t) ≥ 0 provided η < 1, otherwise C (+,0 ;0,0) ≤ 0. On the contrary, if only the imaginary part of α mode is excited, the behaviour will be the opposite.
3. Thirdly, when the eigenvalue α is real or pure imaginary, the GC (or complexity) is a symmetric function about the axis t = π/ω in the first period. This implies

Interestingly, the GC for two single excitation states with conjugated eigenvalues
are symmetric about the axis t = π/ω and hence can be connected by C (α * ;0) (t) = C (α ;0) (2π/ω − t). The relation is valid to the complexity as well because of the above relation for the initial complexity.
6. Finally, from the lower three panels in Fig. 1, we conclude that there exists an identity One may notice that in Fig. 1, the red line describing the GC for the ground thermofield double state is sometimes above the other tinctorial lines. However, this does not mean that the complexity of the ground state is larger than that of an excited state, since our pictures simply show the growth of the complexity, rather than the complexity itself. As a matter of fact, the complexity for an excited state at the initial time has already been large enough, which guarantees that it costs more to prepare the state from the reference state, compared to the ground state. To show this, we plot the difference between the complexity of an excited state and the ground state directly in Fig. 2. In the upper panels, the eigenvalue α is real in the left panel and is pure imaginary in the right panel while in the lower panel, we concentrate on the results for two single excitation states with eigenvalues α and −iα. It is clear that for all these cases, the difference is always positive definite, consistent with our previous argument (108). In surprise, from the figure, we also find a translation formula for the complexity itself: 1. For general α, the complexity for two single excited states with eigenvalues α and ±iα obeys C (±iα ;0) (t) = C (α ;0) (t + π/ω) (recall that the minimal time period for the ground state is π/ω). It also implies that C (±iα * ;0) (t) = C (∓iα ;0) (2π/ω − t) = C (α ;0) (3π/ω − t). To proceed, we would like to further study the effects of the temperature on the GC. For this purpose, we fix the other parameters and let T vary and then compare the complexity for the target states at different temperatures. Some of our numerical results are shown in Fig. 3. From the figure, we conclude as follows.
2. We rediscover that the GC vanishes at t = π/ω for the eigenvalues α = P (1±i) and this point divides the image into two parts. For example, when α = P (1 + i) , P > 0, the extreme value decreases as the temperature increases in the first half period (0, π/ω) while in the other half period it behaves in the opposite way, as shown in the lower left panel.

Two excitations
In this subsubsection, we focus on the two excitation states, which have α = 0 and γ = 0. In the following, we will study three special cases: two equal excitations α = γ, two opposite excitations α = −γ and two conjugated excitations α = γ * respectively, to illustrate some main features about the evolution of complexity for two excitation states. Let us start with two equal excitations. In Fig. 4, we show the evolution of complexity for various eigenvalues α. From the figure, we find some interesting features.
2. From the upper panels, we find that similar to the single excitation case, the complexity for two excited states with conjugated eigenvalues are symmetric about the axis t = π/ω in the first period and hence C (α * ;α * ) (t) = C (α ;α) (2π/ω − t).
3. Secondly, from the lower left panel, we observe that the states with a larger absolute value of the imaginary part | α| have a higher peak value of the complexity. It may suggest that the imaginary part of α or γ gives greater influence in preparing two equal excited states. 4. In particular, from the lower right panel, we find a new interesting relation C min (iα ;iα) = C max (α ;α) , where α is real. The result implies that C (iα ;iα) (t) ≥ C (α ;α) (t ) for any time t , t . This strongly supports the above argument. b) Two opposite excitations (α = −γ).
Next, we turn our attention to the complexity for two opposite excitations. Some interesting results are shown in Fig. 5. From the figure, we conclude as follows. 1. Similar to the two equal excitation case, for general α, there exists a relation C (α ;−α) (t) = C (−α ;α) (t). Moreover, the complexity for two states with conjugated eigenvalues are symmetric about the axis t = π/ω in the first period and hence C (α * ;−α * ) (t) = C (α ;−α) (2π/ω − t).
2. Compared to the two equal excitation case, the two opposite excitation states will have a higher peak value of complexity if they have a larger absolute value of the real parts of the eigenvalues | α|, instead of the imaginary parts, see the lower left panel. In other words, the real parts of the eigenvalues α or γ give greater influence in preparing two opposite excited states. This is very different from the two equal excitation case, which is more affected by the imaginary parts of the eigenvalues. In particular, from the lower right panel, we find that when α is real, C min (α ;−α) = C max (iα ;−iα) , implying that C (α ;−α) (t) ≥ C (iα ;−iα) (t ). This strongly supports the above argument.
3. Furthermore, the last panel also tells us that when the eigenvalues are real or pure imaginary, the complexity for a two opposite excitation is related to that of a two equal excitation. For example, if α is real, then C (α ;−α) (t) = C (iα ;iα) (t + π/ω) and C (iα ;−iα) (t) = C (α ;α) (t + π/ω). The results, together with the relation C min (iα ;iα) = C max for the two equal excitation case, lead to C min c) Two conjugated excitations (α = γ * ).
At last, we close this part with the complexity for two conjugated excited states. Some numerical results are shown in Fig. 6. We conclude as follows.
3. The lower panel tells us that the maximal value of the complexity for two conjugated excited states will be the same if the eigenvalues of the states simply exchange the absolute value of the real and the imaginary parts of the eigenvalues, namely C max (α ;α * ) = C max (−α ;−α * ) = C max (±iα ;∓iα * ) . Compared with the aforementioned two types of excited states, this feature is unique. It implies that the maximal cost to prepare a two when only the real parts of the eigenvalues are excited and when only the imaginary parts are excited, where t and t are arbitrary. Compared to the analogous relations in section 4.2, one can immediately find that they are quite different. It tells us that different choices of the reference state not only affects the complexity for a single target state but also changes the relations between the complexities for different target states.

Comments on general case
One may have noticed that the above results for the several special cases show some common features. This strongly motivates us to further study the complexity for the CT states with general eigenvalues. By scanning the parameters space, we find that there indeed exist some universal relations for the complexity. Without presenting more numerical results, we summarize the relations as follows and try to explain them from a physical or mathematical point of view.
1. The first relation is It is simply a remanent of the fact that interchanging the eigenvalues between the two sides of the system does not change the state under consideration, namely |ψ CT (α , γ ; t) = |ψ CT (γ , α ; t) .

The second is the inversion relation
It strongly implies that the complexity generally depends only on the various quadratic polynomials of the real and imaginary parts of the eigenvalues α and γ. This can be easily understood from our formula (107). The complexity for the ground state C (0) does not depend on the eigenvalues. On the other hand, according to (51) and (101), one finds that the inversion (α , γ) → (−α , −γ) leads to λ → −λ and hence ν → −ν.
However, since the complexity depends quadratically on the vector ν, the result is clearly invariant.
3. The third relation we found is However, it is hard to prove since the rotation generator K is time dependent as well as the complexity for the ground state. Furthermore, we also find Again, we do not know how to prove it, without solving the generator K analytically.
Nevertheless, the above result, together with (107) tells us that under the transformation (α , γ , t) → (α * , γ * , −t), the elements of the vector ν should transform as whereK p ± ,K q ± are 2 × 2 partitioned matrixes. It follows that On the other hand, under the above transformation, we find that λ p → λ p and λ q → −λ q . Thus, symmetry considerations indicate that under the time reversal, the partitioned matrixesK p ± probably transform in an opposite way asK q ∓ , namely To check whether this is the case, we first consider the simplest case: the reference state is the Dirac vacuum. Using the results in sec.4.2, we deducê It is easy to see that under the time reversal, (K p ± ,K q ± ) → ±(K p ± ,K q ± ). This gives us strong confidence that the above argument is reasonable. Moreover, for general Gaussian reference states, we can numerically verify that the partitioned matrixes exactly transform in the same way as the Dirac vacuum case.

Complexity for quantum field theory
In this section, we move to the circuit complexity for generalised coherent states in a (1+1)dimensional free scalar field theory living on a cylinder with circumference L. Based on our experience for the two modes case, the numerical calculations for 2N modes will be straightforward but costly after we built the model of circuit complexity for the generalised coherent states in QFT. Hence, in this section, we would like to focus on analytical treatment of complexity in QFT. Along this line, we start with the Hamilton of the theory, We will regulate the field theory by a lattice model in which the lattice spacing is defined by where N is the site number of the lattice arranged on the spatial circle. For simplicity we redefine the canonical variables as and impose periodic boundary conditions, Q N +1 := Q 1 and P N +1 := P 1 . Then the Hamilton can be rewritten as With the help of the Fourier transformatioñ the Hamilton can be recast into a more compact form: where the frequency is given by The canonical commutation relations are given by In view of the symmetryQ † N −k =Q k ,P † N −k =P k and ω N −k = ω k , we would like to introduce two new canonical quantities which obey the canonical commutation relations Furthermore, in terms of these new canonical variables, the Hamilton simplifies to (152) Therefore, the system can be viewed as a sum of N independent harmonic oscillators with equal mass m = 1/2δ and frequency ω k .
The whole system we are interested in is a double copy of the free scalar theory. It is characterized by the Hamilton H L ⊗H R , as well as the canonical variables The annihilation and creation operators can be defined as They obey the commutation relations As expected, for a free field theory, the modes with different "momenta" k simply commute with each other. From this and the results for a single harmonic oscillator, it will be straightforward to construct some states of interest for the total system. For example, the TFD state can still be defined as where the anti-Hermitian operator U (β) now is given by where θ k (β) = arctanh e −βω k /2 . It turns out that the TFD state can be written as a tensor product as where |TFD k stands for the TFD state for a single harmonic oscillator. Likewise, the coherent thermal (CT) state can still be defined by (43). One finds Of course, the relation can be generalised to the time dependent states directly where |ψ CT (α k , γ k , t) describes the same state (time dependent CT state) defined in Eq.
By definition, it is easy to see that the covariance matrix for such tensor product states can be recast into a corresponding tensor plus form. For the CT state, we have where G T (|ψ CT (α k , γ k , t) ) has the same expression as Eq. (128) except that the parameters (α, γ, ω, θ) now should be replaced by (α k , γ k , ω k , θ k ).
To proceed, we would like to choose the Dirac vacua as the reference state, i.e., the ground state for the Hamilton (142). Its covariance matrix is given by where Following closely the discussions in sec.4.2, we deduce the complexity The result is simply a sum of the complexity for the N-independent harmonic oscillators It is not hard to believe that the relation is valid to more general reference states.
Alternatively, we may take the reference state to be the ground state of the ultralocal Hamilton, which is given by where µ is a constant parameter playing the role of the mass. Using the canonical variables (q k , p k ), we find Notice that the frequency is a same constant ω = 2µ for all the modes. The covariance matrix for this ground state turns out to be with Notice that G k R can be viewed as the covariance matrix (127) for a general Gaussian reference state with the parameter η k = µ/ω k . This implies that the complexity corresponding to this reference state can be obtained by Again, the result is simply a sum of the complexity C k (α k , γ k , ω k , η k ) for the N single harmonic oscillators.

Conclusions
In the inspiring paper [35], the authors first calculated the circuit complexity for timedependent TFD states using the covariance matrix approach. In the present paper, we extended their analysis to consider the complexity of the generalised coherent states in thermal field dynamics. In our case, the target state is not Gaussian anymore and has a nonvanishing one-point function.
We started from the harmonic oscillator system. Based on the construction of the Glauber coherent state and the thermal vacuum state, we introduced Coherent Thermal (CT) state and Thermal Coherent (TC) state, respectively. Also, we found they are related through Eq. (57), and the time-dependent TC state is related to the corresponding timedependent CT state by Eq. (64). Thus we found the complexity for a TC state has a simple connection with the one for the corresponding CT state, as shown in Eq. (65).
As the one-point function of the target state is not vanishing, the covariance matrix approach used in [35] cannot be directly applied in evaluating the circuit complexity for the CT state, as the symmetric two-point function is not enough to determine a CT state.
Nevertheless, by examining the properties of the CT state carefully, we developed the generalized covariance matrix approach and applied it to compute the circuit complexity. We first introduced the one-point function in Eq. (72), and from the transformation law of the 1-point and two-point functions we define the generalized covariance matrix Eq. (81).
The corresponding generators preparing our non-Gaussian state form a group structure R 2N Sp (2N , R). The essential point is that the optimal geodesic is still be generated by the horizontal generator, if the reference state is still Gaussian. After some careful analysis, we derived the circuit complexity for a CT state in Eq. (102). This expression is one of the most important results in this work, so here we write it again, The notable feature of this formula is that there is an extra piece νg R ν T contributing to the complexity.
With the formula (171), we were allowed to study the circuit complexity of the CT state.
We first chose the Dirac vacuum as the reference state, and obtained the complexity which is given in Eq. (102), Based on the formula (171). Surprisingly, we found this result is independent of time. With this formula in hand, we examined some special cases including one single excitation, two equal excitations, two opposite excitations and two conjugated excitations. From these results, we discover that the complexity with two conjugated excitations is always smaller than the one with two equal excitations or two opposite excitations fixing α. Interestingly, the complexity with two equal excitations is more sensitive to the real part of α while the complexity with two opposite excitations relies more on the imaginary part of α. In addition, we also compared the influence of single excitation and two excitations with the same α, we found which complexity in two cases is larger not only depends on the temperature, but also relies on the value of α. This fact implies, the target state with two excitations is highly entangled, and the complexity also reveals the characteristics of entanglement qualitatively.
Furthermore we considered the case that the reference state is a more general Gaussian state instead of the Dirac vacuum. The calculation of the circuit complexity is straightforward, but the details were much more complicated and were asked for numerical method.
Similarly, we looked at four different cases with the single excitation and double excitations including the equal, the opposite and the conjugated excitations. For each case, we found the function of the complexity is time-dependent with a period being 2π/ω rather than π/ω. This is due to the extra piece in (171). We presented some numerical results for the four different cases and read some interesting findings in sec. 4.3. In particular we investigated how various parameters affect the complexity. These parameters include η which reflects the initial state, the temperature T and the level of the excitations P or Q . We studied the symmetry and translation of the complexity in time. Moreover, we gave a study on the exchange symmetry and the parity symmetry of the real part and imaginary part of α = P + iQ. Similar to the discussion on the Dirac vacuum, we compared the complexities of different kinds of excitations with fixed α, especially the one of double excitations. For example, we found that to some extent the imaginary part of α or γ gave more significant influence in preparing two equal excited states with the same |α| or |γ| from the reference state, while the real parts affected more significantly on the two opposite excitations. For two conjugated excitations, there is no obvious dependence on the real parts or imaginary parts of α or γ.
In sec. 5, we briefly showed that the previous analysis can be extended to a free scalar field theory taking the Dirac vacuum as the reference state. Our study was preliminary and could be extended to different directions. We would like to study the circuit complexity in QFT taking a general Gaussian state as the reference state in the future work since there are too many parameters for the generalised coherent states. Besides, in the light of the first law of the complexity proposed in [63] and revisited in [64], it is also interesting to verify whether the first law is valid when we extent to the generalised coherent states. Another interesting topic is to generalize our study to the fermionic case.