With the above motivation, the purpose of this paper is to investigate a possible UV completion of the anomaly-free axion with keV-scale mass and the decay constant at an intermediate scale. In particular, we consider the three Higgs doublet model (3HDM) as such UV completion. The set-up based on the 3HDM was proposed in Refs. [1]}, [2]}, but the detailed analysis of the Higgs sector and various phenomenological bounds
have not been investigated so far.
The 3HDM was first proposed in a context of CP violation in the Higgs sector in the pioneering works [3]}, [4]} where discrete symmetries were imposed for natural flavor conservation [5]}.
In the 3HDM, possible interaction terms and the number of parameters crucially depend on the imposed symmetry.
Possible symmetry groups implemented in 3HDMs are surveyed in Refs. [6]}, [7]}, [8]}, [9]}, [10]}, [11]}, [12]}, [13]}.
Furthermore, phenomenological studies are performed in the 3HDM with CP4 symmetry [14]}, [15]},
\(Z_3\)  symmetry [16]}, [17]}, [18]}, [19]},
\(Z_2\)  symmetries [20]}, [21]}, [22]}, [23]}, [24]}, [25]}, [26]},
Sp(6) symmetry [27]},
U(1)\(\times Z_2\)  flavor symmetry [28]},
and \(S_3\)  symmetry [29]}.
Also, eight extra Higgs bosons are predicted in addition to the SM-like Higgs boson with the mass of 125 GeV.
For this reason, the phenomenology of 3HDM is very rich, and we will focus on the symmetry and matter content that would fit with the anomaly-free axion.
