Theory of dynamical superradiance in organic materials

2.1 Model description

We start by considering a system of N two-level systems coupled to a single cavity mode, as in Figure 1(a)i. In the dipole and rotating-wave approximations, such a system is well described by the TC Hamiltonian [33]. Setting ℏ = 1, the system Hamiltonian is

where a ̂ and a ̂ † are the annihilation and creation operators of the cavity mode of frequency ω _c, σ ̂ i z , ± are Pauli matrices describing the ith emitter satisfying [ σ ̂ i z , σ ̂ j ± ] = ± 2 δ i j σ ̂ i ± , ω ₀ is the energy splitting of the emitters, and g is the individual light–matter coupling (which depends on the mode volume). To include loss and dephasing processes, we employ a GKSL master equation for the density matrix ρ [79],

with L [ X ̂ ] = X ̂ ρ X ̂ † − { X ̂ † X ̂ , ρ } / 2 . The rate κ describes the loss of photons from the cavity, γ individual loss of molecular excitations (e.g., by spontaneous emission into a non-confined radiation mode), and γ _ϕ individual molecular dephasing. The last pure dephasing term can be understood as a Markovian approximation to the effect of vibrational modes coupling to the emitters.

To treat dynamical SR, we choose an initial state ρ(0) = |ψ(θ)⟩⟨ψ(θ)| with

where |↑⟩ _i is an excited (inverted) emitter state in the Hilbert space of the ith two-level system. The initial state has zero photons in the cavity mode, and the emitters are fully inverted and coherently tilted by an angle θ. For θ = 0, the emitters are fully inverted, giving SF, whereas for θ ≠ 0, there is some initial dipole moment present, giving SR [Figure 1(b)]. Note that a rotation of the inverted states about any axis in the xy plane (i.e., exp − i n x σ ̂ i x + n y σ ̂ i y θ / 2 ; n x 2 + n y 2 = 1 ) leads to an emitter state with equivalent initial coherence as in Eq. (4), and we choose here the x-axis for simplicity (i.e., n _x = 1, n _y = 0). Such states of two-level systems are the analog of coherent states in bosonic systems, which are also referred to as spin-coherent states or atomic-coherent states [80], [81], and they have recently been used to study entanglement in Dicke SR [82].

2.2 Numerical method

To solve Eq. (2) numerically, we first make use of the weak permutation symmetry of the emitters [61], [62], [63], [64], [65], [66], [67] by expanding the density matrix as

where |n⟩⟨n′| is the state of the cavity mode with photon numbers n and n′, and O ̂ λ is a permutation-symmetric density matrix depending on the N-emitter state labeled by λ. Since we consider two-level systems, the permutation-symmetric states can be labeled by enumerating the ways of dividing N elements into four bins corresponding to the possible emitter matrix elements of a single two-level system ↑↑, ↓↑, ↑↓, and ↓↓ (↑: excited state, ↓: ground state). Thus, λ describes the 4-element lists m _λ = (m _λ,↑↑, m _λ,↓↑, m _λ,↑↓, m _λ,↓↓) with m _λ,p ≥ 0 and ∑ _p m _λ,p = N. Explicitly, each m _λ corresponds to a matrix O ̂ λ that consists of the sum of all states obtainable by permuting pairs of emitters in a state with m _λ,↑↑ inverted emitters in the left and right states of the density matrix, followed by m _λ,↓↑ ground-state emitters in the left state and inverted emitters in the right state, and so on.

Using the mapping described above, one can, in principle, construct the operators O ̂ λ in the original explicit 2 ^N dimensional two-level-system basis. The permutation approach is, however, based on avoiding the explicit representation and performing all calculations within the compressed space. The number of unique density-matrix elements in the emitter space is given by the number of possible partitions m _λ , which is equal to N + 3 3 .

Since the initial state in Eq. (4) has N excitations in the emitter state and there are no gain terms in Eq. (2), the photon Hilbert space can be truncated at N + 1 without introducing any approximations. The total number of density-matrix elements, therefore, scales as ( N + 1 ) 2 × N + 3 3 , and an exact solution of Eq. (2) becomes possible up to order of N ∼ 30 emitters.

We can further exploit the weak [83], [84] U(1) symmetry of the master equation (2), which states that it is invariant under a ̂ → a ̂ e i ϕ , σ ̂ n − → σ ̂ n − e i ϕ . Let us denote the total number of excitations in the left and right states of a given density-matrix element ρ _λ,n,n′ as ν = n + m _λ,↑↑ + m _λ,↑↓ and ν′ = n′ + m _λ,↑↑ + m _λ,↓↑. As a consequence of the weak U(1) symmetry, only elements with a constant value of ν − ν′ can couple in the dynamics. This motivates the introduction of the block form of the density matrix

where ρ ^{(ν, ν′)} denotes the collection of all elements with left and right excitation numbers equal to ν and ν′, respectively. The most general equations of motion for the individual blocks respecting the weak U(1) symmetry then read

where L k ( ν , ν ′ ) are matrices determined by the master equation. All terms with k < 0 describe gain processes, k = 0 describes the coherent dynamics and dephasing, and k > 0 describes losses. Note that blocks with different values of ν and ν′ do couple, but their difference ν − ν′ is constant. Such a structure of master equation – combining weak U(1) symmetry with weak permutation symmetry – was used in Ref. [68] to discuss steady state superradiance in a system of Rydberg atoms.

Since Eq. (2) features jump operators that change the total excitation number by at most one, and the system we consider has no gain, only the terms k = 0 and k = 1 are nonzero in Eq. (7). Further, we are primarily interested in calculating the expectation value of the number of photons in the cavity mode ⟨ n ̂ ⟩ = ⟨ a ̂ † a ̂ ⟩ , so we can restrict Eq. (7) to only those blocks with ν − ν′ = 0 on which the photon number operator depends. In this case, we arrive at the evolution equation for the blocks ρ ^(ν) ≡ ρ ^(ν,ν)

The matrices L 0 ( ν ) and L 1 ( ν ) can be constructed from Eq. (2). Crucially, since we consider a system without gain, we can solve Eq. (8) sequentially by noting that there exists a block with the highest excitation number ν _max. For this block, the second term on the right-hand side in Eq. (8) vanishes, resulting in a homogeneous system of differential equations for the block ρ ( ν max ) , which can be numerically integrated. The solution ρ ( ν max ) may then be used to solve the block ρ ( ν max − 1 ) . Iterating until ν = 0 yields dynamics for all blocks with ν − ν′ = 0. We note that this method of solution is reminiscent of that used to solve the diagonal elements of the density matrix in Dicke SR [85].

From the initial state in Eq. (4), we see that the maximum excitation number is ν _max = N. The number of elements of the blocks ρ ^(ν) for a fixed N increases monotonically from 1 at ν = 0 to N + 3 3 at ν = N. The largest block must have exactly as many elements as there are emitter states, because for every emitter state, one can always find photon numbers n and n′ such that ν = ν′ = N.

Code for the numerical scheme described is available in the PIBS (Permutational Invariance Block Solver) Python package [86]. This method allows one to solve Eq. (2) with initial conditions given in Eq. (4) up to N = 140 emitters. Currently, only the solution for the diagonal block ν = ν′ is implemented in the PIBS code, but an extension to off-diagonal blocks is straightforward.

2.3 Benchmarking Cumulant approaches

In this section, we compare the exact solution of Eq. (2) based on PIBS to approximate first- and second-order cumulant approaches. In a cumulant approximation of order M, one splits expectation values of M + 1 operators in the Heisenberg equations of motion into expectation values of at most M operators by setting the cumulants of order M + 1 to zero. For example, in a first-order cumulant approximation, the expectation values of any two operators A ̂ and B ̂ are split as ⟨ A ̂ B ̂ ⟩ = ⟨ A ̂ ⟩ ⟨ B ̂ ⟩ , because the second-order cumulant given by ⟨ ⟨ A ̂ B ̂ ⟩ ⟩ = ⟨ A ̂ B ̂ ⟩ − ⟨ A ̂ ⟩ ⟨ B ̂ ⟩ [87] is set to zero. This is equivalent to mean-field theory. Going one order higher, keeping expectations of pairs of operators and setting third-order cumulants to zero, one arrives at the second-order cumulant approximation.

If the initial state respects the U(1) symmetry of the Hamiltonian in Eq. (1), one can simplify the second-order cumulant equations by dropping expectation values of operators that do not conserve the number of excitations, such as ⟨ a ̂ ⟩ and ⟨ σ ̂ i + ⟩ . If these “symmetry-breaking” terms are zero initially, then they remain zero for all times. The resulting equations may be called symmetry-preserving cumulant equations [73].

The initial state in Eq. (4) breaks the U(1) symmetry for θ ≠ 0, since ⟨ σ ̂ i + ⟩ = i ⁡ sin ( θ ) / 2 . As a consequence, for SR initial conditions, one must retain all terms in the cumulant expansion, resulting in symmetry-broken cumulant equations. We use the Julia package QuantumCumulants.jl [88] to calculate and evolve the symmetry-broken cumulant equations.

In Figure 2, we compare the cumulant approaches to the exact solution for different numbers of emitters and initial conditions, by plotting the normalized average number of photons in the cavity ⟨ n ̂ ⟩ / N = ⟨ a ̂ † a ̂ ⟩ / N as a function of time. For brevity, we refer to the first-order cumulant approximation as mean field and to the second-order cumulant approximation simply as cumulants. We do not show a mean-field solution for SF initial conditions, as the mean-field solution in that case is just zero for all times if we start from the symmetric initial state. We also note that in the mean-field treatment, ⟨ n ̂ ⟩ / N = | ⟨ a ̂ ⟩ | 2 / N becomes independent of N.

Figure 2:

Comparison of exact dynamics of the average number of photons in the cavity calculated with PIBS and cumulant approaches. Panels (a) and (c) show the dynamics for SF initial conditions for N = 10 and N = 140. Panels (b) and (d) show the dynamics for SR initial conditions (θ = π/4) for N = 10 and N = 100. Other parameters: g N = 0.4 e V , Δ = ω c − ω 0 = − 0.35 e V , κ = 0.01 e V , γ = 0.001 e V , γ ϕ = 0.0075 e V .

Figure 2(a) and (c) shows the dynamics for N = 10 and N = 140 for SF initial conditions (θ = 0). One can see that cumulants underestimate the dephasing processes, as the amplitude of oscillations is much larger than that in the exact solution. As such, for N = 140, we only see agreement between the cumulant approximation and the exact solution at early times, specifically for t < 10 fs. At later times, there is no agreement over the range of N for which we can obtain exact results. Part of the dephasing (which is different from the dephasing due to the vibrations) causing this mismatch can be explained by the quantum corrections to the energy eigenvalues of the TC Hamiltonian for finite N [34]. As discussed there, the mean-field (classical) solution with periodic oscillations corresponds to the case of equal (harmonic) energy-level spacing. At finite N, however, the level spacing is not perfectly regular. Using the fact that the resulting anharmonicity is known to scale as 1/(ln N)³, and that the cumulant solution for N = 140 matches the exact solution up to t ≈ 10 fs [cf. Figure 2(c)], we can estimate that, e.g., for N = 10⁸, this dephasing would set in after t ≈ 520 fs.

Figure 2(b) and (d) compares the exact solution to cumulants and mean field for an SR initial state with θ = π/4. Generally, the cumulant solution matches the exact solution better than the mean-field solution. Moreover, for increasing N, the early time period where cumulants match the exact solution increases. This indicates that for SR initial states, symmetry-broken cumulants capture the early-time behavior well. We provide quantitative evidence that this is the case in Appendix A, where we show that both mean-field and cumulants converge monotonically toward the exact solution for N up to 100.

As an additional test of mean-field results, we note that for SR initial states, the symmetry-broken cumulant solution converges to the mean-field solution for large enough N, as demonstrated in Figure 3. Since we are ultimately interested in the SR initial state, this leads us to conclude that for large N, it is reasonable to use mean-field theory.

Figure 3:

Comparison of symmetry-broken second-order cumulant approximation to mean-field approximation for (a) N = 100 [same as Figure 2(d)], (b) N = 1,000, and (c) N = 10,000. The SR angle is set to θ = π/4, all other parameters are the same as in Figure 2.

3.1 Model description

In the previous section, we discussed the dissipative TC model with local dephasing, capturing the emitters vibronic coupling within a Markovian approximation. We now introduce a model that includes vibrations explicitly in the system dynamics. As a simplest case, we assume that there is one dominant vibrational mode of frequency ω _ν , which is the same for all molecules. We extend the TC model from Eq. (1) to include this mode,

which is the HTC model [47], [56].

The vibrational mode of the ith molecule is treated as a harmonic oscillator with annihilation and creation operators b ̂ i and b ̂ i † . The coupling strength between the vibrations and the electronic degrees of freedom is set by the Huang–Rhys parameter S. Including the same dissipation processes as before (see Eq. (3) [89]), and adding thermalization processes for the vibrational modes, the master equation is

where D TC was defined in Eq. (3). The thermalization rates are γ _↑ = γ _ν n _B and γ _↓ = γ _ν (n _B + 1), where n B = e ω ν / T − 1 − 1 is the Bose–Einstein distribution at temperature T (we set k _B = 1). The model of thermalization here, which includes a Lamb-shift term H ̂ LS = i ( γ ν / 4 ) ∑ i = 1 N b ̂ i † b ̂ i † − b ̂ i b ̂ i , is one of momentum damping, as discussed in Ref. [90]. The energy of the HTC Hamiltonian (9) depends on the vibrational displacement x ̂ i ∝ b ̂ i † + b ̂ i , and the thermalization drives the system toward x ̂ i such that the energy is minimized. In Figure 1(a)ii, the internal level structure and incoherent processes of the HTC model are shown.

Using Eq. (4), we can write the initial state including vibrations as

where ρ vib , i = 1 − e − ω ν / T e − ω ν b ̂ i † b ̂ i / T is a thermal state of the vibrational mode of molecule i.

The Hilbert space of the HTC model is much larger than the one considered in Section 2. Truncating the vibrational mode to N _ν levels increases the Hilbert space by a factor of N ν N , which makes an exact solution quickly intractable for growing N. However, from Section 2, we have confidence in mean-field theory capturing the essential early-time behavior also for the HTC model and, therefore, solve the dynamics using the mean-field approximation.

In what follows, we compare the dynamics of the model Eq. (10) for two different cases [cf. Figure 1(a)]: (i) S = 0 and variable γ _ϕ , which just reduces to the TC model Eqs. (1) and (2), and (ii) γ _ϕ = 0 and variable S. We do this comparison because both S and γ _ϕ model forms of dephasing by phonons, and we are interested in how these different processes in solid-state environments modify superradiance.

3.2 Results

The results shown here were calculated in the mean-field approximation using the Julia package QuantumCumulants.jl [88]. We consider SR initial conditions with θ = 10⁻³ π throughout, and set the number of molecules N = 10⁸.

Figure 4(a) shows the average number of photons in the cavity mode as a function of time for different values of S. For S = 0, the HTC model reduces to the TC model, resulting in damped Rabi-oscillations (note the logarithmic scale). For S ≤ 0.3, after a build-up time where ⟨ n ̂ ⟩ / N grows from 0 up to roughly 10⁻⁶, there is a region where the photon number increases exponentially with time as ⟨ n ̂ ⟩ ∝ e t / τ . The parameter τ is the characteristic risetime, which is seen to increase with S. Increasing S beyond a value of 0.3 changes the exponential behavior after the build-up time, and a characterization based on a single risetime τ is no longer meaningful. It is, however, clear that a stronger coupling to the vibrational modes inhibits the emission of photons into the cavity mode.

Figure 4:

Mean-field solution of the HTC model starting from SR initial conditions. (a), (b), and (c) Average number of photons in the cavity mode, electronic coherence, and magnitude of the Bloch vector, respectively, as a function of time for different coupling strengths to the vibrational mode. (c), (d), and (f) Same as (a), (b), and (c) but with the vibrational mode replaced by a pure dephasing term with rate γ _ϕ . The curves S = 0 in (a), (b), and (c) are the same as the curves with γ _ϕ = 0 in (d), (e), and (f). Note that the values of S and γ _ϕ are not equidistantly spaced. Parameters: N = 1 0 8 , g N = 0.2 e V , Δ = 0 , κ = 0.01 e V , γ = 1 0 − 6 e V , θ = 1 0 − 3 π , ω ν = 0.15 e V , γ ν = 0.01 e V , T = 0.026 e V . For (a), (b), and (c): γ _ϕ = 0; for (d), (e), and (f): S = 0.

In Figure 4(d), the vibrational mode has been replaced by a pure Markovian dephasing term with rate γ _ϕ . For all values of γ _ϕ , there is a region where the photon number increases exponentially, and the risetime increases with increasing γ _ϕ . Thus, the pure dephasing model can capture the photon dynamics for small S below 0.3 but breaks down for larger values of S.

To characterize the extent to which the vibrational modes destroy superradiance, we additionally calculate the electronic coherence | ⟨ σ ̂ + ⟩ | ≡ | ⟨ σ ̂ i + ⟩ | for the same sets of parameters used in Figure 4(a) and (d), see Figure 4(b) and (e). Since we use SR initial conditions, the coherence at t = 0 is nonzero and given by | ⟨ σ ̂ + ⟩ | = sin ( θ ) / 2 . Figure 4(b) shows that the coherence increases exponentially at early times for small S, similar to the exponential rise in the photon number in Figure 4(a). For S ≳ 0.3, the exponential increase of the coherence is inhibited, and there are regions where the coherence decreases before it has reached its maximum. In Appendix B, we show that increasing the initial amount of coherence (increasing θ) tends to make the initial rise of photon number faster and less dependent on the vibrational coupling, and in Appendix C, we show that the threshold of SR at a certain value of S (here at S ≈ 0.3) emerges as a competition between the light–matter coupling g N and the vibronic coupling ω ν S .

In Figure 4(e), where the vibrational mode has been replaced by a pure dephasing term, we see an exponential increase of the coherence for low dephasing rates and early times, similar to the behavior in Figure 4(b) for small S. Large dephasing rates also destroy superradiance, as can be seen by the initial decrease of coherence for γ _ϕ ≳ 0.06 eV. We note that the long-time behavior of the coherence is quite different in Figure 4(b) and (e). In the former case, the coherence builds up eventually and remains relatively large, whereas in the latter case, the coherence decreases. This shows that, when regarding ⟨ σ ̂ ⁺⟩, the coupling to vibrations does not result in the same kind of decoherence as introduced by a pure dephasing term, even for small S.

Lastly, Figure 4(c) and (f) shows the magnitude of the Bloch vector ⟨ J ̂ 2 ⟩ = ⟨ J ̂ x 2 + J ̂ y 2 + J ̂ z 2 ⟩ , J ̂ α = ∑ i = 1 N σ ̂ i α / 2 , which, within mean-field theory and using N ≫ 1, is given by

The initial state in Eq. (4) is an eigenstate of J ̂ 2 with the maximum possible eigenvalue N/2(N/2 + 1) ≈ N ²/4. In the absence of incoherent processes and vibronic coupling, i.e., for the closed TC model in Eq. (1), the emitter state remains an eigenstate of J ̂ 2 with the same eigenvalue since [ H ̂ TC , J ̂ 2 ] = 0 and, therefore, ⟨ J ̂ 2 ⟩ ≈ N 2 / 4 for all times. This is seen for the curves with S = 0 or γ _ϕ = 0 in Figure 4(c) and (f). We note that the effect of local molecular losses proportional to γ = 10⁻⁶ eV is negligible in the time period considered here, which would otherwise lead to a decrease in ⟨ J ̂ 2 ⟩ .

For nonzero S or γ _ϕ , the magnitude of the Bloch vector does not retain its maximum initial value. In the case of vibrational dressing in Figure 4(c), ⟨ J ̂ 2 ⟩ is subject to oscillations but stays well above zero in the considered time frame. In contrast, for pure dephasing in Figure 4(f), ⟨ J ̂ 2 ⟩ decreases monotonically to zero, with minor bumps for small dephasing rates, showing a stark difference between the two models of treating the vibrational effects.

Such behavior can be related to discussions of how any processes that break spin conservation destroy free-space superradiance [91], [92], [93]. As noted in those works, when the dynamics approaches the equator – specifically when | ⟨ σ ^ z ⟩ | ≲ 1 / N – simple counting arguments show there are many more accessible subradiant states (quantum states with ⟨ J ̂ 2 ⟩ ≪ N ) than superradiant states. As such, any process that can change the modulus of J ̂ will lead to transitions to these subradiant states, reducing ⟨ J ̂ 2 ⟩ and ultimately destroying the superradiance. Notably, these processes only become relevant as one approaches the equator: when ⟨ σ ̂ ^z ⟩≃ 1, scattering to subradiant states is suppressed if the scattering does not directly change ⟨ σ ̂ ^z ⟩. References [91], [92], [93] discuss how direct interactions between different emitters, breaking the permutation symmetry, have such an effect. Here, we see similar behavior driven by individual dephasing or vibrational modes.

Interestingly, for both models, a stronger influence of the vibrational modes (i.e., larger S or larger γ _ϕ ) leads to an increased time period, for which ⟨ J ̂ 2 ⟩ stays close to its maximum initial value. This can be explained by the decreased emission rate of photons for larger S and γ _ϕ , see Figure 4(a) and (d), which means that the excitation state remains close to its initial value of ⟨ σ ̂ z ⟩ ≈ 1 for longer times and thus avoids states near the equator where transitions to subradiant states arise. According to Eq. (12), this results in ⟨ J ̂ 2 ⟩ ≈ N 2 / 4 .

In experiments, the cavity frequency can typically be readily varied; our calculations predict interesting phenomena when varying the cavity detuning Δ = ω _c − ω ₀, with a clear difference between pure dephasing and having coupling to a vibrational mode. The dynamics of the photon mode in the HTC model is plotted for Δ = −0.15 eV in Figure 5(a) and for Δ = −0.3 eV in Figure 5(b). Compared to the case of zero detuning in Figure 4(a), we see that a negative detuning favors the emission of photons for larger values of S. This is seen most clearly in Figure 5(b), where the curves up to S = 0.4 almost collapse into one curve at early times. In the inset of Figure 5(b), we observe that the photon number rises faster for 0.1 ≤ S ≤ 0.4 than for S = 0. This can be explained by first noting that, for ω _ν = 0.15 eV and T = 0.026 eV used here, the phonons are almost fully in the ground state at t = 0. Further, for a negative detuning, the cavity frequency ω _c = ω ₀ + Δ is smaller than the level splitting of the electronic ground and excited ω ₀. Thus, the cavity frequency matches transitions from the excited state manifold with zero vibrational excitations to the ground-state manifold with vibrational excitations. As S increases, the most probable number of vibrational excitations created in an optical transition increases. As such, one sees an optimal value of S in this regime where processes exciting one vibrational mode become dominant. In contrast, if the vibrational mode is treated as a pure dephasing term, such a reordering phenomenon does not occur.

Figure 5:

Average number of photons in the cavity mode as a function of time for different coupling strengths to the vibrational mode for (a) Δ = −0.15 eV and (b) Δ = −0.3 eV. All other parameters are the same as in Figure 4, and we specifically note that the vibrational frequency is set to ω _ν = 0.15 eV.

For a final set of results, we calculate the risetime τ for different values of S, Δ, and γ _ϕ , while making sure to stay in a region where τ is well defined. Figure 6(a) shows τ as a function of S and different values of Δ. For small negative detunings, τ increases monotonically with S, as evident in Figure 4(a). In contrast, large negative detunings Δ ≤ −0.25 eV show regions where increasing S leads to a shorter risetime, which is the behavior observed in Figure 5. This behavior is confirmed in Figure 6(b), which shows τ as a function of Δ for different S. Here, we see in particular how the minimum of τ at zero detuning for the TC model (S = 0) shifts to more negative detunings as S > 0 is increased, and further, the symmetry about the minimum value of τ is lost. This symmetry remains in the case where the vibrational mode is replaced by a pure dephasing rate, irrespective of the strength of the dephasing, see Figure 6(c). In that case, the risetime increases as γ _ϕ increases, and the minimum value is always at zero detuning.

Figure 6:

Photon number risetime extracted from an exponential fit to the linear regions as seen for example in Figure 4(a) and (c). (a) Risetime as a function of S for different detunings. (b) Risetime as a function of detuning for different values of S. (c) Same as (b), but the vibrational mode has been replaced by a pure dephasing term with rate γ _ϕ . Note that the curve with S = 0 in (b) is the same as the curve with γ _ϕ = 0 in (c). All other parameters are the same as in Figure 4.

The risetime behavior for different values of θ is discussed in Appendix B. We find that there is an initial rapid rise to a higher photon number for large θ than for a small one. That is, the photon number rise is faster for a smaller amount of electronic excitation (smaller inversion). This is a curious behavior, opposite, for example, to usual pulsed lasing, where the pulse build-up time becomes shorter for a higher amount of excitation [55]. This is because coherence plays a crucial role in SR, and larger θ (up to π/2) means higher initial coherence.