Cosmologically varying kinetic mixing

The portal connecting the invisible and visible sectors is one of the most natural explanations of the dark world. However, the early-time dark matter production via the portal faces extremely stringent late-time constraints. To solve such tension, we construct the scalar-controlled kinetic mixing varying with the ultralight CP-even scalar’s cosmological evolution. To realize this and eliminate the constant mixing, we couple the ultralight scalar within 10−33eV ≲ m0 ≪ eV with the heavy doubly charged messengers and impose the ℤ2 symmetry under the dark charge conjugation. Via the varying mixing, the keV – MeV dark photon dark matter is produced through the early-time freeze-in when the scalar is misaligned from the origin and free from the late-time exclusions when the scalar does the damped oscillation and dynamically sets the kinetic mixing. We also find that the scalar-photon coupling emerges from the underlying physics, which changes the cosmological history and provides the experimental targets based on the fine-structure constant variation and the equivalence principle violation. To ensure the scalar naturalness, we discretely re-establish the broken shift symmetry by embedding the minimal model into the ℤN-protected model. When N ~ 10, the scalar’s mass quantum correction can be suppressed much below 10−33eV.

We know that dark matter exists, but we do not know the dark matter's particle nature.Even so, we can naturally imagine that the dark matter stays in the invisible sector, and the invisible and visible sectors are connected by the portal.Through the portal, the energy flow from the visible sector to the invisible sector, which is known as freeze-in [1], or vice versa, which is known as freeze-out.Hence, the dark matter reaches Ω DM h 2 0.12, being compatible with the CMB anisotropy [2,3].The kinetic mixing [2], one of the three major portals [2,[4][5][6][7][8][9][10][11][12][13], connects the photon and the dark photon as L ⊃ F µν F µν /2.In the last few decades, the research on the time-independent kinetic mixing has boosted on both the experimental and theoretical sides [14][15][16][17][18][19][20][21][22][23][24][25][26], with only a few discussions on the spacetime-varying scenarios [27][28][29][30].In the meantime, other varying constants are extensively discussed .Moreover, extremely strong tensions exist between the early-time dark matter production through the portal and the late-time constraints on the portal, such as the freeze-in of the dark matter as the dark photon [58][59][60] or the sterile neutrino [61,62], and the freeze-out of the dark matter as the dark Higgs [63].To solve such tension in the simplest way, we allow the portal evolves during the universe's expansion.However, there is no free lunch to evade the constraints without consequences.To be more specific, controlling the portal leaves significant imprints in our universe, which changes the early cosmological history and can be detected by experiments designed for general relativity testing and ultralight dark matter detection.
In this work, we study the scalar-controlled kinetic mixing by meticulously exploring the top-down and bottom-up theories, cosmological history, keV − MeV dark photon dark matter production, and experimental signals from both the dark photon dark matter and the nonrelativistic ultralight scalar relic.To vary the kinetic mixing, we couple the ultralight scalar φ, the CP-even degree of freedom predicted by the string theory [64][65][66][67][68][69][70], to the heavy fermionic messengers doubly charged under the standard model U (1) and dark U (1).Here, the constant kinetic mixing is eliminated when the Z 2 symmetry under the dark charge conjugation is imposed.Given this, in the low energy limit, the varying-mixing operator φF F emerges, along with the scalar-photon coupling, such as φF 2 or φ 2 F 21 .Initially, φ has the early misalignment opening the portal for the dark photon dark matter production with the kinetic mixing FI ∼ 10 −12 , which stems from the early-time Z 2 -breaking of the system.Afterward, φ's damped oscillation gradually and partially closes the portal, which stems from the late-time Z 2 -restoration.Through the evolution during the cosmological expansion, φ sets the benchmark kinetic mixing of the dark photon dark matter, which is free from stringent late-time constraints, such as stellar energy loss [71][72][73][74], direct detection [60,[75][76][77][78][79][80], and late-time decay bounds [58,59,81,82].At the same time, via the scalar-photon coupling, the ultralight scalar as the nonrelativistic Based on this model, the energy flows from the dark sector to the standard model sector in the early universe for the dark matter production with the portal opened.In the late universe, with the portal partially closed, the dark matter is safe from stringent constraints.Right: The cosmologically varying kinetic mixing in the UV and IR theories.In the UV theory, there are heavy messengers Ψ and Ψ , both charged under U (1)Y and U (1) d .Here, Ψ and Ψ carry the same electromagnetic charge but the opposite dark charge.To generate the varying kinetic mixing, Ψ and Ψ are coupled with the scalar Φ via the Yukawa interactions y and y .To eliminate the constant kinetic mixing, the mass degeneracy between Ψ and Ψ is imposed, which is protected by the Z2 symmetry.In the IR theory, the nondynamical heavy messengers are integrated out, which induces the varying kinetic mixing (φF F ) and the scalar-photon coupling (φF 2 or φ 2 F 2 ).In the nonperturbative region, φ should be replaced by f sin(φ/f ) according to Eq. ( 7) and Eq. ( 13) because it is the angular part of Φ.Based on whether the Z2 symmetry is slightly broken by the Yukawas or not, we classify the theory into two types with totally different phenomenologies: The type-A model with the linear scalar-photon coupling (y = 0, y = 0) and the type-B model with the quadratic scalar-photon coupling (y = −y).These scalar-photon interactions lead to nontrivial cosmological history caused by the thermal effect and the signals from the αem variation and the equivalence principle violation.
relic in the mass range 10 −33 eV m 0 eV changes the fine-structure constant, and the scalar as the mediator contributes to the extra forces between two objects.Based on these facts, the experiments such as the equivalence principle (EP) violation test [83][84][85], clock comparison [86][87][88][89][90], resonant-mass detector [91], PTA [92], CMB [93,94], and BBN [93,95,96] can be used to test the scalar-controlled kinetic mixing, and the experimental targets are set by the dark photon freeze-in.In the meantime, the experimental targets for keV − MeV dark photon dark matter detections are set by the ultralight scalar mass.If the signals from the dark photon dark matter and the ultralight scalar experiments appear consistently, we can confidently declare the verification of our model.In addition, given the scalar-photon coupling in the strong region, the scalar's high-temperature evolution is affected by the standard model plasma, which sets the early displacement, modifies the start of oscillation, and even enhances the scalar's signal.To understand the whole setup classified under the exactness of the Z 2 symmetry, one can refer to Fig. 1.
To protect the CP-even scalar's naturalness caused by the heavy messengers and inspired by the former works on the discrete symmetries [97][98][99][100][101][102][103][104], we embed the varying kinetic mixing into N copies of the universes, where the Z N symmetry rebuilds the global U (1) shift symmetry in the discrete form.In such Z N -protected model, the scalar's lowest order mass term becomes Φ N /Λ N −4 , which reveals the exponential suppression of the quantum correction.For example, we only need N ∼ 10 to suppress the scalar mass correction all the way down to 10 −33 eV.Furthermore, to understand the additional cancellation from the exact Z 2 symmetry and gain an accurate analytical result, we expand the Z N Coleman-Weinberg formally calculated to the leading order in [98,102] to the arbitrary orders with two entirely different methods, i.e., the Fourier transformation and the cosine sum rules.Other topics of the varying kinetic mixing from the supersymmetric Dirac gaugino model and the dark matter models via other cosmologically varying portals are preliminarily discussed in our work.
The remainder of this paper is organized as follows.In Sec.II, we build the minimal model for the scalar-controlled kinetic mixing and show that the scalar-photon coupling appears simultaneously.Based on whether the scalarmessenger couplings are Z 2 -invariant, we categorize the theory into the type-A model with the linear scalar-photon coupling and the type-B model with the quadratic scalar-photon coupling.In Sec.III, for the type-A and the type-B models, we do the systematic classification of the scalar evolution jointly affected by the thermal effect, bare potential effect, and cosmological expansion.In Sec.IV, we discuss the dark photon dark matter freeze-in via the varying kinetic mixing.We also discuss the detection of the dark photon dark matter with the experimental targets set by the scalar mass and the experiments of the non-relativistic ultralight scalar relic with targets set by the dark photon dark matter freeze-in.In Sec.V, we build the Z N model to protect the scalar's naturalness from the heavy messengers, discuss the extra cancellation from the exact Z 2 , and calculate the Z N -invariant Coleman-Weinberg potential utilizing the Fourier transformation and the cosine sum rules.In Sec.VI, we generate the varying kinetic mixing and the dark axion portal simultaneously from the Dirac gaugino model.In Sec.VII, we preliminarily explore the dark matter models via other cosmologically varying portals and their experimental signals.Finally, in Sec.VIII we summarize our results.We also provide a series of appendices containing the calculational details, such as the analytical solutions of the high-temperature scalar evolution in Appendix.A, the freeze-in of the dark photon dark matter in Appendix.B, and the exact expansion of the Z N Coleman-Weinberg potential using the Fourier transformation and cosine sum rules in Appendix.C.

II. MINIMAL MODEL
In the beginning, let us recall the well-known constant kinetic mixing of the field strengths of the standard model U (1) Y and the dark U (1), which can be written as Here, the kinetic mixing can be generated at the one-loop level by a pair of vector-like fermions Ψ and Ψ charged as (e, e ) and (e, −e ) under U (1) Y × U (1) d [2].In the low energy limit, the value of the constant kinetic mixing is where M and M are the masses of Ψ and Ψ , respectively.To build the varying kinetic mixing, we eliminate the constant kinetic by imposing the mass degeneracy M = M and promote to a dynamical variable (x) by imposing the Z 2 symmetry under the dark charge conjugation, i.e., By doing so, the constant kinetic mixing is forbidden because F F is odd under C d , whereas the varying kinetic mixing (x)F F is permitted.
To realize this from the top-down theory, we introduce the Lagrangian where Φ is a neutral complex scalar and y ( ) , M , λ, f are chosen to be real parameters for simplicity.In the ground state, the approximate global U (1) symmetry of the potential is spontaneously broken by Φ = f √ 2 e i(φ/f +c) whose angular components include a CP-even pseudo-Goldstone boson φ and an arbitrary phase c.Knowing that the transformation on the real scalar φ → −φ−(2c − π) f is equivalent to Φ → −Φ † , under which the Yukawa interactions in Eq. ( 4) flip signs, we choose c = π/2 to fix the phase factor throughout the paper.Given this, we define the dark charge conjugation of the UV theory as If y = −y and M = M , the Lagrangian is invariant under C d .We will see from the later discussion that as long as φ has a nonzero vacuum expectation value (VEV), the dynamical kinetic mixing can be generated.This is because, in this case, Ψ ( ) 's effective masses become From Eq. ( 6), we know that the degeneracy of the effective masses is broken by φ when y = y .In Eq. ( 6) and the rest part of this paper, we use "( )" to denote the physical quantities in the SM sector (without " ") and the dark sector (with " ").After integrating out Ψ ( ) , the Lagrangian of the effective kinetic mixing becomes When φ f , the kinetic mixing in Eq. ( 7) can be linearized as , where Therefore, in our model, the kinetic mixing varies with φ's cosmological evolution.Since the crucial point of this model is that the kinetic mixing vanishes when φ relaxes to the minimum, we should align the potential's minimum with the kinetic mixing's zero.To realize this, we introduce the C d -even tadpole term which arises as the potential when the phase factor is c = π/2.Conversely, because the system is invariant under C d listed in Eq. ( 5), V 0 has to be the even function of φ, while is the odd function of φ.In this case, the phase difference between and V 0 is (n + 1/2)π, where n = 0, ±1, ±2 • • • .Therefore, under the protection of Z 2 symmetry of the dark charge conjugation C d , the potential's minima are necessarily aligned with the zeros of the kinetic mixing.During the cosmological evolution, φ initially stays at the nonzero displacement, so the portal is opened.As the universe expands, the pseudo-Goldstone boson φ begins the damped oscillation around zero, so the portal is partially closed.In today's universe, φ consists of the nonrelativistic ultralight relic in the mass range 10 −33 eV m 0 eV, so the kinetic mixing has the nonzero residual.More fundamentally, we can interpret the time-varying kinetic mixing from the perspective of symmetry, since the nonzero kinetic mixing results from the Z 2 violation.In the early universe, the Z 2 symmetry is spontaneously broken by φ's nonzero displacement.In the late universe, φ oscillates toward zero, and the Z 2 symmetry is restored.Such early-time breaking and late-time restoration of the global symmetry is regarded as the inverse symmetry breaking in [105][106][107][108] for other motivations.
Let us classify the theory based on whether the Lagrangian is exactly invariant under the dark charge conjugation C d or not.As shown in Eq. ( 5), when y = −y, the Z 2 symmetry of C d is strictly preserved.When y = −y, the Z 2 symmetry is broken, which induces the two-loop constant kinetic mixing and the one-loop tadpole of Φ.Here, the small Yukawas highly suppress the two-loop constant mixing, which is ∼ (M/Λ KM ) 2 /ee .The cancellation of Φ-tadpole needs fine-tuning, which can be entirely solved in the framework of Z N -protected model discussed in Sec.V. Therefore, the y = −y case still has the approximate Z 2 symmetry.Based on the discussion above, we consider the following two types of the models in this work: Due to φ-Ψ ( ) interaction, the scalar-photon coupling emerges from UV physics.At the one-loop level, the coupling between the CP-even scalar φ and the photon can be written as where Utilizing the classification in Eq. ( 11), we have Type-A: where the type-A and type-B models have linear and quadratic scalar-photon couplings, respectively.To compare with the experiments testing the fine-structure constant variation and the equivalence principle violation, we define the dimensionless constants d γ,i (i = 1, 2) through where the indices "i = 1" and "i = 2" denote the type-A and type-B models, respectively.Comparing Eq. ( 15) with Eq. ( 14), we have Some other papers quantify the fine-structure constant variation through the notation ∆α em /α em = φ i /Λ i γ,i (i = 1, 2).Combining this notation with Eq. ( 14), we have Λ γ,i ∼ e Λ KM (i = 1, 2), which means that e determines the hierarchy between the effective energy scale of the varying kinetic mixing and the varying fine structure constant.In the rest part of this paper, we use the first notation, i.e., d γ,i (i = 1, 2).In Eq. ( 16), we find that for the fixed φF F operator inducing the effective kinetic mixing, smaller e leads to larger y ( ) such that the testable signals of the α em variation get stronger.Such small dark gauge coupling can be naturally generated in the large volume scenario in string compactification [18,109,110].
To understand the setup of our model intuitively, one can refer to Fig. 1.The left panel shows that when Z 2 symmetry is imposed, the constant kinetic mixing is canceled, but the scalar-controlled kinetic mixing survives.In this case, as φ evolves during the cosmological evolution, the kinetic mixing becomes time-dependent.This provides a novel mechanism to generate a small but non-zero kinetic mixing.The right panel reveals that the scalar-photon couping is also generated as the byproduct when UV physics is considered.Based on the exactness of the Z 2 symmetry, the theory can be classified as the type-A model with the linear scalar-photon couping and the type-B model with the quadratic scalar-photon couping.We will see from the later discussion that such scalar-photon couplings affect φ's evolution via the thermal effect from the SM plasma.They also change the fine-structure constant and violate the equivalence principle, which provides essential prospects for the experimental tests.

III. COSMOLOGICAL HISTORY
In this section, we jointly discuss the ultralight scalar's cosmological evolution, which is affected by the scalar bare potential, thermal effect, and cosmological expansion.According to [102,111], the lowest order thermal contribution containing α em is at the two-loop level, and the free energy contribution coming from that can be written as In Eq. ( 17), the suppression factor S(T ) is 1 when e ± (as well as other species of heavier particles carrying the U (1) em charges) are relativistic but decreases exponentially after e ± becomes non-relativistic.The factor i q 2 i (q i is the electric charge of the particle "i") quantifies the thermal contribution to the total free energy when U (1) em charged particles are relativistic in the thermal bath.As shown in Eq. ( 13), when φ has nonzero VEV, the fine structure constant is modified, based on which φ's thermal potential can be obtained by replacing α em in Eq. ( 17) by α em (φ).After combining Eq. ( 13) and Eq. ( 17), we have where the definition of r can be found in Eq. ( 6).
Here, we choose T rh to be much smaller than the heavy fermion mass, i.e., M , so Ψ's contribution to the thermal potential is exponentially suppressed.In addition, because we do not include the dark electrons and positrons, the thermal effect from the dark sector is also negligible.Fig. 3 shows how the potential changes as the temperature varies and how φ moves under different circumstances.From lighter (yellow) to darker color (dark red), the ratios of thermal mass over the bare mass are m T /m 0 = 5, 4, 3, 2, 1, 0. For the bare potential, V 0 has the minimum at φ/f = 0.As for the thermal potential, the local minimums of the type-A and type-B models are (Here n = 0, ±1, ±2, • • • ) Type-A: In the following discussion, we focus on the range −2πf ≤ φ ≤ 2πf without loss of generality.For the type-A model, since m T m 0 in the early epoch, the minimum is φ min πf /2.When m T m 0 , the potential's minimum continuously shifts to zero.For the type-B model, when m T > m 0 , within the 2πf periodicity there are three local minimums, i.e., φ min = 0, ±πf .When m T ≤ m 0 , only φ min = 0 is the minimum.We should notice that during the entire process, φ keeps the classical motion without the tunneling: Even though φ is initially located inside the false facuum πf /2 φ rh 3πf /2, because f is large, the vacuum decay from πf to 0 is highly suppressed such that it is much slower than the universe's expansion [112][113][114][115][116][117][118].
Knowing from Eq. ( 18) that m T has the same temperature power law as H when T m e (S 1), we define a dimensionless quantity to classify φ's evolution.Here, r ∼ yf /M , which is defined in Eq. ( 6).The motion of φ is underdamped if η > 1, because the thermal effect dominates over the universe's expansion; In contrast, if η < 1, φ's motion is overdamped under the Hubble friction.η = 1 is the critical value separating the two moving patterns mentioned above.To relate η with experimental observables which are quantified by d γ,i (i = 1, 2), we write η as In the high-T limit where m T dominates over m 0 , φ's thermal evolution can be solved in the linear approximation of the equation of motion, which can be found in Appendix.A. Since we compare m T and m 0 to determine whether the thermal effect dominates over the effect from the bare potential, we define the critical temperature T * at which m T = m 0 .Combining Eq. ( 18) and Eq. ( 20), we have where the concrete formulas of η for the type-A and the type-B models can be found in Eq. (21).Now, let us briefly explain the meaning of Eq. (22).For m T to be smaller than m 0 , we only need one of the following two conditions Here, η ∼ mT /H when T me, which is already defined in Eq. ( 22).Left: η 1.Because the thermal effect from the SM bath is negligible, the ultralight scalar's evolution can be classified as the standard misalignment.In this case, φ stays constant at the early stage and begins the oscillation when H ∼ m0, which is labeled as Tosc.Right: η 1.The thermal effect is important for the scalar's evolution.In the plot, point T * and point "Q" denote the moment when mT = m0 and H ∼ mT , respectively.Here, HQ ∼ 10 −16 / log 2 η eV.When T T * , φ converges to the minimum of the thermal potential obeying the power law |φ − φmin| ∝ T 1/2 where φmin = πf /2 for the type-A model and φmin = ±πf for the type-B model.We merely show the m0 HQ case, because if m0 HQ, φ's movement is nothing more than the late-time standard misalignment with the initial condition fixed by the early time thermal effect.For the type-A model, when T me, φ tracks the potential minimum φmin = f arctan m 2 T /m 2 0 .As the temperature drops far below the electron mass, φ cannot track the minimum and begins to oscillate.For the type-B model, we focus on the πf /2 φ rh 3πf /2 case.When T T * , φ is trapped inside the local minimum φmin = πf .Thereafter, when T T * , φ begins oscillating around the bare minimum φmin = 0.Here we have T * T |3H=m 0 , which means the oscillation is postponed.Therefore, the scalar is in the phase of trapped misalignment.
to be satisfied: The first condition is T m e such that m T is exponentially suppressed.The second condition is that the unsuppressed part of the scalar thermal mass, i.e., m T /S 1/2 , is smaller than m 0 .Since T * is defined as the temperature when m T = m 0 , we use T ≶ T * and m T ≷ m 0 interchangeably.In Fig. 4 showing the typical numerical solutions of the φ evolution, we define the dimensionless temperature T := T /T | 3H=m0 with the extra "∼", then we have the dimensionless critical temperature T * ∼ max[m e /(m 0 m pl ) 1/2 , 1/η 1/2 ] for the moment when m T = m 0 .
Before discussing the scalar's evolution, let us give the following definitions: These definitions are based on the value of η, which quantifies how large the thermal effect is compared with the effect of the bare potential shown in Eq. (10).When η is much smaller than one, the thermal effect is negligible, so it can be classified as the standard misalignment in which one solely needs to compare the effects of the bare potential and the universe expansion.When η is of order unity or even larger, the thermal effect can never be omitted, and one should consider all three effects (thermal effect, bare potential, universe expansion) jointly.Now, let us discuss each situation one by one.
A. Standard Misalignment: η 1 In this case, when 3H m 0 , φ obeys Figure 3.The typical scalar evolution trajectories in the early universe.The lighter (darker) color in these diagrams denotes higher (lower) temperature.From the yellow line to the dark red line, mT /m0 = 5, 4, 3, 2, 1, 0. Here, η quantifies the thermal effect, whose definition can be found in Eq. (20).Upper left and upper right: For both type-A and type-B models, when η 1, the thermal effects are negligible.Here, φ is nearly frozen at the early stage (gray arrowed line from point I to point II).When 3H m0 (point II), φ begins the damped oscillation with the amplitude |φ| ∝ T 3/2 (black arrowed line from point II to point III) and converges to the origin (point III).Lower left: For the type-A model with η 1, in the high-temperature environment, φ converges to πf /2 (point I) following |φ − πf /2| ∝ T 1/2 .As the temperature falls, mT /m0 also decreases.During the stage when T me, the variation of the thermal potential is adiabatic, so φ tracks the potential minimum φmin = f arctan m 2 T /m 2 0 as shown in Eq. ( 19) (gray arrowed line from point I to point II).When T me, since the too fast variation of the thermal potential violates the adiabatic condition, φ can no longer stay in the minimum and starts to oscillate following |φ| ∝ T 3/2 .We do not show this stage of φ's movement in the diagram because the amplitude is too small to visualize compared with f .Lower right: For the type-B model with η 1 and φ initially trapped inside the wrong vacuum, φ does damped oscillation |φ − πf | ∝ T 1/2 when T T * and converges to πf .When T T * , V | φ=πf 0, so φ begins the damped oscillation with |φ| ∝ T 3/2 (black arrowed line from point I to point II).
which means φ is nearly frozen in the early universe.To have an overall understanding of the cosmological history, one could refer to the left panel of Fig. 2: The H-line (orange) is higher than the m T -line (red) during the whole process, meaning that the thermal effect is negligible so one only needs to focus on the H-line (orange) and m 0 -line (blue).After the crossing of H-line and m 0 -line, the scalar begins the damped oscillation whose amplitude obeys the power law |φ| ∝ T 3/2 .The upper left and upper right panels of Fig. 3 show how φ moves for the type-A and type-B models separately when η 1: From point I to point II, φ keeps nearly the same.This means that the initial field displacement determines φ's starting oscillation amplitude |φ| osc .Here, we focus on the model with the natural initial condition, i.e., |φ| osc /f ∼ O(1).In the late time, φ begins to oscillate at the temperature (m 0 10 −28 eV) 2 3 (10 −33 eV m 0 10 −28 eV) β T in Eq. ( 25) and later mentioned β φ in Eq. ( 26) are determined by H's power law of T , which is T 2 at the radiation-dominated universe, and T 3/2 at the matter-dominated universe.In Eq. ( 25), we choose m 0 10 −28 eV as the reference value simply because the scalar oscillation happens at the matter-radiation equality T ∼ eV for such scalar mass.Before or after T ∼ eV, H has different power laws of T .Such oscillation is labeled as black oscillatory lines from point II to point III in Fig. 3 and is also represented as blue lines in Fig. 4. Because the inflation smears out the field anisotropy, φ does spatially homogeneous oscillation, which can be written as From [119,120] one knows that φ satisfies ρ φ ∝ T 3 and w φ = p φ /ρ φ 0, so φ is part of the dark matter with the fraction F = Ω φ /Ω dm .Without loss of generality, we choose F = 10 −3 as a benchmark value, such that astrophysical and cosmological measurements like Lyman-α forest [121][122][123][124][125][126], CMB/LSS [127,128], galaxy rotational curves [129][130][131][132][133], and the observations of the ultra-faint dwarf galaxies [134] do not exclude the parameter space where m 0 10 −19 eV.After substituting Eq. ( 25) into Eq.( 26), one can easily get the starting oscillation amplitude In the case where m 0 10 −25 eV, φ's de Broglie wavelength is much smaller than the scale of the Milky Way halo, so φ behaves like a point-like particle similar to other cold DM particles.For this reason, φ's local density is ρ φ,local Fρ local , where ρ local 0.4GeV/cm 3 is DM's local density near earth.By effectively adding an enhancement factor where ρ average = ρ c Ω DM 1.3 keV/cm 3 on |φ| in Eq. ( 26), we get φ's amplitude today near the earth, which is |φ| 0 (2Fρ local ) 1/2 /m 0 .Here, |φ| 0 denotes φ's local oscillation amplitude today.If m 0 10 −25 eV, oppositely, φ cannot be trapped inside the Milky Way halo's gravitational potential well.In this case, today's φ field is homogeneous, so the enhancement factor in Eq. ( 28) should not be included, from which we have |φ| 0 = (2Fρ average ) 1/2 /m 0 .φ's oscillation amplitude in the middle mass range needs numerical simulation, which is left for future exploration.

B. Thermal Misalignment: η 1
As long as η is of order unity or even larger, the thermal effect from the SM plasma plays a decisive role in φ's evolution in the early universe.To understand the combined effects on φ's movement from the thermal potential, bare potential, and the universe's expansion, one can refer to the right panel of Fig. 2: At the early stage, since m T (red) and H (orange) are both proportional to T 2 , their lines are approximately parallel in the plot with the log-log scale.
Here we neglect the g * variation, which only leads to an O(1) modification.In the plot, we label the cross point of m T -line (red) and m 0 -line (blue) with "T * ", which is already defined in Eq. ( 22).When T T * , the effect from the bare potential dominates over the effect from the thermal potential.When T m e , the m T -line drops fastly T /m 2 0 until the temperature drops below the electron mass.After that, φ cannot track the instant changing of the potential and begins the damped oscillation scaled as |φ| ∝ T 3/2 .In this plot, we have not shown φ's oscillation afterward because its amplitude is too small to be visualized.For the green line (η ∼ 1), when T 1, φ gradually slides to the thermal potential minimum πf /2 following Eq.( 30).When T ∼ 1, φ begins the damped oscillation obeying |φ| ∝ T 3/2 .The blue line (η 1) represents the standard misalignment where the thermal effect is negligible during the process.There, φ const until T 1.After that, φ begins to oscillate following |φ| ∝ T 3/2 .Right: The type-B model (y = −y).The red line represents the situation where η 1 and φ's initial condition is πf /2 φ rh 3πf /2.In this case, when T T * , φ does the damped oscillation |φ − πf | ∝ T 1/2 during which φ converges to πf .When T T * , or equivalently speaking, mT m0, we have V | φ=πf 0, so φ rolls down from πf and oscillates like |φ| ∝ T 3/2 .We can see that the red line's oscillation is postponed compared with the standard misalignment beginning at T 1.For this reason, the red line is in the phase of the trapped misalignment.The green line represents the situation where η ∼ 1 and πf /2 φ rh 3πf /2.In this case, φ gradually slides to the thermal potential minimum πf as shown in Eq. ( 30).At the point T 1, φ begins the damped oscillation scaled as |φ| ∝ T 3/2 .The blue line represents the standard misalignment mechanism where the thermal effect is inferior to the bare potential effect, i.e., η 1.Similar to the discussion of the blue line in the left panel, φ stays constant in the early universe and then begins the oscillation when T 1.Even though the red lines in the left and right panels are qualitatively different, the green and blue lines are quite similar in terms of the oscillation temperature ( T 1).The main difference between the green and blue lines is that the green lines' initial displacements are thermally determined to be |φ| osc |φmin|, but the blue lines' oscillation amplitude depends on their initial conditions.following e −me/T and crosses the H-line.We label the cross point of the m T -line and the H-line as "Q", and we have As we will see in the rest of this section, comparing m 0 and H Q is vital in determining the temperature at which φ starts the late-time oscillation.
Let us first discuss φ's movement in the stage T T * , during which the bare potential can be omitted in the high-temperature environment.In Appendix.A, we solve the scalar evolution for this situation.Because φ − φ min is the linear combination of , to describe the thermal convergence more quantitatively, we need the specific value of η.We recall that the potential minimums are φ min = πf /2 for the type-A model and φ min = 0, ±πf for the type-B model.Given the initial condition φrh 0, we approximately have where η = 1 is the critical value determining how φ evolves towards the local minimum.If η ≤ 1 (but not too small), φ gradually slides to φ min , as shown in the green lines in the left and right panels of Fig. 4. Incidentally, in the limit η 1, we come back to Eq. ( 24) where φ is nearly frozen at the early universe, which is related to the blue lines in both panels of Fig. 4. Here, the term is negligible because it describes the movement with non-zero φrh .If η > 1, φ's movement towards φ min is oscillating because Obeying the power law T 1/2 , φ's oscillation amplitude decreases as the temperature goes down.This can be explained more intuitively: When T m e , the adiabatic condition of the WKB approximation is satisfied ( ṁT /m 2 T ∼ 1/η 1), so there is no particle creation or depletion for φ, or equivalently speaking, its number density is conserved.For this reason, φ's oscillation amplitude obeys |φ − . Eq. ( 30) reveals an interesting phenomenon: As long as there is a hierarchy between T rh and m e , which is natural in most of the inflation models, φ converges to the local minimum of the thermal potential.For example, in η 1 case, given that T rh ∼ 1TeV, when the universe temperature drops to T ∼ m e , the field deviation from the local minimum becomes |φ − φ min |/φ min ∼ 10 −3 .Higher T rh leads to even smaller field displacement from the thermal minimum.Such determination of scalar's misalignment through the thermal effect is named the thermal misalignment in several recent works [102,135,136].Other mechanisms setting scalar's nonzero initial displacement can be found in [137][138][139].As shown in Eq. ( 23), in our paper, we use the term thermal misalignment for the cases satisfying η 1 to distinguish them from the standard misalignment in which the thermal effect does not play a role.We use such a definition because when this condition is satisfied, the thermal effect from the SM bath washes out φ's sensitivity of the initial condition and dynamically sets the field displacement.
When T T * , the bare potential becomes important in φ's evolution.In the m 0 H Q case, the ultralight scalar's evolution is simply the combination of the early time wrong vacuum damped oscillation and the late-time true vacuum damped oscillation around zero with where φ min = πf /2 for the type-A model and φ min = ±πf for the type-B model with φ initially in the wrong vacuum.
We do not put more words on that because this kind of the thermal misalignment can be treated as the standard misalignment with the thermally determined initial condition.Now, we shift our focus to the m 0 H Q case in the rest of this section.Since the total potentials of the type-A and type-B models have different shapes depending on exactness of the Z 2 symmetry, and the scalar evolution depends on the numerical value of η, let us describe the scalar evolution case by case: When T T * , φ does the damped oscillation, which converges to πf /2.Afterward, when T ∼ T * , or equivalently speaking, m T ∼ m 0 , the potential minimum begins to shift towards from πf /2 to 0, obeying φ = f arctan m 2 T /m 2 0 as shown in Eq. (19).When the universe temperature is much higher than m e , the adiabatic condition is satisfied because ṁT /m 2 T ∼ 1/η 1.We show the movement of φ in the lower left panel of Fig. 3 and the red line in the left panel of Fig. 4.However, as the temperature drops below m e , because ṁT /m 2 T ∼ e me/T /η 1, φ is not able to respond to the sudden variation of the potential minimum anymore and begins the oscillation.Here, the oscillation temperature and the starting amplitude can be written as where T Q ∼ m e / log η as defined in Eq. (29).Given f m pl , the hierarchy between T * and T Q strongly suppresses φ's late-time oscillation amplitude, leading to φ's minuscule relic abundance.Even though this phase does not appear in the following context of the dark photon dark matter freeze-in in Sec.IV, we still discuss this phase in this section for completeness.
In the early time when T T | 3H=m0 , φ slides to the thermal minimum πf /2.When T ∼ T | 3H=m0 , the scalar begins the late-time damped oscillation obeying |φ| ∝ T 3/2 with the starting temperature and amplitude In Fig. 3, we do not list such a case because, even though categorized as the thermal misalignment, it can be decomposed into the standard misalignment (upper left panel of Fig. 3) plus the determined initial amplitude.
In the left panel of Fig. 4, the green line describes such situation: When T 1, the green line gradually slides to the πf /2.When T ∼ 1, the green line begins the damped oscillation with the amplitude scaled as |φ| ∝ T 3/2 .
Unlike the type-A model, there is no continuous shift of the potential minimum during the cosmological evolution.Therefore, we only need to focus on the moment when the local minimum flips.Here, we mainly focus on the case in which φ's initial condition satisfies πf /2 φ rh 3πf /2.In this case, because of the thermal effect when T T * , φ converges to the local minimum πf inside the wrong vacuum.When T T * , at the point φ = πf , the second order derivative becomes nonpositive, i.e., V | φ=πf 0, thereafter φ begins the oscillation around zero.From Eq. ( 22), we have According to Eq. ( 22), we recall that For the η 1 case, one could look at the red line in the right panel of Fig. 4. Here, φ oscillates around the thermal minimum πf with the power law |φ − πf | ∝ T 1/2 when T T * .Afterward, when T T * , alternatively speaking, m T m 0 , φ begins the oscillation following |φ| ∝ T 3/2 .The plot shows the apparent postponement of the scalar oscillation for the red line compared with the other two lines.Because the CP-even scalar φ's oscillation is postponed, the evolution of φ can be classified as the trapped misalignment, which is formally investigated in axion models [140,141].In this case, φ's relic abundance is enhanced given the same |φ| osc or f .We can also think about such characteristics inversely: For φ to reach the same abundance quantified by F, one only needs smaller |φ| osc or f .To be more quantitative, one can write the misalignment at the beginning of the oscillation as where |φ| osc,std denotes the starting amplitude for the standard misalignment as shown in Eq. ( 27).From Eq. ( 35), one can see that the necessary early misalignment |φ| osc is rescaled by a factor of (T * /T | 3H=m0 ) 3/2 , which shows that |φ| osc is much smaller compared with |φ| osc,std when F is determined.
For the η ∼ 1 case, one could refer to the green line in the right panel of Fig. 4. In the early stage, i.e., . When T T | 3H=m0 , φ begins the damped oscillation with the power law |φ| ∝ T 3/2 .From the discussion above, we can find that the η ∼ 1 case has the thermal determination of the initial condition but does not have the postponement of the oscillation.Therefore, φ is in the phase of the thermal misalignment but not in the phase of the trapped misalignment.
Finally, let us briefly discuss the case where −πf /2 φ rh πf /2.Here we have that φ oscillates obeying |φ| ∝ T 1/2 when T T * , and then oscillates obeying |φ| ∝ T 3/2 when T T * .Because the late-time amplitude is suppressed by a factor of (T rh /T * ) 1/2 , given the nonnegligible fraction of φ among the dark matter (For example, F ∼ 10 −3 ), f depends on T rh and may be larger than m pl .Therefore, this paper focuses on the case where φ is initially localized inside the wrong vacuum.
Realizing that all the production mechanisms mentioned above do not rely on kinetic mixing, which is indispensable for dark photon detection, we provide a minimal extension of the dark photon dark matter freeze-in where the ultralight scalar's evolution sets the experimental benchmark of kinetic mixing dynamically.In addition, we want to stress the enormous significance of detecting the ultralight scalar in the whole mass range, i.e., 10 −33 eV m 0 eV: Even though it cannot be the main component of dark matter in the mass range 10 −33 eV m 0 10 −17 eV, which is excluded by the fuzzy dark matter constraints [121][122][123][124][125][126]134] and the superradiance constraints from the supermassive black holes [158][159][160][161], such tiny amount of the ultralight scalar relic can still open the gate for the main component dark matter's production.
At the beginning of this section, let us give a brief introduction of the setup: In our model, the dark photon dark matter is produced through the operator L ⊃ φF µν F µν /2Λ KM whose effective kinetic mixing is supported by φ's VEV in the early universe.Thereafter, when T T osc , φ begins the oscillation with its amplitude scaled as |φ| ∝ (T /T osc ) 3/2 .Because most of the existing constraints are imposed when T T osc , the dark photon's parameter space can be vastly extended, and the ratio T 0 /T osc determines today's local kinetic mixing for the future dark photon dark matter detections.In addition, since the φF 2 or φ 2 F 2 operator is induced as the byproduct of the varying kinetic mixing from the UV theory discussed in Sec.II, testing the fine-structure constant variation and the equivalence principle violation through the ground-based experiments (torsion balance, clock comparison, resonant-mass detector, AION, MAGIS) [83,84,[87][88][89][90][91][162][163][164][165][166], satellite-based experiments (MICROSCOPE, Space-Q, AEDGE) [85,[167][168][169], astrophysics (PTA, SKA) [92,170], and cosmology (CMB, BBN, Lyman-α forest) [93][94][95][96]170] open another brand new window for the dark matter experimentalists.

A. Dark Photon Production
Let us begin with the dark photon freeze-in through the varying kinetic mixing.To simplify the discussion, we focus on the case where the dark photon is produced before the ultralight scalar's oscillation.When T T osc , φ is frozen with the nonzero displacement such that the dark photon is produced through γ → A and e − e + → A channels.Our numerical calculation of the required kinetic mixing to reach today's dark matter relic abundance is shown as the yellow line labeled as " FI " in Fig. 5.In the following text, let us briefly introduce how to estimate the dark photon abundance from the freeze-in process semianalytically.
In the case where m A < 2m e , we solve the Boltzmann equation In Eq. ( 36), n γ and n A are the number density of the photon and dark photon, respectively.Γ γ→A is the thermal average of the photon to dark photon transition rate, from which we know that the resonant oscillation, i.e., γ → A , happens when m γ m A .In the mass range 0.1MeV m A 1MeV, when γ → A happens, e ± are relativistic.Therefore, the plasmon mass can be approximately written as m 2 γ 2πα em T 2 /3.Combining m 2 γ and Eq. ( 36), we have Figure 5.The parameter space of the dark photon dark matter (Ω A h 2 0.12) with the constant kinetic mixing.Here, m A is the dark photon mass, and is the constant kinetic mixing.The yellow line labeled with " FI" is the dark photon freeze-in line, i.e., the kinetic mixing necessary for the dark photon to reach today's dark matter relic abundance.In the mass range m A < 2me, the dominant dark photon production channel is photon to dark photon resonant conversion, i.e., γ → A .When m A ≥ 2me, A → e − e + dominates, which explains the slight decreasing of the freeze-in line.In the plot, the purple region represents the constraints from the stellar energy loss of red giants, horizontal branches, and sun [58][59][60][71][72][73][74].The shaded red region labeled by "DD" denotes the constraints from the dark matter direct detection experiments [60,[75][76][77][78][79][80].Up to now, the most stringent constraint within keV − MeV range comes from the XENONnT experiment [80].The dashed red line is the projection of the next-generation LUX-ZEPLIN (LZ) experiment [171].The blue region denotes the constraints from the dark photon dark matter decay, whose dominant channel is A → 3γ when m A < 2me, or A → e − e + when m A ≥ 2me.Since the dark photon decay affects CMB and the late-time photon spectrum, constraints can be imposed on the dark photon parameter space [58][59][60]172].From [60,82,173], we know that the constraints from CMB and the late-time processes are comparable for the constant kinetic mixing model.Therefore, we label the constraint from the dark photon decay with "CMB+Late".In the plot, one can find that the freeze-in line is covered by the imposed constraints, which means that the constant kinetic mixing model is completely ruled out.Such exclusion of the dark photon dark matter freeze-in is formally discussed in [58][59][60].
where T γ→A is the resonant temperature of γ → A , T eq ∼ eV is the temperature of the matter-radiation equality, and Ω A , γ→A is the related dark photon relic abundance.In Eq. (37), Ω A 's dependence on m A is approximately canceled out.Therefore, the freeze-in line in the region where m A < 2m e is approximately a constant.Since Ω A h 2 0.12, from Eq. ( 37) we have which explains the behavior of the yellow line in Fig. 5 in the mass range 0.1MeV m A 1MeV.In the mass range m A 0.1MeV, the resonant oscillation happens when e ± are non-relativistic.To be more quantitative, T res ∼ 0.1m e , which is insensitive to m A .In this case, we have Ω A ∝ 2 m 3 A , which explains the slope of the freeze-in line in Fig. 5 in this mass range.In the discussion above, we focus on the transverse dark photon because the produced longitudinal dark photon is subdominant, which is checked by us numerically and consistent with [72].
In the case where m A 2m e , the inverse decay channel, i.e., e − e + → A , opens up and dominates over γ → A .

Now the Boltzmann equation can be written as ṅA + 3Hn
A n e − n e + σ e − e + →A , where n e − n e + σ e − e + →A ∼ 2 α em m 5/2 Inside the collision term of Eq. ( 39), there is an exponential suppression factor e −m A /T cutting off the dark photon production when the universe cools down.From Eq. ( 39), we estimate the cut-off temperature of e − e + → A near which most of the dark photons are produced and the related dark photon relic abundance as respectively.In Eq. ( 40), one can find that Ω A , e − e + →A is also almost independent of the dark photon mass, which is similar to the result in Eq. ( 37) but different by an α 1/2 em factor.For this reason, in the mass range m A 2m e , the needed kinetic mixing to reach today's relic abundance is slightly smaller than the one in the range m A < 2m e by an O(1) factor, which explains the slight lowering of the freeze-in lines in Fig. 5 and Fig. 6.
Being different from the model with the constant kinetic mixing, our model has extra dark photon production channels where φ is the particle rather than the VEV, such as γ → A φ and e − e + → A φ. Nevertheless, one can see that these channels are subdominant from the discussion below.By estimating the ratio of the A abundance produced through these channels over the A abundance produced through the varying kinetic mixing, we have where Ω γ→A represents the dark photon abundance from the channels in which φ serves as the VEV, and Ω γ→A φ is for the dark photon abundance from the channels where φ is the particle.From Eq. ( 41), one can find that the majority of the dark photon comes from the channels with the φ as the VEV, because the channels with φ as the particle are highly suppressed by the large Λ KM , while the channels with the φ as the VEV are compensated with large |φ| osc .

B. Signatures and Constraints
In our model, the experiments can be categorized into two classes: 1.The detections of the dark photon dark matter.These rely on the portal between the visible and dark sectors, i.e., the varying kinetic mixing through which the energy and momentum flow.Today's local kinetic mixing is set by the scalar mass m 0 .2. The detection of the ultralight scalar φ, i.e., the handle controlling the portal.Such detections are based upon the scalar-photon couplings, which vary the fine-structure constant or violate the equivalence principle.The targets for the experiments detecting the ultralight CP-even scalar are set by the dark photon dark matter freeze-in.

Detection of the Dark Photon Dark Matter
Let us first investigate the phenomenology of the dark photon dark matter established on the nonzero kinetic mixing.Before starting the discussion, one needs to recall that in our model, the kinetic mixing varies during the universe's evolution.Hence, the experimental detectability based on some specific process depends on the universe's epoch when this process happens.To understand this, we start by reviewing the ultralight scalar's evolution.As we know, when T T osc , φ is at rest with the nonzero field displacement.Afterward, when T T osc , φ begins the damped oscillation scaled as |φ| ∝ T 3/2 .Given that the dark photon is the main component of dark matter, we have FI ∼ 10 −12 , based on which today's local kinetic mixing, denoted as 0 , is set by T osc .From Sec.III, we write 0 as Figure 6.The parameter space of the dark photon dark matter (Ω A h 2 0.12) produced from the freeze-in mechanism through the varying kinetic mixing.Here, m A is the dark photon mass, and 0 is today's local kinetic mixing near the earth.One should note that the CMB (light blue) and BBN (green) bounds are imposed on the early universe with larger kinetic mixing.
Based on the evolution of φ, we recast these constraints to the m A − 0 plane.Left: m0 10 −25 eV.In the plot, we choose m0 10 −30 , 10 −28 , 10 −26 eV as the benchmark values.Here, the strongest constraints come from the stellar energy loss (purple), A decay during CMB (light blue) and the late time (dark blue), and the direct detection experiments (red).The dashed red line denotes the projection of LUX-ZEPLIN (LZ).The yellow line labeled by " FI" is the kinetic mixing required to reach the dark matter relic abundance during the freeze-in.If m0 10 −28 eV, one has Tosc TCMB, so CMB FI.For this reason, the dark photon heavier than 0.1MeV is excluded, while the lighter one is not.Right: m0 10 −25 eV.Here we choose m0 10 −25 , 10 −21 , 10 −17 , 10 −13 eV as the benchmark values.The light gray region denotes the parameter space where our calculation breaks up when the oscillation begins before the freeze-in process, i.e., Tosc TFI.Here, the most relevant constraints comes from the A decay during CMB (light blue) and the late time (dark blue), where A → 3γ dominates when m A < 2me, and A → e − e + dominates when m A 2me. Since the kinetic mixing is large during the BBN, the dark photon decay is also constrained by BBN (green).One should notice that the BBN bound is cut off at m A 5 MeV, below which the A decay cannot disintegrate the relevant light elements [174,175].Since the oscillation temperatures for the type-A and the type-B model have some slight quantitative differences in the mass range m0 10 −18 eV, here we mainly use the type-A model as an example.
where the universe temperature today is T 0 ∼ 10 −3 eV, the enhancement factor from the structure formation is E ∼ 600, as shown in Eq. ( 28), and the kinetic mixing during the freeze-in is FI ∼ 10 −12 , as shown in Eq. (38).Based upon Eq. ( 42), we divide the discussion into two parts separately, as shown in Fig. 6 Let us begin with the m 0 10 −25 eV case with the dark photon parameter space shown in the left panel of Fig. 6.One can see that the most relevant constraints come from the dark matter direct detection (red), the stellar energy loss (purple), and the dark photon decay during the CMB (light blue) and the late time (blue).In the plot, the yellow line represents the required kinetic mixing during the dark photon freeze-in to reach the relic abundance Ω A h 2 0.12, from which one can easily find that the constant kinetic mixing scenario is thoroughly excluded.In contrast, the varying kinetic mixing model opens the parameter space and provides the benchmark values determined by m 0 .Here, we choose m 0 10 −30 , 10 −28 , 10 −26 eV to plot the freeze-in lines on the m A − 0 plane.Since T osc ∼ (m 0 m pl ) 1/2 , larger m 0 leads to an earlier start of the oscillation, such that 0 is smaller.Let us briefly introduce the constraints appearing in the left panel of Fig. 6 in the following paragraphs.
In the mass range m A 0.1MeV, the most relevant constraints come from the dark matter direction (red) and the stellar energy loss (purple).The direct detection experiments test scintillation photons from the electron recoil inside the noble liquids or semiconductors due to the absorption of the dark photon [60,[75][76][77][78][79][80].Up to now, the most strict constraint within keV to MeV mass range comes from the XENONnT experiment [80].The LUX-ZEPLIN (LZ) experiment, represented by the dashed red line in the plot, can test smaller kinetic mixing by a half to one order of magnitude [171].Furthermore, we can expect that future dark matter direct detection experiments, such as DarkSide-20k [176] and DARWIN [177], with more giant detectors and lower backgrounds, can have better detection capability.In this mass range, our model is also constrained by the stellar energy loss of the sun, horizontal branch, and red giant via the channel γ → A , which changes the stellar evolutions [58][59][60][71][72][73][74].The direct detection of the non-relativistic dark photons produced by the sun, known as the solar basin, can also impose comparable constraints [178].
In the mass range m A 0.1MeV, our model is constrained by the decay of the dark photon dark matter.When m A < 2m e , since the two-photon channel is forbidden according to the Landau-Yang theorem [179,180], the dark photon decays through A → 3γ induced by the electron loop.When m A 2m e , the dominant channel is A → e − e + .These two channels are constrained by the observations of the CMB and the late-time photon background, which give comparable constraints in the constant kinetic mixing scenario.From [59,60,81,82], we know that Γ A →3γ 10 −9 H 0 and Γ A →e − e + 10 −7 H 0 for m A ∼ MeV.However, in our model, one should realize that these two physical processes happen in different stages of the universe.To be more specific, the stage of CMB is at T CMB ∼ eV while the galactic photon is emitted in today's universe.From the discussion of φ's evolution in Sec.III, we know that the kinetic mixing during CMB is where T osc ∼ (m 0 m pl ) 1/2 .Based on Eq. ( 43), we can recast the constraint on the kinetic mixing from the dark photon dark matter decay during the CMB stage to the m A − 0 diagram.When m 0 10 −28 eV, one knows T osc T CMB , so CMB FI .Given this, one can figure out that the dark photon mass region m A 0.1MeV is excluded by CMB, while the smaller dark photon mass region is not.This explains CMB bound's cutting off at m A 0.1MeV in the left panel of Fig. 6.When m 0 10 −28 eV, one has T osc T CMB , so the kinetic mixing during the stage of CMB is CMB ∼ 0 (T CMB /T 0 ) 3/2 , which explains the strengthening of the CMB bound at the range m A 0.1MeV compared with the constraints from the late-time photon.Now let us discuss the m 0 10 −25 eV case appearing on the right panel of Fig. 6.Here, the most relevant constraints come from the decay of the dark photon dark matter during CMB (light blue) and BBN (green).To recast the constraints in the early universe to the m A − 0 plane, we use the formula From the discussion in Subsec.III B, we know that thermal effect may change T osc for the type-B model when the scalar is heavier than 10 −16 / log 2 η eV.Even though the type-B model does not cause the qualitative differences, to simplify the discussion, we only discuss the type-A model as an example such that the oscillation temperature still satisfies T osc ∼ (m 0 m pl ) 1/2 in the experimentally allowed region.Similar to the discussion above, we recast the CMB bound to the m A − 0 plot.During the BBN, BBN ranges from 10 −12 to 10 −14 , which depends on the concrete value of m 0 .In this range, the dark photon decay through A → e − e + leads to the light elements' disintegration, based on which the BBN constraint is imposed [174,175].However, as long as m A 5MeV, the dark photon decay cannot change the light element abundance since the injected energy is smaller than the deuterium binding energy, which is the smallest among all the relevant light elements (except 7 Be whose abundance does not affect the main BBN observables).In the right panel of Fig. 6, there is a gray region in the lower left corner, which represents the parameter space where the calculations of freeze-in lines break because the freeze-in happens after the start of the oscillation.At the end of this paragraph, we point out one nontrivial character: Through the varying kinetic mixing, the dark photon heavier than 1MeV can be produced and avoid the constraints of the decay A → e − e + because the kinetic mixing portal fastly closes right after the dark photon production.
Before ending this subsection, we want to point out that since the dark photon is produced through γ → A with the initial momentum p A ∼ T , it washes out the dark matter substructure in the late universe, which is the character of the warm dark matter [181][182][183][184].To avoid this, we need m A few × 10 keV.Since imposing such constraints needs detailed analyses based on the dark matter phase space distribution, which deviates from the main point of this paper, we leave this to future work. .The parameter space of the scalar-photon couplings.The nonrelativistic scalar relic's fraction among the dark matter is F 10 −3 .The kinetic mixing before the oscillation is FI 2×10 −12 , which is the benchmark value for the freeze-in production of the dark photon dark matter satisfying Ω A h 2 0.12.Left: The type-A model (y = 0).We choose e = 1, 10 −3 , 10 −6 to demonstrate the experimental targets for the dark photon freeze-in through the varying kinetic mixing, which appear as the black lines.The thick dashed magenta line is the η 1 contour, upon which the scalar's early displacement is set to be πf /2 by the thermal misalignment.The gray region denotes the constraints from the equivalence principle tests (gray: EP tests) from the MICROSCOPE satellite [85] and the torsion balance experiments [83,84].The blue region represents current constraints from the clock comparisons [86-88, 90, 185].We can see, for the type-A model, most of the model's parameter space (e 1) is inside the projections of the proposed experiments, such as Lyman-α UVES (dashed purple) [170], the clock comparison (dashed pink: optical-optical, dashed red: optical-nuclear) [86], the cold-atom interferometer (dashed dark green: AEDGE, dashed green: AION-km, dashed light green: MAGIS-km) [165,166,168,186], and the resonant-mass oscillator (dashed orange: DUAL) [91].Right: The type-B model (y = −y) beginning with the wrong vacuum (πf /2 φ rh 3πf /2).We choose e = 10 −6 , 10 −8 , 10 −10 as the benchmark values to plot the black lines.The thick magenta line is the η 1 contour.In the region beyond such a magenta line, the scalar acquires the field displacement πf via the thermal misalignment.The thick dashed cyan line, i.e., the contour of mT |3H=m 0 m0, envelopes the region where the thermal effect postpones the oscillation starting, i.e., Tosc < T |3H=m 0 .This is the region where the ultralight scalar does the trapped misalignment, which is previously investigated in [140,141] for the axion models.Given the scalar abundance, since the trapped misalignment needs smaller oscillation amplitude, as discussed in Eq. ( 35), the power law of the black lines is changed within the upper right corner enveloped by the thick dashed cyan line.In the plot, one can find that the constraints from the cosmology (light orange: CMB, purple: BBN) are comparable with the current constraints from clock comparison experiments (blue) [86-88, 90, 162, 185] and the equivalence principle tests (gray: EP tests) [83][84][85].For the type-B model, one needs e 10 −6 for the dark photon dark matter freeze-in to be detected by future experiments.

Detection of the Ultralight Scalar
From the former discussion in Sec.II, we know that the scalar-photon coupling appears along with the varying kinetic mixing from the UV physics, so our model can be tested from another entirely different perspective, i.e., the phenomenology of the ultralight CP-even scalar.Depending on whether the Z 2 symmetry is exactly preserved or not, the scalar-photon coupling is φF 2 (Type-A), or φ 2 F 2 (Type-B).From Eq. ( 16), we know that the dimensionless coefficients d γ,i (i = 1, 2) contain the production of the e and Λ KM , and Λ KM can be written as Λ KM |φ| osc / FI .We also know that in the freeze-in model of the dark photon dark matter, there is FI ∼ 10 −12 .Therefore, the targets of the scalar experiments are set by the freeze-in of the dark photon dark matter.To be more quantitative, we have where "i = 1" and "i = 2" denote the type-A and type-B models, respectively.From Eq. ( 45), we know that with FI determined, the smaller e is, the larger d γ,i (i = 1, 2) are, so our model can be more easily to be detected.The reason is that from the UV model in Sec.II, to maintain a determined kinetic mixing, y ( ) should be turned up when e is turned down.Therefore, the coefficients in front of the scalar-photon couplings are turned up.Since the phenomenologies for the linear (type-A) and quadratic (type-B) scalar-photon couplings are quite different, we discuss these two models individually in the rest of this subsection.a. Type-A Model Let us first discuss the type-A model shown in the left panel of Fig. 7. To illustrate how the experimental signals vary with the dark gauge coupling, we plot the d γ,1 lines (black) choosing e = 1, 10 −3 , 10 −6 as the benchmark values.From Eq. ( 25) we know that T osc ∝ m 2/3 0 when m 0 10 −28 eV, and T osc ∝ m 1/2 0 when m 0 10 −28 eV.This explains the shapes of the black lines, which are represented as (m 0 10 −28 eV) const (10 −33 eV m 0 10 −28 eV) (46) In Fig. 7, the thick dashed magenta line is the contour of η 1.In the upper right corner beyond this contour, η is larger but still of the order of one to ten, so the "arctan" suppression in the type-A η 1 case mentioned in Sec.III B does not appear.According to our former discussion in Sec.III B, for the type-A model with η ∼ 1, the ultralight scalar's early field displacement is set to be πf /2 by the thermal misalignment.Since |φ| osc /f ∼ O(1), from Eq. ( 21) one can know the power law of the thick dashed magenta line satisfies d γ,1 ∝ m −1/4 0 .In the region below this thick dashed magenta line, we have η 1, so the ultralight scalar does the standard misalignment.One of the strongest constraints for the type-A model comes from the equivalence principle experiments testing the acceleration difference of two objects made by different materials attracted by the same heavy object.In the left panel of Fig. 7, such a constraint is shown as the shaded gray region.Until now, the most stringent constraint is imposed by the space mission named MICROSCOPE [85] and the torsion balance experiments named Eöt-Wash [83,84], which gives d γ,1 10 −4 .Since such acceleration difference is caused by the Yukawa interaction mediated by the ultralight scalar, these constraints are independent of the nonrelativistic scalar relic's fraction among the dark matter, i.e., F. The clock comparisons of Dy/Dy [87], Rb/Cs [88], and Al + /Hg + [90] also give stringent constraints based on testing the time-varying α em .These constraints are shown as the shaded blue region.
To go beyond the current constraints, there are several experiments proposed.For the future proposed clock comparison experiments [86], the projection of the improved optical-optical clock comparison is shown as the pink line, and the projection of the optical-nuclear clock comparison is shown as the red line.For these projections, the dashed parts denote the projection of the α em oscillation testing, and the dotted parts denote the projection of the α em drift testing.According to [86], the projection has the optimistic assumption that the measurement takes place when the scalar is swiping through the zero such that αem is independent of m 0 .Following this, we extrapolate the projections of the optical-optical and optical-nuclear experiments to 10 −32 eV considering the homogeneity of the ultralight scalar when the scalar's de Broglie wavelength is much larger than the size of the Milky Way halo.The cold-atom interferometer experiments such as AEDGE [168], AION [165], and MAGIS [166] have strong detection capability in the mass range 10 −19 eV m 0 10 −12 eV.In the plot, the projections of AEDGE (broadband, resonant mode), AION-km, and MAGIS-km are shown in the dashed dark green, dashed green, and dashed light green lines, respectively.The region of the thermal misalignment located in the upper right corner of the left panel of Fig. 7 can be tested by the proposed resonant-mass detectors, such as DUAL shown in the dashed orange line [91].The CP-even scalar in the mass range 10 −32 eV m 0 10 −28 eV can be tested via the proposed Lyman-α UVES observation [170] shown in the dashed purple line.
Since all the constraints and projections except the EP tests rely on the α em variation, and ∆α em ∝ d γ,1 |φ|, they are all scaled by F −1/2 for different ultralight scalar fractions among the dark matter.We know from Eq. ( 46) that, as F varies, the relative position between the experimental targets (black lines) and the constraints/projections from the non-EP experiments keeps the same.Knowing that e cannot be much larger than O(1) (Otherwise, the Landau pole is shifted to the low energy and the theory's perturbativity is broken), one can find from the left panel of Fig. 7 that most of the theory space is within the detection capabilities of the proposed experiments.
b. Type-B Model Now let us discuss the type-B model whose parameter space is shown in the right panel of Fig. 7. Here, we make φ to have the initial wrong vacuum (πf /2 φ rh 3πf /2).In the plot, the thick dashed magenta line is the contour obeying η 1.In the region upon (below) this line, the ultralight scalar does the thermal (standard) misalignment.From Eq. ( 21) we know that η ∼ (α em d γ,2 ) 1/2 .Thus, the thick dashed magenta line approximately obeys d γ,2 ∼ 10 2 .As discussed in Subsec.III B, the thermal effect makes φ converge to πf in the early universe in the parameter space upon the thick dashed magenta line.The thick dashed cyan line denotes the contour satisfying T osc T | 3H=m0 , or equivalently speaking, m T | 3H=m0 m 0 .In the region enveloped by the thick dashed cyan line, we have T osc < T | 3H=m0 meaning that the thermal effect postpones φ's oscillation, which represents the phase of trapped misalignment.From the plot, one can find that the thick dashed cyan line is approximately flat when m 0 10 −16 eV.Such flatness appears because m T m 0 happens when e ± are relativistic, which makes m T /H ∼ η.As long as η 1, the starting of the oscillation is postponed by the thermal effect 4 .One can also find that the left side of the thick dashed cyan line becomes vertical, because m T m 0 happens when e ± are non-relativistic.In this case, the condition for T osc < T | 3H=m0 becomes m 0 m 2 e /m pl .The experimental benchmarks are shown as the black lines in m 0 − d γ,2 plane.In the region of the standard misalignment (η 1), the black lines obey where m * denotes the abscissa of the cross point of the black line and the thick dashed cyan line.To understand the power law of d γ,2 in Eq. ( 47), one can refer to Eq. ( 25) which tells that T osc ∝ m 2/3 0 when m 0 10 −28 eV, and when m 0 10 −28 eV.After plugging T osc into Eq.( 45), one can have Eq.( 47).Now let us discuss the behaviors of the black lines in the scalar mass range m 0 m * .If the black line intersects with the thick dashed cyan line on its vertical side, i.e., m * 10 −16 eV, one can have T osc ∼ m e .If the black line crosses with the thick dashed cyan line on its horizontal side, i.e., m * 10 −16 eV, we have Substituting the equations of T osc into Eq.( 45), we have which explains the black lines' tilting up when entering the region of the trapped misalignment, i.e., the area enclosed by the thick dashed cyan line.Such enhancement of the signal comes from the postponement of the oscillation which is the major feature of the trapped misalignment.Similar to the type-A model, the type-B model can be tested through terrestrial experiments.In the plot, the current constraints based on clock comparison [87,88,90,162,163] are shown as the shaded blue region, and the constraints from the equivalence principle tests [85,162,163] are shown as the shaded gray region.The projections of the proposed optical-optical and optical-nuclear clock comparisons are shown in orange and red, respectively [86,163].Here, the dashed and dotted parts of the projections denote the detection capabilities of the α em oscillation and the α em drift, respectively.In addition, our model can also be tested by the cold-atom interferometers.The projection of AEDGE [168] and AION-km [165] are shown in the dashed dark green and dashed green lines, respectively.One may find that the constraints and projections from the ground-based and low-altitude experiments get weakened or cut off in the strong coupling region.This is caused by the scalar's matter effect sourced by the earth, one of the characteristics of the quadratic scalar-SM coupling.Following [162,163], we have where d γ,2,crit is the scalar's critical value for the matter effect to appear, R ∼ 6 × 10 3 km is the earth's radius, M ∼ 6 × 10 24 kg is the earth's mass, and Q γ, ∼ 10 −3 is the earth's dilaton charge.When d γ,2 d γ,2,crit , the scalar's Compton wavelength is smaller than the earth's radius.In this situation, the scalar easily overcomes the spatial gradient and is pulled toward the origin, so its near-ground and underground oscillation amplitudes are highly 4 Here, one may be curious about the nonoverlapping between the thick dashed cyan and thick dashed magenta lines when m 0 10 −16 eV.
To understand such splitting, one needs to review the exact definition of η in Eq. ( 20), which is η := 2m T /H (Here S 1).When 3H = m 0 , to make m T > m 0 , one needs η > 6.Since η ∝ d suppressed.Even so, if the experiments are carried on in space with altitudes comparable with the earth's radius (for example, AEDGE), the screen effect sourced by the earth can be largely alleviated.
Unlike the type-A model, the type-B model has strong constraints from the early universe processes, such as CMB and BBN, which are comparable with the terrestrial constraints.The reason for this character is that when tracing back to the early universe with the φ's abundance determined, ∆α em increases much more for the type-B model than the type-A model.For the former analysis of the CMB and BBN bounds on the quadratic scalar-photon coupling, one can refer to [93,95,96] 5 .Considering the thermal effect in Subsec.III B, we impose the constraints on the type-B scalar's parameter space using the current cosmological constraints on ∆α em /α em , which are [93,94] 6 In the right panel of Fig. 7, the CMB constraint is shown as the shaded orange region, and the BBN constraint is shown as the shaded purple region.
Based on the left panel of Fig. 7, let us discuss the detectability of the dark photon dark matter freeze-in via the varying kinetic mixing.Given the scalar fraction F 10 −3 , for the model to be tested by the proposed experiments but not excluded by CMB/BBN, the dark gauge coupling needs to be within the range 10 −10 e 10 −6 .When F varies, the relative position between the black lines and the non-EP experiments does not change qualitatively because all these constraints essentially come from ∆α em ∝ d γ,2 |φ| 2 .In this case, their positions are rescaled by F −1 in the weak coupling region.Differently, what the equivalence principle experiments are testing is the acceleration difference of the objects A and B, which φ|∇φ| according to [162,163].From this, we know that the constraints on d γ,2 from the equivalence principle test should be rescaled by F −1/2 in the weak coupling region.When F 10 −8 , the equivalence principle constraints are stronger than the current CMB, BBN, and clock constraints.We leave detailed discussions on the F-dependence of the constraints and theoretical benchmarks to our future work.
At the end of this section, we want to comment on the constraints imposed by the black hole superradiance.For the F 1 case, the current constraints from the supermassive black holes exclude the region 10 −21 eV m 0 10 −17 eV [158][159][160][161]187], and the constraints from the solar mass black holes exclude the region 10 −13 eV m 0 10 −11 eV [188].However, since φ's self-interaction is λ φ ∼ m 2 0 /f 2 , smaller F leads to smaller f , which increases λ φ .From [161,188] we know that, for the scalar as the subfraction of the dark matter, the superradiance constraints are alleviated by the scalar's large attractive self-interaction.

V. ZN -PROTECTED SCALAR NATURALNESS
Since the Yukawa interaction in Eq. ( 4) breaks Φ's global U (1) symmetry, the ultralight CP-even scalar φ suffers from large quantum correction, which is a common problem in many ultralight scalar models.Taking the type-B model in Sec.II as an example, we have the mass correction ∆m φ ∼ yM , and its benchmark value which is much larger than the scalar's bare mass.For the type-A model, the Coleman-Weinberg potential is V cw ∼ yM 3 f sin(φ/f ).Given that r 1, the type-A model's mass quantum correction ∆m φ ∼ yM/r 1/2 is even much larger.Because V cw 's minimum of the type-A model is not at zero, one need to introduce additional tadpole counterterm to forbid φ to have the late-time nonzero displacement.
Although such a naturalness problem exists in the minimal model, by imposing an extra Z N symmetry [97,98], the global U (1) symmetry can be approximately restored, so the quantum correction is exponentially suppressed.In such setup, for the type-A and type-B models, the mass fine-tuning is exponentially solved.In addition, for the type-A are all coupled with the complex scalar Φ, where "k" labels the kth world and k = 0 is our world which experienced the reheating.Because of the suppression from y ( ) , kth and jth universes (k = j) do not talk with each other.Here, the ZN rotation Φ → Φ exp(i 2π/N ), Ψ k+1 is labeled by the gray arrowed circles.The Z2 symmetry between Ψ and Ψ k , which is exact for the type-B model (y = −y) and approximate for the type-A model (y = 0, y = 0), eliminates the time-independent part of the kinetic mixing.model, the one-loop tadpole is also negligible.To realize this, we introduce N copies of the worlds containing the standard model sector (SM k ), the dark sector (DS k ), and the portal (O P k ) where k = 0, 1, . . ., N − 1.Because of the Z N symmetry, the system is invariant under the transformation so the lowest order effective operator is which is invariant under the Z N transformation but not invariant under the U (1) transformation.Such dimension-N operator shows that even though the global U (1) symmetry is broken, as long as the symmetry under the discrete subgroup Z N still exists, the quantum correction ∆m . As long as ∆m 2 φ m 2 0 , the scalar's bare mass m 0 softly breaking the Z N symmetry can be naturally small.

A. ZN -Protected Model
To build a concrete varying kinetic mixing model without the naturalness problem, we embed the simple model described by Eq. ( 4) into the Lagrangian being invariant under the Here, Ψ k and Ψ k , two doubly charged fermions in the k-th universe, carry the same kth-hypercharge but the opposite kth-dark hypercharge.Following the phase choice in Sec.II, we choose c = π/2, Φ = if √ 2 e iφ/f and introduce the Z 2 -invariant bare potential in Eq. ( 10) to ensure the late-time closing of the kinetic mixing portal.When the ultralight scalar φ is nonzero, Ψ To gain more intuition of the setup of the Z N -invariant model, our readers can refer to Fig. 8, the schematic diagram of the Z N UV Lagrangian described by Eq. ( 54).In this figure, Φ is the axle attached to all N copies of the universes, and "Z 2 " refers to Ψ ( ) k 's mass degeneracy which is exact when y = −y (Type-B) or approximate when y y (Type-A).
In the low energy limit, the heavy fermions Ψ ( ) k are integrated out, so the IR Lagrangian becomes where and From Eq. ( 56), Eq. ( 57) and Eq. ( 58), one can easily find that the IR Lagrangian is invariant under the k+1 .Here, (SM + O P + DS) 0 is our universe which experiences the reheating, but the other universes are not reheated at all.When m 0 H, φ begins the damped oscillation, so the kinetic mixing 0 between SM 0 and DS 0 gradually decreases towards zero as described in Sec.II and Sec.III.

B. Quantum Correction of φ
For the ultralight scalar φ in Eq. ( 54), the leading order quantum correction can be described by the one-loop Coleman-Weinberg potential As discussed in [103], the Fourier series of any scalar potential respecting the Z N symmetry only receives contributions on lN th modes (l is a positive integer number) so that the Coleman-Weinberg potential can be written as In Eq. ( 60), V cw ( N ) denotes the Fourier coefficient of the single-world Coleman-Weinberg potential, the prefactor N comes from the contribution from N worlds which have the equal contribution to the Fourier coefficient, and cos(lN θ) comes from the effective operator L ⊃ const × (Φ lN + Φ † lN )/Λ lN −4 l .Now, let us calculate the Fourier coefficient in (60) through the equation Expanding the polynomial (1 − r cos θ) 4 and carrying the integration over θ, the Fourier coefficient when N > 4 can be exactly written as 9 In the Appendix.C, we list two ways to derive a compact formula for the Z N Coleman-Weinberg potential precisely expanded to all orders of r, which is obtained for the first time as far as we have known.After omitting r's higher order terms in the exact Fourier coefficient (63), we express the Coleman-Weinberg potential as where "• • • "s in the brackets of (64) denote r ( ) 's higher order terms with the same parity as the lowest order terms.In (64), the function G(n) is Alternatively, Eq. ( 64) can be derived from Eq. ( 59) by applying the cosine function sum rules According to Eq. ( 55), M k is a cosine function of θ, so we can expand the Coleman-Weinberg potential in Eq. ( 59) as a polynomial function of cos(θ + 2πk/N ).Afterward, by employing these cosine sum rules, we can also write the Coleman-Weinberg potential in the form as shown in Eq. (64).The concrete formulas of the constants C lmN and D m can be found in Appendix.C, Eq. (C5).Eq. (66) shows that only when the cosine function's power in the effective potential is greater than N , the non-constant terms emerge.Since the cosine function appearing in the potential is always accompanied by r, based on the cosine sum rules, one can find that the lowest order Z N potential is proportional to r N if there is no further cancellation coming from the exact Z 2 symmetry. 9In [189] we find the integral Given that r 1, the factor r N G(N ) in Eq. ( 64) indicates that the quantum correction of m 2 φ is suppressed by the (r/2) N −5/2 factor in the type-A model, and (r/2) N −2 factor in the type-B model.For the freeze-in of dark photon dark matter, r's benchmark value is r ∼ 10 −10 In such a case, for e ∼ O(1), as long as N 7, φ's naturalness problem can be solved entirely.For smaller e , r is larger, so N needs to be larger to solve the naturalness.However, as long as e 10 −10 which is in the permitted region of the current constraints, r 1, the Z N scenario protecting the pseudo Goldstone boson φ's naturalness always works.
Since V cw receives contribution from both r and r , N 's parity matters.For even N , the first term in Eq. ( 64) shows that both r and r provide a positive contribution, so the leading order correction starts from r N .In contrast, if N is odd, r and r have opposite signs, so the quantum correction from r N + r N is reduced.When r = −r (Type-B model), the (r N + r N ) cos[N (φ/f + π/2)] term has the exact cancellation, so the leading order correction starts from the (r 2N + r 2N ) cos[2N (φ/f + π/2)] term.We want to remind our readers that such cancellation happens in all orders because the front factor of the higher order terms cos[N (φ/f + π/2)] in Eq. (63) or Eq.(C1) are r ( )N +2j which has the same parity as r ( )N .This is consistent with the EFT analysis at the beginning of this section.
Before closing the discussion in this section, we want to make some extra comments besides solving the ultralight CP-even scalar's naturalness problem in our concrete situation: Our model, i.e., Z N -protected varying kinetic mixing, provides a systematic approach to solve the scalar's naturalness problem for both the linear (Type-A) and quadratic (Type-B) scalar-photon couplings.This is the extension and unification of [102] and [163], which focus on the linear and quadratic scalar-photon coupling, respectively.Furthermore, armed with two different methods, i.e., the Fourier transformation and the cosine sum rules discussed in this section and in Appendix.C, for the first time, one can completely and precisely calculate the pseudo-Goldstone boson's Z N -invariant Coleman-Weinberg potential previously appearing in [97,98,102] to arbitrary orders.We wish such calculations could shed light on the formalism of the Z N -naturalness itself.

VI. VARYING KINETIC MIXING FROM DIRAC GAUGINO
Motivated by stabilizing the hierarchy between the light scalars and the heavy fermions, we discuss one of the possible supersymmetric extensions of the varying kinetic mixing in the Dirac gaugino model [190,191] with the superpotential In Eq. ( 68), W j is the gauge field strength of the SM gauge group G SM,j where the label j = 1, 2, 3 denotes the SM gauge groups U (1) Y , SU (2) L and SU (3) c , respectively, W is the gauge field strength of U (1) d , and A j is the chiral multiplet.The operator in Eq. ( 68) in hidden sector models was firstly introduced by [192] and further understood by [190] afterward as the supersoft operator such that it provides the Dirac gaugino masses and does not give logarithmic divergent radiative contributions to other soft parameters.Writing A j , W j and W in terms of the Taylor expansion of the Grassmann variable θ, we have A j ⊃ (S j + iP j ) / √ 2 + √ 2θã j , W α j ⊃ λ α j + F µν j (σ µν θ) α + D θ α , and W α ⊃ F µν (σ µν θ) α + D θ α .Plugging these expanded fields into Eq.( 68), and doing the integration over θ 2 firstly and then the auxiliary fields, we have the effective Lagrangian where S j (P j ) is the CP-even (CP-odd) scalar, λ j is the gaugino, and m D,j = D /Λ D,j is the gaugino's Dirac mass given by the scale of SUSY breaking.One should note that in Eq. ( 69), S j 's mass term L ⊃ −2m 2 D,j S 2 j comes from integrating out the D-term (On the contrary, P j has no extra mass contribution).It is because the supersymmetry is protected that the mass of S j is correlated with the mass of the gaugino λ j .Consequently, S j is pushed to the heavier mass range given that the gaugino mass is highly constrained: According to [193][194][195][196], the LHC has already excluded the electroweakinos and the gluinos masses below O(100GeV) and O(TeV) respectively.Unlike the previous discussion, in the Dirac gaugino model, S j cannot be the ultralight scalar where the misalignment mechanism provides a natural way to open the portal in the early time but gradually close it in the late time.Even so, it is still possible to realize the temporary period of S j = 0 in the early universe through the two-step phase transition, also referred to as the VEV Flip-Flop in some specific dark matter models [39][40][41].The concrete model building and the phenomenology are beyond the scope of this paper.
In the j = 1 case, the first term in Eq. ( 69) containing the CP-even scalar S 1 corresponds to the kinetic mixing portal between U (1) Y and U (1) d which is determined by S 1 's VEV.The second term in Eq. ( 69) containing P 1 (axion) leads to the dark axion portal which is investigated in [197][198][199][200][201][202][203][204][205][206][207].There are several major differences between the kinetic mixing portal L ⊃ −S 1 F µν F µν /Λ D,1 and the dark axion portal L ⊃ −P 1 F µν F µν /Λ D,1 : 1.In the Dirac gaugino model, S 1 's mass is correlated with the λ α 1 mass which is pushed to O(100GeV) scale by LHC, while P 1 's shift symmetry protests its arbitrarily small bare mass.2. The VEV of P 1 does not play a direct physical role because it only contributes to the total derivative term in the Lagrangian (P 1 's time or spatial derivative still has nontrivial physical effects nonetheless).
In the j = 2, 3 cases, if S a j = 0, the first term in Eq. ( 69) would mix the non-Abelian gauge field with the dark photon such that the non-Abelian gauge symmetry is broken.Being referred to as the non-Abelian kinetic mixing [21,[208][209][210][211][212][213], the constant mixing models are highly constrained by the collider experiments.However, in the high-temperature environment of the early Universe, the large non-Abelian kinetic mixing can possibly be realized for the non-Abelian vector dark matter production and other intriguing phenomena.We leave the detailed discussion in future work.

VII. OTHER COSMOLOGICALLY VARYING PORTALS
Let us begin with the general form of the cosmologically varying portals through which the dark and the visible sectors are connected.To be more generic, we write them as In Eq. ( 70), O SM and O DS are the operators of the visible sector and the dark sector, respectively, d SM and d DS denote the dimensions of these two operators, and d is the dimension of the time-varying portal.To simplify the notation of Eq. ( 70), we drop the (spacetime, spin, flavor, . . . ) indices of O SM and O DS whose contraction makes the varying portal to be a singlet.For simplicity, we only keep the linear form of the CP-even scalar φ, even though in the UV theory, the non-linearity may appear, as we have seen in Eq. (7).Based on the EFT, we know that when the effective operator Eq. ( 70) is introduced, the operators merely containing φ and O SM also appear because the symmetry does not forbid them.The co-appearance of the effective operator shown in Eq. ( 70) and the scalar-SM coupling provides an excellent chance to test these kinds of models from the experiments detecting the portal itself and the ones measuring the scalar-SM coupling.In the rest of this section, we will give some specific examples of the varying portals and show how these minimal extensions illuminate the dark matter model building.
Let us briefly review the varying kinetic mixing portal in the EFT language.After choosing Eq. ( 70) goes back to the operator L ⊃ φ F µν F µν /Λ discussed before.Through this operator, the dark photon dark matter can be produced without violating the stringent constraints as shown in Sec.IV.Since the spacetime indices of F µν need to be contracted, the lowest order operator of the scalar-SM coupling is φ F µν F µν .If there is an exact Z 2 symmetry invariant under the dark charge conjugation φ → −φ, F µν → −F µν , φF µν F µν is forbidden, so the lowest order operator of the scalar-SM coupling becomes φ where s is a scalar singlet in the dark sector.Here, L ⊃ λ s O DS O SM = λ s s 2 |H| 2 is well-known as the singlet-scalar Higgs portal (SHP) through which the dark matter s reaches today's relic abundance (The dominant channels are ss → f − f + , W − W + , ZZ, hh, • • • .Here, f refers to the SM fermions.)[10][11][12].Besides, because this Lagrangian is invariant under the Z 2 transformation s → −s, s is stable.Although SHP provides such a simple and effective way to generate the scalar dark matter in accordance with the CMB measurement (Ω s h 2 0.12), in the mass range m s 1TeV, most of its parameter space except the narrow window of the resonance (m s m h /2) is excluded by the Higgs invisible decay h → ss [214,215], the dark matter direct detection experiments (LUX, XENON1T, PandaX-II) and the indirect detections (AMS, Fermi) [63,[216][217][218].By introducing the time-varying SHP as the minimum extension, the parameter space is widely extended.There, φ's misalignment supports s's freezeout at the early universe and then starts the damped oscillation such that σv ss ∝ (T /T osc ) 3 when T T osc .For this model, there are two types of experiments: 1.The future direct detection experiments (XENONnT, LUX-ZEPLIN, DarkSide-20k, DARWIN) which rely on partially opened SHP today.2. The table-top experiments and the astrophysical observations testing the φ|H| 2 or φ 2 |H| 2 operator (See [136,186,219,220] for the detections of the φ|H| 2 operator).
Another example of the singlet O SM is that In this model, s is the scalar mediator which interacts with the dark matter particle χ via the CP-odd coupling L ⊃ iy χ s χγ 5 χ (Here we denote the scalar mediator by "∼" to distinguish it from the scalar dark matter discussed in the previous paragraph).Given the constant portal L ⊃ O DS O SM /Λ = y f s f f where y f = v h / √ 2Λ, χ reaches today's relic abundance through the freezeout channel χ − χ + → f − f + whose cross section is σv ss ∼ y 2 f y 2 χ m 2 χ /m 4 s .Here, v h is the Higgs VEV, and Λ is the effective scale of the constant portal.Since χ − χ + → f − f + is s-wave (The leading order of σv χ − χ + is independent of dark matter relative velocity), the region m χ 10 GeV is excluded by the measurement of CMB anisotropy [3].To produce lighter but CMB-friendly dark matter, we introduce the operator as the minimal extension.Through the φ-dependent Yukawa coupling L ⊃ y f (φ)s f f where y f (φ) = φv/ √ 2Λ 2 and the aforementioned CP-odd Yukawa coupling L ⊃ iy χ s χγ 5 χ, the dark matter lighter than 10 GeV can reach today's relic abundance without violating CMB annihilation bound as long as φ's damped oscillation begins earlier than the CMB stage (T CMB ∼ eV).For this model, there are two kinds of experiments: 1.The next-generation CMB experiments (such as the CMB-S4 [221][222][223][224][225]) and the future direct and indirect detection experiments.2. The table-top experiments and astrophysical observations testing the variation of the SM fermion masses (φ f f or φ 2 f f operator) [86,88,91,92,162,163,185,186].

VIII. CONCLUSION
In this work, we study the scalar-controlled kinetic mixing varying with the ultralight CP-even scalar's evolution during the cosmological evolution.Via such a dynamical portal, keV − MeV dark photon dark matter is frozenin, free from the late-universe exclusions, and has the kinetic mixing set by the scalar mass.More generally, this model can serve as the minimal solution for the tension between the portal dark matter's early-time production and late-time constraints.Furthermore, the scalar-photon coupling inevitably emerges from UV physics, modifying the universe's thermal history, varying the fine-structure constant, and violating the equivalence principle.These phenomenons provide excellent targets to test our model via the ultralight scalar experiments in the mass range 10 −33 eV m 0 eV.In the meantime, the scalar mass also determines the target values of the kinetic mixing for keV − MeV dark photon dark matter experiments.In the rest of the paper, we study the Z N -protected naturalness of the ultralight CP-even scalar, the varying kinetic mixing from the Dirac gaugino model, and the dark matter via other cosmologically varying portals.
We begin our study by constructing the minimal model of the scalar-controlled kinetic mixing arising from the heavy doubly charged messengers.In our model, the Z 2 symmetry under the dark charge conjugation is imposed to elimite the constant kinetic mixing but preserve the varying kinetic mixing.From the perspective of symmetry, the early-time openning and the late-time closing of the kinetc mixing portal results from the system's inverse Z 2breaking.Importantly, in such a setup, the scalar-photon coupling also occurs from the top-down theory.To classify the model, we name the theory with the approximate Z 2 symmetry as the type-A model and the theory with the exact Z 2 symmetry as the type-B model.Given this, the type-A model has linear scalar-photon coupling, and the type-B model has quadratic scalar-photon coupling.Then we systematically categorize the scalar's evolution under the combined effects of the thermal correction, bare potential, and the universe expansion.For the type-A and type-B models, there exist the experimentally allowed regions for the thermal misalignment, where the scalar acquires the nonzero displacement insensitive to the initial condition.We also find the thermal postponement of the oscillation starting in the type-B model initially staying in the wrong vacuum, which is the major feature of the trapped misalignment.Such delaying of the oscillation enhances the experimental signals.Afterward, we study keV − MeV dark photon dark matter production via the varying kinetic mixing, where the kinetic mixing is FI ∼ 10 −12 for the dark photon freeze-in when T ∼ m A and decreases obeying the power law T 3/2 when the scalar does the lateuniverse damped oscillation.Therefore, the late-time kinetic mixing set by the scalar mass is free from constraints from stellar energy loss, dark matter direct detection, and the dark photon dark matter decay.In addition, the kinetic mixing during freeze-in provides the perfect experimental benchmarks in the ultralight scalar experiments in the mass range 10 −33 eV m 0 eV, including the test of the fine-structure constant variation and the equivalence principle violation.Since smaller dark gauge coupling leads to stronger signals, most of the e space in the type-A model and the region e 10 −6 in the type-B model can be tested by proposed experiments such as clock comparison and cold-atom interferometry.For the smaller scalar relic fraction among dark matter, the constraints from the finestructure constant variation keep the same, while the constraints from the equivalence principle experiments become more stringent.
We also detailedly study the Z N -protection of the ultralight CP-even scalar's naturalness in the varying kinetic mixing model.There, we embed the aforementioned minimal model into the model containing N worlds, such that the broken U (1) shift symmetry is approximately restored.Generally, the lowest order potential in this setup comes from the Φ N /Λ N −4 operator, which exponentially suppresses the scalar's quantum correction from heavy messengers.As long as N ∼ 10, the quantum correction of the scalar mass can be much lighter than 10 −33 eV.For the type-B model with odd N , the lowest order operator is Φ 2N /Λ 2N −4 because of the cancellation from the exact Z 2 .Moreover, we expand the Z N Coleman-Weinberg potential to all orders utilizing the Fourier transformation and the cosine sum rules.In the paper's final part, we discuss the Dirac gaugino model simultaneously inducing the varying kinetic mixing and dark axion portal.We also briefly discuss the dark matter models via other cosmologically varying portals.More specifically, we discuss the freeze-out of the dark Higgs dark matter lighter than 1TeV via the varying Higgs portal and the CMB-friendly fermionic dark matter lighter than 10GeV via the varying fermion portal. .This plot describes the ultralight scalar's evolution in the high-temperature environment and compares the numerical (solid lines) and analytical (dotted dark lines) results listed in Eq. (A4), Eq. (A5), and Eq.(A6).In the plot, similar to Fig. 4, T is defined as T := T /T |3H=m 0 .|δφ| represents the absolute value of the deviation from the thermal potential's minimum.Without loss of generality, we choose the type-A model as an example in this plot, and the discussion of the type-B model is quite similar.Here, the red, green, and blue colors represent the η > 1, η = 1, and η < 1 cases, respectively, and one can find that the numerical and analytical results for each case match magnificently well.We find that the red and green lines have approximately the same power law, i.e., |δφ| ∝ T 1/2 , whose slight deviation can be explained by the log term in Eq. (A5).We can also find that the red line's oscillation period appears the same in the plot with the log-scaled 1/ T axis.This is because the only temperature-dependence in Eq. (A6) comes from log T , or, more intuitively, is caused by the same temperature power law of mT and H. small η, the latter solution contributes to the nonzero initial velocity, while the former one does not.The η = 1 case is the critical point in the parameter space which separates the sliding and oscillating phases.In this case, the scalar moves toward the thermal minimum obeying the relation δφ ∝ T 1/2 approximately.When η goes beyond the critical value, i.e., in the case where η > 1, the scalar's evolution is in the oscillating phase whose amplitude dwindles like |δφ| ∝ T 1/2 as the universe expands.From Eq. (A6), we can find that the front coefficient of the cosine function dominates over the sine function as the result of the initial condition φrh = 0. Another feature of the oscillating solution for the η > 1 case worthwhile to be discussed is that the T -dependent part is the log function.Such a log term comes from the imaginary part of T 's power, or, in other words, is the feature of the homogeneous linear equation, which takes its form because m T and H have the same temperature power law.
Those who are curious about comparing the analytical and numerical solutions can refer to Fig. 9, where the solid lines represent the numerical solutions, and the dotted dark lines denote the analytical solutions, as shown in Eq. (A4), Eq. (A5), and Eq.(A6), in the linear approximation.The red, green, and blue colors are related to the η > 1, η = 1, and η < 1 cases, respectively.Since the type-A and the type-B models have similar behaviors when moving toward the thermal minimum in the high-temperature universe, we take the type-A model as an example to draw the plot, and the numerical-analytical comparison for the type-B model can be identically transplanted.In Fig. 9, one can easily see that the analytical solutions listed in Eq. (A4), Eq. (A5), and Eq.(A6) are perfectly consistent with the numerical results.production is dominated by the resonant transition γ → A .When m A 2m e , the dark photon mainly comes from the inverse decay e − e + → A .where "res" means all the temperature-dependent quantities in Eq. (B2) are chosen to be the ones when the plasmon mass and the dark photon mass match with each other.In Eq. (B2), s 0 is today's universe entropy, ρ c is the critical density, and d log m In Eq. (B2), m 3 A and T 3 res cancel with each other, so for fixed dark photon relic abundance, is the constant of m A .If all the dark matter is comprised of the dark photon, the needed kinetic mixing is

(B6)
For m A 10 −2 eV, γ → A takes place at the temperature T 10keV.Since the symmetric e − e + annihilates away, the plasmon mass comes from the asymmetric e − whose number density is conserved.In this case, the plasmon mass is scaled as m γ ∝ n e ∝ T 3/2 , and the resonant temperature is scaled as T res ∝ m (B10) 10 Here Tres is T γ→A in Eq. (37).We use the lower index "res" to emphasize the character of the resonance transition for γ → A . 11Even though the warm dark matter bounds exclude the 100% dark photon dark matter via freeze-in when m A few × 10keV, the subcomponent dark photon dark matter can still be produced.In this case, one can rescale the value of FI by Ω substituting it into Eq.(B10), using the Maxwell-Boltzmann distribution to approximate the electron/positron phase space, and doing the phase space integration as shown in [226], we have From Eq. (B13), we know that the to reach dark matter's relic abundance when m A 2m e is nearly the constant of m A , which is similar as Eq.(B4) but a bit smaller.

Figure 1 .
Figure 1.Left: The schematic diagram of the cosmologically varying kinetic mixing.The dark and standard model sectors are connected through the kinetic mixing controlled by the CP-even scalar φ, the subcomponent of the dark matter in the late universe.Based on this model, the energy flows from the dark sector to the standard model sector in the early universe for the dark matter production with the portal opened.In the late universe, with the portal partially closed, the dark matter is safe from stringent constraints.Right: The cosmologically varying kinetic mixing in the UV and IR theories.In the UV theory, there are heavy messengers Ψ and Ψ , both charged under U (1)Y and U (1) d .Here, Ψ and Ψ carry the same electromagnetic charge but the opposite dark charge.To generate the varying kinetic mixing, Ψ and Ψ are coupled with the scalar Φ via the Yukawa interactions y and y .To eliminate the constant kinetic mixing, the mass degeneracy between Ψ and Ψ is imposed, which is protected by the Z2 symmetry.In the IR theory, the nondynamical heavy messengers are integrated out, which induces the varying kinetic mixing (φF F ) and the scalar-photon coupling (φF 2 or φ 2 F 2 ).In the nonperturbative region, φ should be replaced by f sin(φ/f ) according to Eq. (7) and Eq.(13) because it is the angular part of Φ.Based on whether the Z2 symmetry is slightly broken by the Yukawas or not, we classify the theory into two types with totally different phenomenologies: The type-A model with the linear scalar-photon coupling (y = 0, y = 0) and the type-B model with the quadratic scalar-photon coupling (y = −y).These scalar-photon interactions lead to nontrivial cosmological history caused by the thermal effect and the signals from the αem variation and the equivalence principle violation.

Figure 2 .
Figure2.The evolution of mT (red), m0 (blue) and H (orange) as the universe expands.Because H and mT are both proportional to T 2 when T me, the orange and red lines are approximately parallel with each other in the log-log diagram until T ∼ me.Here, η ∼ mT /H when T me, which is already defined in Eq. (22).Left: η 1.Because the thermal effect from the SM bath is negligible, the ultralight scalar's evolution can be classified as the standard misalignment.In this case, φ stays constant at the early stage and begins the oscillation when H ∼ m0, which is labeled as Tosc.Right: η 1.The thermal effect is important for the scalar's evolution.In the plot, point T * and point "Q" denote the moment when mT = m0 and H ∼ mT , respectively.Here, HQ ∼ 10 −16 / log 2 η eV.When T T * , φ converges to the minimum of the thermal potential obeying the power law |φ − φmin| ∝ T 1/2 where φmin = πf /2 for the type-A model and φmin = ±πf for the type-B model.We merely show the m0 HQ case, because if m0 HQ, φ's movement is nothing more than the late-time standard misalignment with the initial condition fixed by the early time thermal effect.For the type-A model, when T me, φ tracks the potential minimum φmin = f arctan m 2 T /m 2 0 .As the temperature drops far below the electron mass, φ cannot track the minimum and begins to oscillate.For the type-B model, we focus on the πf /2 φ rh 3πf /2 case.When T T * , φ is trapped inside the local minimum φmin = πf .Thereafter, when T T * , φ begins oscillating around the bare minimum φmin = 0.Here we have T * T |3H=m 0 , which means the oscillation is postponed.Therefore, the scalar is in the phase of trapped misalignment.

Figure 4 .
Figure 4. φ's typical numerical solutions.In these two plots, we define the dimensionless quantity T := T /T |3H=m 0 , so at point T = 1 we have 3H = m0.The red, green, and blue lines represent the η 1, η ∼ 1, and η 1 cases, respectively.Left: The type-A model (y = 0).For the red line (η 1), when T T * , φ does the damped oscillation scaled as |φ − πf /2| ∝ T 1/2 , so it converges to πf /2.Remembering the T * 's definition in Eq. (22), we label the point where T T * as "mT m0".As the temperature decreases, φ tracks potential's minimum φmin = f arctan m 2T /m 2 0 until the temperature drops below the electron mass.After that, φ cannot track the instant changing of the potential and begins the damped oscillation scaled as |φ| ∝ T 3/2 .In this plot, we have not shown φ's oscillation afterward because its amplitude is too small to be visualized.For the green line (η ∼ 1), when T 1, φ gradually slides to the thermal potential minimum πf /2 following Eq.(30).When T ∼ 1, φ begins the damped oscillation obeying |φ| ∝ T 3/2 .The blue line (η 1) represents the standard misalignment where the thermal effect is negligible during the process.There, φ const until T 1.After that, φ begins to oscillate following |φ| ∝ T 3/2 .Right: The type-B model (y = −y).The red line represents the situation where η 1 and φ's initial condition is πf /2 φ rh 3πf /2.In this case, when T T * , φ does the damped oscillation |φ − πf | ∝ T 1/2 during which φ converges to πf .When T T * , or equivalently speaking, mT m0, we have V | φ=πf 0, so φ rolls down from πf and oscillates like |φ| ∝ T 3/2 .We can see that the red line's oscillation is postponed compared with the standard misalignment beginning at T 1.For this reason, the red line is in the phase of the trapped misalignment.The green line represents the situation where η ∼ 1 and πf /2 φ rh 3πf /2.In this case, φ gradually slides to the thermal potential minimum πf as shown in Eq. (30).At the point T 1, φ begins the damped oscillation scaled as |φ| ∝ T 3/2 .The blue line represents the standard misalignment mechanism where the thermal effect is inferior to the bare potential effect, i.e., η 1.Similar to the discussion of the blue line in the left panel, φ stays constant in the early universe and then begins the oscillation when T 1.Even though the red lines in the left and right panels are qualitatively different, the green and blue lines are quite similar in terms of the oscillation temperature ( T 1).The main difference between the green and blue lines is that the green lines' initial displacements are thermally determined to be |φ| osc |φmin|, but the blue lines' oscillation amplitude depends on their initial conditions.
: The left panel denoting the case m 0 10 −25 eV where ρ φ, local Fρ average , and the right panel representing the case m 0 10 −25 eV where ρ φ, local Fρ local .

Figure 7
Figure 7.The parameter space of the scalar-photon couplings.The nonrelativistic scalar relic's fraction among the dark matter is F 10 −3 .The kinetic mixing before the oscillation is FI 2×10−12 , which is the benchmark value for the freeze-in production of the dark photon dark matter satisfying Ω A h 2 0.12.Left: The type-A model (y = 0).We choose e = 1, 10 −3 , 10 −6 to demonstrate the experimental targets for the dark photon freeze-in through the varying kinetic mixing, which appear as the black lines.The thick dashed magenta line is the η 1 contour, upon which the scalar's early displacement is set to be πf /2 by the thermal misalignment.The gray region denotes the constraints from the equivalence principle tests (gray: EP tests) from the MICROSCOPE satellite[85] and the torsion balance experiments[83,84].The blue region represents current constraints from the clock comparisons[86-88, 90, 185].We can see, for the type-A model, most of the model's parameter space (e 1) is inside the projections of the proposed experiments, such as Lyman-α UVES (dashed purple)[170], the clock comparison (dashed pink: optical-optical, dashed red: optical-nuclear)[86], the cold-atom interferometer (dashed dark green: AEDGE, dashed green: AION-km, dashed light green: MAGIS-km)[165,166,168,186], and the resonant-mass oscillator (dashed orange: DUAL)[91].Right: The type-B model (y = −y) beginning with the wrong vacuum (πf /2 φ rh 3πf /2).We choose e = 10 −6 , 10 −8 , 10 −10 as the benchmark values to plot the black lines.The thick magenta line is the η 1 contour.In the region beyond such a magenta line, the scalar acquires the field displacement πf via the thermal misalignment.The thick dashed cyan line, i.e., the contour of mT |3H=m 0 m0, envelopes the region where the thermal effect postpones the oscillation starting, i.e., Tosc < T |3H=m 0 .This is the region where the ultralight scalar does the trapped misalignment, which is previously investigated in[140,141] for the axion models.Given the scalar abundance, since the trapped misalignment needs smaller oscillation amplitude, as discussed in Eq.(35), the power law of the black lines is changed within the upper right corner enveloped by the thick dashed cyan line.In the plot, one can find that the constraints from the cosmology (light orange: CMB, purple: BBN) are comparable with the current constraints from clock comparison experiments (blue)[86-88, 90, 162, 185]  and the equivalence principle tests (gray: EP tests)[83][84][85].For the type-B model, one needs e 10 −6 for the dark photon dark matter freeze-in to be detected by future experiments.

Figure 8 .
Figure 8.The schematic diagram of the ZN invariant UV model.The messenger particles Ψ ( )

Figure 9
Figure 9.This plot describes the ultralight scalar's evolution in the high-temperature environment and compares the numerical (solid lines) and analytical (dotted dark lines) results listed in Eq. (A4), Eq. (A5), and Eq.(A6).In the plot, similar to Fig.4, T is defined as T := T /T |3H=m 0 .|δφ| represents the absolute value of the deviation from the thermal potential's minimum.Without loss of generality, we choose the type-A model as an example in this plot, and the discussion of the type-B model is quite similar.Here, the red, green, and blue colors represent the η > 1, η = 1, and η < 1 cases, respectively, and one can find that the numerical and analytical results for each case match magnificently well.We find that the red and green lines have approximately the same power law, i.e., |δφ| ∝ T 1/2 , whose slight deviation can be explained by the log term in Eq. (A5).We can also find that the red line's oscillation period appears the same in the plot with the log-scaled 1/ T axis.This is because the only temperature-dependence in Eq. (A6) comes from log T , or, more intuitively, is caused by the same temperature power law of mT and H.

20 Figure 10 . 3 A
Figure 10.The plot showing how Tres and d log m 2 γ /d log T res change in terms of m A .Here, Tres represented by the red line is the temperature in which the resonant transition γ → A happens.Its numerical value can be read from the red vertical axis on the left-hand side of the plot.The dimensionless quantity d log m 2 γ /d log T res represents how fast the plasmon mass changes in terms of the temperature at which the resonant transition happens.Readers can refer to the blue axis on the right-hand side for its numerical value.In the mass range m A 10 5 eV, Tres ∝ m A and d log m 2 γ /d log T res 2 because m 2 γ ∝ T 2 when electrons are relativistic.When m A 10 −2 eV, Tres ∝ m 2/3 A and d log m 2 γ /d log T res 3, as a result of the plasmon mass power law m 2 γ ∝ n e − ∝ T 3 when the symmetric e − e + annihilate away but the asymmetric e − remains.In the middle mass range where 10 −2 eV m A 10 5 eV, the γ → A resonance happens when the plasmon mass is experiencing the exponential drop as the universe temperature goes below the electron mass.In this case, Tres is in the range 10 −2 me − me and insensitive to the value of m A .The semi-analytical estimations of Tres and d log m 2 γ /d log T res in such a case can be found in Eq. (B5).

1 . 2 A m 2 γA 3 2 • 5 3/ 2 2 3 • π 5 / 2 • g * ,S g 1 / 2 * 2 m 3 A m pl s 0 T 3 d
m A < 2me Knowing that the dominant contribution comes from the resonant transition of the transverse photon when m A < 2m e , we write down the Boltzmann equation ṅA + 3Hn A n γ Γ γ→A , average of γ → A transition probability and n γ = 2ζ(3)T 3 /π 2 is the transverse photon number density.The delta function in Γ γ→A reveals the feature of the resonant production: Most of the dark photons below the electron mass are produced when m .Doing the time integration of the Boltzmann equation, we get today's dark photon relic abundance Ω log m 2 γ /d log T ρ c res , (B2)

2 γ e 2 T 2 / 6 , 10 T res 8m A and d log m 2 γ
2 γ /d log T res , an order 1 − 10 factor, represents how fast m 2 γ passes through the vicinity of m 2 A .The numerical values of T res and d log m 2 γ /d log T res are shown in Fig. 10.For m A 10 5 eV, the resonant transition happens when electrons are relativistic.In this time, the plasmon mass is m based on which we have

FI ∼ 10 − 12 Ω A h 2 For 10 −2 eV m A 10 5 1 /2 π 3/ 2 e 2 m e m A − 1 and d log m 2 γ− 3 / 2 Ω A h 2
eV, the resonant transition happens when 10keV T m e .In this epoch, the plasmon mass is m 2 γ e 2 n e /m e , where n e g e (m e T /2π) 3/2 e −me/T .By solving the resonant condition, we have T res ∼ m e log 2 T res is insensitive to m A 's changing, given the dark photon relic abundance, ∝ m −3/2 A .More quantitative representation for the − m A relation is given as FI ∼ 10 −11 m A 1keV

2 / 3 A 3 and d log m 2 γ 11 FI ∼ 10 2 .) 3 d 3
. The concrete formulas for T res and d log m 2 γ /d log T res are T res 15keV × m A 10 −2 eV 2/photon relic abundance, the kinetic mixing is m A ≥ 2me For the heavy dark photon which satisfies m A ≥ 2m e , the inverse decay channel e − e + → A opens up and dominates over the γ → A channel.In such a case, the Boltzmann equation is ṅA + 3Hn A n e − n e + σ e − e + →A v , (B9) where the collision term is n e − n e + σ e − e + →A v g e − g e + d 3 p e − (2πp e + (2π) 3 f e − f e + σ e − e + →A v Møl .

n) 2 8π 3 m 3 A T 1 2 (
e − n e + σ e − e + →A v ( eusing the approximation K 1 (m A /T ) (πT /2m A )1/2 e −m A /T in the low temperature limit T m A and solving the Boltzmann equation Eq. (B9), we get the dark photon relic abundanceΩ A 3 4 • 5 5/2 2 17/2 • π 11/ 2F µν F µν .The table-top experiments testing the α em variation and the equivalence principle violation can be used to test this model because of the existence of φF µν F µν or φ 2 F µν F µν .In other situations where O SM is invariant under arbitrary global and gauge transformations, there are no more indices to contract, so the lowest order operator of the scalar-SM coupling is φ O SM or φ 2 O SM depending on whether the exact Z 2 symmetry, i.e., the invariance under φ → −φ, O DS → −O DS transformation, exists or not.One typical example is that