On polarization parameters of spin-$1$ particles and anomalous couplings in $e^+e^-\to ZZ/Z\gamma$

We study the anomalous trilinear gauge couplings of $Z$ and $\gamma$ using a complete set of polarization asymmetries for $Z$ boson in $e^+e^-\to ZZ/Z\gamma$ processes with unpolarized initial beams. We use these polarization asymmetries, along with the cross-section, to obtain a simultaneous limit on all the anomalous couplings using Markov Chain Monte Carlo (MCMC) method. For an $e^+e^-$ collider running at $500$ GeV center-of-mass energy and $100$ fb$^{-1}$ of integrated luminosity the simultaneous limits on the anomalous couplings are $1\sim3\times 10^{-3}$.


Introduction
The Standard Model (SM) of particle physics is a well established theory and its particle spectrum has been completed with the discovery of the Higgs boson [1] at the Large Hadron Collider (LHC). The predictions of the SM are being confirmed from time to time in various colliders with great success and yet phenomena such as CP-violation, neutrino oscillation, baryogenesis, dark matter etc., require one to look beyond the SM. Most of the beyond SM (BSM) models needs either new particles or new couplings among SM particles or both. Often this leads to a modified electroweak sector with modified couplings. To test the SM (or BSM) predictions for the electro-weak symmetry breaking (EWSB) mechanism one needs precise measurements of the strength and tensorial structure of the Higgs (H) coupling with all other gauge boson (W ± , γ, Z), Higgs self couplings and couplings among gauge boson themselves.
In this work, we focus on the precise measurement of trilinear gauge boson couplings, in a model independent way, at proposed International Linear Collider (ILC) [2,3]. The possible trilinear gauge boson interactions in electro-weak a email:rr13rs033@iiserkol.ac.in b email:ritesh.singh@iiserkol.ac.in (EW) theory are WW Z, WW γ, ZZγ, ZZZ, γγZ and γγγ, out of which the SM, after EWSB, provides only WW Z and WW γ self couplings. Other interactions among neutral gauge boson are not possible at the tree level in the SM and hence they are anomalous. Thus any deviation from the SM prediction, either in strength or the tensorial structure, would be a signal of the BSM physics. There are two ways to study these anomalous couplings in a model independent way: First way is to write down an effective vertex using the most general set of tensorial structures for it satisfying Lorentz invariant, U(1) em invariant and Bose symmetry weighted by corresponding CP-even/odd form-factors [4][5][6]. There is, a priory, no relation between various form factors. The other way to study anomalous couplings is to add a set of higher dimension operators, invariant under (say) SU(2) L ⊗ U(1) Y [7], to the SM Lagrangian and then obtain the effective vertices with anomalous couplings after EWSB. This method has been used to study anomalous triple neutral gauge boson coupling with the SM gauge group [8][9][10].
On the experimental side, the anomalous Lagrangian in Eq. (1) have been explored at the LEP [32][33][34][35][36], the Tevatron [37][38][39] and the LHC [40][41][42][43][44]. The tightest bounds on f V i (i = 4, 5) [40] and on h V j ( j = 3, 4) [44] comes from the CMS collaboration (see Table 1). For the ZZ process the total rate has been used [40] while for the Zγ process both the cross-section and p T distribution of γ has been used [41][42][43][44] for obtaining the limits. All these analyses vary one parameter at a time to find the 95% confidence limits on the formfactors. For the Zγ process the limits on the CP-odd formfactors, h V 1 , h V 2 , are comparable to the limits on the CP-even form-factors, h V 3 , h V 4 , respectively. To put simultaneous limits on all the form-factors one would need as many observables as possible, like differential rates, kinematic asymmetries etc. A set of asymmetries with respect to the initial beam polarizations has been considered [11, 13-16, 21, 23-25] at the e + e − and eγ colliders. Refs. [11,15,16,23] also include CP-odd asymmetries made Table 1 List of tightest limits on anomalous couplings of Eq. (1) available in literature.

Coupling
Limits Remark Ref. [44] out of the initial beam polarizations alone, which are instrumental in putting strong limits on CP-odd form-factors. Additionally, Ref. [11] includes the polarization of produced Z as well for forming the asymmetries. The latter does not necessarily require initial beams to be polarized and can be generalized to the hadron colliders as well where such beam polarization may not be available.
To this end, we discuss angular asymmetries in colliders corresponding to different polarizations of Z boson in particular and any spin-1 particle in general. For a spin-s particle, the polarization density matrix is a (2s+1)×(2s+1) hermitian, unit-trace matrix that can be parametrized by 4s(s + 1) real parameters. These parameters are different kind of polarization. For example, a spin-1/2 fermion has 3 polarization parameters called longitudinal, transverse and normal polarizations (see for example [45,46]). Similarly, for a spin-1 particle we have a total of 8 such parameters. Three of them are vectorial like in the spin-1/2 case and other five are tensorial [46,49] as will be described in Section 2 for completeness. Eight polarization parameters for a massive spin-1 particle has been discussed earlier in the context of anomalous trilinear gauge coupling [47,48], for the spin measurements studies [46] and to study processes involving W ± bosons [49]. In this work we will investigate all anomalous couplings (upto dimension-6 operators) of Eq. (1) in the processes e + e − → ZZ/Zγ with the help of total cross section and 8-polarization asymmetry of final state Z boson.
The plan of this paper is follows. In Section 2 we discuss the polarization observables of a spin-1 boson in detail using the language of polarization density matrices. Section 3 has brief discussion about the anomalous Lagrangian and corresponding off-shell vertices and required on-shell vertices. In Section 4 we study the sensitivity of the polarization asymmetries to the anomalous couplings and in Section 5 we perform a likelihood mapping of the full coupling space using Markov Chain Monte Carlo (MCMC) method for two benchmark points along with a likelihood ratio based hypothesis testing to resolve the two benchmark points. We conclude in Section 6.

Polarization observables for spin-1 Boson
In this section, we discuss the complete set of polarization observables of a spin-1 particle in a generic production process and the method to extract them from the distribution of its decay products.  Fig. 1 Feynman diagram for a general process for production of spin-1 particle (Z boson) and its decay to a pair of fermions.

The production and decay density matrices
We consider a generic reaction A B → V C D E ..., where particle V is massive and has spin-1, see Fig. 1. The production density matrix for particle V can be written as, Here, M (λ ) is the helicity amplitude for the production of V with helicity λ ∈ {−1, 0, 1}, while the helicities of all the other particles (A, B,C, D, E ...) are suppressed. The differential cross section for the production of V would be given by, Here Ω V is the solid angle of the particle V , while the phase space corresponding to all other final state particle is integrated out. The production density matrix can be written in terms of polarization density matrix P(λ , λ ) as [46], Here the matrix P(λ , λ ) is a 3 × 3 hermitian matrix that can be parametrized in terms of a vector p = (p x , p y , p z ) and a symmetric, traceless, second-rank tensor T i j as, Here θ , φ are the polar and the azimuthal orientation of the final state fermion f , in the rest frame of V with its would be momentum along z-direction. The parameters α and δ are given by, δ = with Combining the production and decay density matrices, the total rate for the process shown in Fig. 1, with particle V being on-shell, is given by [46], where σ = σ V BR(V → ff ) is the total cross section for production of V including its decay and s = 1 being spin of the particle. Using Eqs. (5) and (6) in Eq. (9), the angular distribution for the fermion f becomes: Using partial integration of this differential distribution of the final state f one can construct several asymmetries to probe various polarization parameters in Eq. (5) which will be discussed in the next section.

Estimation of Polarization parameters
For the process of interest, one might want to calculate the polarization parameters of the particle V . This can be achieved at two levels: At production process level and at the level of decay products. At the production level we first calculate the production density matrix, Eq. (2), using helicity amplitudes of the production process and then calculate polarization matrix, Eq. (5). The polarization parameters can be extracted from polarization matrix elements as:  All other polarization parameters can be obtained in similar manner, using Eq. (10), as: We note that pure azimuthal asymmetries (A x , A y , A xy , A x 2 −y 2 ) listed above are already discussed in Ref. [46] however a complete set of asymmetries for W ± -boson (δ = 0, α = −1) is given in Ref. [49]. While extracting the polarization asymmetries in the collider/event generator we have to make sure that the analysis is done in the rest frame of V . The initial beam defines the z-axis in the lab while the production plane of V defines the xz plane, i.e. φ = 0 plane. While boosting to the rest frame of V we keep the xz plane invariant. The polar and the azimuthal angle of the decay products of V are measured with respect to the would be momentum of V .
As a demonstration of the two methods mentioned above, we look at two processes: e + e − → ZZ and e + e − → Zγ. The polarization parameters are constructed both at the production level, using Eq. (11), and at the decay level, using Eqs. (12) and (13). We observe that out of eight polarization asymmetries only three, A x , A x 2 −y 2 and A zz , are nonzero in the SM. The asymmetries A x 2 −y 2 and A zz are calculated analytically for the production part and shown as a function of beam energy in Fig. 2 with solid lines. For the same processes with Z → ff decay, we generate events using MadGraph5 [50] with different values of beam energies. The polarization asymmetries were constructed from these events and are shown in Fig. 2 with points. The statistical error bars shown corresponds to 10 4 events. We observe an agreement between the asymmetries calculated at the production level (analytically) and the decay level (using event generator). Any possible new physics in the production process of Z boson is expected to change the cross-section, kinematical distributions and the values of the polarization parameters/asymmetries. We intend to use these asymmetries to probe the standard and beyond standard physics.

Anomalous Lagrangian and their probes
The effective Lagrangian for the anomalous trilinear gauge boson interactions in the neutral sector is given in Eq. (1), which include both dimension-6 and dimension-8 operators as found in the literature. For the present work we restrict our analysis to dimension-6 operators only. Thus, the anoma-lous Lagrangian of our interest is: This yields anomalous vertices ZZZ through f Z 4,5 couplings, γZZ through f couplings. There is no γγγ vertex in the above Lagrangian. We use FeynRules [51] to obtain the vertex tensors and they Fig. 3 Feynman diagram for a general anomalous triple gauge boson vertex with V = Z/γ. are given as: Notations for momentum and Lorentz indices are shown in the Fig. 3. We are interested in possible trilinear gauge boson vertices appearing in the processes e + e − → ZZ and e + e − → Zγ with final state gauge bosons being on-shell. For the process e + e − → ZZ, the vertices γ ZZ and Z ZZ appear with on-shell conditions k 2 1 = k 2 2 = M 2 Z . The terms proportional to k ν 1 and k σ 2 in Eqs. (15) and (16) vanish due to the transversity of the polarization states. Thus, in the onshell case, vertices for e + e − → ZZ reduce to: For the process e + e − → Zγ the vertices γZ Z and γ γZ appear with corresponding on-shell and polarization transversity conditions. Putting these conditions in Eqs. (15) and (17) and some relabelling of momenta etc. in Eq. (15) the relevant vertices Z γZ and γ γZ can be represented together as, The off-shell V is the propagator in our processes and couples to the massless electron current, as shown in Fig. 4(c) and Fig. 4(d). After above said reduction of the vertices, there were some terms proportional to q µ that yields zero upon contraction with the electron current, hence dropped from the above expressions. We note that although h Z and f γ appears together in the off-shell vertex of γZZ in Eq. (15), they decouple after choosing separate processes; the f V appear only in e + e − → ZZ while h V appear only in e + e − → Zγ. This decoupling simplifies our analysis as we can study two processes independent of each other when we perform a global fit to the parameters in Section 5.

Asymmetries, limits and sensitivity to anomalous couplings
In this section we thoroughly investigate the effects of anomalous couplings in the processes e + e − → ZZ and e + e − → Zγ. We use tree level SM interactions along with anomalous couplings shown in Eq. (14) for our analysis. The Feynman diagrams for these processes are given in the Fig. 4 where the anomalous vertices are shown as big blobs. The helicity amplitudes for anomalous part together with SM contributions for both these processes are given in Appendix A. These are then used to calculate the polarization observables and the cross section which are given in the Appendix B.

Parametric dependence of asymmetries on anomalous couplings
The dependence of the observables on the anomalous couplings for the ZZ and Zγ processes are given in Table 2 and Table 3, respectively. In the SM, the helicity amplitudes are real, thus the production density matrix elements in Eq. (2) are all real. This implies A y , A xy and A yz are all zero in the SM, see Eq. (11). The asymmetries A z and A xz are also zero for SM couplings due to the forward-backward symmetry of the Z boson in the c.m. frame, owing to the presence of both tand u-channel Table 3 Dependence of polarization observables on the anomalous coupling in Zγ final state.

Observables
Linear terms Quadratic terms diagrams and unpolarized initial beams. After including anomalous couplings A y and A xy receives non-zero contribution, while A z , A xz and A yz remain zero for the unpolarized initial beams.
From the list of non-vanishing asymmetries, only A y and A xy are CP-odd, while others are CP-even. All the CP-odd observables are linearly dependent upon the CP-odd couplings, like f V 4 and h V 1 , while all the CP-even observables have only quadratic dependence on the CP-odd couplings. In the SM, the Z boson's couplings respect CP symmetry thus A y and A xy vanish. Hence, any significant deviation of A y and A xy from zero at the collider will indicate a clear sign of CPviolating new physics interactions. Observables, that have only linear dependence on the anomalous couplings yield a single interval limits on these couplings. On the other hand, any quadratic appearance (like ( f V 5 ) 2 in σ ) may yield more than one intervals of the couplings while putting limits. For the case of ZZ process , A x and A y do not have any quadratic dependence, hence they yield the cleanest limits on the CPeven and -odd parameters, respectively. Similarly, for the Zγ process we have A y , A xy and A x 2 −y 2 that have only linear dependence and provide clean limits. These clean limits, however, may not be the strongest limits as we will see in the following sections.

Sensitivity and limits on anomalous couplings
Sensitivity of an observable O dependent on parameter f is defined as, , the error is given by, where, N + + N − = N T = Lσ , L being the integrated luminosity of the collider. While the error in cross section σ will Fig. 4 Feynman diagrams for the production of ZZ or Zγ at e + e − collider .
be given by, Here ε A and ε are the systematic fractional error in A and σ respectively, while remaining are statistical error. For numerical calculations, we choose ILC running at c.m. energy √ŝ = 500 GeV and integrated luminosity L = 100 fb −1 . We use ε A = ε = 0 for most of our analysis however we do discuss the impact of systematic error on our results. With this choice the sensitivity of all the polarization asymmetries, Eq. (13), and the cross section has been calculated varying one parameter at a time. These sensitivities are shown in Fig. 5 and Fig. 6 for the ZZ and Zγ processes, respectively, for each observable. In e + e − → Zγ process we have taken cut-off on polar angle, 10 • ≤ θ γ ≤ 170 • to keep way from the beam pipe. For these limits the analytical expressions, shown in Appendix B , are used.
We see that in the ZZ process the tightest constraint on f γ 4 at 1σ level comes from the asymmetry A y owing to its linear and strong dependence on the coupling. For f γ 5 , both A x and the cross-section σ ZZ , give comparable limits at 1σ but σ ZZ gives tighter limit at higher values of sensitivity. This is because the quadratic term in the σ ZZ comes with higher power of energy/momenta and hence a larger sensitivity. Similarly, the strongest limit on f Z 4 and f Z 5 as well comes from σ ZZ . Though the cross-section gives the tightest constrain on most of the coupling in ZZ process, our polarization asymmetries also provide comparable limits. Another noticeable fact is that σ ZZ has linear as well as quadratic dependence on f Z 5 and the sensitivity curve is symmetric about a point larger than zero. Thus, when we do a parameter estimation exercise, we will always have a bias towards a +ve value of f Z 5 . The same is the case with the coupling f γ 5 as well, but the strength of the linear term is small and the sensitivity plot with σ ZZ looks almost symmetric about f γ 5 = 0.
In the Zγ process, the tightest constraint on h γ 1 comes from A xy , on h γ 3 comes from σ Zγ , on h Z 1 comes from A y and on h Z 3 comes from A x . Cross section σ Zγ and A zz has linear as well as quadratic dependence on h γ 3 and σ Zγ gives two intervals at 1σ level. Other observables can help resolve the degeneracy when we use more than one observables at a time. Still, the cross-section prefers a −ve value of h γ 3 and it will be seen again in a multi-variate analysis. The coupling h Z 3 also has quadratic appearance in the cross-section and yields a bias towards −ve values of h Z 3 . The tightest limits on the anomalous couplings (at 1σ level), obtained using one observable at a time and varying one coupling at a time, are listed in Table 4 along with the corresponding observables. A comparison between Table 1 and Table 4 shows that an e + e − collider running at 500 GeV and 100 fb −1 provides better limits on the anomalous coupling ( f V i ) in the ZZ process than the 7 TeV LHC at For the Zγ process the experimental limits are available from 8 TeV LHC with 19.6 fb −1 luminosity (Table 1) and they are comparable to single observable limits shown in Table 4. These limits can be further improved if we use all the observables in a χ 2 kind of analysis.
We can further see that the sensitivity curves for CPodd observables, A y and A xy , has zero or a very mild dependence on the CP-even couplings. The mild dependence comes through the cross-section σ , that is sitting in the denominator of the asymmetries. We see that CP-even observables provide tight constraint on CP-even couplings and CPodd observables provide tight constraint on the CP-odd couplings. Thus, not only we can study the two processes independently, it is possible to study the CP-even and CP-odd couplings almost independent of each other. To this end, we shall perform a two-parameter sensitivity analysis in the next subsection.
A note of systematic error is in order. The sensitivity of an observable is inversely proportional to the size of its estimated error, Eq. (20). Including systematic error will increase the size of estimated error and hence decrease the sensitivity. For example, including ε A = 1% to L = 100 fb −1 increases δ A by a factor of 1.3 and dilute the sensitivity by the same factor. This modifies the best limit on | f γ 4 |, coming from A y , to 2.97 × 10 −3 (dilution by a factor of 1.3), see Table 4. For cross-section, adding ε = 2% systematic error increases δ σ by a factor of 1.5. The best limit on | f Z 4 |, coming from cross-section, changes to 5.35 × 10 −3 , a dilution by a factor of 1.2. Since inclusion of the above systematic errors modifies the limits on the couplings only by 20% to 30%, we shall restrict ourself to the statistical error for simplicity for rest of the analysis.

Two parameter sensitivity analysis
We vary two couplings at a time, for each observable, and plot the S = 1 (or ∆ χ 2 = 1) contours in the corresponding parameter plane. These contours are shown in Fig. 7 and Fig. 8 for ZZ and Zγ processes, respectively. Asterisk ( ) marks in the middle of these plots denote the SM value, i.e., (0, 0) point. Each panel corresponds to two couplings that are varied and all others are kept at zero. The shapes of contours, for a given observable, are a reflection of its dependence on the couplings as shown in Table 2 and Table 3. For example, let us look at the middle-top panel of Fig. 7 The contours corresponding to the cross-section (solid/red) and A zz (short-dash-dotted/orange) are circular in shape due to their quadratic dependence on these two couplings with same sign. The small linear dependence on f γ 5 makes these circles move towards a small +ve value, as already observed in the one-parameter analysis above. The A y contour (short-dash/blue) depends only on f γ 4 in the numerator and a mild dependence on f γ 5 enters through the cross-section that sits in the denominator of the asymmetries. The role of two couplings are exchanged for the A x contour (big-dash/black). The A xy contour (dotted/magenta) is hyperbolic in shape indicating a dependence on the product f  Table 2 and the expressions in appendix B. Finally, the shape of the A x 2 −y 2 contour (big-dash-dotted/cyan) indicates a quadratic dependence on two couplings with opposite sign. Similarly, all other panels can be read as well. Note that taking any one of the coupling to zero in these panels gives us the 1σ limit on the other coupling as found in the one parameter analysis above.
In the contours for the Zγ process, Fig. 8, one new kind of shape appears: the annular ring corresponding to σ Zγ in middle-top, left-bottom and right-bottom panels. This shape corresponds to a large linear dependence of the cross-section on h γ 3 along with the quadratic dependence. By putting the other couplings to zero in above said panels one obtains two disjoint internals for h  depend upon both the couplings linearly and hence the slanted (almost) parallel lines. Rest of the panels can be read in the same way. Till here we have used only one observable at a time for finding the limits. A combination of all the observables would provide a much tighter limit on the couplings than provided by any one of them alone. Also, the shape, position and the orientation of the allowed region would change if the other two parameters were set to some value other than zero. A more comprehensive analysis requires varying all the parameters and using all the observables to find the parameter region of low χ 2 or high likelihood. The likelihood mapping of the parameter space is performed using MCMC method in the next section.

Likelihood mapping of parameter space
In this section we perform a mock analysis of parameter estimation of anomalous coupling using pseudo data generated by MadGraph5. We choose two benchmark points for coupling parameters as follows: For each of these benchmark points we generate events in MadGraph5 for pseudo data corresponding to ILC running at 500 GeV and integrated luminosity of L = 100 fb −1 . The likelihood of a given point f in the parameter space is defined as: where f 0 defines the benchmark point. The product runs over the list of observable we have: cross-section and five nonzero asymmetries. We use MCMC method to map the likelihood of the parameter space for each of the benchmark point and for both the processes. The 1-dimensional marginalized distributions and the 2-dimensional contours on the anomalous couplings are drawn from the Markov chains using GetDist package [52].

MCMC analysis for e + e − → ZZ
Here we look the process e + e − → ZZ followed by the decays Z → l + l − and Z → qq, with l − = e − , µ − in the MadGraph5 simulations. The total cross-section for this whole process would be, The theoretical values of σ (e + e − → ZZ) and all the asymmetries are obtained using expressions given in Appendix B and shown in the second column of Tables 5 and 6 for bechmark points SM and aTGC, respectively. The MadGraph5 simulated values for these observables are given in the third column of the two tables mentioned for two benchmark points. Using these simulated values as pseudo data we perform the likelihood mapping of the parameter space and obtain the posterior distributions for parameters and the observables. Last two columns of Tables 5 and 6 show the 68% and 95% Bayesian confidence interval (BCI) of the observables used. One naively expects 68% BCI to roughly have same size as 1σ error in the pseudo data. However, we note that the 68% BCI for all the asymmetries is much narrower than expected, for both the benchmark points. This can be understood from the fact that the cross-section provides the strongest limit on any parameter, as noticed in earlier section, thus limiting the range of values for the asymmetries. However, this must allow 68% BCI of cross-section to match with the expectation. This indeed happens for the aTGC case (Table 6), but for the SM case even the cross-section is narrowly constrained compared to naive expectation. The reason for this can be found in the dependence of the crosssection on parameters. For most of the parameter space, the cross-section is larger than the SM prediction and only for a small range of parameter space it can be smaller. This was already pointed out while discussing multi-valued sensitivity in Fig. 5. We found the lowest possible value of crosssection to be 37.77 fb obtained for f γ,Z Thus, for most of the parameter space the anomalous couplings cannot emulate the negative statistical fluctuations in the cross-section making the likelihood function to be , effectively, a one-sided Gaussian function. This forces the mean of posterior distribution to a higher value. We also note that the upper bound of the 68% BCI for cross-section (38.92 fb) is comparable to the expected 1σ upper bound (38.78 fb). Thus we have an overall narrowing of range of posterior distribution of cross-section values. This, in turns, leads to a narrow range of parameters allowed and hence narrow ranges for the asymmetries in the case of SM benchmark point. For the aTGC benchmark point, it is possible to emulate the negative fluctuations in cross-section by varying the parameters, thus corresponding posterior distributions compare with the expected 1σ fluctuations. The narrow ranges for posterior distribution for all the asymmetries are due to the tighter constraints on parameters coming from cross-section and correlation between the observables.
We are using a total of six observables, five asymmetries and one cross-section for our analysis of two benchmark points, however we have only four free parameters. This invariably leads to some correlations between the observables apart from the expected correlations between parameters and observables. Fig. 9 shows most prominently correlated observable for each of the parameters. The CP nature of observables is reflected in the parameter it is strongly  Table 6 List of observables shown for the process e + e − → ZZ for the benchmark point aTGC with √ŝ = 500 GeV. Rest of the details are same as in Table 5.

Observables
Theoretical (      correlated with. We have A y and A xy linearly dependent upon both f γ 4 and f Z 4 , however A y is more sensitive to f γ 4 as shown in Fig. 5 as well. Similarly, for other asymmetries and parameters one can see correlation which is consistent with sensitivity plots in Fig. 5. The strong (and −ve) correlation between A zz and σ shown in Fig. 10 indicates that any one of them is sufficient for the analysis, in principle. However, in practice cross-section puts much strong limit than A zz , which explains the much narrower BCI for it as compared to the 1σ expectation.
Finally, we come to the discussion of parameter estimation. The marginalized 1-dimensional posterior distributions for the parameters of ZZ production process are shown in Fig. 11, while the corresponding BCI along with best fit points are listed in Table 7 for both the benchmark points. The vertical lines near zero corresponds to the true value of parameters for SM and the other vertical line corresponds to aTGC. The best fit points are very close to the true values except for f Z 5 in the aTGC benchmark point due to multi-valuedness of cross-section. The 95% BCI of the parameters for two benchmark points overlap and it appears as if they cannot be resolved. To see the resolution better we plot 2dimensional posteriors in Fig. 12, with the benchmark points shown with an asterisk. Again we see that the 95% contours do overlap as these contours are obtained after marginalizing over unshown parameters in each panel. Any higher dimensional representation is not possible on paper, but we have checked 3-dimensional scatter plot of points on the Markov chains and conclude that the shape of the good likelihood region is ellipsoidal for the SM point with the true value at its center. The corresponding 3-dimensional shape for the aTGC point is like a part of an ellipsoidal shell. Thus in full 4-dimension there will not be any overlap (see Section 5.3) and we can distinguish the two chosen benchmark points as it is quite obvious from the corresponding cross-sections. However, left to only the cross-section we would have the entire ellipsoidal shell as possible range of parameters for the aTGC case. The presence of asymmetries in our analysis helps narrow down to a part of the ellipsoid and hence aids parameter estimation for the ZZ production process.

MCMC analysis for e + e − → Zγ
Next we look at the process e + e − → Zγ and Z → l + l − with l − = e − , µ − in the MadGraph5 simulations. The total crosssection for this whole process is given by, The theoretical values of cross-section and asymmetries (using expression in Appendix B) are given in second column of Tables 8 and 9 for SM and aTGC points, respectively. The tables contain the MadGraph5 simulated data for L = 100 fb −1 along with 68% and 95% BCI for the observables obtained from the MCMC analysis. For the SM point, Table 8, we notice that the 68% BCI for all the observables are narrower than the 1σ range of the pseudo data from MadGraph5. This is again related to the correlations between observables and the fact that the cross-section has a lower bound of about 111 fb obtained for h γ 3 ∼ −4.2 × 10 −3 with other parameters close to zero. This lower bound of crosssection leads to narrowing of 68% BCI for σ and hence for other asymmetries too, as observed in the ZZ production process. The 68% BCI for A x 2 −y 2 and A zz are particularly narrow. For A zz , this is related to the strong correlation between (σ − A zz ), while for A x 2 −y 2 the slower dependence on h γ 3 along with strong dependence of σ on h γ 3 is the cause of a narrow 68% BCI.
For the aTGC point, there is enough room for the −ve fluctuation in the cross-section and hence no narrowing of the 68% BCI is observed for it, see Table 9. The 68% BCI for A x and A y are comparable to the corresponding 1σ intervals while the 68% BCI for other three asymmetries are certainly narrower than 1σ intervals. This narrowing, as discussed earlier, is due to the parametric dependence of the observables and their correlations. Each of the parameters has strong correlation with one of the asymmetries as shown in Fig. 13. The narrow contours indicate that if one can improve the errors on asymmetries, it will improve the parameter extraction. Steeper is the slope of the narrow contour, larger will be its improvement. We note that A x and A y have steep dependence on the corresponding parameters, thus even small variations in parameters leads to large variations in the asymmetries. While for A xy and A x 2 −y 2 the parametric dependence is weaker leading to their smaller variation with parameters and hence narrower 68% BCI.
For parameter extraction we look at their 1-dimensional marginalized posterior distribution function shown in Fig. 14 for the two benchmark points. The best fit points along with 68% and 95% BCI is listed in Table 10. The best fit points are very close to the true values of the parameters and so are the mean of the BCI for all parameters except h γ 3 . For it there is a downward movement in the value owing to the multi-valuedness of the cross-section. Also, we note that the 95% BCI for the two benchmark points largely overlap making them seemingly un-distinguishable at the level of 1dimensional BCIs. To highlight the difference between two benchmark points, we look at 2-dimensional BC contours as shown in Fig. 15. The 68% BC contours (dark shades) can be roughly compare with the contours of Fig. 8. The difference being, Fig. 15 has all four parameters varying and all six observables are used simultaneously. The 95% BC contours for two benchmark points overlap despite the fact that cross-section can distinguish them very clearly. In full 4dimensional parameter space the two contours do not overlap and in next section we try to establish that.

Separability of benchmark points
To depict the separability of the two benchmark points pictorially, we vary all four parameters for a chosen process as a linear function of one parameter, t, as: such that f(0) = f SM is the coupling for the SM benchmark point and f(1) = f aTGC is the coupling for the aTGC point.
In Fig. 16 we show the normalized likelihood for point f(t) assuming the SM pseudo data, L (f(t)|SM), in solid/green line and assuming the aTGC pseudo data, L (f(t)|aTGC), in dashed/blue line. The left panel is for the ZZ production process and the right panel is for Zγ process. The horizontal line corresponds to the normalized likelihood being e − 1 2 , while the full vertical lines corresponds to the maximum value, Table 8 List of observables shown for the process e + e − → Zγ for the SM point with √ŝ = 500 GeV and L = 100 fb −1 . Rest of the details are same as in Table 5.

Observables
Theoretical ( Table 9 List of observables shown for the process e + e − → Zγ for the aTGC point with √ŝ = 500 GeV and L = 100 fb −1 . Rest of the details are same as in Table 5.

Observables
Theoretical (     which is normalized to 1. It is clearly visible that the two benchmark points are quite well separated in terms of the likelihood ratios. We have L (f aTGC |SM) ∼ 8.8 × 10 −19 for the ZZ process and it means that the relative likelihood for the SM pseudo data being generated by the aTGC parameter value is 8.8 × 10 −19 , i.e. negligibly small. Comparing the likelihood ratio to e −n 2 /2 we can say that the data is nσ away from the model point. In this case, SM pseudo data is 9.1σ away from the aTGC point for the ZZ process. Similarly we have L (f SM |aTGC) ∼ 1.7 × 10 −17 , i.e. the aTGC pseudo data is 8.8σ away from the SM point for the ZZ process. For the Zγ process we have L (f aTGC |SM) ∼ 1.7 × 10 −24 (10.5σ ) and L (f SM |aTGC) ∼ 1.8×10 −25 (10.7σ ). In all cases the two benchmark points are well separable as clearly seen in the Fig. 16.

Conclusions
There are angular asymmetries in collider that can be constructed to probe all 8 polarization parameters of a massive spin-1 particle. Three of them, A y , A xy and A yz , are CPodd and can be used to measure CP-violation in the production process. On the other hand A z , A xz and A yz are Podd observables, while A x , A x 2 −y 2 and A zz are CP-and Peven. The anomalous trilinear gauge coupling in the neutral sector, Eq. (14), is studied using these asymmetries along with cross-section. The one and two parameter sensitivity of these asymmetries, together with cross-section, are explored and the one parameter limit using one observable is listed in Table 4 for an unpolarized e + e − collider. For finding the best and simultaneous limit on anomalous couplings, we have performed likelihood mapping using MCMC method and the obtained limits are listed in Tables 7 and 10 for ZZ and Zγ processes, respectively. For the ZZ process, ILC ( √ŝ = 500 GeV, L = 100 fb −1 ) limits are tighter than the available LHC ( √ŝ = 7 TeV, L = 5 fb −1 ) limits [40], while ILC limits on the Zγ anomalous couplings are slightly weaker than the available LHC ( √ŝ = 8 TeV, L = 19.6 fb −1 ) limits [44]. The LHC probes the interactions at large energies and transverse momentum, where the sensitivity to the anomalous couplings are enhanced. While we perform our analysis at √ŝ = 500 GeV leading to a weaker, though comparable, limits on h V i s in the Zγ process. With polarized initial beams, the P-odd observables, A z , A xz and A yz , will be non-vanishing for both the process with appropriate kinematical cuts. This gives us three more observables to add in the likelihood analysis, which can lead to better limits. At LHC, one does not have the possibility of initial beam polarization, however, W + W − and ZW ± processes effectively have initial beam polarization due to chiral couplings of W ± . Study of W and Z processes at LHC and polarized e + e − colliders are under way and will be presented elsewhere.